Kilonova Parameter Estimation with LSST at Vera C. Rubin Observatory

The upcoming Vera Rubin Observatory’s Legacy Survey of Space and Time (LSST) opens a new opportunity to rapidly survey the southern sky at optical wavelengths (i.e., ugrizy bands). In this study, we aim to test the possibility of using LSST observations to constrain the mass and velocity of different kilonova (KN) ejecta components from the observation of a combined set of light curves from afterglows of γ-ray bursts and KNe. We used a sample of simulated light curves from the aforementioned events as they would have been seen during the LSST survey to study how the choice of observing strategies impacts the parameter estimation. We found that the design of observing strategy that is the best compromise between light-curve coverage, observed filters, and reliability of the fit involves a high number of visits with long-gap pairs of about 4 hr every two nights in the same or different filters. The features of the observing strategy will allow us to recognize the different stages of the evolution of the light curve and gather observations in at least three filters.


Introduction
The detection of GW170817, a binary neutron star (BNS) merger, using both gravitational waves (GWs) and photons, marked a groundbreaking milestone in multimessenger astronomy.Initially, GW170817 was identified solely by its gravitational-wave signal (Abbott et al. 2017a(Abbott et al. , 2017b)); subsequently, an array of electromagnetic (EM) signals from ground-based and space-borne telescopes covering the entire EM spectrum confirmed the presence of a luminous electromagnetic counterpart to the event (Abbott et al. 2017b).In particular, approximately 11 hr after the GW detection, the search for the EM signal of GW170817 led to the discovery of an electromagnetic counterpart named AT2017gfo associated with the GW signal (Cannon et al. 2012;Abbott et al. 2017b;Andreoni et al. 2017;Arcavi et al. 2017;Coulter et al. 2017;Díaz et al. 2017;Drout et al. 2017;Evans et al. 2017;Hu et al. 2017;Kilpatrick et al. 2017;Lipunov et al. 2017;McCully et al. 2017;Pian et al. 2017;Smartt et al. 2017;Soares-Santos et al. 2017;Tanvir et al. 2017;Troja et al. 2017;Utsumi et al. 2017;Valenti et al. 2017;Buckley et al. 2018).
This discovery played a crucial role in addressing numerous issues in high-energy astrophysics and fundamental physics.It significantly contributed to resolving the origins of short γ-ray bursts (sGRBs), the existence of kilonovae (KNe), and the processes behind heavy element synthesis.Additionally, it offered valuable independent constraints on two key aspects of astrophysics.First, it provided insights into the previously unknown equation of state of neutron stars (NSs), as discussed in Abbott et al. (2018).It also helped refine the understanding of the Hubble constant (Abbott et al. 2017c(Abbott et al. , 2021;;Cantiello et al. 2018;Fishbach et al. 2019;Hotokezaka et al. 2019;Kashyap et al. 2019;Coughlin et al. 2020aCoughlin et al. , 2020b;;Dietrich et al. 2020;Doctor 2020).These findings represent significant steps forward in our comprehension of the cosmos and have opened new avenues for future research in these fields.
The concept of KNe as transient phenomena powered by the radioactive decay of synthesized heavy r-process elements resulting from the ejection of neutron star matter during compact mergers was first highlighted by Lattimer & Schramm (1974).Several subsequent studies have contributed to our understanding of this event, including Li & Paczyński (1998), Freiburghaus et al. (1999), Lattimer & Prakash (2000), Metzger et al. (2010), Roberts et al. (2011), Tanaka & Hotokezaka (2013), Grossman et al. (2014), Metzger & Fernández (2014), Kasen et al. (2015), and Barnes et al. (2016).Along with the ejection of neutron-rich material, a relativistic jet is also produced.The jet, moving close to the speed of light, emits a powerful beam of γ-ray radiation, leading to a so-called sGRB prompt emission.As the jet interacts with the interstellar Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.medium, it decelerates and produces a detectable afterglow powered by synchrotron emission, observable from X-rays to radio frequencies.Before the detection of GRB 170817A, the connection between sGRBs and compact object mergers had been supported only by indirect evidence (Tanvir et al. 2013;Fong et al. 2015;Abbott et al. 2017b;Goldstein et al. 2017;Savchenko et al. 2017).However, the simultaneous detection of GWs and γ-rays demonstrated that at least a portion of sGRBs are indeed associated with the merging of BNSs.Currently, the conventional scheme for the progenitors of GRBs is a subject of debate, as counterexamples have emerged in recent years (Ahumada et al. 2021;Mei et al. 2022;Rastinejad et al. 2022;Troja et al. 2022;Yang et al. 2022).To complicate matters further, indirect detections of KN emission have been proposed, supported by the identification of optical and near-infrared (NIR) excesses in the flux of some GRB afterglows (Tanvir et al. 2013;Troja et al. 2019;Jin et al. 2020;Rossi et al. 2020;Rastinejad et al. 2022).
Constraining the properties of these particular events is crucial in resolving the ongoing debate regarding the dominant site for the production of r-process nuclei in the Universe.Some studies (e.g., Kasen et al. 2017;Anand et al. 2023) argue that BNS mergers are the primary source, while others (Siegel 2019) suggest the collapse of massive stars as the main contributor.
sGRB detection rates range between 10 and 40 per year for the GRB instruments on board the Neil Gehrels Swift Observatory (Gehrels et al. 2004) and the Fermi satellite respectively (Abdo et al. 2008).However, the optical counterparts for these bursts have proven to be elusive, mainly because the localization of Fermi sGRBs typically spans hundreds of square degrees (e.g., Mong et al. 2021;Ahumada et al. 2022).The follow-up of BNS and neutron star-black hole (NSBH) mergers detected by the International Gravitational-Wave Network, consisting of Advanced LIGO, Advanced Virgo, and KAGRA (LVK), during the third observing run (O3) has not been fruitful, possibly due to the fact that the GW sky maps are similarly large (Andreoni et al. 2019(Andreoni et al. , 2020a;;Coughlin et al. 2019;Goldstein et al. 2019;Gompertz et al. 2020;Kasliwal et al. 2020;Chang et al. 2021;Petrov et al. 2022).The lack of counterpart detections can therefore be explained by the fast fading nature of both KNe and afterglows, the large sky maps to observe, and the low local rate of compact binary mergers (Dichiara et al. 2020).
Empirical constraints on KN rates by optical surveys set an upper limit of R < 900 Gpc -3 yr -1 (Andreoni et al. 2020b;Andreoni et al. 2021) for KNe similar to AT2017gfo.Moreover, only a fraction of those will be detectable, as they could be beyond the detection limit of available telescopes.Obscuration and absorption by the Galactic plane is also a significant limitation to the detection of counterparts.In light of these constraints, the expected BNS detection rate of 4-80 events per year for the LVK network after 2020 (Kagra Collaboration et al. 2018b;Petrov et al. 2022), based only on GW searches, will likely provide only a few tens of detections throughout the next decade.
The new Legacy Survey of Space and Time (LSST; LSST Collaboration et al. 2009) is expected to be a game-changing facility in astrophysics.Time-domain astronomy will particularly benefit from the large ∼10 deg 2 field of view of the camera combined with the depth achievable with the 8.4 m diameter primary mirror, with an effective aperture of 6.423 m.
Depending on the choice of the LSST cadence, the project could unveil a large number of KNe and other types of fast fading transients (e.g., Andreoni et al. 2022).
The current estimates of KNe rates tell us that LSST will be able to detect ≈10 2 -10 3 events within z = 0.25 during the entire survey (Della Valle et al. 2018).Moreover, Andreoni et al. (2022) demonstrated that LSST is expected to find more than 300 KNe out to ≈1400 Mpc over a ten-year survey.Among those, we expect about 3-32 KNe recognizable as fastevolving transients similar to the one associated to GW170817.Furthermore, KNe have been analyzed only in association with other events such as GRB or GW detections, thus the possibility to detect and recognize such events is strictly related to the ability to survey as fast as possible the wide error boxes from GW signals and, once located, to promptly analyze their EM emission.
In spite of the technological and instrumental advances across multiple wavelengths, the fast-evolving nature of KNe will likely impede their spectral analysis.For this reason, this paper aims to analyze multiple observational strategies that only rely on photometry to derive physical parameters of KN sources without using spectra.Throughout this paper we assume that we know the location and the energy of the merger from other messengers (GW and GRB).
We consider that KNe can potentially be detected as a possible additional component to the optical and NIR afterglow of short GRBs in the temporal window that lasts from few hours to a few weeks after the onset of the burst (e.g., Kasen et al. 2015;Barnes et al. 2016;Fernández & Metzger 2016;Metzger 2017).This assumption follows the findings of Rossi et al. (2020) where the authors were able to isolate a golden sample of GRB afterglows, the behavior of which indicates the presence of a KN component in the afterglow light curves.However, they stated that strong constraints on redshift or NIR observation are needed to be able to find a KN contribution in the afterglow.In this paper we used a similar method to Rossi et al. (2020) for comparing GRB afterglow light curves with KNe in all the LSST observable bands to define an optimal observing strategy that can enhance the ability to detect KNe and characterize their sources, even using their indirect observations through GRB afterglows.
The paper is organized as follows: in Section 2 is described the impact of parameter estimation to understand the physics behind the KN explosion; in Section 3 the methodology of the simulations and analysis of the extracted data are described; in Section 4 the observing strategies are summarized; Sections 5 and 6 are dedicated to the analysis of the model and features of the observing strategies that impact the parameter estimation; Sections 7 and 8 are dedicated to the description of the obtained results; finally Section 9 is dedicated to the discussion of the results and the conclusion.

The Need for Parameter Estimation
KNe observations encode information on both ejecta properties and ejection mass processes that happen during the merger and post-merger (Metzger 2017;Villar et al. 2017;Coughlin et al. 2018).Part of the NS (whether the system is a BNS or an NSBH) can be expelled and become unbound (Davies et al. 1994;Rosswog et al. 1999).Tidal forces right before the merger can cause partial disruption of the NSs (or of the single NS in the case of the NSBH) with material launched at mildly relativistic velocities in the orbital plane of the system.Once the accretion disk is formed, neutrino radiation and nuclear recombination with magnetodynamic viscosity can drive mass outflow from the disk (see Ascenzi et al. 2021, and references therein).The velocity, the mass, and the geometry of the ejecta strictly depend on the properties of the system involved, thus a statistical study of a population of KNe would help us in understanding the distribution of those parameters.
Several studies of GW170817 attempted to infer the amount of material ejected (e.g., Alexander et al. 2017;Chornock et al. 2017;Cowperthwaite et al. 2017;Metzger 2017;Smartt et al. 2017;Tanvir et al. 2017;Coughlin et al. 2018Coughlin et al. , 2019;;Breschi et al. 2021;Heinzel et al. 2021;Ristic et al. 2022;Collins et al. 2023).However, the extent and the properties of the ejected material from this event remain uncertain and in considerable tension with theoretical expectations for the amount of each type of ejecta component (Radice et al. 2018;Korobkin et al. 2021;Nedora et al. 2021;Collins et al. 2023), likely in part because of underestimated uncertainties in these theoretical evaluations of ejecta (Henkel et al. 2023).KNe are rare and faint events compared to supernovae and other common classes of extragalactic transients, so they are hard to detect.Observational constraints, due to their intrinsic peculiarity, affect our understanding of the processes the source undergoes.This implies that we are not able to grasp the complexity of the compact object merger without making some important simplification, such as assuming a simplified treatment of the radiative transport that overlooks detailed three-dimensional anisotropic considerations.Similarly, opacity is considered in a simplified manner, without taking into account sophisticated nuclear reaction networks and composition variations (e.g., Wollaeger et al. 2018Wollaeger et al. , 2021)).Recent calculations using improved anisotropic radiative transfer and opacity calculations still arrive at similar conclusions for the description of the AT2017gfo event (Metzger 2017;Almualla et al. 2021;Heinzel et al. 2021;Ristic et al. 2022;Collins et al. 2023): a relatively high mass of the blue component is expelled from the poles while a significant mass of the red component is preferentially ejected toward the equator.
In particular, Ascenzi et al. (2019) show the posterior distribution of the KN parameters (ejecta mass, velocity, and lanthanide content) extracted from the observed multiwavelength afterglow light curves, assuming that the excess flux of the GRB afterglow light curve can be related to the radioactive decay from the KN ejecta.To estimate the parameters, the authors used a joint model of KN and GRB afterglow to reproduce the observed data; they produced a statistical distribution of the value from a sample of observed GRB afterglow light curves they claimed to be associated with KN events.The uncertainties on the derived parameters are considerable, primarily because of the sparsely populated light curves and the incomplete understanding of the nuclear reaction network responsible for producing the KN component.As a consequence of these uncertainties, certain cases result in parameter distributions that appear nearly uniform, lacking distinct patterns or trends.
Well constrained parameters will allow us to: (i) add features to the theoretical description of the event or (ii) break degeneracy between models (e.g., data from AT2017gfo agree with both two-and three-component ejecta models as shown in Cowperthwaite et al. 2017).The lack of both photometric and spectroscopic data poses a limitation, compelling us to address the challenge of constraining models solely based on photometric data.This situation highlights the opportunity to enhance observing strategies to optimize the chances of verifying theoretical predictions concerning KN events.By improving observational techniques and data collection, we can better assess and validate the theoretical models of these events.

Detection Rates in O4
Rubin Observatory's primary survey, named Wide Fast Deep (WFD), covers an area of 18,000 deg 2 through its "universal cadence," while approximately 10% of the observing time is reserved for other programs, including intensive observation of Deep Drilling Fields (DDFs).Compared to typical points on the sky, the DDFs will receive deeper coverage and more frequent temporal sampling in at least some of the LSST camera's ugrizy filters.
To estimate the detection rate of EM counterparts potentially detectable by LSST, we constructed a population of KNe and GRB afterglow light curves, starting from a realistic population of BNSs, following the method described in Colombo et al. (2022Colombo et al. ( , 2023)), summarized in Appendix A. We considered two sets of limiting magnitudes: one shallow and related to the expected depth of visits of LSST for WFD and one deep related to DDFs.
Assuming a network of LIGO, Virgo, and KAGRA detectors with the projected O4 sensitivities and a 70% uncorrelated duty cycle for each detector, we find a GW detection rate of 7.4 5.5 11.3 -+ yr −1 , 77% of which can produce a KN and 53% a relativistic jet.The shallow limiting magnitudes for WFD are sufficient to the detect the majority of KNe in the y and z bands and all the KNe in the i, r, g, and u bands, with a corresponding detection rate of 5.7 4.2 8.7 -+ yr −1 .In contrast, the fraction of events with a detectable otpical GRB afterglow is between 2% and 5%, with a maximum detection rate of 0.37 0.27 0.56 -+ yr −1 .This is due to the large abundance of off-axis jets within the estimated GW horizon, with a corresponding faintness of these counterparts.In Table 1 we report all the detection rates and the assumed limiting magnitudes.

Method
Motivated by the discovery of AT2017gfo and its luminous blue emission, several attempts were made to find similar cases in archival short-GRB observations (e.g., Troja et al. 2019;Rossi et al. 2020).For instance, Troja et al. (2019) found that some nearby events have optical luminosities comparable to AT2017gfo.In particular, they showed that the sGRB 150101B was a likely analog to GW170817, characterized by a latepeaking afterglow and a luminous optical KN emission, dominating at early times.This finding suggests that KNe similar to AT2017gfo could have been detectable in the optical spectrum even though they might not have been explicitly identified before the discovery of GW170817.Driven by this motive, we aim to study the performance of the KNe parameter estimation to better understand the physical properties of the KN ejecta and their evolution using LSST observational strategies.The work is divided into three major parts: 1. simulation of a sample of KN + GRB light curves (Section 3.1); 2. simulation of the observed light curves using a realistic cadence strategy (Section 4); 3. estimation of the parameters' variance using a Bayesian fitting algorithm to retrieve the posterior distribution (see Section 5).
The parameter estimation of a transient light curve is something that is usually done after the follow-up, so the ansatz here are that we already took care about distance estimation, contaminants, and candidate selections, thus we are only interested in analyzing how the search design impacts the parameter estimation.

Light-curve Simulations
Considering the combinations of KN and GRB events, simulated light curves allow us to build up a science case for KNe that are well localized and that have a constrained estimate for the distance.To approach the complexity of the KN models (Pang et al. 2020) we use nuclear-multimessengerastronomy algorithm, nmma13 (Pang et al. 2023).The software gives us the possibility to generate a distribution of realistic ejecta masses described by a population of BNS mergers, using the procedure from Dietrich et al. (2020), where they developed a framework to combine multiple constraints on the masses and radii of NSs, including data from GWs, EM observations, and theoretical nuclear physics calculations.The simulations to derive the prior distribution of the ejecta mass employ quasiequilibrium circular initial data in the constant rotational velocity approach, i.e., they are consistent with Einstein equations and in hydrodynamical equilibrium.The model assumes the SFHo (Steiner-Fischer-Hempel baseline model in Steiner et al. 2013) equation of state (EOS), which satisfies the current astrophysical constraints (e.g., Miller et al. 2019).
In our study, we employed the KN model introduced by Perego et al. (2017).This model takes into account the radiation produced by two distinct components: dynamical ejecta and disk ejecta.The disk ejecta can be further divided into two parts.The first part is wind ejecta (Ruffert et al. 1997;Kiuchi et al. 2015;Fernández et al. 2017), which is propelled in directions close to the polar axis by the neutrino flux originating from the hotter regions of the disk during the neutrino-dominated phase.The second part of the disk ejecta is known as secular ejecta (Fernández & Metzger 2013;Radice et al. 2018) and arises from viscous angular momentum transport.
By analyzing the distribution of ejecta masses, we were able to derive the distribution of ejecta velocities.This velocity distribution is influenced by the explosion energy, and for more details on this aspect, one can refer to Metzger (2017).
To take into account all the possible parameter correlations we construct distribution of priors for the KN parameters shown in Figure 1.The distribution of the model parameters we injected in nmma to simulate KN light curves, the ejecta mass, and the velocity distribution is modeled by the binary mass distribution in Dietrich et al. (2020); the open angle for the wind ejecta component is drawn from a uniform distribution ( ) 6, 4  p p .Then drawing randomly from these distributions we characterize the entire sample of simulated sources.
To create a science case to frame the experiment we made some assumptions: 1. we have the information on the distance to the event thanks to a GW trigger; 2. we have the information on the energy of the explosion because we detect a correlated GRB; 3. if there is an associated GRB we can assume the localization of the KN as known.We combine the simulated population with a single afterglow model for each viewing angle, so that the differences in the resulting light curves are due to the KN contribution (the model parameters for the afterglow and the KN are listed in Tables 2 and 3).This choice is driven by the possibility to detect KNe as flux excess in the afterglow evolution (Rossi et al. 2020), and because the frequency of observed afterglows exceeds the rate of KNe (≈5 KNe yr −1 versus ≈100 afterglows yr −1 ; see LSST Collaboration et al. 2009, as reference for the reported rate of afterglows) we expect to recognize KNe using well sampled afterglow light curves.All sources are simulated assuming three reference distances along the line of sight: 42, 100, and 300 Mpc.Eventually, the effect of the viewing angle cannot be neglected, thus we consider also three reference viewing angles of 0, π/4, and π/2 rad.The particular choice was made so that the ability to retrieve the physical parameters from the light curves would not be influenced by effects on the light curve due to the distance (e.g., selection effects due to the limiting depth of the survey or Malmquist bias) or viewing angle.The lower distance correspond to the distance value of AT2017gfo as a reference, while the median and higher distances were set based on considerations about the detectable sGRB rate; indeed Dichiara et al. (2020) -sGRBs within 200 Mpc are detectable.Scaling this value to greater distances, we estimated N sGRB (d L < 350 Mpc) ≈ 1.2 yr −1 , considering the lower bound of the uncertainty range; hence we set the higher reference distance to 300 Mpc.

Kilonova Model
nmma uses fitting formulae based on numerical simulation of the merger and post-merger dynamics to compute the ejecta properties (Radice et al. 2018;Krüger & Foucart 2020) as a function of the binary parameters (namely the component masses and the EOS).The procedure used is presented in Dietrich et al. (2020), where they survey 5000 EOSs that provide possible descriptions of the structure of NSs, recovering those that reproduce astrophysical constraints, such as NS maximum mass.For more details see the Supplementary Material of the referenced paper.
We then evaluate the accretion disk mass using the fitting formula from Barbieri et al. (2020), whose predictions are consistent with numerical simulations of both symmetric and The computation is based on a semianalytical model in which axisymmetry relative to the direction of the binary angular momentum is assumed.The ejecta, assumed to be in homologous expansion, are divided into polar angle bins, and thermal emission at the photosphere of each angular bin along radial rays is computed following Grossman et al. (2014) and Martin et al. (2015), taking into account the projection of the photosphere in each bin.See Table 2 for the distribution of reference parameters for the light-curve simulations.

Afterglow Model
GRBs associated with gravitational-wave events are, and will likely continue to be, viewed at a larger inclination than GRBs without detections of gravitational waves.As demonstrated by the afterglow of GW170817A, this requires an extension of the common GRB afterglow models, which typically assume emission from an on-axis top-hat jet.We used the Python package afterglowpy (Ryan et al. 2020), which characterizes the afterglows arising from structured jets, providing a framework covering both successful and choked jets.The temporal slope before the jet break is found to be a simple function of the ratio between the viewing angle and effective opening angle of the jet.
To accommodate an initial structure profile E(θ) in afterglowpy we consider the flux as a function of the polar angle θ.This assumes that each annulus of constant θ evolves independently, as an equivalent top hat of initial width θ j = θ.This is a very good approximation when transverse velocities are low: when the jet is ultrarelativistic and has not begun to spread, and when the jet is nonrelativistic and the spreading has ceased (van Eerten et al. 2010).The model allows for several angular structures of the GRB jet.In our exercise we use one GRB model for the afterglow since we are interested in the ability to infer the KN parameters.The parameters used for the GRB model are shown in Table 3, and we assumed a Gaussian jet structure so as not to correlate effects on the viewing angle or beam direction with effects coming from the peculiar jetenvironment interaction due to the particular geometry.The parameters are set to produce simulations that are as realistic as possible; for this reason we assumed the parameters from the Notes.The luminosity distance , D L , is needed to generate the model, as well as the angle for the polar emission of the wind ejecta, θ w , the ejecta velocity, v ej , the ejecta masses, M ej,dym and M ej,wind , the extinction, Ebv, and the exponent β of the relation , where M is the total mass, v is velocity of the mass envelope, and v 0 is the average minimum velocity of the ejecta.The last equation is used to reproduce the structure of the matter within the moving ejecta (see Perego et al. 2017, for details of the KN model).(2014).Thus the reference afterglow template can be considered representative of the short-GRB population.We also considered different reference values for short GRBs from Fong et al. (2015), but different assumptions appear not to dramatically impact the results, thus we considered the values in Table 3.

MAF and OpSim
The comprehensive discussion of the software made available by the Rubin Observatory for community contribution to the survey design is not within the scope of this paper.Interested readers are referred to the opening paper (Bianco et al. 2021) and its references for a full examination of the software's workings (Delgado et al. 2014;Delgado & Reuter 2016;Yoachim et al. 2016;Naghib et al. 2019) and for more details and information on this topic.
The Operations-Simulator software (OpSim14 ; Delgado et al. 2014) generates a simulated strategy based on a set of criteria, such as total number of images per field per filter, including simulated weather, telescope downtimes, and other occasional interruptions.The survey requirements (survey strategy) are the input to an OpSim run and the output is a database of observations with associated attributes (e.g., image 5σ depth) that specify a succession of simulated observations for the 10 yr survey.Since its creation, the Rubin OpSim has gone through various revisions, the main differences in which are the methods used to optimize the pointing sequences and filters to achieve the desired survey features (Bianco et al. 2021).
The Metric Analysis Framework (MAF15 ) API is a software package created by the Rubin Observatory (Jones et al. 2014) to evaluate how various simulated Legacy Survey of Space and Time (LSST) observation strategies impact different specific science goals.The MAF has been made public upon its creation to facilitate community input in the strategy design and it enables interaction with OpSim primarily by SQL, allowing the user to select filters or time ranges (e.g., the first year of the survey).Further, the choice of slicers allows the user to group observations.For example, one may "slice" the survey by equal-area spatial regions, using the HEALPIX scheme of Górski et al. (2005).Throughout, we choose a Healpix-elSlicer with resolution parameter NSIDE = 16, corresponding to a pixel area of 13.4 deg 2 (and thus the choice that most closely matches the size of the Rubin LSST field of view; see Ivezić et al. 2019, for reference to Rubin LSST characteristics).Thus, to pass from the simulated theoretical light curves produced according the procedure described in Section 3 to the simulated observed light curves we used the MAF.In this way we are able to apply parameter estimation tools to the simulated observation to analyze the impact of the observing strategies on the ability to retrieve parameters injected in nmma to produce the theoretical simulations (more details in Section 7).

Impacts of Observing Strategy on Parameter Estimation
The ability to populate light curves is very limited, as the number of filters and the number of detections in each filter vary depending on the intrinsic properties of the event (see Figure 2).This could impact our capability to infer the KN parameters.To evaluate the performance of nmma in estimating the KN parameters and to reproduce the injected light curve, we used the posterior's variance as a metric for the performance of the fitting procedure.We analyze three features that typically impact the performance of a sampler: 1. the number of detected points on the light curve; 2. the number of available filters; 3. the peak magnitude of the light curve.
In the upper left panel of Figure 3 we analyzed the lightcurve sampling by changing the time resolution of the template.
When we refer to the general trend, the figure shows that for light-curve data spaced more than ≈6 hr apart, dynamical ejecta is more poorly constrained than the other parameters, with the exception of two cases.However, in those two cases the values of ejecta mass and velocity are closer to the others, suggesting as filled colored points we used the baseline_ v2.0_10yrs strategy design.The triangles represent nondetections.The way the strategy plans to look at the footprint implies that we will be able to detect the events but we will not have the same ability to characterize because we will lose information about the color and the morphology of the light curve.Indeed, we simulate the light curve in all the bands but only the z-band appears to be detectable in this region and at the survey time.
that the effect could be related to a small fluctuation around the best-fit parameter configuration.This is interpreted as an effect of the fitting procedure: because of the time gap between the detections we miss the possibility of getting the maximum.However, constraining the rise and the fall of the light curve would help in constraining the KN model parameters.Specifically, the sampler cannot constrain the global minimum of the cost function (in our case the likelihood of the detections, see Appendix B) in the M dyn -v dyn -M disk -θ w hypercube when the light curve is not populated.This is because there is a large degeneracy of states that reproduce the same collection of fluxes (top panels in Figure 4).When the sampler can constrain the values of the parameters in the hypercube, the performance appear to be better (bottom panels in Figure 4).This happens for the other panels too; however, when analyzing the performance of the sampler with respect the number of filters (upper right panel in Figure 3) we see that this is not a very important observational feature in constraining the posterior's variance.For the experiment shown in upper right panel of Figure 3 we set the detection cadence at 9 hr.
Eventually the peak magnitudes have a similar impact to the upper left panel on the performance as shown in the bottom panel of Figure 3. Thus, this behavior can be interpreted similarly to the case shown in the upper left panel, considering that the closer to the limiting magnitude the peak is, the fewer detections of the light-curve features we get.The main problem is that the degeneracy in the parameter space produces a different minimum in the cost function because of the systematics we analyzed here.Thus, the area surveyed to find the global minimum is larger.Eventually we are able to infer the values of the injected parameter within some confidence level with the drawback of losing precision.

Impact of the Model's Description on Parameter Estimation
The assumptions on the radiation transport and on the nuclear network, which are at the foundation of a model that attempts to describe the observed event, can influence the ability of the cost function to find the global minimum of the parameter space that gives the best configuration of the parameters to describe them.Due to this connection between the model and the cost function, we analyze whether the behavior we highlighted in the previous paragraph and in Figure 3. Analysis of the fitting performance for a single light curve.Each panel shows the logarithm of the posterior's variance as a function of (i) the time gap between two consecutive points on the light curve (upper left panel), (ii) the number of available filters (upper right panel), and (iii) the maximum magnitude of the light curve (bottom panel).The plot can be read as follows: lower values on the y-axis represent better performances.Interpreting the plots, we find that the most information is gained when we populate the light curve with points that are very close in time, and the event is bright (upper left and bottom panels).There is a net improvement in the performance when considering all the filters together; however, in the other cases, the difference in performance is negligible.For each event, the signal-to-noise ratio (S/N) of the GW strain has been evaluated for the LVK detectors.For events above the S/N threshold in the population, M ej and v ej are estimated using Equations ( 18) and ( 22) in Radice et al. (2018) and considering the SFHo EOS: Taking into account that both dynamical and disk ejecta contribute to the mass of the ejecta, the total m ej is estimated as m ej = M dyn + M disk .We aim to analyze the impact of observations on constraining the mass ratio, thus with X = [m ej , v ej ] we estimated the uncertainties on the KN parameters: where we consider the uncertainties on the two merging NS masses to be equal in the last equation.We assume the uncertainties on the masses as Finally, we obtain the following: 1 -which is the ratio between the variance squared of the mass ratio and the ejecta mass or velocity-can be interpreted as the sensitivity to the system's photometric observation.Table 4 shows the form for ( ) A M q X 1 for q < 1 and q ≈ 1 related to ejecta mass or ejecta velocity estimations (see Figure 5. The inference of model parameters from a physical model assumes that the value inferred represents the value from the event's underlying model.However, due to oversimplifications in the theoretical treatment or technological limitation, this is not always true. Figure 5 shows how the model impacts on the parameter values (i.e., ejecta mass and ejecta velocity) when we try to infer them assuming we can measure other observables-the total binary mass and the mass ratio.If the physical model used to describe the event is highly degenerate for those parameters, we cannot distinguish between sets of parameters that produce the same set of observables, i.e., the uncertainties on the parameters are so high that the range of possible inferences related to that measure is very broad.Our case shows this is the case when we infer the ejecta velocity.Indeed, Figure 5(b) shows that if we survey the parameter space following v ej the fact that we have high uncertainty translates into a broader region of local minima in the cost function.Thus, with every inference we make to look for the global minimum we are likely to end up in a very similar state to where we started; this is because the uncertainty distribution is almost uniform, which means that whatever the true ejecta velocity is the chance of being in any other region close to that value is the same, meaning that we are likely to miss the global minimum.Conversely, Figure 5(a) shows the case in which we survey the parameter space in the direction of the ejecta mass m ej .This direction of the parameter space appears to be very helpful in constraining the specific parameter, showing that the uncertainty on the inferred parameter can change by an order of magnitude.Hance, the possibility of matching the global minimum can be higher if we survey the m ej direction of the hypercube.Eventually, we can conclude that we expect to have a much better constrained inference for ejecta mass than for ejecta velocity.

Observing Constraints from the Simulations
In Section 3 we described how we simulated light curves, with a combination of KN and afterglow emissions from the same source.These simulations are then used as reference templates to produce mock observed KN light curves during the operation time of LSST.To tackle the problem, we associate to the templates of KN + afterglow light curves an explosion time, uniformly chosen, within the 10 years of the survey over the whole observed sky.Eventually, the ability to discover fast and faint transients, such those we simulated, largely depends on the area observed-which in our experiment is the field of view of the pointings-, the depth of those observations, the cadence, and the filters adopted by the survey (all the simulated survey strategies are listed in Table 5).
Using the baseline as a test for our machinery, we applied the observational constraints from this OpSim to simulate what the observed light curves look like.We simulated for each reference distance and viewing angle a set of 100,000 KNe + afterglow events, with a total of 900,000 sources.
Figure 6 shows that the contours of the detectability regions change dramatically for further events, and that the best filters to follow up KN + afterglow events are g + z bands, due to the possibility to follow the events for more time up to 100 Mpc.Whereas optical filters are seen to perform well on the bluer side of the spectrum with closer events, the relative importance of bluer and infrared filters appears to be unchanged as the source's distance increases.Analyzing the detected eventsdefined as all the events with a light curve that has more than two detected observations in any filter-we find that as the source is more distant, events will be observable for a duration that is almost equal to that between two consecutive observations.Light curves with this duration above the S/N will have just two detected points at the specific wavelength, thus a multiwavelength analysis is mandatory to be able to use these data, otherwise no parameter estimation will be possible.Figure 6 does not change dramatically under the assumption of different viewing angles, thus for all the cases this figure appears to be a good reference for the description of the results.
Because of the nominal limiting magnitudes (LSST Collaboration et al. 2009) we have a very small time window in which to catch farther KNe with a maximum duration of ∼6 days in NIR bands for the simulated sources.The nonuniformity of the filter coverage through the entire survey impacts on the possibility of characterizing the explosions.Moreover, even though simulations are produced in all the six LSST bands ugrizy, light curves simulated in Figure 2 show only the filter that produced a detection for the specific event in that time window (MJD 60000-62000).Late-time evolution of the light curve will not be detectable for sources close to 300 Mpc; thus, to be able to constrain model parameters a higher priority would rather be given to closer sources if and when detected, because they will be characterized by a well populated light curve.The drawback that has to be stressed is the very small time window in which the light curve is detectable, which for closer sources (i.e., 42 Mpc in our simulations) is from 5 to 10 days and for farther sources (i.e., 300 Mpc in our simulations) from 1 to 5 days from the explosion.For a simulated event at 300 Mpc, the main features that are affected are the peak magnitude and the duration of the event above the limiting magnitude.This is because we get less flux from farther events and thus we reach the limiting magnitude of the survey earlier when observing the light curve's evolution.Similarly, this impacts on the fall rate distribution, which cannot be accurately measured in all the cases because we lose information on the late-time morphology of the light curve.

Constraints on Kilonova Model
Follow-up of the events in one or more filters will allow us to infer the model's parameters within the 3σ uncertainties.We perform the fit using dynesty (Koposov et al. 2022), a Python package to estimate Bayesian posteriors and evidences (marginal likelihoods) using dynamic nested sampling methods.By adaptively allocating samples based on posterior structure, dynamic nested sampling has the benefits of Markov Chain Monte Carlo (MCMC) algorithms that focus exclusively Table 4 The Form for ( ) A M q X 1 for q < 1 and q ≈ 1 Related to Ejecta Mass or Ejecta Velocity Estimations Notes.For an NS the assumption M M 0 i i * -» has been considered.The tidal deformability parameter, L, is modeled from Chatziioannou (2020).A relation between these variables can be found in Timmes et al. (1996): on posterior estimation while retaining nested samplingʼs ability to estimate evidence and sample from complex, multimodal distributions.Nested sampling is a method for estimating Bayesian evidence that was first proposed and developed by Skilling (2006).The basic idea is to approximate the evidence by integrating the prior in nested "shells" of constant likelihood.Unlike MCMC methods, which can only generate samples proportional to the posterior, nested sampling simultaneously estimates both the evidence and the posterior (see Appendix B). Figure 5 shows that the model itself acts as a source of uncertainty on the estimation of KN parameters, with a greater ratio of uncertainties on the ejecta mass and velocity as the total mass of the progenitor system grows.The behavior of the curve suggests that the uncertainties of the parameters tend to be lower when the system is more massive.However, when we consider a fixed mass for the progenitor system we measure smaller uncertainties on the ejecta mass as the mass ratio increases, while we note the opposite behavior for the ejecta velocity.The trend of the uncertainties with respect to the total mass of the system implies that there is high degeneracy on the value of the parameters when looking for the best configuration to replicate the observed light curves.This is because the model reproduces similar light curves for all the combination of parameters within the uncertainties' range.Figure 7, indeed, shows that for the majority of the events the uncertainty on all parameters is 100%.However, for simulated events with higher ejecta mass and smaller velocities, there is a tail of ≈10% of events with better constrained parameters.This is important because to constrain the EOS models we need very accurate measurement of KN parameters, so the 30% accuracy, even if good, is not sufficient; thus a dedicated target of opportunity (ToO) appears to be a necessity to follow up KNe from an external trigger if the baseline is the chosen strategy for the survey.To support this we compare, among the OpSims, the mean number of detection per filter (i.e., the higher this number, the better the accuracy; see Figure 8).
The results show that the OpSims that allow revisits within the same night have higher a fraction of well populated light curves, which implies a greater ability to constrain the uncertainties of the model parameters.The takeaway message from this work is that from survey observations we can expect to improve our detection ability, because we observe deeper and wider, changing the configuration of filters from time to time.However, to efficiently constrain the model's parameters we need to maximize the information content of our observations, improving the number of filters and the number of detections we have to observe the evolution of the event.

Discussion and Conclusions
This work aims to understand whether the LSST observing strategy can help collect data that will improve our understanding of KN sources.Sagués Carracedo et al. (2021) analyzed how to optimize the strategy for the distribution of filters and survey depth to boost the detection efficiency for these faint and fast-evolving transients.They explored the dependence on the mass of the ejecta, the geometry, the viewing angle, the wavelength coverage, and the source distance.Eventually they claim that the detection efficiency has a strong dependence on viewing angle, especially for filters blueward of the i band.This loss of sensitivity can be mitigated by early, deep observations.Efficient searches for the gri counterpart of KNe at ∼200 Mpc would require reaching a limiting magnitude mag 23 lim = mag within 5 days from explosion, to ensure good sensitivity over a wide range of the model phase space.Toward this end, Andreoni et al. (2022) analyze different choices of filter setting and exposure time, and they find that observations in redder izy bands are crucial for identification of nearby (within 300 Mpc) KNe that could be spectroscopically classified more easily than more distant sources.LSSTʼs potential for serendipitous KN discovery could be improved by increasing the efficiency with the use of individual 30 s exposures (as opposed to 2 × 15 s snap pairs), with the addition of red-band observations coupled with samenight observations in the g or r bands, and possibly with the further development of a new rolling-cadence strategy.However, even if detected, the KNe are not sampled enough to allow parameter estimation without ancillary data.
This work showed how to constrain the parameters of the KNe model using data from the LSST-Vera Rubin Observatory.This new facility is expected to push forward our knowledge of the physics of compact objects and improve the statistics of unique transient events such as KNe (Andreoni & Kool 2020b).However, to be able to deeply comprehend the compact objects' EOS and thermalization processes (Korobkin et al. 2012;Barnes et al. 2016), together with the energy-dependent photon opacities in r-process matter (Even et al. 2020;Tanaka et al. 2020), it is essential to constrain the uncertainties originating from various assumptions in the modeling.This is due to the complexity of the underlying physics, which is affected by diverse interactions and scales (see Metzger 2017, and references therein).
The possibility to extract information about the source of a KN event depends on a number of assumptions, including: 1. the KN model; 2. the available filters for the observation; 3. the distance; 4. the time window in which the event is above the observation limit.
The search for KNe light curves can be achieved through discovery of a transient during searches on the entire probability distribution map from the GW trigger or as a targeted search in a small number of specific and very limited regions of the sky.Below we refer to those two scenarios as "All-sky search" and "Targeted search." We analyzed the possibility to constrain the ejecta mass, velocity, and opacity from the photometric multiwavelength search for KN events assuming we know the distance and position of the KN from other messengers (i.e., GW, GRB).
There is a tradeoff between the fitting performance and the features of the observed light curves.In order to be able to extract precise fits of a model's parameters some conditions of measurability need to be satisfied.They can be summarized as follows: 1. recognize the different phases of the light-curve evolution; 2. gather observations in at least three filters.
The criteria listed above can be translated to survey strategy design: more revisits (with the best strategy considering longgap pairs of about 2 hr every two nights) in the same or different filters will weight differently in the observing strategy in the two scenarios (All-sky search and Targeted search).We use as a proxy of the performance the average number of detections per filter (see Appendix C for details on the figure of merit, FoM).We then normalize the FoM with respect to the FoM for the baseline, N balseline = 7 for both "All-sky search" and "Targeted search." For targeted searches, the possibility of considering at least one of the two criteria performs better with respect to the baseline, as is visible in Figure 9. Eventually the configuration of a ToO in the main strategy will improve our ability to extract precise information using only the photometric light curve to constrain the source's ejecta parameters.
As shown in the middle panel of Figure 9, when a targeted search is considered presto_gap (IDs = 63, 68, 74) and long_gaps_np (IDs = 5,12,19,37,40) OpSims almost double the performance of the baseline in constraining the model parameters for farther sources; this is due to the color information these strategies allow one to obtain.Indeed, the presto_gap_half adds a third visit within the same night for half the nights of the survey, with variations on the time interval between the first pair of visits (standard separation of 33 minutes) and the third visit.Among this family the best strategy is presto gap3.5, which consider triples spaced 3.5 hr apart (g + r, r + i, i + z are the initial pairs).long_gaps_np similarly extends the gap between the pair of visits, modifying it to a variable time period of between 2 and 7 hr.The pair of visits are both in the same filter, in any of griz (g + r, r + i, or i + z pairs).In some of the simulations, these long-gap visits are obtained throughout the survey, while for other simulations the longer time separations do not start until year 5.Among this family of OpSims the best performing strategy is long_gaps_nightsoff0, which considers long-gap pairs every night.When we analyze farther sources it appears that having long-gap pairs every 4-7 nights allows one to better constrain the model parameters, with long_gaps_np_ nightsoff7 being the best performing OpSim for this family in the case of a source at 300 Mpc.Overall the best performing OpSim is vary_gp_gpfrac0.30(ID = 104).The vary_gp family is a set of simulations that investigate the effect of varying the amount of survey time spent on covering the background (non-WFD-level) Galactic plane area.The combination of image quality, set of filters per observation, and cadence allow the best coverage for our simulated light curves.
Across the different panels in Figure 9 the difference in performances when the population of events is considered at different viewing angles (from top to bottom, / / 0, 4, 2 v q p p = ) can be analyzed.The general discussion still holds; however, we see that a worsening of the performance is evident.The long_gaps_np family appear to outshine the baseline in all cases.This indicates the importance of revisiting in the same night with different filter pairs to have a well constrained characterization of the events.
When an All-sky search (i.e., a search on the entire probability distribution map from the GW trigger) is considered, the main criterion for a better performance appears to be the homogeneity of the filter coverage (see Appendix C), meaning that strategies that respect the criteria in a higher number of regions perform better than the baseline.From .Comparison plot of all the v2.0 OpSims for an All-sky search, the performance is normalized with respect to the baseline.The metric counts the median number of detections per filter, which is used as a proxy to evaluate the strategy that will allow the most accurate parameter estimation, as described in Section 5. OpSim indexes shown on the x-axis are described in Table 5.

Figure 9.
Comparison plot of all the v2.0 OpSims for a targeted search in a fixed pointing in the sky.The performance is normalized with respect to the baseline, N balseline = 7.Each plot considers a population of KN + afterglow simulated with a fixed viewing angle; from top to bottom the viewing angle is / / 0, 4, 2 v q p p = .The metric counts the median number of detections per filter, which is used as a proxy to evaluate the strategy that will allow the most accurate parameter estimation, as described in Section 5. OpSim indexes shown on the x-axis are described in Table 5.
the results (see Figure 8) it is shown that the baseline performs better than almost all the OpSims; those with better performance differ from this strategy in varying the exposure time or the number of images per exposure to force the limiting magnitude to be homogeneous over the whole sky (see the vary family details, LSST Collaboration et al. 2009). 16ndeed, the best performing OpSim is multi_short, which takes four short (5 s) visits per filter in a row, and it stops after 12 short visits per filter in a year and it achieves ∼700 visits per pointing.
In short, the baseline is a great compromise among all the strategies for KNe science, and in future an improvement of the ability to constrain parameters of serendipitously discovered KN events is also foreseen.However, small changes of this strategy oriented to adding a third image for color information within the 4 hr gap or ad hoc ToO strategies to follow up the evolution of the light curve will enhance by a factor of 2 the ability of the targeted search to describe the KNe events with the most reliable KN models known up to now.

Figure 2 .
Figure 2.An example of a single template (thick lines) observed at (R.A., decl.)=(197.45,−23.38) (the position of NGC 4993, the host galaxy of AT2017gfo), at the three reference distances [42, 100, 300] Mpc during three time-windows through the 10 years of the survey.To reproduce the observed detections shown in the panels as filled colored points we used the baseline_ v2.0_10yrs strategy design.The triangles represent nondetections.The way the strategy plans to look at the footprint implies that we will be able to detect the events but we will not have the same ability to characterize because we will lose information about the color and the morphology of the light curve.Indeed, we simulate the light curve in all the bands but only the z-band appears to be detectable in this region and at the survey time.

Figure 5
Figure 5 is somehow related to the model used in the fitting procedure.We consider a BNS population fromDietrich et al. (2020) as pointed out in Section 3 and fit to currently available observational constraints from both GW-detected and Galactic BNS binaries as described in Appendix A inColombo et al. (2022).The merger rate is obtained by convolving the delay time distribution (represented as the time gap between the formation of the binary system and its merger) with the cosmic star formation rate(Madau & Dickinson 2014) normalized to the local rate density R 347 Gpc yr 0 256 536 3 1

Figure 4 .
Figure 4. Corner plots with an example for a worst-case scenario (top row) and best-case scenario (bottom row) when applying the fitting procedure.

Figure 5 .
Figure5.The sensitivity plot described in detail in Section 6.The panels show the sensitivity of the parameters to the variation of the BNS represented by the primary mass and the mass ratio, q.

Figure 6 .
Figure6.The distribution of the features we can extract from the observed light curves.Each panel represents the distribution of two features: peak magnitude and duration of the light curve above the limiting magnitude; different lines represent a distance at which the event is simulated.Details in Section 7.

Figure 7 .
Figure 7. Distribution of the uncertainties of the model parameters from the fit.Top panels represent the PDF of the parameters' relative error; bottom panels represent the cumulative density function (CDF) of the same variable.

Figure 8
Figure8.Comparison plot of all the v2.0 OpSims for an All-sky search, the performance is normalized with respect to the baseline.The metric counts the median number of detections per filter, which is used as a proxy to evaluate the strategy that will allow the most accurate parameter estimation, as described in Section 5. OpSim indexes shown on the x-axis are described in Table5.

Table 1
Estimation of the Expected Observation Rate of KNe in Association with a GW and GRB Afterglow Notes.Details of how the single values in the table have been estimated are given in Appendix A. Below each rate, we report in parentheses the fraction of the total O4 BNS GW rate (HLVK O4).The GW detection limits refer to the S/N net threshold.Limiting magnitudes of LSST filters are in the AB system(LSST Collaboration  et al. 2009); detection rates are in units per year.The reported errors, given at the 90% credible level, stem from the uncertainty of the overall merger rate, while systematic errors are not included.These results and the underlying methodology are described in Appendix A.
Figure 1.The distribution of the model parameters we injected in nmma to simulate KN light curves, the ejecta mass and the velocity distribution is modeled by the binary mass distribution in Dietrich et al. (2020); the open angle for the wind ejecta component is drawn from a uniform

Table 2
The KN Parameters and Their Probability Density Functions (PDFs) Used to Create Simulated Light Curves for the KN Components with nmma

Table 3
van Eerten et al. (2010), D'Avanzo et al. (2014that Were Used to Create Simulated Light Curves for the GRBs The luminosity distance, D L , is needed to generate the model, as well as the viewing angle, θ v , the half-opening angle, f c , the outer truncation angle, θ w , the isotropic-equivalent energy, E 0 , the circumburst density, n 0 , the electron energy distribution index, p, and the fraction of energy imparted both to the electrons, ò e , and to the magnetic field, ò B , by the shock.usualvaluesfrom both observed and modeled afterglow populations as expressed invan Eerten et al. (2010), D'Avanzo et al. (2014), and van Eerten (2018).However specific, the parameters selected for the reference afterglow template are the median values of an observed short-GRB population taken as reference from D'Avanzo et al.

Table 5
Opsim's Names and IDs, as Plotted in Figures9 and 8