Updated observing scenarios and multi-messenger implications for the International Gravitational-wave Network’s O4 and O5

Advanced LIGO and Virgo’s third observing run brought another binary neutron star merger (BNS) and the ﬁrst neutron-star black-hole (NSBH) mergers. While no conﬁrmed kilonovae (KNe) was identiﬁed in conjunction with any of these events, continued improvements of analyses surrounding GW170817 allow us to project constraints on the Hubble Constant ( H 0 ), the Galactic enrichment from r -process nucleosynthesis, and ultra-dense matter possible from forthcoming events. Here, we describe the expected constraints based on the latest expected event rates from the international gravitational-wave network (IGWN) and analyses of GW170817. We show the expected detection rate of gravitational waves and their counterparts, as well as how sensitive potential constraints are to the observed numbers of counterparts. We intend this analysis as support for the community when creating scientiﬁcally-driven electromagnetic follow-up proposals. During the next observing run O4, we predict an annual detection rate of electromagnetic counterparts from BNS of 0 . 43 +0 . 58 − 0 . 26 (1 . 97 +2 . 68 − 1 . 2 ) for the Zwicky Transient Facility (Rubin Observatory).

However, triggering target-of-opportunity (ToO) observations on gravitational-wave (GW) events comes at the cost of precious telescope time that could otherwise be employed for alternative science cases.Therefore, to make the most of available ToO time, we must understand how the potential targeted observations contribute to our specific scientific goals.To support this effort, "observing scenarios" are produced to simulate the detection and localization of GW events, (e.g., Singer et al. 2014;Abbott et al. 2016).Petrov et al. (2022) recently produced observing scenarios tuned to open public alerts from the third observing run (O3), improving the consistency of the localization performance in O3 by simulating the actual GW signal-to-noise ratio (S/N) threshold used during O3 and allowing for the inclusion of single detector searches (Godwin et al. 2020;Nitz et al. 2020).Studies like Nissanke et al. (2013b), Petrov et al. (2022) and Colombo et al. (2022), which provide a set of simulated merger signals detected by the international gravitational-wave network (IGWN) during each observing run, lend the ability to realistically predict how well we can address specific scientific questions pertaining to the nature of compact objects, r-process nucleosynthesis, and the expansion rate of the Universe within the next decade.
In this paper, we describe the simulations produced for the observing scenarios currently available to the user community in the IGWN User Guide 1 , as well as simulate potential science constraints based on self-consistent counterpart search simulations.In Sec. 2, we summarize the simulations of the observation scenarios expected during the next observation campaigns.In Sec. 3, we report the results of our simulations on the tracking of 1 https://emfollow.docs.ligo.org/userguide/capabilities.htmlEM counterparts of GWs by optical telescopes, notably the Zwicky Transient Facility (ZTF; Bellm et al. 2019;Graham et al. 2019;Masci et al. 2019;Dekany et al. 2020), a time-domain optical survey with a very wide field of view (FOV) of 47 deg 2 mounted on the Samuel Oschin 48-inch (1.2 m) Telescope at the Palomar Mountain, and the future Vera C. Rubin Observatory's Legacy Survey of Space and Time (Rubin Observatory; Ivezić et al. 2019), a large (8.4 m), wide-field (9.6 deg 2 FOV) ground-based telescope designed to conduct deep 10 yr survey of the Southern sky.In Sec.3.4, we present an estimate of how future multimessenger observations during O4 and O5 will help us to measure the Hubble constant H 0 .We present our conclusions in Sec. 4.

Overview
Here, we perform detailed simulations for the upcoming fourth and fifth LIGO-Virgo-KAGRA (LIGO;LSC et al. 2015, Virgo;Acernese et al. 2014 and KAGRA;Akutsu et al. 2021), observing runs (O4 and O5, respectively) and present the multi-messenger constraints on the Hubble constant that can be derived from future events.The continuous improvement of GW detector sensitivities allows probing farther into the Universe to detect more compact binary coalescences (CBCs).Furthermore, since KAGRA joins the LIGO and Virgo detectors for the next observing runs, we will have a total of 4 detectors online, which might result in an increased detection probability depending on KAGRA's sensitivity.Following Petrov et al. (2022), we simulate realistic astrophysical distributions of the mass, spin, and sky locations of CBCs by assessing the likelihood of detection for the networks considered.We can estimate the distributions of sky-localization areas and distances that we expect for detected events, as well as the rate of GW event detection.We will use two characteristic surveys to assess counterpart detection chances, i.e., ZTF, which will be sensitive to BNS mergers similar to GW170817/AT2017gfo up to ∼ 300 Mpc (Coughlin et al. 2019b;Anand et al. 2020;Kasliwal et al. 2020), andRubin Observatory (Andreoni et al. 2022), which will observe well beyond the IGWN horizon of current GW detectors.
Figure 1 provides a flowchart for the observing scenarios pipeline and summarizes the overall process.Each step of the workflow will be further described in the following subsections.

