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

R. Weizmann Kiendrebeogo, Amanda M Farah, Emily M Foley, Abigail Gray, Nina Kunert, Anna Puecher, Andrew Toivonen, R. Oliver Vandenberg, Shreya Anand, Tomás Ahumada, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

R. Weizmann Kiendrebeogo, Amanda M Farah, Emily M Foley, Abigail Gray, Nina Kunert, et al. Up- dated observing scenarios and multi-messenger implications for the International Gravitational-wave Network's O4 and O5. The Astrophysical Journal, 2023, 958 (2), pp.158. 10.3847/1538-4357/acfcb1 . hal-04152506

HAL Id: hal-04152506

https://hal.science/hal-04152506v1

Submitted on 15 Apr 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Updated Observing Scenarios and Multimessenger Implications for the International Gravitational-wave Networks O4 and O5

R. Weizmann Kiendrebeogo<sup>1,2,3</sup>, Amanda M. Farah<sup>4</sup>, Emily M. Foley<sup>3,5</sup>, Abigail Gray<sup>3</sup>, Nina Kunert<sup>6</sup>,
Anna Puecher<sup>7,8</sup>, Andrew Toivonen<sup>3</sup>, R. Oliver VandenBerg<sup>3</sup>, Shreya Anand<sup>9</sup>, Tomás Ahumada<sup>9</sup>,
Viraj Karambelkar<sup>9</sup>, Michael W. Coughlin<sup>3</sup>, Tim Dietrich<sup>6,10</sup>, S. Zacharie Kam<sup>1</sup>, Peter T. H. Pang<sup>7,8</sup>,
Leo P. Singer<sup>11</sup>, and Niharika Sravan<sup>12</sup>

*Received 2023 July 13; revised 2023 September 13; accepted 2023 September 24; published 2023 November 21*

## Abstract

An advanced LIGO and Virgo's third observing run brought another binary neutron star merger (BNS) and the first neutron-star black hole mergers. While no confirmed kilonovae were identified 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 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 scientifically 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.26}^{+0.58}$ ($1.97_{-1.2}^{+2.68}$) for the Zwicky Transient Facility (Rubin Observatory).

*Unified Astronomy Thesaurus concepts: Gravitational wave sources (677); Gravitational wave astronomy (675); Gravitational wave detectors (676); High energy astrophysics (739); Hubble constant (758); Astrophysical black holes (98); Neutron stars (1108); Gravitational waves (678); Cosmology (343)*

## 1. Introduction

After the detection of AT2017gfo (Abbott et al. 2017a; Coulter et al. 2017; Smartt et al. 2017), associated with the binary neutron star (BNS) merger GW170817 (Abbott et al. 2017b) and the short gamma-ray burst GRB170817A (Goldstein et al. 2017; Savchenko et al. 2017; Abbott et al. 2017c), there have been significant electromagnetic (EM) follow-up efforts of both further BNS detections such as GW190425 (Coughlin et al. 2019b; Hosseinzadeh et al. 2019; Kuin & Swift Team 2019; Song et al. 2019; Abbott et al. 2020a; Antier et al. 2020a; Gompertz et al. 2020; Saleem et al. 2020) and the neutron-star black hole (NSBH) coalescences GW200105 and GW200115 (Antier et al. 2020b; Abbott et al. 2021; Anand et al. 2021; Chattopadhyay et al. 2022; Wang & Zhao 2022).

Many science cases motivate the follow-up of neutron star mergers; these science cases include constraints on the neutron star equation of state (Bauswein et al. 2017; Margalit & Metzger 2017; Coughlin et al. 2018a, 2019a, 2020c;

Annala et al. 2018; Most et al. 2018; Radice et al. 2018; Lai et al. 2019; Dietrich et al. 2020; Huth et al. 2022), the Hubble constant (Abbott et al. 2017d; Hotokezaka et al. 2019; Coughlin et al. 2020a, 2020b; Dietrich et al. 2020), and r-process nucleosynthesis (Chornock et al. 2017; Coulter et al. 2017; Cowperthwaite et al. 2017; Pian et al. 2017; Rosswog et al. 2017; Smartt et al. 2017; Kasliwal et al. 2019; Watson et al. 2019). Even in the absence of further counterpart detections (e.g., Andreoni et al. 2019, 2020; Coughlin et al. 2019b; Goldstein et al. 2019; Gomez et al. 2019; Lundquist et al. 2019; Ackley et al. 2020; Antier et al. 2020b; Gompertz et al. 2020; Kasliwal et al. 2020; Anand et al. 2021), these "upper limit" observations place constraints on possible kilonovae (KNe) counterparts and their potential progenitor parameters (e.g., Coughlin et al. 2020c; Kasliwal et al. 2020).

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.,

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1"><sup>1</sup> Laboratoire de Physique et de Chimie de l'Environnement, Université Joseph KI-ZERBO, Ouagadougou, Burkina Faso; weizmann.kiendrebeogo@oca.eu</footnote>

<footnote n="2"><sup>2</sup> Artemis, Observatoire de la Côte d'Azur, Université Côte d'Azur, Boulevard de l'Observatoire, F-06304 Nice, France</footnote>

<footnote n="3"><sup>3</sup> School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA</footnote>

<footnote n="4"><sup>4</sup> Department of Physics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="5"><sup>5</sup> B.S. Physics, B.A. Mathematics, Wake Forest University, USA</footnote>

<footnote n="6"><sup>6</sup> Institute of Physics and Astronomy, Theoretical Astrophysics, University Potsdam, Haus 28, Karl-Liebknecht-Str. 24/25, D-14476, Potsdam, Germany</footnote>

<footnote n="7"><sup>7</sup> Nikhef—National Institute for Subatomic Physics, Science Park 105, 1098 XG Amsterdam, The Netherlands</footnote>

<footnote n="8"><sup>8</sup> Institute for Gravitational and Subatomic Physics (GRASP), Utrecht University, Princetonplein 1, 3584 CC Utrecht, The Netherlands</footnote>

<footnote n="9"><sup>9</sup> Cahill Center for Astrophysics, California Institute of Technology, Pasadena CA 91125, USA</footnote>

<footnote n="10"><sup>10</sup> Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, Potsdam D-14476, Germany</footnote>

<footnote n="11"><sup>11</sup> Astroparticle Physics Laboratory, NASA Goddard Space Flight Center, Mail Code 661, Greenbelt, MD 20771, USA</footnote>

