Predictions for Electromagnetic Counterparts to Neutron Star Mergers Discovered during LIGO-Virgo-KAGRA Observing Runs 4 and 5

We present a comprehensive, configurable open-source framework for estimating the rate of electromagnetic detection of kilonovae (KNe) associated with gravitational wave detections of binary neutron star (BNS) mergers. We simulate the current LIGO-Virgo-KAGRA (LVK) observing run (O4) using up-to-date sensitivity and up-time values as well as the next observing run (O5) using predicted sensitivities. We find the number of discoverable kilonovae during LVK O4 to be ${ 1}_{- 1}^{+ 4}$ or ${ 2 }_{- 2 }^{+ 3 }$, (at 90% confidence) depending on the distribution of NS masses in coalescing binaries, with the number increasing by an order of magnitude during O5 to ${ 19 }_{- 11 }^{+ 24 }$. Regardless of mass model, we predict at most five detectable KNe (at 95% confidence) in O4. We also produce optical and near-infrared light curves that correspond to the physical properties of each merging system. We have collated important information for allocating observing resources and directing search and follow-up observations including distributions of peak magnitudes in several broad bands and timescales for which specific facilities can detect each KN. The framework is easily adaptable, and new simulations can quickly be produced as input information such as merger rates and NS mass distributions are refined. Finally, we compare our suite of simulations to the thus-far completed portion of O4 (as of October 14, 2023), finding a median number of discoverable KNe of 0 and a 95-percentile upper limit of 2, consistent with no detection so far in O4.