Population Models
In this section, we outline the processes used to generate two distributions of CBCs which we will compare Each distribution consists of three astrophysical subpopulations of CBCs: BNS, NSBH, and BBH.These subpopulations are separated based on the masses of their components: m 1 is the mass of the primary component, and m 2 is the mass of the secondary, with m 2 ≤ m 1 by definition.The mass of a nonrotating neutron star cannot exceed the Tolman-Oppenheimer-Volkoff (TOV) limit M max, TOV ≈ 2 − 2.5 M ⊙ (Legred et al. 2021;Ruiz et al. 2018;Shibata et al. 2019;Margalit & Metzger 2017b;Rezzolla et al. 2018;Dietrich et al. 2020b).However, rotating NSs can exceed this limit (Baumgarte et al. 2000;Stergioulas 2003).Ad-ditionally, a population analysis of all CBCs detected by the IGWN finds a sharp feature in the mass distribution of compact objects at 2.4 +0.5 −0.5 M ⊙ (90% credible interval) and interprets this feature as a delineation between NSs and BHs, due to its proximity to the TOV limit (Farah et al. 2022;Abbott et al. 2023).In order to provide a conservative upper bound on the number of NS-containing events so that follow-up programs can make optimistic plans for observing EM counterparts, we take the 95% upper bound on the location of the feature found by Farah et al. (2022) and choose to delineate between NSs and BHs at 3 M ⊙ .This high value comes at the potential expense of contaminating the BNS and NSBH samples with a few low-mass BBHs, but we find this preferable to the possibility of wrongly classifying an event that could result in a bright EM counterpart as a BBH.The choice of choosing the subpopulation boundary to be 3 M ⊙ also maintains a consistency with previous analyses (Abbott et al. 2020b).
We follow Farah et al. (2022), who proposed a resolution to thisNS-BH discrepancy by using all publicly available CBCs in the GWTC-2.1 catalog in a single population analysis, thereby foregoing the need for a priori classifications and instead allowing the population fit to pick out distinct subpopulations of CBCs.We use a similar procedure.Figure 2 shows all publicly available CBCs in the GWTC-3 catalog (Abbott et al. 2023).The ambiguously classified event GW190814 is shown in dark purple.The gray band indicates the approximate location of the purported lower-mass gap.GW190814 is the only event within this region at more than 90% credibility LRR distribution.Here, as in Petrov et al. (2022) and Abbott et al. (2020b), we use separate models to describe each CBC subpopulation (BNS, NSBH, and BBH).For the BNS population, we draw from a truncated Gaussian mass distribution centered at 1.33 M ⊙ ( Özel & Freire 2016) from the interval [1, 3] M ⊙ with a standard deviation of 0.09 M ⊙ ( p(m) ∼ N (µ, σ 2 ) ) and spin uniformly distributed magnitudes in the interval [0, 0.05].For the BBHs, sampling is performed from [3, 50] M ⊙ (Abbott et al. 2019a) using a truncated power-law distribution p(m) ∝ m a with a = −2.3(Salpeter 1955).Both masses are independently drawn from this distribution and paired according to the rule m 2 ≤ m 1 .The spins are also uniformly distributed, with magnitudes smaller than 0.99.In both cases, the spins are either aligned or antialigned, i.e., we neglect the possibility of misaligned spin and precessing systems in this work.Lastly, the NSBH population mass and spin distributions are described by drawing one component each from the BNS and BBH distributions above.For each of the three subpopulations, we draw 10 6 samples.
GWTC-3 distribution.This distribution is drawn from a model that describes the full population as a continuous function, foregoing the need to specify different models for each individual subpopulation (Fishbach et al. 2020).The mass and spin distributions are described by the PDB model from Farah et al. (2022), andAbbott et al. (2023).This model consists of a broken power law with a notch filter n(m|M gap low , M gap high , A) that suppresses the merger rate between NSs and BHs (M gap low and M gap high ) by a factor of A (Fishbach et al. 2020;Farah et al. 2022;Abbott et al. 2023), allowing for a potential lower-mass gap in that region ( Özel et al. 2010;Farr et al. 2011).It additionally includes a low-pass filter at the upper end of masses of black holes to take into account a possible tapering of the mass distribution at these locations.The component mass distribution is then as follows: where h(m|m min , η min ) and l(m|m max , η max ) are the lowmass and high-mass tapering functions, respectively.The 1D mass distribution, p(m|λ), is shown in Figure 3 for a specific choice of λ (λ represents the 12 parameters of the model; see Tab. 7 in Appendix A.1).
The 2D mass distribution is constructed by assuming that both the primary and secondary masses are drawn from p(m|λ) and related via a "pairing function" (Fishbach & Holz 2020;Doctor et al. 2020).As defined in Fishbach & Holz (2020), the pairing assumed here is a power law in the mass ratio, q ≡ m 1 /m 2 .Explicitly, The values of the hyperparameters Λ = {λ, β} are listed in Appendix A.1 and were chosen by fitting this model to all CBCs in GWTC-3 and choosing the maximum a posteriori value for Λ.The effects of neglecting the hyperparameter uncertainty are estimated in Appendix A.1.
The PDB model assumes a spin distribution with isotropically oriented component spins and uniform component spin magnitudes.The spin magnitude distribution for objects with masses less than (m < 2.5 M ⊙ ) is defined in the range of [0, 0.4], and that for objects with masses larger than 2.5 M ⊙ is defined in the range [0, 1].
A set of 106 CBCs were drawn from the PDB model, constituting the GWTC-3 distribution.These samples were then split into the three subpopulations by defining neutron stars as objects with masses below 3 M ⊙ and black holes as objects with masses above 3 M ⊙ .This yields 892, 762 BNSs,35,962 NSBHs,and 71,276 BBHs.One resulting difference between the LRR and GWTC-3 distributions is that the LRR distribution is drawn from a model defined only below 50 M ⊙ whereas the GWTC-3 distribution is drawn from a model defined up to 100 M ⊙ , allowing for higher-mass black holes in the latter case (though the tapering of the PDB mass distribution above M max = 54.38 does somewhat limit the number of high-mass black holes).
Figure 4 shows the mass distributions of the CBCs subpopulations.Figure 5 shows the simulated mass dis- tributions for each model that survives the S/N cut (see Section 2.3).