<footnote n="12"><sup>12</sup> Department of Physics, Drexel University, Philadelphia, PA 19104, USA</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Flowchart of observing scenarios process. Here, $\mu = 1.33 M_\odot$, and $\sigma = 0.09 M_\odot$; then, in LRR case, $m$ represents the primary mass $m_1$ or secondary mass $m_2$, since they are drawn in the same way with $m_2 \le m_1$.</figcaption></figure>

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 Petrov et al. (2022), 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,<fnref n="13" /> as well as simulate potential science constraints based on self-consistent counterpart search simulations. In Section 2, we summarize the simulations of the observation scenarios expected during the next observation campaigns. In Section 3, we report the results of our simulations on the tracking of EM 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² 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² FOV) ground-based telescope designed to conduct deep 10 yr survey of the Southern sky. In Section 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 Section 4.

## 2. Observing Scenarios

### 2.1. 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; KAGRA; Akutsu et al. 2021), observing runs (O4 and O5, respectively) and present the multimessenger 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. 2021; Kasliwal et al. 2020), and Rubin 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.

<footnote n="13">¹³ https://emfollow/docs.ligo.org/userguide/capabilities.html</footnote>
## 2.2. Population Models

In this section, we outline the processes used to generate two distributions of CBCs, which we will compare in this work. The first distribution, LRR (Living Reviews in Relativity), is drawn from the population model outlined in Abbott et al. (2020b), Petrov et al. (2022) and used for previous IGWN observing scenarios. This distribution consists of the normal distribution (NS) for neutron star masses and the power law for black hole masses. The second distribution, PBD/GWTC-3 (GWTC-3), is derived from the population model described in Farah et al. (2022), and Abbott et al. (2023) is used for the current observing scenarios. Specifically, we fit the POWER LAW + DIP + BREAK (PDB) model used in Farah et al. (2022) to all CBCs in GWTC-3 from Abbott et al. (2023) and then use the maximum a posteriori value of the resulting fit.

Each distribution consists of three astrophysical subpopulations of CBCs: binary neutron star merger (BNS), neutron-star black hole merger (NSBH), and binary black hole merger (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 \le m_1$ by definition. The mass of a nonrotating neutron star cannot exceed the Tolman–Oppenheimer–Volkoff (TOV) limit $M_{max,TOV} \approx 2-2.5 M_\odot$ (Margalit & Metzger 2017; Rezzolla et al. 2018; Ruiz et al. 2018; Shibata et al. 2019; Dietrich et al. 2020; Legred et al. 2021). However, rotating NSs can exceed this limit (Baumgarte et al. 2000; Stergioulas 2003). Additionally, 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_\odot$ (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_\odot$. 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_\odot$ also maintains a consistency with previous analyses (Abbott et al. 2020b).

We follow Farah et al. (2022), who proposed a resolution to this NS–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).

LRR distribution. Here, as in Petrov et al. (2022), 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_\odot$ (Özel & Freire 2016) from the interval $[1, 3] M_\odot$ with a standard deviation of $0.09 M_\odot$ ($p(m) \sim N(\mu, \sigma^2)$) and spin uniformly distributed magnitudes in the interval $[0, 0.05]$. For the BBHs, sampling is performed from $[3, 50] M_\odot$ (Abbott et al. 2019a) using a truncated power-law distribution $p(m) \propto m^a$ with $a = -2.3$ (Salpeter 1955). Both masses are independently drawn from this distribution

and paired according to the rule $m_2 \le 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 model from Farah et al. (2022), Abbott et al. (2023). This model consists of a broken power law with a notch filter $n(m|M_{low}^{gap}, M_{high}^{gap}, A)$ that suppresses the merger rate between NSs and BHs ($M_{low}^{gap}$ and $M_{high}^{gap}$) 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:

$$p(m|\lambda) \propto n(m|M_{low}^{gap}, M_{high}^{gap}, \eta_{low}, \eta_{high}, A) \times h(m|M_{min}, \eta_{min}) \times l(m|M_{max}, \eta_{max}) \times \begin{cases} (m/M_{high}^{gap})^{\alpha_1} & \text{if } m < M_{high}^{gap} \\ (m/M_{high}^{gap})^{\alpha_2} & \text{if } m \ge M_{high}^{gap} \end{cases} \quad (1)$$

With $1 \le m/M_\odot \le 100$;

here, $n(m|M_{low}^{gap}, M_{high}^{gap}, \eta_{low}, \eta_{high}, A) = (1 - A h(m|M_{low}^{gap}, \eta_{low})) l(m|M_{high}^{gap}, \eta_{high}))$;

where

$$h(m|M_{min}, \eta_{min}) = 1 - l(m|M_{min}, \eta_{min})$$

$$l(m|M_{max}, \eta_{max}) = (1 + (m/M_{max})^{\eta_{max}})^{-1}.$$

$h(m|m_{min}, \eta_{min})$ and $l(m|m_{max}, \eta_{max})$ are the low-mass and high-mass tapering functions, respectively.

The 1D mass distribution, $p(m|\lambda)$, is shown in Figure 3 for a specific choice of $\lambda$ ($\lambda$ represents the 12 parameters of the model; see Table 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|\lambda)$ and related via a “pairing function” (Doctor et al. 2020; Fishbach & Holz 2020). As defined in Fishbach & Holz (2020), the pairing assumed here is a power law in the mass ratio, $q \equiv m_1/m_2$. Explicitly,

$$p(m_1, m_2|\Lambda) \propto p(m=m_2|\Lambda)p(m=m_1|\Lambda) \left(\frac{m_2}{m_1}\right)^{\beta} \quad (2)$$

The values of the hyperparameters $\Lambda = \{\lambda, \beta\}$ 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 $\Lambda$. The effects of neglecting the hyperparameter uncertainty are estimated in Appendix A.1.