INTRODUCTION
Currently, observable gravitational waves are primarily produced by the coalescence of binary compact objects (Abbott et al. 2016;Abbott et al. 2019bAbbott et al. , 2021b;; The LIGO Scientific Collaboration et al. 2021).Specifically, binary neutron star (BNS) mergers, like GW170817 (Abbott et al. 2017a), are of interest as these events can yield a post-merger, electromagnetic counterpart known as a kilonova (Abbott et al. 2017c).These transient events are fueled by the radioactive decay of heavy nuclei which are synthesized through r-process nucleosynthesis reactions possible given the neutron rich environment.Under certain conditions, black hole-neutron star mergers can produce kilonovae as well; however, it is much less likely (Fragione 2021), so we focus on BNS mergers here.
As the two neutron stars inspiral, they become tidally disrupted, causing neutron-rich material to be ejected from the system.The amount of material ejected depends, among other factors, on the ★ E-mail: vedgs2@illinois.eduequation of state (EOS) (Sekiguchi et al. 2015;Lattimer & Prakash 2016) being "stiff" or "soft" (Lattimer & Prakash 2016;Shibata 2016).A neutron star with a stiff EOS exhibits greater pressure, for a given density, and has larger radius causing it to experience greater tidal forces from its companion.In this work we use the SFHo EOS (Steiner et al. 2013a) used in Setzer et al. (2023) for modeling the kilonova population.Several spectral-energy density (SED) models exist that are parameterized by, for example, the mass and velocity of the ejecta, electron fraction, or opacity (e.g., Barnes & Kasen 2013;Kasen et al. 2015;Metzger 2017;Radice et al. 2018b).For this work we use the bns_m3_3comp model grid developed in Bulla (2019); Dietrich et al. (2020), henceforth referred to as the Bulla (2019) grid since it has consideration for observing constraints like viewing angles in its parameter space.
SSS2017a or AT2017gfo is the first optically confirmed kilonova from a binary neutron star merger (Coulter et al. 2017;Lipunov et al. 2017;Soares-Santos et al. 2017;Valenti et al. 2017), which was detected in conjunction with the gravitationalwave event GW170817 (Abbott et al. 2017b) and the gamma ray burst GRB170717A (Abbott et al. 2017d).This was a landmark discovery for the field of multi-messenger astronomy (MMA) since it was the first detection of a cosmic event via gravitational-waves, a kilonova, and gamma rays.
However, GW170817 remains the only such KN discovery to-date.This is in part due to the current limitations in GW event localization, the coordination required to perform proper follow-up, and the expected rarity of such events.Nonetheless, these events promise many scientific opportunities, such as studying the neutron stars and their EOS (Kilpatrick et al. 2017;Margalit & Metzger 2017a;Siebert et al. 2017;Radice et al. 2018a;Coughlin et al. 2019;Dietrich et al. 2020), understanding r-process nucleosynthesis (Chornock et al. 2017;Cowperthwaite et al. 2017a;Drout et al. 2017;Murguia-Berthier et al. 2017;Shappee et al. 2017), and measuring the expansion of the Universe (Abbott et al. 2018(Abbott et al. , 2019a;;Coughlin et al. 2020;Dietrich et al. 2020).But, to capitalize on these scientific promises, observers must be prepared to discover and follow-up future BNS events.Understanding the number of observable kilonovae expected during gravitational-wave observing runs would provide critical input for the follow-up efforts within the MMA community.
To address this need, we present here a new methodology to quantify the rates of observable kilonovae during the LVK's ongoing and future observing runs, complimenting similar analysis done recently (Colombo et al. 2022;Frostig et al. 2022;Weizmann Kiendrebeogo et al. 2023).We base our calculation on a number of factors to obtain realistic estimates of these rates (summarized in Figure 1).First, we sample from the appropriate distributions of BNS masses, astronomical extinction, merger rates, and distances adopted from the literature.We use these sampled parameters, either directly or as inputs to compute flux parameters, to perform interpolation on radiative transfer SED models that we then use to determine the likelihood of electromagnetic counterpart detection.We implement Monte Carlo trials to sample from the distributions in our parameter space to get the distributions of discovery and peak magnitudes, the distances of detected events, and the number of counterpart detections expected in the LVK O4 and O5 observing runs.The framework is also expandable and can support new parameter models, telescopes, and PSDs from future observing runs can be added as they become available.
The paper is outlined as follows.In Section 2, we detail our usage of existing SED models to build synthetic photometry.In Section 3, we describe the BNS parameter distributions that we use in our analysis.In Section 4, we explain our process of sampling from these distributions while accounting for instrumental downtime and other observational constraints.We present the resulting kilonova detection rates in Section 5.

SED APPROXIMATION
Running comprehensive, independent simulations (Kasen et al. 2017;Bulla 2019) to produce SEDs for each merger over all trials is computationally unfeasible.Thus, we use interpolation methods over existing SED grids to approximate the EM radiation.Bulla (2019) produced a model (Possis) for a grid of kilonovae SEDs simulated using three-dimensional Monte Carlo radiative transfer.These models are parameterised by two different components of the ejecta matter,  total ej : the lanthanide-rich dynamical component,  dyn ej which is released during the merger and the typically larger, lanthanide-free wind component,  wind ej released after the merger as a result of unbinding disk matter.Another parameter is the half-opening angle of the lanthanide-rich component of the dynamical ejecta, Φ, and the model has a dependence on the cosine of the observing angle, cos Θ. Dietrich et al. (2020) further improved the model to account for thermal efficiencies and time dependence for the temperature.We choose to use this model in our work.Thus, our SED model is parameterised by {m dyn ej , m wind ej , Φ, cosΘ}.It is important to note that our SED model does not have spin parameter.Although high spin values will have an effect on the resultant kilonova (Raaĳmakers et al. 2021), the vast majority of Milky Way neutron stars have very low spin (Zhu et al. 2018), suggesting that high-spin systems are uncommon.

Interpolation method
All SEDs from the above mentioned model, except those with Φ = 0 and Φ = 90, were used to create a 4 dimensional grid.These two Φ values were excluded as they either lack SEDs for different observing angles or are not available for all permutations of the  wind ej and  dyn ej .Table 1 describes the discrete points at which SEDs are computed using radiative transfer.While computing the SED for parameters not on the model grid, linear interpolation was used via a regular grid interpolator.The motivation here was that since the flux at every wavelength is known at several finely spaced points in our parameter space through robust simulations, it is reasonable to interpolate between two known points.We use linear interpolation since it is very fast to generate new SEDs on the fly which eliminates the need for pre-computing them, however other interpolation methods with different speed trade-offs also exist within packages like Nmma (Pang et al. 2022).
Distributions of the ejecta masses computed during our trials (Figure 2) indicate that a non-negligible fraction of binary neutron star mergers will produce  dyn ej that is greater than the maximum value on the grid (0.02 M ⊙ ) or  wind ej that is less than the minimum value on the grid (0.01 M ⊙ ), when sampling component BNS masses from realistic distributions.This necessitates some method for estimating SEDs when the m ej parameters fall outside the grid range.
Given the linear relationship between energy radiated and ejecta mass (Section 3.1 Barnes (2020); Equation 4Li & Paczyński (1998)), we have computed scaling laws for the total energy radiated for each cos Θ and Φ pair.If the  total ej from our BNS merger exceeds the grid limit, we use these linear laws to scale the closest grid SED.
If our  total ej is lower than the minimum  total ej value on the grid, we scale down the closest grid SED using a power law fit since it has the additional benefit of predicting zero flux when the  total ej = 0, according to: where  is the scaling factor and SED nn is the nearest neighboring SED: where  and  denote the slope and intercept for the linear fit respectively and  denotes the exponent for the power law fit.Note that all the best fit scaling parameters (i.e., , and ) were pre-computed for every pair of (Φ, cos Θ). Figure D1 shows the best fits for the linear and power scaling laws along with the relative errors for some pairs of (Φ, cos Θ).Table D2 and table D3 document all the parameters for linear and power law scaling respectively.We use this  total ej -dependent interpolation scheme instead of linear extrapolation beyond the regular grid range since a small negative slope over a large extrapolated grid range eventually result in negative fluxes at many wavelengths.These extrapolation artifacts result in non-physical SEDs.
Since we want to sample the (Φ, cos Θ) parameters from a continuous range rather than the discrete points computed above to avoid quantization, we fit a spline function to our data for all three scaling parameters (namely , , and ) using the smooth bivariate spline.Figure 3 shows the spline surfaces fit to the discrete points.We compute our scaling parameters from this surface for all our Monte Carlo trials.The sum of residuals from the surfaces are provided in Table D1.
Using this piecewise extrapolation method ensures that we can always get SEDs that have reliable total flux since we are scaling the SEDs based on the ejecta mass.However, this method fails to take into account any changes in color as a function of the  total ej since there was no obvious statistical trend for how the spectrum shifted.Doing this correctly will require updating the original radiative transfer simulations for a larger range of  ej values which is outside the scope of this paper.

Redshift and Extinction
Since most binary neutron star mergers are expected to be of extragalactic origin, we treat our SEDs for both host and Milky Way extinction.For the host galaxy, considered to be at rest, we use the CCM89Dust effect based on work from Cardelli et al. (1989).The  host − is computed for each SED using the   sampled from the distribution described in Equation 10 and   = 3.1 using the following equation For the Milky Way galaxy, considered to be the observing frame, we use the F99Dust effect to redden the SED based on work from Fitzpatrick (1999).Here, we use the SFD dust map Schlegel et al. (1998) to find the  MW − based on the RA and Dec of each event.Finally, we redshift the SED based on the luminosity distance of our kilonova.Details about how this distance is sampled are provided in the parameter distribution section (Section 3).All effects are applied to the SED within Sncosmo (Barbary et al. 2016).

PARAMETER DISTRIBUTION
Our Monte Carlo simulations use the aforementioned pipeline to estimate SEDs and produce the associated light curves.We sample the component neutron star masses in order to determine if the merger can be detected via gravitational waves.The component masses are also used to compute the  wind ej and  dyn ej which, in conjunction with the sampled Φ and cos Θ, are used to estimate the SEDs.Finally we sample   and the event coordinates to treat the SEDs and generate the synthetic observables.The distributions for all of these inputs are described below.Table 3 summarizes these distributions.

BNS mass distribution and computing ejecta mass
For our BNS pairs, we consider the standard formation scenario where binaries consist of a first-born recycled neutron star sped up from accretion (with mass  recycled ) and a second-born slow neutron star (with mass  slow ).
Analyzing the the mass distributions of these two distinct NS populations has been the subject of numerous studies.Farrow et al. ( 2019) used a two-peak Gaussian for the recycled NS (Table 2) and a flat distribution with the range [1.16, 1.42] M ⊙ for the slow NS, with further analysis done by Golomb & Talbot (2022).Galaudage et al. (2021) used a two-peak Gaussian to describe both the slow and recycled NS mass distributions (Table 2).The motivation behind the two peak Gaussian model for slow neutron stars is to reconcile the disagreement between the empirical galactic data on BNS pairs from radio sources (Farrow et al. 2019) and the distribution that would be required to explain gravitational wave events like GW190425 (Abbott et al. 2020).It is worth noting that both results use the same two peak Gaussian to explain the recycled NS distribution.
Our simulations support three different BNS mass models: the Farrow et al. ( 2019) and Galaudage et al. (2021) models discussed above (Figure 4) and a flat distribution with some astrophysical priors from the Kilopop package (Setzer et al. 2023).While presenting the results, we only use the Farrow et al. (2019) and Galaudage et al. (2021) models since the uniform distribution is not empirically motivated.Ultimately, we find that while the choice of mass distribution changes the distribution of ejecta properties, it makes little difference in the number of dicoverable KNe estimated by our simulations (Table 5).

Discussion of EOS
In addition to the components masses, the equation of state employed will also have an impact on the both the  TOV (Oppenheimer & Volkoff 1939) and the ejected matter.Consequently, the post merger remnant and the resulting kilonova also depend on the assumed EOS.A neutron star with a stiff EOS exhibits greater pressure, for a given density, and has larger radius causing it to experience greater tidal forces from its companion.This results in greater dynamical ejecta.The post-merger remnant influences the opacity of the ejecta (Kasen et al. 2015;Radice et al. 2018b).A long-lived neutron star remnant will emit neutrinos, causing the electron fraction to increase and thus reduce the production of lanthanides.In the case of GW170817, Margalit & Metzger (2017b) were able to constrain the nature of the merger remnant to be a short-lived hyper-massive neutron star through the gravitational wave and electromagnetic signals.
While there is still no clear EOS that is most favorable, significant work concerned with fitting functions (Setzer et al. 2023) to ejecta properties is done using the SFHo EOS (Steiner et al. 2013b).Since we make use of these fitting functions in our work, we opt to use the same EOS.It is important to note that the equation of state used in this work places constraints on both minimum and maximum masses for neutron stars.For this reason, we cut off the tails of our mass distributions at 1M ⊙ and 2.05M ⊙ respectively.
There is some evidence that the SFHo EOS results in low  TOV compared to empirical data (Foley et al. 2020) and thus may not capture the population diversity.However, recomputing the ejecta fits for different EOS models is outside the scope of this paper.

Computing ejecta masses
While the two NS masses are not themselves parameters for the chosen SED model, we use the masses to compute the ejecta parameters.Although we are not introducing any new ejecta fits in this work, we use this section to recapitulate the methodology for computing these ejecta parameters.Setzer et al. (2023) used the fitting function introduced by Coughlin et al. ( 2019) and data from 259 Numerical Relativity (NR) simulations (Radice et al. 2018c) to come up with the final fitting function for the dynamical ejecta mass which depends on the masses (M 1,2 ) and compactness (C 1,2 ) of the merging neutron stars and is given by log 10 ( where  = −0.0719, = 0.2116,  = −2.42,and  = −2.905,and [1 ↔ 2] refers to repetition of the preceding fit with the indices interchanged.The compactness of the neutron star is given by is given by Due to the logarithmic nature of our fit, two heavy neutron stars (∼ 2M ⊙ each) result in dynamical ejecta mass that can exceed 1M ⊙ .Metzger (2019) (Section 3.1.1)and work referenced within find that the total dynamical ejecta from BNS mergers lie in the range 10 −4 M ⊙ − 10 −2 M ⊙ .Thus, we limit the maximum dynamical ejecta from our mergers to 0.09M ⊙ which means our final The other ejecta parameter,  wind ej , is some fraction of the disk mass,  disk , where the two are related by and  is the unbinding efficiency which is sampled uniformly from the range 10-40% and  disk is computed as follows ( where  = −31.335, = −0.9760, = 1.0474, and  = 0.05957.Finally,  thr , the mass threshold for prompt black hole collapse, is computed using the following (Bauswein et al. 2013), If  1 +  2 ≥  thr , then we set both dynamical and wind ejecta to zero, ensuring there is no luminous remnant.

𝐴 𝑉 distribution
Modeling the extinction of kilonovae host galaxies based on empirical data is unlikely to yield accurate results given the single data point -host NGC 4993 for GW170817 (Pan et al. 2017).For this reason, we sample from an extinction distribution for galaxies known to host supernovae.Kessler et al. (2009) (Equation 18) inferred the mean reddening parameter for host galaxies using a SDSS-II sample of Type 1a supernovae and found that their extinction can be well-explained by an exponential function of the form where   = 0.334 ± 0.088.We sample from this distribution with a fixed   = 0.334.

Spatial and distance distributions
For each set of trials, we first define a cube of length  within which we simulate the events.We sample , , and  coordinates from a uniform, random distribution in the range of [− l 2 , + l 2 ] where  is specific to the simulation and detailed in Section 4. The Cartesian coordinates are converted to spherical coordinates to compute the event RA and Dec values.The euclidean distance is computed using the following: Phase (Days) The additional 0.05   term is added to ensure a minimum distance for events.We use these distances to compute the redshifts, assuming a flat Lambda-CDM cosmology 1 .Since we are only simulating mergers ≤ 500 Mpc, sampling distances is roughly equivalent to sampling redshifts.

Φ distribution
Φ describes the half-opening angle at which the lanthanide-rich component of the dynamical ejecta is distributed close to the equatorial plane during the merger giving rise to the red kilonova.The remainder of the ejecta is expelled further from the equator resulting in the blue kilonova.The lack of comprehensive data on kilonovae means that we do not have reliable distributions for Φ.For this reason, we choose to sample from a uniform continuous distribution of values for Φ in the range [15, 75], the entire range for values for which complete simulations exist (See Section 2.1).

cos Θ distribution
Given the random distribution of BNS mergers in space, it is safe to assume that the cosine of the observing angle will be uniformlydistributed.Specifically, we sample from a continuous distribution of cos Θ values in the range [0, 1].We also use the observing angle to compute the inclination, defined as the angle between the line of sight and the total angular momentum, (Chen et al. 2019) (Ω) which is used as a parameter in the GW waveform generation (Section 4.2).

MONTE CARLO TRIALS
For each trial, we sample the BNS component masses, distances, coordinates,   , Φ, and cos Θ from the distributions mentioned above.
The criteria used for determining gravitational wave and electromagnetic detections are specified below.We found that the expected numbers of discoverable kilonovae begin to converge after a few hundred iterations of our simulations for both the O4 and O5 observing runs.For this reason, we run all our simulations for 1000 iterations since we do not expect the results to change significantly with further increase in the number of trials.

Finding the number of events
For each trial, we first need to find the number of events, henceforth called  events where: The time is computed from the overlap duration of our chosen optical survey and the LVK observing run.Based on Table C2, we know that the maximum distances for BNS GW detections is ∼ 253 Mpc and ∼ 449 Mpc for LVK observing runs O4 and O5, respectively.Thus, the lengths of our event cube for simulations, , are set to 510 Mpc and 910 Mpc respectively.This ensures that the limiting factor for kilonova discovery is always either the sensitivity of the LVK detectors or the limiting magnitude of our survey.
The merger rate model is another configurable option in our simulations.Considerable work has been done to understand the frequency of BNS mergers (Nitz et al. 2021;Abbott et al. 2023).For this work, we set the BNS merger rate to 210 +240 −120 Gpc −3 yr −1 with log-normal uncertainties, in accordance with the LVK user guide2 (Figure 6) (Abbott et al. 2023).

Detecting gravitational waves from mergers
The maximum distance values at which a BNS merger is detectable is a function of the component masses of our binary system and its inclination (Chen et al. 2021).
Once we sample our mass distributions to find the component masses, we compute its gravitational waveform which acts as our signal.This waveform is parameterized by the masses of our coalescing binary and the inclination angle, Ω.Since we are dealing with binary neutron stars, we use the TaylorF2 waveform (Buonanno et al. 2009;Messina et al. 2019), which assumes neutron stars are non-spinning point masses and has been used for BNS merger rate modeling before (Nitz et al. 2021).Given that we only use the waveform to determine the maximum distance at which a detector would be able to discover a merger, the TaylorF2 waveform is sufficiently  accurate (compared to a more comprehensive model like IMRPhe-nomPv2_NRTidal) and much faster to compute, a key advantage for the speed of our Monte Carlo trials.
For each instrument, we use the PSD 2 which describes the noise at a given frequency.We then use a signal to noise ratio threshold of 8 (to remain consistent with LVK 2 ) (Petrov et al. 2022) to determine the maximum distance at which such a merger would produce a detection.Since we already know the luminosity distances for each of our mergers, we can determine if the event would produce a detection for each of the instruments.
Another aspect to take into account is the correlation in uptimes between the LVK detectors and their respective duty cycles.We used data on the correlation between different detectors from the LIGO O3a run3 to create a correlation matrix where the rows and columns are ordered by LIGO Hanford, LIGO Livingston, Virgo, and KAGRA respectively.Since we do not have duty cycle correlation data for KAGRA during O3, we assume the same ∼ 56 − 58% correlation as Virgo.These values are consistent with the current ∼ 58% duty cycle correlation between LIGO -Livingston and LIGO -Hanford reported for LVK O4 during the September 21, 2023 LVEM call4 .
1.0 0.56 0.56 0.56 0.56 1.0 0.58 0.58 0.56 0.58 1.0 0.56 0.56 0.58 0.56 1.0 We use this to create a detector uptime correlation matrix of dimension 4 × 4. For each trial we also create a matrix of random numbers in the range [0,1] of dimension  events × 4. We multiply the random numbers with the correlation matrix and scale all values to the range [0,1], using a min-max scaler, resulting in a matrix of dimensions  events × 4. Each column of this matrix is a series that represents the probability the detector is active when each merger in the trial is taking place.Since we already know the duty cycles for each detector (Table 4), we can set the detector to observe during an event if this probability during the event is less than the value of the duty cycle for that instrument.We can repeat this for all four   detectors and determine which instruments would be on during each of our mergers in a given trial.If a detector is on and observing during a merger and the merger is within the detection range for the component masses at their inclination, then we consider the event to be detected since we are not modeling any coincidental terrestrial noise, antenna patterns on the detector, or accounting for software failures.We then determine how many instruments will detect the merger by checking each detector's status and its detection capability.
It is important to note that non-detections, in the context of GW events, can encode vital information that can help improve sky map localization.For instance, the non-detection of GW170817 by Virgo, despite observing during the event, helped narrow down the localization (Abbott et al. 2017a).However, since we are not generating sky maps for this work, we choose to ignore antenna patterns while considering detections.
A duty cycle of 70% was used based on the observing capabilities as reported in the LVK Userguide 2 .However, both Virgo and KA-GRA will not be operating for the entirety of the 18 month period, so we encode this information into their duty cycles.duty cycle = 0.7 operating months 18 (15) Based on the latest observing plan 5 available at this time (dated October 14, 2023), we assume that Virgo will operating at optimum sensitivity for 12 months and KAGRA for 7 months.Table 4 details the duty cycles for all 4 detectors in the LVK network used in our simulations.
At this stage of the pipeline, we already have data on the BNS mergers that were detected in our simulation.Fig 7 shows the properties of these mergers.We find that the median number of GW detections for merging neutron stars over LVK O4 is ∼ 3 − 4, depending on the mass model used.

Detecting EM counterparts from mergers
As mentioned in the Section 3, we can find the  wind ej and  dyn ej for each merger.If the  total ej > 0, then we conclude that the merger has left behind a kilonova.Thus, we use the SED approximation method described in Section 2 to produce synthetic light curves for all mergers that have non-zero  total ej .We can use the synthetic photometry and the survey's detection thresholds to find the discovery magnitude and peak magnitude of each detectable kilonova.
Next, we need to account for the fraction of events that would be lost to light from the sun, called  sun loss .This value changes for different surveys and is a configurable option.We uniformly sample a number between 0 and 1 for each of the  events ; if this number is greater than  sun loss , then the event is not lost to the sun.Finally, we label an event as  good if it was detected by  GW detector(s), has non-zero ejecta, has a peak magnitude that can be detected by the survey, and is not lost to the sun since these filters provide all the necessary, but not sufficient, conditions for the discovery of the kilonova.Based on our simulations, typically ∼ 1 − 3% of the BNS merger have a prompt collapse to black holes, resulting in zero ejecta.
Since there are 4 detectors in the LVK network we will be doing our analysis for n = 1, 2, 3, and 4. It is worth noting that a single instrument detection will typically yield very poorly localized sky maps (on the order of ten thousand sq degrees) and lower network SNR which makes targeted search for kilonovae difficult.
This concludes the entire methodology we use in order to estimate the rate of discoverable BNS kilonovae.

RESULTS
In this section we will discuss the results of our Monte Carlo trials.It is important to note that these results indicate the best case scenario for kilonova detection since they do not account for observing inefficiencies (like weather, tiling etc.) or poor gravitational wave skymap localizations.

LVK O4 observing run
For this simulation we use both the Galaudage et al. (2021) and Farrow et al. (2019) mass distribution, a detection threshold of 23 magnitude, the DECam r passband for detections, the LVK User guide rates for BNS mergers, and a sun loss fraction of 0.5.In reality, the sun loss fraction is dependent on the specific follow-up survey we consider, the site(s) for ground-based facilities, and the location of the KN on the sky (and in particular, for space-based facilities the overlap between the GW localization and the allowed viewing area).Determining if an event will be lost to light from the sun must be done on an event by event basis, taking into account the survey strategy.Additionally, given that the search for real KNe will be done by a network of both public and private telescopes, which may elect to not share information about a counterpart for several hours  after discovery, makes modeling this effect infeasible without several assumptions.For the sake of simplicity, we choose to encode this information using a constant 0.5 fraction.As evident from Table 5 and Figure 9, the median number of mergers with detectable electromagnetic counterpart over all of the LVK O4 is ∼ 1 − 2, depending on the mass model used.The figure also shows the expected values and distributions for event distances and magnitudes while Table 5 breaks down the events by the number coincidental GW detections.
Additionally, we found the distribution for the discovery window, defined as the time for which the kilonova is brighter than the limiting magnitude of the survey, for 1, 2, and 3-detector events (Table 8).Even though the times shown are shorter than the 10+ days of observations obtained for GW170817 (Figure 5), we should still, on average, have several days to get detections of the EM counterpart.

5.82
Table 8.Discovery windows (in days) for KNe during LVK O4 for N = 1, 2, and 3 detector events in DECam r-band.The missing statistics for the 4detector events is due to a negligible number of mergers having coincidental detections on 4 instruments.All these discovery windows are significantly shorter than the 10+ days for which SSS17a was discoverable.

Looking ahead -LVK O5
LVK O5 presents the next opportunity for finding kilonovae, post-GW trigger.With the observing run slated to begin at the end of 2026 with a proposed end in the middle of 2029, we expect the Vera Rubin Observatory to be operational for the entirety of O5.This section aims to paint a picture of what the next ∼ 5 years of kilonova discovery could look like.Once again, we use the most updated PSDs for the O5 run from LVK which, notably, represent the high end targets of BNS ranges for LIGO and KAGRA and the low end for Virgo.The predictions presented in this section must be assessed with the added context that the sensitivities used are the targeted, optimistic values and real PSDs during O5 might not achieve these goals.For this reason, another analysis for LVK O5 will likely be required once the PSDs and observing plans are defined more concretely.Regardless, we present tentative numbers here since they will be useful for planning and forecasting.
With this caveat, we predict the median number of mergers with detectable electromagnetic counterpart over all of the LVK O5 to be ∼ 19.As evident from Figure 10 and Table 5, an updated LVK network during O5 with significantly improved sensitivities may present the first opportunity to discover a small sample of kilonovae which would enable exciting new population studies furthering both transient astronomy and cosmology.

Comparison with current LVK O4 results
To test the consistency of our pipeline with empirical observations, we simulate a partial LVK O4 run to compare with the actual ongoing (2021) mass model.Left: Distribution of the number of EM detectable events for 1, 2, and 3 GW detector events.Solid black line shows the distribution of the total number of discoverable kilonovae.Center: Distribution of distances of EM detectable events for 1, 2, and 3 GW detector events.Right: Distribution of peak and discovery magnitudes of EM detectable events for 1, 2, and 3 GW detector events.Solid lines represent distributions of peak magnitude, while dotted lines represent distributions of discovery magnitude.Table 5, Table 6, and Table 7 summarizes the results.
LVK O4 run.We adjust the simulation parameters to match the current O4 run by setting: • the duration to be ∼ 4.67 months, reflecting the current O4 period of May 24, 2023 to October 14, 2023; • the duty cycles for Virgo and KAGRA to zero, since they are not observing during the current period; • the duty cycles of the two LIGO detectors to 70% 2 .
All other configurable parameters mirror the full O4 simulation discussed in section 5.1.
This procedure allows us to assess the validity of our model, given that we have not detected any BNS mergers during the first ∼ 4.67 months of LVK O4 (as of October 14, 2023).Using the Galaudage et al. (2021) mass model, we found that the median number of disoverable kilonovae in the first ∼ 4.67 months of LVK O4, to be 0 +2 −0 with 1 +3 −1 BNS merger detections.

Comparison with complimentary work
Complimentary work has been done in the past to understand the detection rates of KN for surveys like the Zwicky Transient Facility (ZTF), the Wide-Field Infrared Transient Explorer (WINTER), and the Vera Rubin Observatory (Colombo et al. 2022;Frostig et al. 2022;Weizmann Kiendrebeogo et al. 2023).Table 9 summarizes these results while Table 5 reports the results from this work.Our independent analysis with the ZTF -r band, a limiting magnitude  2021) mass model predicts the number of discoverable kilonovae to be 1 +3 −1 over the 18 month LVK O4.

Retrospective analysis for LVC O2 and O3
Since our flexible framework can easily adopt different PSDs over arbitrary observing durations to simulate discovery rates, we have performed our analysis for the LVC O2 and O3 observing campaigns.We note that this since the BNS merger rates used in this work are themselves inferred from the results of the O2 and O3 runs, these simulations are somewhat self fulfilling.However, to the extent that this model accurately represent the true BNS merger rate, we can (2021) mass model.Left: Distribution of the number of EM detectable events for 1, 2, 3, and 4 GW detector events.Solid black line shows the distribution of the total number of discoverable kilonovae.Center: Distribution of distances of EM detectable events for 1, 2, 3, and 4 GW detector events.Right: Distribution of peak and discovery magnitudes of EM detectable events for 1, 2, 3, and 4 GW detector events.Solid lines represent distributions of peak magnitude, while dotted lines represent distributions of discovery magnitude.Table 5, Table 6, and Table 7 summarizes the results.We predict the number of detectable BNS mergers over the O2 run to be 1 +2 −1 with 0 +2 −0 discoverable KNe (90% credible).Only ∼32% of our trials had ≥ 1 discoverable KNe during this period, consistent with the single discovery of GW170817 (Abbott et al. 2017a) and SSS17a (Coulter et al. 2017) during O2.
For the O3 run, we predict the number of detectable BNS mergers to be 1 +3 −1 with 0 +2 −0 discoverable KNe (90% credible).Despite the increased sensitivity from O2, only ∼49% of our trials had ≥1 discoverable KNe during this period.Nevertheless, these results are also consistent with the single discovery of GW190425 (Abbott et al. 2020) and no KN discovery during O3.

CONCLUSION
These results paint kilonovae like SSS2017a with confirmed GW detections, gamma ray bursts, counterpart photometry, and spectroscopy as incredibly rare.The low intrinsic rate of detectable kilonovae is compounded by difficulties in locating the counterparts in poorly localized skymaps that, at this time, routinely have 90% confidence intervals that span on the order of 10 3 sq deg.Such large localization are simply impractical to probe efficiently without wide field of view surveys (Weizmann Kiendrebeogo et al. 2023).These inefficiencies are difficult to model accurately and thus the numbers presented here are upper limits.
Moreover, Table 8 illustrates how our window of opportunity for finding future kilonovae will likely be significantly shorter than GW170817's counterpart.These factors demonstrate the need for improved tooling, infrastructure, and search strategies for future KN discovery, such as Teglon8 and the systems described in Almualla et al. (2021), Bom et al. (2023) and Chatterjee et al. (2022).
Finally, a prompt chirp mass estimate from LVK, even if provided only to a tenth of a solar mass or with a small random offset applied, would allow forecasting of the electromagnetic signal.This, in turn, would enable observers to prioritize and coordinate follow-up resources more effectively, improving the yield of counterpart discoveries.Our synthetic photometry pipeline can be integrated into alerts systems, like the one described in Section A, to inform discovery and follow up strategies.

Figure 3 .
Figure 3. Left: Spline surfaces for the slopes () of the linear scaling laws for KN SEDs.Middle: Spline surfaces for the intercept () of the linear scaling laws.Right: Spline surfaces for the the exponents () of the power scaling laws.These laws are used to scale SEDs in cases where the ejecta masses exceed the grid limits of the SED model.

Figure 6 .
Figure 6.Distribution of BNS merger rates approximates the log normal distribution mentioned in the LVK user guide.

Figure 7 .
Figure 7. Predictions for BNS mergers detected by LVK during O4 using either the Farrow et al. (2019) (blue curves) or Galaudage et al. (2021) (orange curves) mass models.Top: Distribution of the number of GW-detected BNS mergers for our Monte-Carlo trials.The average number of mergers is represented by a vertical dashed line.Middle: Distribution of chirp massesfor all BNS mergers with GW detections.Bottom: Kernel density estimate contours (corresponding to 20%, 50%, and 80% of the probability mass) for the luminosity distance as a function of chirp mass for all BNS mergers with GW detections.GW170817 and GW190425 are represented by red and black points (dashed vertical lines), respectively, in the bottom (middle) panel.

Figure 8 .
Figure 8. Distribution of apparent peak magnitude for all kilonovae simulated, including those without an EM or GW detection (up to 40 magnitudes) for LVK O4.

Figure 9 .
Figure 9. LVK O4 Top pane shows results from using the Farrow et al. (2019) mass model while the bottom pane shows results from using the Galaudage et al.(2021) mass model.Left: Distribution of the number of EM detectable events for 1, 2, and 3 GW detector events.Solid black line shows the distribution of the total number of discoverable kilonovae.Center: Distribution of distances of EM detectable events for 1, 2, and 3 GW detector events.Right: Distribution of peak and discovery magnitudes of EM detectable events for 1, 2, and 3 GW detector events.Solid lines represent distributions of peak magnitude, while dotted lines represent distributions of discovery magnitude.Table5, Table6, and Table7summarizes the results.

Figure 10 .
Figure 10.LVK O5 Top pane shows results from using the Farrow et al. (2019) mass model while the bottom pane shows results from using the Galaudage et al.(2021) mass model.Left: Distribution of the number of EM detectable events for 1, 2, 3, and 4 GW detector events.Solid black line shows the distribution of the total number of discoverable kilonovae.Center: Distribution of distances of EM detectable events for 1, 2, 3, and 4 GW detector events.Right: Distribution of peak and discovery magnitudes of EM detectable events for 1, 2, 3, and 4 GW detector events.Solid lines represent distributions of peak magnitude, while dotted lines represent distributions of discovery magnitude.Table5, Table6, and Table7summarizes the results.

Figure D1 .
Figure D1.Scaling laws for all Φ and cos Θ = 0 or 0.1 pairs along with the ΔF/F errors.All parameters for the linear scaling laws and the power scaling laws are provided in table D2 and D3 respectively.

Table 1 .
List of Bulla SED grid values for each of the 4 parameters.

Table 2 .
BNS population parameters. defines the fraction of binaries in the low mass peak.
Galaudage et al. (2021) masses of BNS populations for LVK O4 and O5 simulation.The uniform distribution for the slow NS (blue) is fromFarrow et al. (2019).The two peak Gaussian for the slow NS (orange) is fromGalaudage et al. (2021).The distribution for the recycled NS (green) is shared by both models.
Figure 5.Light curves for GW170817 constructed using interpolated spectra with fit parameters computed by Dietrich et al. (2020), overplotted with real photometry.Appendix B provides a complete list of the sources for all the photometry used in this plot.

Table 3 .
Table of input parameters distributions used for LVK O4.

Table 4 .
Duty cycles for detectors used for the LVK O4 Monte Carlo simulations The average number of mergers is represented by a vertical dashed line.Middle: Distribution of chirp masses for all BNS mergers with GW detections.Bottom: Kernel density estimate contours (corresponding to 20%, 50%, and 80% of the probability mass) for the luminosity distance as a function of chirp mass for all BNS mergers with GW detections.GW170817 and GW190425 are represented by red and black points (dashed vertical lines), respectively, in the bottom (middle) panel.

Table 5 .
Summary of kilonova discovery estimates over LVK O4 and O5 from this work.Last column presents the expected number for KNe.All other columns in this table represent the middle 90% credible intervals with the median, 5  ℎ , and 95  ℎ percentile numbers being reported.

Table 6 .
Summary of mean luminosity distances (in Mpc) of discoverable KNe for LVK O4 and O5 based on simulations with 1, 2, 3, or 4 coincidental GW detection(s) and an EM detection.The missing statistics for the 4-detector events during LVK O4 is due to a negligible number of mergers having coincidental detections on 4 instruments.

Table 9 .
Comparison of results from analysis for kilonova detection rates from complimentary work for LVK O4. of 21.4, and the Galaudage et al. (

Table D3 .
table D2 and D3 respectively.Parameters for power scaling