Simulation Campaign
We use the public software suite ligo.skymap2, which provides tools to read, write, generate, and visualize GW sky maps from the IGWN.After having drawn the intrinsic parameters, masses and spins of the CBCs from each of our distributions (LRR and PBD/GWTC-3 (GWTC-3)), we distribute all samples uniformly in comoving volume and isotropically in both sky location and orbital orientation.This choice reflects our expectation that GW sources are not spatially clustered or preferentially facing toward or away from Earth (Ade et al. 2016).We add Gaussian noise to the simulation of the GW signals of our CBCs.All the source code to reproduce these simulations, from the drawing of intrinsic parameters3 to the statistical production of sky maps4 , passing successively by the filtering of CBC events that pass the S/N cut5 , as well as their location in the sky 6 Following Petrov et al. (2022), we apply an S/N threshold of 8 for the entire GWTC-3 distribution, and the BNSs and NSBHs populations of the LRR distribution.This S/N threshold is set to 9 for the BBH population of the LRR distribution, consistent with the localization area and distance distributions of O3 alerts (Petrov et al. 2022).This simulation set yields estimates of the GW sky-localization area for all subpopulations, the luminosity distance, and the comoving volume.We provide a 90% credible prediction of the comoving region and volume, containing the total posterior probability.As in Abbott et al. (2020b) and Petrov et al. (2022), the localization of the sky area is provided by Bayestar, the rapid localization code used in production IGWN alerts (Singer & Price 2016).
According to the IGWN, four detectors, namely LIGO Hanford, LIGO Livingston, Virgo, and KAGRA (Ab-bott et al. 2018), will be engaged during the next two observing campaigns, O4 and O5.In our simulation, we adopt this configuration, along with the assumption that the four detectors each have a 70% operating cycle, independently of each other.However, the recent update states that KAGRA will start the run with LIGO Hanford, LIGO Livingston and Virgo, then return to extended commissioning to rejoin with greater sensitivity late in O4.The noise power spectral density (PSD), also known as sensitivity curves, is applied to each observation run and for each detector 9 .We use the publicly available noise curves released in LIGO-T2200043-v3 10 .For O4, we used aligo O4high.txt,avirgo O4high NEW.txt, kagra 10Mpc.txtrespectively for LIGO (LHO, LLO), Virgo, and KAGRA respectively, while, for O5, we used AplusDesign.txt,avirgo O5low NEW.txt, and kagra 128Mpc.txt.
In order to measure the performance of the different interferometers, we used those sensitivities to calculate the BNS inspiral range of 1.4 M ⊙ binary system detected with S/N = 8, during the next observation O4 and O5.The distances (in megaparsecs) from the BNS inspiral range are recorded in Table 1.
For the simulations, we must assume an astrophysical merger rate (taken to be in a frame that is comoving with the Hubble flow).As in Petrov et al. (2022), it is averaged over an assumed nonevolving mass and spin distribution.In this set, we use the merger rate per unit comoving volume per unit proper time provided by the PDB (pair) model in the first row of Table II in Abbott et al. (2023) to standardize our merger rates.PDB (pair) model uses the mass and spin distribution that is the closest match to GWTC-3 distribution.
We normalize the initial distribution of the GWTC-3 with the total rate density of mergers, integrated across all masses and spins, taken to be fixed at 240 +270 −140 Gpc −3 yr −1 (which is in the first row and last column on this table).For the LRR distributions, we also used PDB (pair) model, by taking the numbers from the first line and the columns (1), (2), and (3) of the same Table II Abbott et al. (2023).We reproduce these astrophysical density rates in Table 2.

Results
We make the results of the observing scenarios, including the sky-localization FITS files, available on for PBD/GWTC-3 (GWTC-3) (Coughlin et al. 2022) and LRR (Kiendrebeogo et al. 2023) in separate repositories.These simulations allow us to estimate the rates of detection by the IGWN over O4 and O5; in the following, we will focus on O4 as an example for follow-up simulations by ZTF and Rubin Observatory (although, in practice, Rubin Observatory is not expected to contribute significantly until O5).However, we reproduce some of the analyses for O5, and report them in Appendix A.

O5
In Figure 6, we summarize the detection results of the simulation set.We provide predictions for the annual detection rates of CBCs in O4 and O5 for both the LRR and GWTC-3 distributions in Table 3.The confidence interval combines both the log-normal distribution of the merger rate and uncertainties from the Poisson counting statistics.The low number of NSBHs predicted by the PDB model is due to the existence of a nearly empty mass gap in that model, combined with a pairing function (Fishbach & Holz 2020) that favors equal-mass binaries.NSBHs must straddle the mass gap, with one component on each side.This leads to asymmetric mass ratios, which are in turn disfavored by the model fit, as most binaries in the population are consistent with being equal mass.A version of the PDB model with a partially filled mass gap would predict more NSBH events relative to the other types of CBCs.
In Table 4, we also provide statistics regarding the GW signal sky-localization area, luminosity distance, and comoving volume.Sky-localization area (volume) is given as the 90% credible region, defined as the smallest  area (volume) enclosing 90% of the total posterior probability.This corresponds to the area (volume) of the sky that must be covered to have a 90% chance of including the source.We have adopted the same statistical treat-  ment process as the one used in Petrov et al. (2022).
The results from the simulation of the GWTC-3 distribution are also available in the IGWN Public Alerts User Guide (see footnote 1).
There are notable differences that arise due to the improved mass distributions measured in GWTC-3, which were derived from the maximum a posteriori fit to all compact binaries detected so far; this is vastly different from the normal distribution (NS) and power law (BH) for masses assumed by LRR.This difference, coupled with the inclusion of single-detector triggers, yields large differences from previous reports.This GWTC-3 distribution accounts for an increase in the predicted number of detected events by a factor ∼ (0.83%/0.22%) = 3.772 for O4, then ∼ (1.22%/0.48%)= 2.542 for O5.Breaking down by population type, the estimated annual detec-tion rate is a factor of ∼2 (5) times higher in BNS (BBH) subpopulations for GWTC-3 but drops to ∼ 0.6 for NSBH.The median luminosity distance predicted by GWTC-3 is about ∼ 1.14 (1.19), 1.36 (1.31), and 2.44 (2.36) larger than the value for LRR for BNS, NSBH, and BBH events respectively during O4 (O5), with the median of sensitive volume similarly increasing by ∼ 1.7 (1.6), 1.9 (2), and 8.14 (6.65).For sky-localizations, the results are broadly similar to previous results.For example, during O4, GWTC-3 predicts ∼ 11% smaller sky-localization area for BNS subpopulation than that from LRR, while in the BBH case GWTC-3 predicts about ∼ 52% larger than the value for LRR.