The model assumes a spin distribution with isotropically oriented component spins and uniform component spin magnitudes. The spin magnitude distribution for objects with
<figure><img src="image_5.png" /><figcaption>**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 \text{ 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.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** The 1D mass distribution, $p(m|\lambda)$ on the interval $[1, 100] M_\odot$ for a specific choice of hyperparameters $\lambda$. See Table 7 in Appendix A.1, for the other parameters of the mass distribution.</figcaption></figure>

masses less than ($m < 2.5 M_\odot$) is defined in the range of $[0, 0.4]$, and that for objects with masses larger than $2.5 M_\odot$ is defined in the range $[0, 1]$.

A set of $10^6$ CBCs were drawn from the 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_\odot$ and black holes as objects with masses above $3 M_\odot$. 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_\odot$ whereas the GWTC-3 distribution is drawn from a model defined up to $100 M_\odot$,

<figure><img src="image_7.png" /><figcaption>**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_\odot$, but are allowed to be as massive as $100 M_\odot$ in the GWTC-3 distribution.</figcaption></figure>

allowing for higher-mass black holes in the latter case (although the tapering of the mass distribution above $M_{\text{max}} = 54.38$ does somewhat limit the number of high-mass black holes).

Figure 4 shows the mass distributions of the CBCs sub-populations. Figure 5 shows the simulated mass distributions for each model that survives the S/N cut (see Section 2.3).

## 2.3. Simulation Campaign

We use the public software suite ligo.skymap,<fnref n="14" /> 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 (Planck Collaboration 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 parameters<fnref n="15" /> to the statistical production of sky maps,<fnref n="16" /> passing successively by the filtering of CBC events that pass the S/N cut,<fnref n="17" /> as well as their locations in the sky (see footnote 17) are publicly accessible on GitHub<fnref n="18" /> for LRR distribution, and on GitHub<fnref n="19" /> (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 area and distance distributions of O3 alerts (Petrov et al. 2022). This simulation set yields estimates of the

<footnote n="14">¹⁴ https://git.ligo.org/leo-singer/ligo.skymap</footnote>

<footnote n="15">¹⁵ https://lscsoft/docs.ligo.org/ligo.skymap/tool/bayestar_inject.html</footnote>

<footnote n="16">¹⁶ https://lscsoft/docs.ligo.org/ligo.skymap/tool/ligo_skymap_stats.html</footnote>

<footnote n="17">¹⁷ https://lscsoft/docs.ligo.org/ligo.skymap/tool/bayestar_localize_coincs.html</footnote>

<footnote n="18">¹⁸ https://github.com/lpsinger/observing-scenarios-simulations/tree/v1</footnote>

<footnote n="19">¹⁹ https://github.com/lpsinger/observing-scenarios-simulations/tree/v2</footnote>
<figure><img src="image_8.png" /><figcaption>**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 Appendix A.1.1, Figure 14.</figcaption></figure>

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), 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 (Abbott 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.<fnref n="20" /> We use the publicly available noise curves released in LIGO-T2200043-v3.<fnref n="21" /> For O4, we used aligo_O4high.txt, avirgo_O4high_NEW.txt, kagra_10Mpc.txt respectively for LIGO (LIGO Hanford Observatory (LHO), LIGO Livingston Observatory (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

<footnote n="20">²⁰ https://observing/docs.ligo.org/plan/index.html</footnote>

<footnote n="21">²¹ https://dcc.ligo.org/T2200043-v3/public</footnote>
<table><caption>**Table 1**
The Cosmology-corrected Inspiral Sensitive Distance (in Mpc) from a GW Strain PSD</caption><thead><tr><th>Run</th><th>LIGO (L1—H1)</th><th>Virgo (V1)</th><th>KAGRA (K1)</th></tr><tr><td></td><td colspan="3">BNS Inspiral Range of the in Mpc</td></tr></thead><tbody><tr><td>O4</td><td>224</td><td>145</td><td>37</td></tr><tr><td>O5</td><td>494</td><td>183</td><td>390</td></tr></tbody></table>

<table><caption>**Table 2**
The Merger Rate per Unit Comoving Volume Used for LRR and PBD/GWTC-3 (GWTC-3) Distributions</caption><thead><tr><th>Distribution</th><th>BNS</th><th>NSBH</th><th>BBH</th></tr><tr><td></td><td colspan="3">Merger Rate Density (Gpc<sup>-3</sup> yr<sup>-1</sup>)</td></tr></thead><tbody><tr><td>LRR</td><td>170<sup>+270</sup><sub>-120</sub></td><td>27<sup>+31</sup><sub>-17</sub></td><td>25<sup>+10</sup><sub>-7</sub></td></tr><tr><td>GWTC-3</td><td>210<sup>+240</sup><sub>-120</sub></td><td>8.6<sup>+9.7</sup><sub>-5</sub></td><td>17.1<sup>+19.2</sup><sub>-10</sub></td></tr></tbody></table>

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 2 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_{-140}^{+270}$ 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 2 (Abbott et al. 2023). We reproduce these astrophysical density rates in Table 2.

## 2.4. Results

We make the results of the observing scenarios, including the sky-localization FITS files, available on for PBD/GWTC-3 (GWTC-3) (Coughlin & Farah 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.

### 2.4.1. Detection Rates and Summary Statistics for O4 and 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 model is due to the existence of a nearly empty mass gap in that model, combined with a pairing

<figure><img src="image_9.png" /><figcaption>**Figure 6.** 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, Figure 15.</figcaption></figure>

<table><caption>**Table 3**
Annual Detection Rates for Compact Binary Coalescences that We Expect for the Runs O4 and O5</caption><thead><tr><th>Run</th><th>Distribution</th><th>BNS</th><th>NSBH</th><th>BBH</th></tr><tr><td></td><td></td><td colspan="3">Annual Number of Detections</td></tr></thead><tbody><tr><td rowspan="2">O4</td><td>LRR</td><td>17<sup>+35</sup><sub>-13</sub></td><td>10<sup>+18</sup><sub>-8</sub></td><td>46<sup>+23</sup><sub>-17</sub></td></tr><tr><td>GWTC-3</td><td>36<sup>+49</sup><sub>-22</sub></td><td>6<sup>+11</sup><sub>-5</sub></td><td>260<sup>+330</sup><sub>+150</sub></td></tr><tr><td rowspan="2">O5</td><td>LRR</td><td>86<sup>+171</sup><sub>-59</sub></td><td>48<sup>+71</sup><sub>-30</sub></td><td>190<sup>+80</sup><sub>-58</sub></td></tr><tr><td>GWTC-3</td><td>180<sup>+220</sup><sub>-100</sub></td><td>31<sup>+42</sup><sub>-20</sub></td><td>870<sup>+1100</sup><sub>-480</sub></td></tr></tbody></table>

**Note.** 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.

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 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. The 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 treatment process
<table><caption>**Table 4**
Summary Statistics for O4 and O5</caption><thead><tr><th>Run</th><th>Dist.</th><th>BNS</th><th>NSBH</th><th>BBH</th></tr></thead><tbody><tr><td colspan="5">Median 90% Credible Area (deg<sup>2</sup>)</td></tr><tr><td rowspan="2">O4</td><td>LRR</td><td>2100<sup>+150</sup><sub>-220</sub></td><td>2090<sup>+130</sup><sub>-130</sub></td><td>653<sup>+53</sup><sub>-36</sub></td></tr><tr><td>GWTC-3</td><td>1860<sup>+250</sup><sub>-170</sub></td><td>2140<sup>+480</sup><sub>-530</sub></td><td>1428<sup>+60</sup><sub>-55</sub></td></tr><tr><td rowspan="2">O5</td><td>LRR</td><td>2050<sup>+100</sup><sub>-160</sub></td><td>2110<sup>+100</sup><sub>-100</sub></td><td>682<sup>+25</sup><sub>-30</sub></td></tr><tr><td>GWTC-3</td><td>2050<sup>+120</sup><sub>-120</sub></td><td>2000<sup>+350</sup><sub>-220</sub></td><td>1256<sup>+48</sup><sub>-53</sub></td></tr><tr><td colspan="5">Median 90% Credible Comoving Volume (10<sup>6</sup>Mpc<sup>3</sup>)</td></tr><tr><td rowspan="2">O4</td><td>LRR</td><td>46.5<sup>+6.6</sup><sub>-7.0</sub></td><td>159<sup>+26</sup><sub>-16</sub></td><td>207<sup>+21</sup><sub>-20</sub></td></tr><tr><td>GWTC-3</td><td>67.9<sup>+11.3</sup><sub>-9.9</sub></td><td>232<sup>+101</sup><sub>-50</sub></td><td>3400<sup>+310</sup><sub>-240</sub></td></tr><tr><td rowspan="2">O5</td><td>LRR</td><td>240<sup>+29</sup><sub>-26</sub></td><td>785<sup>+68</sup><sub>-62</sub></td><td>857<sup>+63</sup><sub>-60</sub></td></tr><tr><td>GWTC-3</td><td>376<sup>+36</sup><sub>-40</sub></td><td>1350<sup>+290</sup><sub>-300</sub></td><td>8580<sup>+600</sup><sub>-550</sub></td></tr><tr><td colspan="5">Median Luminosity Distance (Mpc)</td></tr><tr><td rowspan="2">O4</td><td>LRR</td><td>349<sup>+12</sup><sub>-14</sub></td><td>564<sup>+15</sup><sub>-13</sub></td><td>1102<sup>+33</sup><sub>-32</sub></td></tr><tr><td>GWTC-3</td><td>398<sup>+15</sup><sub>-14</sub></td><td>770<sup>+67</sup><sub>-70</sub></td><td>2685<sup>+53</sup><sub>-40</sub></td></tr><tr><td rowspan="2">O5</td><td>LRR</td><td>619<sup>+15</sup><sub>-19</sub></td><td>1007<sup>+20</sup><sub>-22</sub></td><td>1948<sup>+34</sup><sub>-24</sub></td></tr><tr><td>GWTC-3</td><td>738<sup>+30</sup><sub>-25</sub></td><td>1318<sup>+71</sup><sub>-100</sub></td><td>4607<sup>+77</sup><sub>-82</sub></td></tr><tr><td colspan="5">Sensitive Volume: Detection Rate and/or Merger Rate (Gpc<sup>3</sup>)</td></tr><tr><td colspan="5">Sensitive volume = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">detection rate</div><div style="padding: 0 5px;">merger rate</div></div></td></tr><tr><td rowspan="2">O4</td><td>LRR</td><td>0.1011<sup>+0.0066</sup><sub>-0.0064</sub></td><td>0.403<sup>+0.021</sup><sub>-0.020</sub></td><td>1.861<sup>+0.077</sup><sub>-0.074</sub></td></tr><tr><td>GWTC-3</td><td>0.172<sup>+0.013</sup><sub>-0.012</sub></td><td>0.78<sup>+0.14</sup><sub>-0.13</sub></td><td>15.15<sup>+0.42</sup><sub>-0.41</sub></td></tr><tr><td rowspan="2">O5</td><td>LRR</td><td>0.507<sup>+0.027</sup><sub>-0.026</sub></td><td>1.809<sup>+0.070</sup><sub>-0.068</sub></td><td>7.62<sup>+0.19</sup><sub>-0.19</sub></td></tr><tr><td>GWTC-3</td><td>0.827<sup>+0.044</sup><sub>-0.042</sub></td><td>3.65<sup>+0.47</sup><sub>-0.43</sub></td><td>50.7<sup>+1.2</sup><sub>-1.2</sub></td></tr></tbody></table>

Note. 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.

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 13).

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 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 detection 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.

### 3. Predictions for Detection Rates and Science with Gravitational-wave Counterparts

In this section, we use the results from Section 2 to make predictions for potential counterpart detections during O4 and O5.

#### 3.1. Light Curves

To simulate the light curves corresponding to the BNSs and NSBHs from Section 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. (2020) 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, 2023). For the BNS case, we use the same geometry and lanthanide fractions for each component as presented in Dietrich et al. (2020), shown to yield good fits to GW170817. For the NSBH case, we use the same configurations as presented in Anand et al. (2021). This model takes as input the dynamical ejecta $M_{ej}^{dyn}$, disk wind ejecta $M_{ej}^{wind}$, half opening angle $\Phi$, and the observation angle $\Theta_{obs}$ (see Dietrich et al. 2020 for details). We take the inclination measurement from the simulated value in the observing scenarios, and we draw $\Phi$ uniformly between 15° and 75°. To interpolate the grid-based model, we use a Gaussian process regression method (Coughlin et al. 2018a, 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 GitHub.<fnref n="22" />

#### 3.2. 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 OPTimization<fnref n="23" /> (gwemopt; Coughlin et al. 2018b, 2019c). We take ZTF to have a base sensitivity in g, r, and i bands with $g \sim 21.7$ mag, $r \sim 21.4$ mag, $i \sim 20.9$ mag (Almualla et al. 2021), and Rubin Observatory with u, g, r, i, z, and y bands with $u = 23.9$ mag, $g = 25.0$ mag, $r = 24.7$ mag, $i = 24.0$ mag, $z = 23.3$ mag, 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.2$ mag, $H = 22.6$ mag, and $K = 22.6$ mag. Figure 8 shows a 1D 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

<footnote n="22">²² https://github.com/nuclear-multimessenger-astronomy/nmma</footnote>

<footnote n="23">²³ https://github.com/mcoughlin/gwemopt</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** 2D 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.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

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
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The plots at the top are BNS distributions of the missed and found injections (on the left, BNS light curves observed by ZTF, and the right one shows the same using Rubin Observatory), as a function of the peak *r*-band absolute magnitude and 2D probability enclosed. We encode in the color bar and the transparency of the points the 3D probability of finding the transient as measured by simulated injections. The 1D histograms are marginalized versions of the 2D scatter plot, with the detected set in green and missed set in red. The plots at the bottom are the same for the NSBH light curves. Those plots concern the run O4. In Appendix B.3, Figure 15 shows the same distribution for the run O5.</figcaption></figure>

<figure><img src="image_13.png" /></figure>

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² 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 photometry uncertainties. Finally, the forced photometry upper limits are estimated using KDE fits and used to reject 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², 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 (Table 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.

### 3.3. *KNe Sample Constraints*

For the objects detected by this process, we perform parameter estimation of the resulting EM light curves. To do so, we use the *NMMA* framework (Pang et al. 2022), designed to
<table><caption>**Table 5**
Annual Detection Rate of EM Counterparts that We Expect for ZTF and Rubin Observatory during the Runs O4 and O5</caption><thead><tr><th>Run</th><th>Telescope</th><th>BNS</th><th>NSBH</th></tr></thead><tbody><tr><td></td><td colspan="3">EM Annual Number of Detections</td></tr><tr><td rowspan="2">O4</td><td>ZTF</td><td>0.43<sup>+0.58</sup><sub>-0.26</sub></td><td>0.13<sup>+0.24</sup><sub>-0.11</sub></td></tr><tr><td>Rubin</td><td>1.97<sup>+2.68</sup><sub>-1.2</sub></td><td>0.03<sup>+0.06</sup><sub>-0.03</sub></td></tr><tr><td rowspan="2">O5</td><td>ZTF</td><td>0.43<sup>+0.44</sup><sub>-0.2</sub></td><td>0.09<sup>+0.12</sup><sub>-0.06</sub></td></tr><tr><td>Rubin</td><td>5.39<sup>+5.59</sup><sub>-2.99</sub></td><td>0.43<sup>+0.59</sup><sub>-0.28</sub></td></tr></tbody></table>

perform Bayesian inference of multimessenger signals, incorporating all available data on the neutron star EoS in the process (Dietrich et al. 2020; 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(\theta|x, M)$, for model source parameters $\theta$ under the hypothesis or model $M$ with mock-up data $x$ as

$$p(\theta|x, M) = \frac{p(x|\theta, M) p(\theta|M)}{p(x|M)}, \quad (3)$$

where $p(x|\theta, M)$, $p(\theta|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. 2020) and the detection of the shortest long γ-ray burst ever confirmed (Ahumada et al. 2021).

### 3.3.1. Ejecta Constraints

We begin by evaluating 90% upper limits possible from the sample considered here (KNe detections for these objects), then constrain them for the ejecta model parameters $M_{ej}^{dyn}$ and $M_{ej}^{wind}$, as well as a systematic contribution to the dynamical ejecta $\alpha$ and the fraction of the disk mass contributing $\zeta$. This enables us to make an empirical constraint of the fraction of the disk contributing to KNe. We constrain $M_{ej}^{dyn}$ to 10%–40% and $M_{ej}^{wind}$ 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_{ej}^{wind} = [-2.8 to 2.0]$, resulting in strong constraints on the presence of such a disk.

### 3.3.2. ZTF Proposal for GW Follow-up and Triggering Criteria

The properties of a GW event released in low-latency via the General Coordinates Network notice can prove to be useful in determining whether or not to trigger telescope time on a given GW event. Information such as the false 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”

<figure><img src="image_14.png" /><figcaption>**Figure 10.** 90% upper limit on the $M_{ej}^{wind}$ measurements for the prompt collapse simulation set, where the $M_{ej}^{wind}$ contribution is set to 0.</figcaption></figure>

<table><caption>**Table 6**
Triggering Based on GW Candidate Event Properties</caption><thead><tr><th>Parameter</th><th>Go-deep</th><th>Go-wide</th><th>Deliberate</th><th>No Go</th></tr></thead><tbody><tr><td>Strategy</td><td>300 s<br>Push<br>distance</td><td>30 s<br>Push<br>localization</td><td colspan="2">Action item: human<br>interaction</td></tr><tr><td>Frequency of<br>triggers</td><td>1 per<br>month<br>3 nights</td><td>2 per month<br>5 nights</td><td colspan="2"></td></tr><tr><td>FAR min(FAR)<br>– “Best”</td><td>&lt; 1 per<br>century<br>Any<br>pipeline</td><td>&lt; 1 per<br>decade<br>Any pipeline</td><td>1 per year<br>–century</td><td>&gt; 1 per year<br>All pipelines</td></tr><tr><td>max(p-astro)</td><td>&gt; 0.9</td><td>&gt; 0.9</td><td>0.1–0.9</td><td>&lt; 0.1</td></tr><tr><td>HasNS</td><td>&gt; 0.9</td><td>&gt; 0.9</td><td>0.1–0.9</td><td>&lt; 0.1</td></tr><tr><td>log(BCI)</td><td>&gt; 4</td><td>&gt; 4</td><td>-1 to 4</td><td>&lt; -1</td></tr><tr><td>HasRemnant?</td><td>&gt; 0</td><td>&gt; 0</td><td>...</td><td>= 0</td></tr><tr><td>pBNS/pNSBH</td><td>&gt; 0</td><td>&gt; 0</td><td>...</td><td>= 0</td></tr></tbody></table>

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 \text{ day}^{-1}$ (user guide$^{24}$; 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} \sim -16$ mag at peak, and taking to account ZTF’s limiting magnitude in 300 s exposures ($m_{AB} \approx 22$ mag), we estimate that ZTF could detect KNe falling within 400 Mpc.

<footnote>$^{24}$ https://emfollow/docs.ligo.org/userguide/analysis/index.html</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 11.** 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.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**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²), we could fully probe the localization for nearly 40% of all events.</figcaption></figure>

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 Section 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², 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.

### 3.4. $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 IMRPHE-NOMD_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 & Harsh (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 (2019, 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) high-mass 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

$$cz \approx H_0 d_L, \quad (4)$$

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_0^4$ prior factor, we correct for this in our study similar to Abbott et al. (2017d). 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 distance-inclination angle degeneracy by including additional information on the binary's orbital inclination from other 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 2(a) 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 counterpart 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.74 \text{ km s}^{-1} \text{ Mpc}^{-1}$; this is mainly caused by the large uncertainty of our measurement $H_0 = 60.68^{+9.24}_{-7.47} \text{ km s}^{-1} \text{ 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} \text{ km s}^{-1} \text{ 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} \text{ km s}^{-1} \text{ Mpc}^{-1}$. We do a similar study for O5, combining 23 events, and estimate $H_0 = 66.74^{+0.39}_{-0.33} \text{ km s}^{-1} \text{ Mpc}^{-1}$. Both estimates recover the injected value of $H_0$.

In Figure 13 (bottom panels), we show our $H_0$-results 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} \text{ km s}^{-1} \text{ Mpc}^{-1}$ and approximately 5 combined observations in O5 with $H_0 = 65.30^{+2.31}_{-2.99} \text{ km s}^{-1} \text{ 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} \text{ km s}^{-1} \text{ Mpc}^{-1}$, which is close to the injection value. For O5, we recover $H_0 = 66.23^{+0.39}_{-0.33} \text{ km s}^{-1} \text{ 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.

## 4. 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
<figure><img src="image_17.png" /><figcaption>**Figure 13.** *H*<sub>0</sub>-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 *H*<sub>0</sub>-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 (Planck Collaboration et al. 2016; Planck, violet), to the Hubble measurement via type-Ia supernovae (Riess et al. 2016; SHOES, light blue), and to the *H*<sub>0</sub>-measurement of superluminal motion of the jet in GW170817 (Hotokezaka et al. 2019; superluminal, gray). All uncertainties are reported at 90% credible interval.</figcaption></figure>

can summarize our conclusions and implications for future GW
follow-ups 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 nearby Chen & Holz (2016). Proximity is an especially important consideration for 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*₀, most events contribute equally, especially as relative uncertainty due to virial velocities decreases as distance increases.

Acknowledgments

M.W.C. acknowledges support from the National Science
Foundation with grant Nos. PHY-2010970 and OAC-2117997.

M.W.C. and R.W.K. were supported by the Preparing for Astrophysics with LSST Program, funded by the Heising Simons Foundation through grant 2021-2975, and administered by Las Cumbres Observatory. A.M.F. is supported by the National Science Foundation Graduate Research Fellowship Program under grant No. DGE-1746045.

The authors are very grateful to Sharan Banagiri for his
review on behalf of Publication and Presentation (P&P), which
has been useful in the improving of this paper. We would like
to thank the Cosmology GW group and in particular Rachel
Gray and Thomas Dent for the suggestions and comments that
allowed us to improve our paper.

The authors acknowledge the International Gravitational-
wave Network Computing Grid (CIT, LHO, LLO) for
providing resources to realize our simulations. We acknowl-
edge the Minnesota Supercomputing Institute (MSI; http://
msi.umn.edu) at the University of Minnesota for providing
resources that contributed to the research results reported
within this paper under project "Identification of Variable
Objects in the Zwicky Transient Facility." This research used
resources of the National Energy Research Scientific Comput-
ing Center (NERSC), a U.S. Department of Energy Office of
Science User Facility operated under contract No. DE-AC02-
05CH11231 under project "Toward a complete catalog of
variable sources to support efficient searches for compact
binary mergers and their products.” This work used the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS). This work used the Extreme Science and Engineering Discovery Environment (XSEDE) COMET at SDSU through allocation AST200016 and AST200029. This material is based upon work supported by NSF’s LIGO Laboratory, which is a major facility fully funded by the National Science Foundation.

# Appendix A
## Observing Scenarios

### A.1. Values of Hyperparameters

We fix the parameters in both of our assumed models for the astrophysical distribution of compact binary mergers. These population model parameters (termed “hyperparameters” when used in a hierarchical Bayesian analysis) are listed in Table 7 for the PBD/GWTC-3 model and were chosen because they correspond to the maximum value of the hyperposterior obtained by fitting the model to GWTC-3 (Abbott et al. 2023). We are therefore not considering the full posterior uncertainty in these hyperparameters. To estimate the effect of neglecting this uncertainty on the number of detected events in each subpopulation, we compare the uncertainty astrophysical merger rate predicted the full PDB analysis presented in Abbott et al. (2023) to those assumed in this work. The percent uncertainty (90% CI) in Table 2 for the PDB/GWTC-3 model is

171% each for BNS, NSBH, and BBH. The percent error reported by Abbott et al. (2023) is 229%, 178%, and 68% for BNS, NSBH, and BBH events, respectively. We therefore conclude that we underestimate the uncertainty in merger rate BNS and NSBH by factors of 1.34 and 1.04, respectively, and we overestimate the errors for the BBH merger rate by a factor of 2.5. This is likely because the shape of the BBH mass spectrum is well-measured compared to that of the NS-containing events, making the extra uncertainty introduced by fitting the hyperparameters relatively small. The opposite is true for NSBHs and BNSs since only ~4 events are used to constrain the shape of the low-mass end of the distribution.

Naively translating the additional uncertainty in the astrophysical merger rate to an uncertainty in the annual number of detected events in O4 yields an increase in uncertainty of 12.2 events for BNS, an increase in uncertainty of 0.24 events for NSBH, and a decrease in uncertainty of 104 events for BBH, assuming symmetric errors. These corrections are approximate and are meant to give an estimate of the effect of neglecting the uncertainty in population hyperparameters under the PBD/GWTC-3 model. A similar estimate is not available for the LRR model, as the authors do not know of a population fit performed using that framework on GWTC-3 data.

The difference between the maximum a posteriori underlying population and the full underlying population hyperposterior can be seen in Figure 5 of Farah et al. (2022).

<table><caption>**Table 7**
Summary of Population Model Parameters</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>α<sub>1</sub></td><td>Spectral index for the power law of the mass distribution at low mass.</td><td>-2.16</td></tr><tr><td>α<sub>2</sub></td><td>Spectral index for the power law of the mass distribution at high mass.</td><td>-1.46</td></tr><tr><td>A</td><td>Lower-mass gap depth.</td><td>0.97</td></tr><tr><td>M<sub>gap low</sub></td><td>Location of lower end of the mass gap.</td><td>2.72 M<sub>☉</sub></td></tr><tr><td>M<sub>gap high</sub></td><td>Location of upper end of the mass gap.</td><td>6.13 M<sub>☉</sub></td></tr><tr><td>η<sub>how</sub></td><td>Parameter controlling how the rate tapers at the low end of the mass gap.</td><td>50</td></tr><tr><td>η<sub>high</sub></td><td>Parameter controlling how the rate tapers at the high end of the mass gap.</td><td>50</td></tr><tr><td>η<sub>min</sub></td><td>Parameter controlling tapering of the power law at low mass.</td><td>50</td></tr><tr><td>η<sub>max</sub></td><td>Parameter controlling tapering of the power law at high mass.</td><td>4.91</td></tr><tr><td>β</td><td>Spectral index for the power-law-in-mass-ratio pairing function.</td><td>1.89</td></tr><tr><td>M<sub>min</sub></td><td>Minimum mass of the mass distribution.</td><td>1.16 M<sub>☉</sub></td></tr><tr><td>M<sub>max</sub></td><td>Onset location of high-mass tapering.</td><td>54.38 M<sub>☉</sub></td></tr><tr><td>a<sub>max,NS</sub></td><td>Maximum allowed component spin for objects with mass &lt; 2.5 M<sub>☉</sub>.</td><td>0.4</td></tr><tr><td>a<sub>max,BH</sub></td><td>Maximum allowed component spin for objects with mass ≥ 2.5 M<sub>☉</sub>.</td><td>1</td></tr></tbody></table>

**Note.** The first several entries describe the rate and mass distribution parameters, and the last two entries describe the spin distribution parameters.
A.1.1. Simulated Mass Distributions for O5

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

<figure><img src="image_18.png" /><figcaption>**Figure 14.** 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.</figcaption></figure>
A.1.2. Histogram of Simulated Detections for O5

In Figure 15, we summarize the O5 detection results of the
simulation set.

<figure><img src="image_19.png" /><figcaption>**Figure 15.** 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 in our simulation.</figcaption></figure>
# Appendix B
## Predictions for Detection Rates and Science with Gravitational-wave Counterparts

### B.1. 2D Histogram of the O5 Simulated BNS and NSBH Light Curves

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.

<figure><img src="image_20.png" /><figcaption>**Figure 16.** 2D histograms of simulated BNS (left) and NSBH (right) light curves for O5. 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 dashed 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.</figcaption></figure>
## B.2. Histograms of Peak Magnitudes

In Figure 17, we 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><img src="image_21.png" /><figcaption>**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 runs 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.</figcaption></figure>
### 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><img src="image_22.png" /><figcaption>**Figure 18.** The plots at the top are BNS distributions of the missed and found injections (on the left, BNS light curves observed by ZTF, and the right one shows the same observation in Rubin Observatory), as a function of the peak *r*-band absolute magnitude and 2D probability enclosed. We encode in the color bar and the transparency of the points the 3D probability of finding the transient as measured by simulated injections. The 1D histograms are marginalized versions of the 2D scatter plot, with the detected set in green, and missed set in red. The plots at the bottom are the same for the NSBH light curves. Those plots concern the run O5.</figcaption></figure>
ORCID iDs

R. Weizmann Kiendrebeogo https://orcid.org/0000-0002-9108-5059

Amanda M. Farah https://orcid.org/0000-0002-6121-0285

Emily M. Foley https://orcid.org/0000-0003-1955-2983

Abigail Gray https://orcid.org/0000-0001-9763-2351

Nina Kunert https://orcid.org/0000-0002-1275-530X

Anna Puecher https://orcid.org/0000-0003-1357-4348

Andrew Toivonen https://orcid.org/0009-0008-9546-2035

R. Oliver VandenBerg https://orcid.org/0000-0002-7826-6269

Shreya Anand https://orcid.org/0000-0003-3768-7515

Tomás Ahumada https://orcid.org/0000-0002-2184-6430

Viraj Karambelkar https://orcid.org/0000-0003-2758-159X

Michael W. Coughlin https://orcid.org/0000-0002-8262-2924

Tim Dietrich https://orcid.org/0000-0003-2374-307X

S. Zacharie Kam https://orcid.org/0000-0002-8146-0177

Peter T. H. Pang https://orcid.org/0000-0001-7041-3239

Leo P. Singer https://orcid.org/0000-0001-9898-5597

References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, LRR, 19, 1

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, ApJL, 848, L12

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, PhRvL, 119, 161101

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017c, ApJL, 848, L13

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017d, Natur, 551, 85

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, LRR, 21, 3

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019a, PhRvX, 9, 031040

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019b, PhRvX, 9, 011001

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019c, PhRvX, 9, 031040

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020a, ApJL, 892, L3

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020b, LRR, 23, 3

Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, ApJL, 915, L5

Abbott, R., Abbott, T. D., Acernese, F., et al. 2023, PhRvX, 13, 011048

Acernese, F., Agathos, M., Agatsuma, O., et al. 2014, CQGra, 32, 024001

Ackley, K., Amati, L., Barbieri, C., et al. 2020, A&A, 643, A113

Ahumada, T., Singer, L. P., Anand, S., et al. 2021, NatAs, 5, 917

Akutsu, T., Ando, M., Arai, K., et al. 2021, PTEP, 2021, 05A102

Almualla, M., Anand, S., Coughlin, M. W., et al. 2021, MNRAS, 504, 2822

Anand, S., Coughlin, M. W., Kasliwal, M. M., et al. 2021, NatAs, 5, 46

Andreoni, I., Coughlin, M. W., Kool, E. C., et al. 2021, ApJ, 918, 63

Andreoni, I., Goldstein, D. A., Anand, S., et al. 2019, ApJL, 881, L16

Andreoni, I., Goldstein, D. A., Kasliwal, M. M., et al. 2020, ApJ, 890, 131

Andreoni, I., Margutti, R., Salafia, O. S., et al. 2022, ApJS, 260, 18

Annala, E., Gorda, T., Kurkela, A., & Vuorinen, A. 2018, PhRvL, 120, 172703

Antier, S., Agayeva, S., Aivazyan, V., et al. 2020a, MNRAS, 492, 3904

Antier, S., Agayeva, S., Almualla, M., et al. 2020b, MNRAS, 497, 5518

Ashton, G., Hübner, M., Lasky, P. D., et al. 2019, ApJS, 241, 27

Baumgarte, T. W., Shapiro, S. L., & Shibata, M. 2000, ApJL, 528, L29

Bauswein, A., Just, O., Janka, H.-T., & Stergioulas, N. 2017, ApJL, 850, L34

Bellm, E. C., Kulkarni, S. R., Barlow, T., et al. 2019, PASP, 131, 068003

Bulla, M. 2019, MNRAS, 489, 5037

Bulla, M. 2023, MNRAS, 520, 2558

Chattopadhyay, D., Stevenson, S., Broekgaarden, F., Antonini, F., & Belczynski, K. 2022, MNRAS, 513, 5780

Chen, H.-Y., & Holz, D. E. 2016, arXiv:1612.01471

Chornock, R., Berger, E., Kasen, D., et al. 2017, ApJL, 848, L19

Colombo, A., Salafia, O. S., Gabrielli, F., et al. 2022, ApJ, 937, 79

Coughlin, M. W., Ahumada, T., Anand, S., et al. 2019b, ApJL, 885, L19

Coughlin, M. W., Antier, S., Corre, D., et al. 2019c, MNRAS, 489, 5775

Coughlin, M. W., Antier, S., Dietrich, T., et al. 2020b, NatCo, 11, 4129

Coughlin, M. W., Dietrich, T., Antier, S., et al. 2020c, MNRAS, 492, 863

Coughlin, M. W., Dietrich, T., Doctor, Z., et al. 2018a, MNRAS, 480, 3871

Coughlin, M. W., Dietrich, T., Heinzel, J., et al. 2020a, PhRvR, 2, 022006

Coughlin, M. W., Dietrich, T., Margalit, B., & Metzger, B. D. 2019a, MNRAS Lett., 489, L91

Coughlin, M. W., Farah, A. M., Foley, E. M., et al. 2022, LIGO/Virgo/KAGRA Observing Capabilities: Simulated Detections and Localization for

O4 and O5 (October 2022 edition), v1, Zonodo, doi:10.5281/zenodo.7026209

Coughlin, M. W., Tao, D., Chan, M. L., et al. 2018b, MNRAS, 478, 692

Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Sci, 358, 1556

Cowperthwaite, P. S., Berger, E., Villar, V. A., et al. 2017, ApJL, 848, L17

Dekany, R., Smith, R. M., Riddle, R., et al. 2020, PASP, 132, 038001

Dietrich, T., Coughlin, M. W., Pang, P. T. H., et al. 2020, Sci, 370, 1450

Dietrich, T., Samajdar, A., Khan, S., et al. 2019, PhRvD, 100, 044003

Doctor, Z., Wysocki, D., O'Shaughnessy, R., Holz, D. E., & Farr, B. 2020, ApJ, 893, 35

Farah, A., Fishbach, M., Essick, R., Holz, D. E., & Galaudage, S. 2022, ApJ, 931, 108

Farr, W. M., Sravan, N., Cantrell, A., et al. 2011, ApJ, 741, 103

Fishbach, M., Essick, R., & Holz, D. E. 2020, ApJL, 899, L8

Fishbach, M., & Holz, D. E. 2020, ApJL, 891, L27

Foucart, F., Hinderer, T., & Nissanke, S. 2018, PhRvD, 98, 081501

Godwin, P., Essick, R., Hanna, C., et al. 2020, arXiv:2010.15282

Goldstein, A., Veres, P., Burns, E., et al. 2017, ApJL, 848, L14

Goldstein, D. A., Andreoni, I., Nugent, P. E., et al. 2019, ApJL, 881, L7

Gomez, S., Hosseinzadeh, G., Cowperthwaite, P. S., et al. 2019, ApJL, 884, L55

Gompertz, B. P., Cutter, R., Steeghs, D., et al. 2020, MNRAS, 497, 726

Graham, M. J., Kulkarni, S. R., Bellm, E. C., et al. 2019, PASP, 131, 078001

Hosseinzadeh, G., Cowperthwaite, P. S., Gomez, S., et al. 2019, ApJL, 880, L4

Hotokezaka, K., Nakar, E., Gottlieb, O., et al. 2019, NatAs, 3, 940

Huth, S., Pang, P. T. H., Tews, I., et al. 2022, Natur, 606, 276

Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111

Janquart, J., & Harsh, N., 2022 RelativeBilbying: a package for relative binning with bilby, https://github.com/lemmis12/relativebilbying

Kasliwal, M. M., Anand, S., Ahumada, T., et al. 2020, ApJ, 905, 145

Kasliwal, M. M., Kasen, D., Lau, R. M., et al. 2019, MNRAS Lett., 510, L1

Kiendrebeogo, R. W., Foley, E. M., Coughlin, M. W., & Singer, L. P., 2023 Updated simulation of the next O4 and O5 observation and detection scenarios of the LIGO-Virgo-KAGRA collaboration using the LRR distribution described in Petrov et al.2022 (September 2022 edition), v1, Zenodo, doi:10.5281/zenodo.7623166

Koposov, S., Speagle, J., Barbary, K., et al., 2022 joshspeagle/dynesty: v1.2.2, v1.2.2, Zenodo, doi:10.5281/zenodo.6456387

Krüger, C. J., & Foucart, F. 2020, PhRvD, 101, 103002

Kuin, N. P. M. & Swift Team 2019, GCN, 24767, https://gcn.nasa.gov/circulars/24767.txt

Lai, X., Zhou, E., & Xu, R. 2019, EPJA, 55, 60

Legred, I., Chatziioannou, K., Essick, R., Han, S., & Landry, P. 2021, PhRvD, 104, 063003

Leslie, N., Dai, L., & Pratten, G. 2021, PhRvD, 104, 123030

LSC, Aasi, J., Abbott, B. P., et al. 2015, CQGra, 32, 074001

Lundquist, M. J., Paterson, K., Fong, W., et al. 2019, ApJL, 881, L26

Margalit, B., & Metzger, B. D. 2017, ApJL, 850, L19

Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, PASP, 131, 018003

Most, E. R., Weih, L. R., Rezzolla, L., & Schaffner-Bielich, J. 2018, PhRvL, 120, 261103

Nitz, A. H., Dent, T., Davies, G. S., & Harry, I. 2020, ApJ, 897, 169

Özel, F., & Freire, P. 2016, ARA&A, 54, 401

Özel, F., Psaltis, D., Narayan, R., & McClintock, J. E. 2010, ApJ, 725, 1918

Pang, P. T. H., Dietrich, T., Coughlin, M. W., et al. 2022, arXiv:2205.08513

Pang, P. T. H., Tews, I., Coughlin, M. W., et al. 2021, ApJ, 922, 14

Petrov, P., Singer, L. P., Coughlin, M. W., et al. 2022, ApJ, 924, 84

Pian, E., D'Avanzo, P., Benetti, S., et al. 2017, Natur, 551, 67

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13

Radice, D., Perego, A., Zappa, F., & Bernuzzi, S. 2018, ApJL, 852, L29

Rezzolla, L., Most, E. R., & Weih, L. R. 2018, ApJL, 852, L25

Riess, A. G., et al. 2016, ApJ, 826, 56

Romero-Shaw, I. M., et al. 2020, MNRAS, 499, 3295

Rosswog, S., Feindt, U., Korobkin, O., et al. 2017, CQGra, 34, 104001

Ruiz, M., Shapiro, S. L., & Tsokaros, A. 2018, PhRvD, 97, 021501

Saleem, M., Resmi, L., Arun, K. G., & Mohan, S. 2020, ApJ, 891, 130

Salpeter, E. E. 1955, ApJ, 121, 161

Savchenko, V., Ferrigno, C., Kuulkers, E., et al. 2017, ApJL, 848, L15

Shibata, M., Zhou, E., Kiuchi, K., & Fujibayashi, S. 2019, PhRvD, 100, 023015

Singer, L., Kiendrébéogo, W., & Tnarikawa, 2022 Ipsinger/observing-scenarios-simulations: Version 2, v2, Zenodo, doi:10.5281/zenodo.7305534

Singer, L. P., & Price, L. R. 2016, PhRvD, 93, 024013

Singer, L. P., Price, L. R., Farr, B., et al. 2014, ApJ, 795, 105
Skilling, J. 2006, *BayAn*, 1, 833

Smartt, S. J., Chen, T.-W., Jerkstrand, A., et al. 2017, *Natur*, 551, 75

Song, H.-R., Ai, S.-K., Wang, M.-H., et al. 2019, *ApJL*, 881, L40

Speagle, J. S. 2020, *MNRAS*, 493, 3132

Stergioulas, N. 2003, *LRR*, 6, 3

Tews, I., Pang, P. T. H., Dietrich, T., et al. 2021, *ApJL*, 908, L1

Veitch, J., & Vecchio, A. 2010, *PhRvD*, 81, 062003

Wang, S., & Zhao, Z.-C. 2022, *EPJC*, 82, 9

Watson, D., Hansen, C. J., Selsing, J., et al. 2019, *Natur*, 574, 497

Zackay, B., Dai, L., & Venumadhav, T. 2018, arXiv:1806.08792