PREDICTIONS FOR DETECTION RATES AND SCIENCE WITH GRAVITATIONAL-WAVE COUNTERPARTS
In this section, we use the results from Sec. 2 to make predictions for potential counterpart detections during O4 and O5.

Light Curves
To simulate the light curves corresponding to the BNSs and NSBHs from Sec. 2, we first predict the dynamical and disk wind ejecta produced for each simulated object.For the BNS case, we use the fit from Coughlin et al. (2018a) for the dynamical and Dietrich et al. (2020a) for the disk wind.For the NSBH case, we use the fit from Foucart et al. (2018) for the dynamical and Krüger & Foucart (2020) for the disk wind.We note that the NSBH case often results in no KNe due to large mass ratios resulting in a direct plunge of the neutron star.To translate these ejecta quantities into light curves, we use the multidimensional Monte Carlo radiative transfer code POSSIS (Bulla 2019(Bulla , 2023)).For the BNS case, we use the same geometry and lanthanide fractions for each component as presented in Dietrich et al. (2020a), shown to yield good fits to GW170817.For the NSBH case, we use the same configurations as presented in Anand et al. (2020).This model takes as input the dynamical ejecta M dyn ej , disk wind ejecta M wind ej , half opening angle Φ, and the observation angle Θ obs (see Ref. Dietrich et al. (2020a) for details).We take the inclination measurement from the simulated value in the observing scenarios, and we draw Φ uniformly between 15 • and 75 • .To interpolate the grid-based model, we use a Gaussian process regression method (Coughlin et al. 2018a(Coughlin et al. , 2019a)). Figure 7 shows a histogram of the light curves in the optical to near-infrared bands for the O4 simulation set, simulated with NMMA, publicly available on GitHub11 .

Simulated follow-up
We use these light curves to perform a simulation campaign for optical observations of GW counterparts with Rubin Observatory and ZTF.To simulate ZTF and Rubin Observatory observations, we use the Gravitational-Wave ElectroMagnetic OPTimization12 (gwemopt; Coughlin et al. 2018bCoughlin et al. , 2019d)).We take ZTF to have a base sensitivity in g, r, and i bands with g ∼ 21.7 mag, r ∼ 21.4 mag, i ∼ 20.9 mag (Almualla et al. 2021), and Rubin Observatory with u, g, r, i, z, and y bands with u = 23.9mag, g = 25.0 mag, r = 24.7 mag, i = 24.0mag, z = 23.3mag, and y = 22.1 mag.Then, in order to obtain the full range of bands used in Figure 7, we incorporate Gemini Observatory, which uses sensitivity tuned u, g, r, i, z, y, J, H, and K bands with u = 24.1 mag, g = 25.0 mag, r = 25.0 mag, i = 25.3 mag, z = 24.5 mag, y = 23.0 mag, J = 23.2mag, H = 22.6 mag, and K = 22.6 mag. Figure 8 shows a one dimensional histogram of prediction of the magnitude peak in all ZTF and Rubin observatory bands during the forthcoming observing runs O4 and O5.gwemopt is the package that was used to create follow-up schedules for a number of facilities during O3, including ZTF and the Dark Energy Camera, and therefore serves as reasonably realistic software to use for this purpose for O4 and O5.We inject simulated KNe consistent with the GW localization and simulate follow-up observations, taking into account the sensitivity and FOV of the telescopes; this yields the expected fraction of detections for KNe within the simulation set.Here, a detection means at least one photometry in the light curve.
In our search algorithms, we sample the light curves according to ZTF Phase-II public cadence and private i band survey cadence by drawing revisit cadences from kernel density estimate (KDE) fits to the same.For every visit, we assign a random night's observing filter sequence.We also add 300 s ToO observations in g and r bands during the first day or first two days after the trigger, for GW localizations of greater or lesser than 1000 deg 2 , respectively, by following the observations taken by ZTF during O3 (Kasliwal et al. 2020).These observations are executed as part of the GW EM search program within the ZTF collaboration (Anand et al. 2021).Given mean magnitudes, we simulate magnitude uncertainties using a skew normal fit to forced photome-  measurements that are fainter.We require that the light curves meet the trigger criterion of S/N > 3 (Andreoni et al. 2021).For Rubin observatory, we used ToO observations based on the strategy presented in Andreoni et al. (2022).
For this analysis, we inject the 1004 ( 2003) BNS and 184 (356) NSBH of the PBD/GWTC-3 (GWTC-3) distribution that have passed the S/N cutoff in O4 (O5). Figure 9 shows the distribution of the missed and found events for both the BNS and NSBH cases for observing run O4.We compare these to the peak absolute magnitude in r band and the 2D probability enclosed by the simulated observations for ZTF and Rubin Observatory.There are many more BNS than NSBH detections due to the much smaller subset of NSBH injections with nonzero ejecta masses.For comparison, we show marginalized 1D histograms for both the detected and missed sets of objects, which show distinct differences in both.The detected set shows a distinct preference toward brighter objects.It also shows a preference for both higher 2D and 3D probability coverage.Despite ZTF's wide FOV of 47 deg 2 , its relatively lower sensitivity limits it to ∼ 12 (4) and 4 (1) detections of the injections respectively for BNS and NSBH events during the observing run O4 (O5).Rubin Observatory instead finds ∼ 55 (60) and 1 (5) of the injections for BNS and NSBH events, respectively during the next observing run O4 (O5).By combining the EM detection fraction and the GWTC-3 CBCs annual detection rates (Tab.3 ), in Table.5, we provide the predictions we expect for the detection rates of EM counterparts during the forthcoming run under the specific assumptions described above.For the objects detected by this process, we perform parameter estimation of the resulting EM light curves.

KNe sample constraints
To do so, we use the NMMA framework (Pang et al. 2022), designed to perform Bayesian inference of multimessenger signals, incorporating all available data on the neutron star EoS in the process (Dietrich et al. 2020a;Pang et al. 2021;Tews et al. 2021).This is an efficient package to evaluate Bayes' theorem in order to obtain posterior probability distributions, p( ⃗ θ|x, M), for model source parameters ⃗ θ under the hypothesis or model M with mock-up data x as where p(x| ⃗ θ, M), p( ⃗ θ|M), and p(x|M) are the likelihood, prior, and evidence, respectively.This framework has been used in the measurement of the NS EoS and H 0 using GW170817 (Dietrich et al. 2020a) and the detection of the shortest long γ-ray burst ever confirmed (Ahumada et al. 2021).

Ejecta constraints
We begin by evaluating 90% upper limits possible from the sample considered here (KNe for these objects), then constrain them for the ejecta model parameters M dyn ej and M wind ej , as well as a systematic contribution to the dynamical ejecta α and the fraction of the disk mass contributing ζ.This enables us to make an empirical constraint of the fraction of the disk contributing to KNe.We constrain M dyn ej to 10-40% and M wind ej to 10-30%.
We also desire to differentiate between the prompt collapse and formation of hypermassive and/or supramassive neutron star.In the prompt-collapse scenario, where we assume that no disk wind is produced, we show in Figure 10 the 90% upper limit on the disk wind contribution for our sample.We derive limits ranging from log 10 M wind ej = [−2.8− 2.0], resulting in strong constraints on the presence of such a disk.

ZTF proposal for GW follow-up and triggering criteria
The properties of a GW event released in low-latency via the General Coordinates Network (GCN) notice can prove to be useful in determining whether or not to trigger telescope time on a given GW event.Information such as the alse alarm rate (FAR), probability that the source is of astrophysical origin (p-astro), probability of the GW merger containing a neutron star (HasNS), and probability of the GW merger leaving behind a remnant (HasRemnant), and the Bayes factor of coherence between multiple detectors (log(BCI)) can provide indirect clues about whether a GW event is astrophysical, significant, and could harbor an EM counterpart.Thus, in Table 6, we define triggering criteria for ZTF based on these low-latency properties.An event satisfying all of the "Go-deep" requirements merits triggering ToO observations, provided the localization and distance are accessible for ZTF.A "Go-wide" event prompts reweighting public ZTF survey fields to observe the localization in the nominal 30 s exposures.Events for which any properties fall within the "Deliberate" or the "No Go" categories merit human interaction and discussion to decide whether to trigger.Since the threshold for LIGO-Virgo-KAGRA public alerts release has been lowered to 2 day −1 ( user guide13 , Abbott et al. 2020b) in O4, having such triggering criteria will be paramount for wisely allocating existing telescope resources.
In Figure 11, we show the expected annual number of triggers within 400 Mpc during the O4 and O5 run, based on the predictions in Section 2.4.Assuming a AT2017gfo-like KNe with M abs ∼ −16 mag at peak, and taking to account ZTF's limiting magnitude in 300 s exposures (m AB ≈ 22 mag), we estimate that ZTF could detect KNe falling within 400 Mpc.Due to the overall low numbers of NSBH mergers, we only expect 0-2 NSBH mergers within 400 Mpc during the O4 and O5 runs, which is consistent with the KNe detection prospects discussed in Sec. 3.However, our calculations yield < N > = 13 (< N > = 18) BNS mergers within 400 Mpc during O4 (O5), providing ample opportunities to conduct sensitive searches for counterparts to BNS mergers.
Finally, we assess the distribution of GW events as a function of sky area.In Figure 12, we show the fraction of O4/O5 triggers whose 90% confidence region. of the GW localization falls within a given sky area threshold.Assuming a typical 8 hr night and ZTF's footprint of ∼50 deg 2 , and a three-filter tiling strategy (i.e., g−r−g) we find that, with ZTF Partnership time alone (comprising 30% of the night), we can fully tile the localization for ∼30% of GW alerts.With the addition of the private Caltech allocation (comprising ∼50% of the night), we can probe the localization for nearly 40%, and by using the public survey allocation as well (100% of the night), we can fully tile the localization for 50% of all events.These figures of merit can be easily computed from the data set presented in this paper to estimate the number of ToO triggers and time request needed for a successful GW follow-up campaign with wide-field telescopes.In particular, these calculations will prove especially useful for informing Rubin's triggering strategy once it comes online.

H 0 constraints
Given the high computational costs for GW analysis runs, we reduce the sample size of events for O4 and O5, and study only the loudest 30 sources in terms of network S/N.This selection of the 30 events was made on the CBCs of PBD/GWTC-3 (GWTC-3) distribution that passed the S/N threshold.Therefore, we note that this number is not a function of the expected EM detection rate during the next campaigns O4 and O5.For the GW simulations, we employ the IMRPhenomD NRTidalv2 model of Dietrich et al. (2019) and the EOS with maximum likelihood inferred in Huth et al. (2022).We inject signals in Gaussian noise, using the detectors' PSD predicted for O4 and O5, as explained above.We perform parameter estimation analyses, where the GW likelihood landscape is explored employing the BILBY library (Ashton et al. 2019;Romero-Shaw et al. 2020) and the included nested sampling (Skilling 2006;Veitch & Vecchio 2010) algorithm DYNESTY (Speagle 2020;Koposov et al. 2022).To reduce the computational cost of the analysis for injected GW signals, we use the relative binning method (Zackay et al. 2018;Leslie et al. 2021), following the implementation in Janquart (2022).We analyze signals starting from 20 Hz, and we employ the usual uniform in comoving volume prior to the luminosity distance used in GW parameter estimation.The same prior setting is used for KNe inferences detailed in the following.With regard to potentially associated KNe signals, we use the NMMA framework to infer posterior distributions on KNe properties as well as the luminosity distance using the KNe model of Bulla (2019Bulla ( , 2023)).Moreover, we perform two sets of EM parameter estimates, for ZTF and Rubin Observatory observations.In order to connect the binary masses in the remaining sample of 30 BNS systems to ejecta material masses powering the KNe, we use the phenomenological relations established in Dietrich et al. (2020).In this process, there were nine (seven) highmass BNS systems in the O4 (O5) sample, which would directly form a BH leaving no ejected material powering an EM counterpart.
For the remaining 21 ( 23) O4 (O5) BNS mergers with an EM counterpart, we use the inferred luminosity distance posterior distributions of both GW simulations and EM simulations to determine the posterior distribution for H 0 .Moreover, we include inferred luminosity distance posteriors of the GW signal GW170817 from Abbott et al. (2019b), Abbott et al. (2019c) and the KN signal AT2017gfo, which we inferred with the same KN model yielding a total number of 22 (24) O4 (O5) BNS coalescences.In this context, we assume the standard cosmology model.Since all of our BNS systems lie within a range of 300 Mpc or within a small redshift regime, we assume the appropriate linear Hubble relation for nearby events in which d L and z are the luminosity distance and the redshift to the source, respectively, and c is the speed of light, and H 0 is the Hubble constant.Since a volumetric prior on the luminosity distance inherits a 1/H 4 0 prior factor, we correct for this in our study similar to Abbott et al. (2017).Concerning the underlying distance and/or redshift distribution of our injected O4 (O5) samples, we use the injected distance and the injected Hubble constant to calculate the corresponding injected redshift, which is the mean of the Gaussian distribution with a 1% relative uncertainty set as standard deviation.In this study, we do not break the distanceinclination angle degeneracy by including additional information on the binary's orbital inclination from other > 0.9 > 0.9 0.1-0.9< 0.1 HasNS > 0.9 > 0.9 0.1-0.9< 0.1 potential observations, such as GRBs.Nevertheless, we are able to improve the Hubble constant measurements through the combined distance measurements from GW and KNe observations; see Figure (2a) for a similar approach taken in Dietrich et al. (2020).We point out that, for a more detailed study with a larger number of considered detections, further selection effects would need to be considered for our population analysis.In particular, we would need to correct for our choice that we consider only the loudest 30 GW signals for which we simulated potential KNe observations.Apart from a selection effect that would result from the redshift measurement of a potential host galaxy, we would also need to consider the dependency of H 0 on the BNS component masses as high-mass BNS systems will quickly form as BHs leaving no traceable EM counter-  part and, consequently, were not considered for our H 0 projections.We will leave a detailed analysis of selection biases for future work at the present moment.However, although additional selection bias corrections might be needed, our analysis does not show biases in our recovery, which overall shows the robustness of our study, but we expect this situation to change when the population sample sizes increase or observations at higher redshifts are included.
In Figure 13 (top panels), we show our H 0 -results for the ZTF scenario obtained from GW and EM parameter estimations on the luminosity distance and, in addition, the results when combining GW+EM information and contrast these to state-of-the-art measurements.The uncertainties of our results are reported at 90% credible interval.We highlight our estimated annual BNS detection rate from Table 5 as gray dashed line and mark the corresponding 90% credible interval as gray regions.With only one joint GW+EM observation, as estimated for ZTF-O4, we are not able to provide strong constraints on the Hubble constant, i.e., while we recover the injection value of H 0 = 67.74km s −1 Mpc −1 ; this is mainly caused by the large uncertainty of our measurement H 0 = 60.68 +9.24  −7.47 km s −1 Mpc −1 .Most notably, this shows how unlikely it is that within O4 we will be able to break the Hubble tension.Similarly, one GW+EM observation in O5 provides an estimate of H 0 = 61.26+17.73  −18.97 km s −1 Mpc −1 .To show what might be possible with several more GW events that have a corresponding KNe and to understand whether there is a systematic bias being introduced with our methodology, we combine 22 BNS events in O4 and estimate H 0 = 66.37 +0.58  −0.95 km s −1 Mpc −1 .We do a similar study for O5, combining 23 events, and estimate H 0 = 66.74 +0.39  −0.33 km s −1 Mpc −1 .Both estimates recover the injected value of H 0 .
In Figure 13 (bottom panels), we show our H 0results for the Rubin Observatory scenario.We find that roughly two joint observations as estimated for O4 with H 0 = 62.56 +5.27  −4.70 km s −1 Mpc −1 and approximately 5 combined observations in O5 with H 0 = 65.30+2.31  −2.99 km s −1 Mpc −1 can recover the injection.For a total number of 22 combined events in O4, we estimate H 0 = 67.01+0.43  −0.53 km s −1 Mpc −1 which is close to the injection value.For O5, we recover H 0 = 66.23 +0.39  −0.33 km s −1 Mpc −1 .Overall, our study shows that subpercent level measurements will be possible through a combination of GW+EM (KNe) information with a sufficient number of KNe detections.Furthermore, we emphasize that, while we include only the distance measurement from the KNe, further information, e.g., through the GRB afterglow might further help to reduce uncertainties by breaking degeneracies between the distance and the inclination.

CONCLUSION
In this paper, we have performed an end-to-end GW survey simulation for O4 and O5 follow-up of BNS and NSBH mergers.We simulate event follow-up for the ZTF and Rubin Observatory, showing the impact on their detection capability.Based on the GW and EM posteriors from these analyses, we simulate the potential H 0 using a distance threshold of 300 Mpc.With the simulation study performed above, we can summarize our conclusions and implications for future GW followups as follows: (i) During O4 and O5, there will be many more events than are reasonable to follow up.Unsurprisingly, the objects that are best to follow-up are the best localized and near-by Chen & Holz (2016).Proximity is an especially important consideration for Figure 13.H0-estimates for the ZTF (top row) and Rubin Observatory (bottom row) observation scenarios for O4 (left column) and O5 BNS samples (right column).We show H0-estimates from our GW simulations (in orange), EM simulations (in violet), and for the combined GW+EM results (in blue) in the top panels, whereas relative errors are shown in the bottom panels.We indicate the expected O4 and O5 detection rates in alignment with Table 5 as gray dashed line and show the 90% credible interval as gray shaded regions.In the bottom panel, we contrast our results to the Planck measurement of the cosmic microwave background (Ade et al. 2016;Planck, violet), to the Hubble measurement via type-Ia supernovae (Riess et al. 2016;SHOES, light blue), and to the H0-measurement of superluminal motion of the jet in GW170817 (Hotokezaka et al. 2019;superluminal, gray).All uncertainties are reported at 90% credible interval.
the NSBH mergers, given their large expected distances and faint intrinsic luminosities.
(ii) During the next O4 and O5 runs, in contrast to Rubin Observatory, ZTF, due to shallower limits in its bands, will have difficulty detecting the EM counterparts, especially at larger distances (see Figure 9).
(iii) The GW contribution dominates under the assumptions in this document, but again, the high S/N, nearby events dominate the sensitivity improvement.A few well-localized nearby events are worth much more than many far-away, poorly localized events.For H 0 , most events contribute equally, especially as relative uncertainty due to virial velocities decreases as distance increases.Here, Figure 16 shows a histogram of the light curves in the optical to near-infrared bands for the O5 simulation set, simulated with NMMA and publicly available on https://github.com/nuclear-multimessenger-astronomy/nmma.

B.3. Distribution of missed and found BNS and NSBH injections by ZTF and Rubin Observatory
In Figure 18, we show the distribution of the missed and found Kilonovae (KNe) events for both the BNS and NSBH cases during observing run O5 with ZTF and Rubin Observatory observations, using Gravitational-Wave ElectroMagnetic OPTimization (https://github.com/mcoughlin/gwemopt).

Figure 1 .
Figure 1.Flowchart of observing scenarios process.Here, µ = 1.33M⊙, and σ = 0.09M⊙, then in LRR case, m represents the primary mass m1 or secondary mass m2 , since they are drawn in the same way with m2 ≤ m1.

Figure 2 .
Figure 2. 90% posterior credible intervals for the component masses for all CBCs in the GWTC-3 catalog (Abbott et al. 2023) study assuming uniform priors in detector-frame masses and fixed FAR about 0.25 year −1 (Abbott et al. 2023).Events classified by the LIGO-Virgo-KAGRA collaboration as BNSs, NSBHs, and BBHs are shown in dark violet, olive (dark yellowish-green), and orange, respectively.The ambiguously classified event GW190814 is shown in dark purple.The gray band indicates the approximate location of the purported lower-mass gap.GW190814 is the only event within this region at more than 90% credibility

Figure 3 .
Figure 3.The 1D PDB mass distribution, p(m|λ) on the interval [1, 100] M⊙ for a specific choice of hyperparameters λ.See Tab. 7 in Appendix A.1, for the other parameters of the mass distribution.

Figure 4 .
Figure 4. Cutoffs between subpopulations of compact binary coalescences for both the GWTC-3 and LRR samples.It should be noted that the BHs in the LRR field are limited to be below 50 M⊙, but are allowed to be as massive as 100 M⊙ in the GWTC-3 distribution.

Figure 5 .
Figure 5. Simulated mass distributions for O4.The left panel shows draws from the LRR distribution, and the right panel shows draws from the PBD/GWTC-3 (GWTC-3) distribution.The upper panels are the 2D mass distributions of the components of each CBC in the context of the detector, and the lower panels are the 2D primary mass and distance distributions.All axes are shown on a logarithmic scale.The color base shows the number of CBC events per pixel.For O5 results, see in Appendix A.1.1, Figure 14. are publicly accessible on GitHub 7 for LRR distribution, and on GitHub 8 / (Singer et al. 2022) for GWTC-3.Following Petrov et al. (2022), we apply an S/N threshold of 8 for the entire GWTC-3 distribution, and the BNSs and NSBHs populations of the LRR distribution.This S/N threshold is set to 9 for the BBH population of the LRR distribution, consistent with the localization

Figure 6 .
Figure6.This figure shows the simulated detections for O4.The upper panel shows the number of injections (CBCs from each population drawn from both distributions) in colored dots and the bar chart representing the number of events passing the S/N cut.The bottom panel shows the percentage of detection relative to the number injections for the two distributions.The colored dots represent the percentages of each population that passed the S/N, while the blue and green lines with respectively 0.83% and 0.22% are successively percentages of detection of all the events (BNS + NSBH + BBH) of the PBD/GWTC-3 (GWTC-3) and LRR distributions injected in our simulation.For O5, see Appendix A.1.2,Figure15

Figure 7
Figure 7. 2-D histograms of simulated BNS (left) and NSBH (right) light curves for O4.We show light curves in u-, g-, r-, i-, z-, y-, J-, H-bands in order to include all the bands used by surveys considered in this work (ZTF, Rubin and Gemini).In each panel, we plot three dash lines; on top the 10th percentile, at middle the 50th percentile and on bottom the 90th percentile.The color bar shows the number of detections in the different bands.For the observing run O5, see Figure.16, in Appendix B.1.

Figure 8 .Figure 9 .
Figure 8. 1D histograms of the peak magnitudes in the ZTF bands (in green) and the Rubin observatory bands (in red).On the left, we show the BNS peak mag in O4; on the right, the same plot for the run O5.See Figure 17, in Appendix B.2 for the peaks related to each band of each telescope.tryuncertainties.Finally, the forced photometry upper limits are estimated using KDE fits and used to reject

Figure 10 .
Figure 10.90% upper limit on the M wind ej measurements for the prompt collapse simulation set, where the M wind ej contribution is set to 0.

Figure 11 .
Figure11.Cumulative histograms of 100,000 realizations of the number of BNS and NSBH mergers predicted to be detected during LIGO-Virgo-KAGRA O4 and O5 within 400 Mpc based on the observing scenarios predictions in this paper.The mean number of expected detections is quoted for each merger type.

Figure 12 .
Figure 12.Cumulative histogram of simulated GW skymaps for O4 that satisfy our triggering criteria as a function of the 90% credible region.With both Caltech and Partnership time spanning ∼50% of the night (corresponding to a maximum area of 800 deg 2 ), we could fully probe the localization for nearly 40% of events.

Figure 14
Figure14describes the simulated O5 mass distributions for each model that meets the signal-to-noise (S/N) threshold outlined in Section 2.3.

Figure 14 .
Figure14.Here we show all the O5 CBCs from all populations.On the left the distribution of LRR, on the right that of PBD/GWTC-3 (GWTC-3).The high panels are the mass distributions of the components of each CBC in the context of the detector and the low panels are the primary mass distributions as a function of the collision distance.Mass and distance distributions are shown on a logarithmic scale.The color base shows the number of CBC events per pixel.

Figure 15
Figure 15 we summarize the O5 detection results of the simulation set.
B. PREDICTIONS FOR DETECTION RATES AND SCIENCE WITH GRAVITATIONAL-WAVE COUNTERPARTS B.1.2D histogram of the O5 BNS and NSBH light curves simulated.

Figure 17
Figure17we present a one-dimensional histogram of predictions for the peaks related to each band of ZTF and the Rubin Observatory during the forthcoming observing runs O4 and O5.

Figure 15 .
Figure15.This figure shows the simulated detections for O5.The upper panel shows the number of injections (CBCs from each population drawn from both distributions) in colored dots and the bar chart represents the number of events passing the cutoff point.The bottom panel shows the percentage of detections relative to the number injections for the two distributions.The colored dots represent the percentages of each population that passed the S/N, while the blue and green lines with respectively 1.22% and 0.48% are successively percentages of detection of all the events (BNS + NSBH + BBH) of the PBD/GWTC-3 (GWTC-3) and LRR distributions injected inn our simulation.

Figure 17 .
Figure 17.Those plots are the 1D histograms of the peak magnitudes in the ZTF and the Rubin observatory bands.On the left we show ZTF BNS peak mag in the run O4 (on top) and O5( at bottom), then on the right the same plot for the case of Rubin Observatory.The black line shows the peak mag in all the bands of each telescope detected during the run O4 and O5 simulation.)

Table 1 .
The Cosmology-corrected Inspiral Sensitive Distance (in Mpc) from a GW Strain PSD

Table 3 .
Annual Detection Rates for Compact Binary Coalescences That we expect for the Runs O4 and O5.These uncertainties do not incorporate the Monte Carlo method, but only combine both the log-normal distribution of the merger rate and the Poisson counting statistics

Table 4 .
Summary Statistics for O4 and O5.These recorded values are given as 90% credible interval calculated with the 5% and 95% quantile.Those uncertainties have been described by Monte Carlo sampling.

Table 5 .
Annual detection rate of EM counterparts that we expect for ZTF and Rubin Observatory during the Run O4 and O5.

Table 6 .
Triggering Based on GW Candidate Event Properties

Table 7 .
Summary of Population Model parameters.The first several entries describe the rate and mass distribution parameters, and the last two entries describe the spin distribution parameters.Parameter controlling how the rate tapers at the low end of the mass gap.50 η high Parameter controlling how the rate tapers at the high end of the mass gap.50 ηmin Parameter controlling tapering of the power law at low mass.50 ηmax Parameter controlling tapering of the power law at high mass.4.91 β Spectral index for the power-law-in-mass-ratio pairing function.1.89