## 1. Introduction

Multimessenger gravitational-wave astronomy is a new field that already boasts a wealth of new scientific achievements, despite only having one detected event to date. The future of this field relies on continued improvement in both gravitational-wave and traditional astronomical instruments. Third-generation gravitational-wave observatories such as the Einstein Telescope (Punturo et al. Reference Punturo2010) and Cosmic Explorer (Abbott et al. Reference Abbott2017a) are billion-dollar scale instruments with approximately a factor ten improvement in sensitivity over current, second-generation observatories like Advanced LIGO (Aasi et al. Reference Aasi2015), Advanced Virgo (Acernese et al. Reference Acernese2015), and KAGRA (Akutsu et al. Reference Akutsu, Ando, Arai, Arai, Araki and Araya2019). Slated for the mid-2030s, these third-generation observatories still require significant technology development to see them realised. To this end, so-called 2.5-generation observatories have been proposed that serve as both dedicated astronomy observatories as well as technology drivers and facilitators for third-generation observatories. These 2.5-generation observatories include broadband interferometers such as A+ (Miller et al. Reference Miller, Barsotti, Vitale, Fritschel, Evans and Sigg2015; Abbott et al. Reference Abbott2018) and Voyager (Adhikari et al. Reference Adhikari2020), as well as dedicated kilohertz observatories (Martynov et al. Reference Martynov2019) such as the newly proposedNeutron Star Extreme Matter Observatory, NEMO (Ackleyet al. Reference Ackley2020).

The NEMO design proposal held neutron star physics from binary neutron star mergers as the primary science driver. It proposed comparable high-frequency ( ${\gtrsim}1 \mathrm{kHz}$ ) sensitivity to Einstein Telescope and Cosmic Explorer, sacrificing low-frequency sensitivity ( ${\lesssim}500 \mathrm{Hz}$ ) to significantly reduce costs and the requirement for new technologies. The original proposal argued that NEMO must co-exist with 2.5-generation instruments to maximise scientific impact—those instruments would provide source localisation with their superior low-frequency sensitivity. However, during the proposed period of operation (late 2020 s/early 2030 s), there may be significant periods of time when NEMO must operate without a global network of gravitational-wave observatories. In this work, we show that NEMO-like observatories can operate in isolation without a heterogeneous global array because of the presence of electromagnetic telescopes with survey and/or all-sky capabilities that will make routine detections of the electromagnetic counterparts of binary neutron star mergers out to relevant distances.

The Vera C. Rubin observatory (Ivezić et al. Reference Ivezić2019) is an8.4-m optical/UV telescope capable of detecting kilonovae associated with binary neutron star mergers. Providing they are in the correct hemisphere, Vera Rubin will detect 100% of kilonovae using only $30 \mathrm{s}$ exposures out to ${\sim}450\,\mathrm{Mpc}$ (Chen et al. Reference Chen, Cowperthwaite, Metzger and Berger2021), approximately equivalent to the maximum horizon distance of NEMO (Ackley et al. Reference Ackley2020). In optical/UV, Vera Rubin is the tip of the iceberg; other likely facilities with even better sensitivity include the James Webb Space Telescope (Gardner et al. Reference Gardner2006) and Keck Wide-Field imager (Gillingham et al. Reference Gillingham, Cooke, Glazebrook, Mould, Smith and Steidel2020) as well as the thirty-metre class of ground-based facilities such as the Giant Magellan Telescope (GMT; Johns et al. Reference Johns, Stepp, Gilmozzi and Hall2012), and the Thirty-Metre Telescope (TMT; Skidmore, TMT International Science Development Teams, & Science Advisory Committee 2015). These will be complemented with all-sky survey facilities including Zwicky Transient Facility (ZTF; Bellm et al. Reference Bellm2019) and Pan-STARRS (Chambers et al. Reference Chambers2016), capable of providing nightly all-sky surveys to observe transient kilonovae and gamma-ray burst afterglows. Observatories like Evryscope may also provide continuous all-sky coverage (Law et al. Reference Law2015).

By the late 2020 s, there will likely be multiple successor x-ray and gamma-ray telescopes to the Neil Gehrel’s *Swift* Observatory (Gehrels et al. Reference Gehrels2004) and Chandra (Weisskopf et al. Reference Weisskopf, Tananbaum, Van Speybroeck, O’Dell, Truemper and Aschenbach2000), including proposed instruments such as ATHENA (Barcons et al. Reference Barcons2017) and Theseus (Amati et al. Reference Amati2018), and currently operational or upcoming facilities like SVOM (Wei et al. Reference Wei2016) and eROSITA (Merloni et al. Reference Merloni2012). Current instruments have the capability to detect on-axis bursts well beyond the horizon of NEMO (e.g., Howell et al. Reference Howell, Ackley, Rowlinson and Coward2019), and proposed missions will be more sensitive with equivalent or better sky coverage, and faster slew times. These instruments even offer the possibility of detecting binary neutron star precursor emission (e.g., Most & Philippov Reference Most and Philippov2020; Sridhar et al. Reference Sridhar, Zrake, Metzger, Sironi and Giannios2021; Ascenzi et al. Reference Ascenzi, Oganesyan, Branchesi and Ciolfi2021), as well as potential emission from the central engine following the merger (e.g., Sarin & Lasky Reference Sarin and Lasky2021). Afterglows in radio will be detectable by instruments like the Square Kilometre Array beyond the horizon distance of NEMO (Dobie et al. Reference Dobie, Murphy, Kaplan, Hotokezaka, Bonilla Ataides, Mahony and Sadler2021), which will facilitate very-long baseline interferometry follow-up to determine jet structure and orientation of some short gamma-ray bursts from binary neutron star mergers (Duque, Daigne, & Mochkovitch Reference Duque, Daigne and Mochkovitch2019). Observatories like the Cherenkov Telescope Array will also provide all-sky high energy gamma-ray survey capabilities (Cherenkov Telescope Array Consortium et al. 2019).

With all of these detections across multiple wavebands and with high-frequency gravitational waves, the problem becomes: *how confidently can we associate any two measurements as coming from the same source?* Without this association, we cannot fulfil the rich promise that multimessenger gravitational-wave science brings. In this paper, we calculate predictions for the rate of binary neutron star mergers NEMO alone can detect and that can be confidently associated with their various electromagnetic counterparts. We show the gamma-ray prompt emission will be confidently associated with the gravitational-wave emission at an overall rate of
$2^{+10}_{-1} \mathrm{yr^{-1}}{}$
, where the uncertainties are 90% confidence intervals arising from uncertainty in current merger-rate estimates. Gravitational waves and kilonovae will be associated at
$13^{+23}_{-10}\mathrm{yr^{-1}}{}$
with broadband afterglows at
$4^{+18}_{-3}\mathrm{yr^{-1}}{}$
. Combined, this implies a coincident detection rate of
$14^{+25}_{-11}\mathrm{yr^{-1}}{}$
out to
$300\,\mathrm{Mpc}$
.

In this work, we discuss only neutron star binary mergers. The aLIGO/aVirgo network have already detected $2^{+2}_{-2}$ (100% confidence interval) neutron star-black hole binaries (Abbott et al. Reference Abbott2020; Abbott et al. Reference Abbott2021b; Abbott et al. Reference Abbott2021a; Abbott et al. Reference Abbott2021c); however, we note that each of these purported mergers have black hole masses and spins such that the neutron stars are not expected to be tidally disrupted, and therefore no electromagnetic counterparts are expected (e.g., Foucart Reference Foucart2020). The rates of neutron star-black hole binary mergers with tidal disruption are therefore unknown (e.g., Sarin et al. Reference Sarin, Lasky, Vivanco, Stevenson, Chattopadhyay, Smith and Thrane2022), and we do not speculate further on the expected science outcomes of those systems in this work. We further note that NEMO-like observatories are potentially good for detecting gravitational-wave signals from supernovae (e.g., Powell & Müller Reference Powell and Müller2020; Mezzacappa et al. Reference Mezzacappa2020; Pan et al. Reference Pan, Liebendörfer, Couch and Thielemann2021), however the rates of gravitational-wave only, and coincident gravitational-wave and electromagnetic detections, are highly uncertain. Finally, NEMO-like observatories may also be able to detect nearly monochromatic signals from millisecond pulsars; however, the gravitational-wave amplitude is highly uncertain (e.g., see Lasky Reference Lasky2015; Haskell & Schwenzer Reference Haskell and Schwenzer2021, and references therein; although see Woan et al. Reference Woan, Pitkin, Haskell, Jones and Lasky2018 for evidence of a minimum ellipticity).

This paper is set out as follows. In Section 2 we detail the salient aspects of future gravitational-wave and electromagnetic facilities for this study and describe the modelling of the gravitational-wave and electromagnetic signal. In Section 3 we provide the main results of the paper, using a Bayesian method for determining coincidences between multiple different observations which is presented in Appendix A. Section 3 focuses first on a hypothetical ‘jackpot’ event (Section 3.1) that is both on axis and nearby at 40 Mpc—the distance of the first gravitational-wave multimessenger event GW170817. In Sections 3.2 and 3.3 we provide more realistic scenarios by looking at populations of events with varying distances and inclination angles to the observer’s line of sight, respectively. In Section 4 we conclude by discussing the promising multimessenger science that can be learned from these future coincident detections.

## 2. Observations in the 2030s

### 2.1. Gravitational-wave observations

In this work, we assume only a single ground-based gravitational-wave interferometer is operating around the globe with NEMO-like sensitivity (Ackley et al. Reference Ackley2020). While we hope this situation never eventuates, we are mindful that there may be a data gap between the end-of-life stage of 2.5-generation interferometers such as A+ and Virgo+ (Abbott et al. Reference Abbott2018), and the beginning of third-generation observatories such as Einstein Telescope and Cosmic Explorer (Punturo et al. Reference Punturo2010; Abbott et al. Reference Abbott2017a).

We simulate a realistic population of binary neutron star gravitational-wave signals at fixed distances but with the inclination angle $\iota$ drawn from a uniform in $\cos\iota$ distribution. We inject all signals with masses drawn from a Gaussian mixture model motivated by analysis of the mass distribution of all neutron stars in our Galaxy (Alsing, Silva, & Berti Reference Alsing, Silva and Berti2018). This mass distribution mixture model is

where $\mathcal{N}(\mu,\,\sigma)$ denotes a normal distribution of mean $\mu$ and standard deviation $\sigma$ , $\mu_1=1.32\,{\rm M}_\odot$ and $\sigma_1 = 0.11$ , $\mu_2=1.80\,{\rm M}_\odot$ , $\sigma_2=0.21\,{\rm M}_\odot$ , and mixing fraction $\epsilon=0.35$ . We inject these signals isotropically over the sky.

We inject these signals into a NEMO detector located in Gingin, near Perth, Australia, to determine the fraction of events detectable at a given distance. We define a gravitational-wave detection as having an optimal matched-filter signal-to-noise ratio greater than eight. This threshold is somewhat arbitrarily defined—it is approximately the single-detector threshold currently used for the LIGO/Virgo Collaborations—however a more nuanced signal-to-noise threshold would require a detailed understanding of the instrument’s noise properties. Regardless, the precise numerical values for the fraction of detectable events is not important, but the overall approximate numbers are useful to understand for the context of this paper. We note that the number of events is $\propto 1/\textrm{SNR}^3$ .

In the top right panel of Figure 1, we plot the total number of detectable events in a single NEMO as a function of distance as the solid coloured band. We use the current best estimate for the binary neutron star merger event rate of $320^{+490}_{-240} \mathrm{Gpc^{-3} yr^{-1}}$ (Abbott et al. Reference Abbott2021b), where the uncertainties are the 90% confidence intervals, which are shown as the dashed curves in each panel of Figure 1. This figure indicates that NEMO will detect 95%, 65% and 5% of binary neutron star mergers at 40, 100, and $300\,\mathrm{Mpc}$ , respectively.

### 2.2. Electromagnetic observations

The 2030s will see the full promise of several current and proposed electromagnetic observatories with all-sky or survey capabilities. These observatories will enable routine and confident detections of kilonovae, short gamma-ray bursts, and their afterglows beyond the horizon distance of NEMO (e.g., Amati et al. Reference Amati2018; Ivezić et al. Reference Ivezić2019). These observatories may also enable detection of precursor electromagnetic emission from binary neutron star mergers as well as their post-merger central engines out to distances relevant to NEMO (see Ascenzi et al. Reference Ascenzi, Oganesyan, Branchesi and Ciolfi2021; Sarin & Lasky Reference Sarin and Lasky2021, and references therein).

In this work, we focus on whether the electromagnetic counterpart of binary neutron star mergers are detectable, and how confidently such signals can be associated with their gravitational-wave counterpart. We consider an electromagnetic signal to be detectable if it is above a pre-defined threshold for each instrument (described below). For gamma-ray burst prompt emission, their afterglows, and kilonovae, these thresholds are relatively well estimated for various electromagnetic telescopes (e.g., Coward et al. Reference Coward2011; Kanner et al. Reference Kanner, Camp, Racusin, Gehrels and White2012; Siellez, Boër, & Gendre Reference Siellez, Boër and Gendre2014). However, the thresholds, rates, and rate of false positives for precursor emission and post-merger emission from the central engine are unknown (e.g., Ascenzi et al. Reference Ascenzi, Oganesyan, Branchesi and Ciolfi2021; Sridhar et al. Reference Sridhar, Zrake, Metzger, Sironi and Giannios2021). We therefore only focus on the relatively well-observed prompt gamma-ray emission, broadband afterglow, and kilonova.

The prompt flash of gamma rays in short gamma-ray bursts is typically followed by a long-lasting, broadband afterglow. While the mechanism responsible for producing the prompt gamma-ray emission is ill understood (e.g., Kumar & Zhang Reference Kumar and Zhang2015), the broadband afterglow is relatively simple to understand. Most afterglows are believed to be from the interaction of a relativistic jet with the surrounding interstellar medium (Sari, Piran, & Narayan Reference Sari, Piran and Narayan1998). There are additional types of afterglows that we do not consider in detail here, such as those suggesting an additional contribution from a neutron star (e.g., Zhang & Mészáros Reference Zhang and Mészáros2001; Rowlinson et al. Reference Rowlinson, O’Brien, Metzger, Tanvir and Levan2013; Lü et al. Reference Lü, Zhang, Lei, Li and Lasky2015; Sarin, Lasky, & Ashton Reference Sarin, Lasky and Ashton2019).

The multi-wavelength observations of GRB170817A confirmed the long-held suspicion that gamma-ray burst jets are structured, that is, the energy and Lorentz factor have a non-trivial angular dependence (Troja et al. Reference Troja2017; Alexander et al. Reference Alexander2018; Mooley et al. Reference Mooley2018a, b). While the true jet structure is unknown, several phenomenological models have been developed and fit to gamma-ray burst observations. One such commonly used model is the Gaussian structured jet (e.g., Lamb, Mandel, & Resmi Reference Lamb, Mandel and Resmi2018; Ryan et al. Reference Ryan, van Eerten, Piro and Troja2019)

Here, $E_{0}$ is the isotropic equivalent kinetic energy in the afterglow, $\theta$ is the angle from the jet axis, and $\theta_{c}$ is the half-opening angle of the ultra-relativistic core or, equivalently, the opening angle of the jet. None of these parameters are well constrained empirically, but multi-wavelength observations of GRB170817A and other gamma-ray bursts offer some clues.

The ultra-relativistic jet emanating from a neutron star merger interacts with the surrounding interstellar medium accelerating a fraction of electrons, $\xi_{n}$ , with some fraction of the total energy of the jet, $\epsilon_{e}$ , and some fraction of the energy in the magnetic field, $\epsilon_b$ . The radiation produced by these electrons is responsible for the observed broadband afterglow. This interaction is relatively easy to model, but computationally expensive. Fortunately, semi-analytic approximations make modelling these broadband afterglows tractable. We use afterglowpy (Ryan et al. Reference Ryan, van Eerten, Piro and Troja2019) with the above jet structure profile to model our afterglow lightcurves.

The prompt emission is even less understood than the broadband afterglow. There is no robust generative model that allows predictions or modelling of the prompt emission energetics. Observations of the off-axis prompt emission from GRB170817A (e.g., Matsumoto, Nakar, & Piran Reference Matsumoto, Nakar and Piran2019; Ioka & Nakamura Reference Ioka and Nakamura2019) further highlighted our lack of understanding of the angular structure (e.g., Howell et al. Reference Howell, Ackley, Rowlinson and Coward2019; Ioka & Nakamura Reference Ioka and Nakamura2019). By the late 2020 s/early 2030 s, we will undoubtedly have more insight into prompt energetics and angular structure. We therefore plough forward in our endeavour, acknowledging these known unknowns.

In light of GRB170817A, one common approach is to compute the peak flux at a given observer viewing angle using a similar structured jet profile used to describe the afterglow (e.g., Howell et al. Reference Howell, Ackley, Rowlinson and Coward2019; Salafia et al. Reference Salafia, Ghirlanda, Ascenzi and Ghisellini2019; Amati et al. Reference Amati2018). Such analyses make two assumptions: (1) the jet structure in the afterglow and prompt emission phase are the same, and (2) the jet is transparent to the gamma rays at all viewing angles. Both these assumptions could be incorrect. Alternatively, one may use a range of estimated gamma-ray efficiencies for prompt emission inside the opening angle of the jet (e.g., Fong et al. Reference Fong, Berger, Margutti and Zauderer2015), and model the emission outside the jet as a cocoon-shock breakout (e.g., Gottlieb et al. Reference Gottlieb, Nakar, Piran and Hotokezaka2018). This gives a qualitatively similar answer for prompt emission energetics outside the jet opening angle as the approach of Howell et al. (Reference Howell, Ackley, Rowlinson and Coward2019), Salafia et al. (Reference Salafia, Ghirlanda, Ascenzi and Ghisellini2019). We therefore adopt the latter approach. By the 2030s, we may better understand the energetics of off-axis emission if we observe more multimessenger events like GW170817 (Biscoveanu, Thrane, & Vitale Reference Biscoveanu, Thrane and Vitale2020).

All of the above ignores complications surrounding jet launching itself. Numerical simulations suggest that no neutron star remnant can launch a jet that can produce a short gamma-ray burst (Murguia-Berthier et al. Reference Murguia-Berthier2017). However, this is in tension with x-ray afterglow observations of many short gamma-ray bursts (e.g., Rowlinson et al. Reference Rowlinson, O’Brien, Metzger, Tanvir and Levan2013). In this work, we stick with numerical results, placing a constraint that a gamma-ray burst and subsequently the afterglow emission is only possible when the remnant mass is above the threshold to form a black hole, that is, $M_{\rm {rem}} \gtrsim 1.2 M_{\rm {TOV}}$ , where $M_{\rm {TOV}}$ is the maximum non-rotating neutron star mass; a property of the nuclear equation of state. For each merger, we estimate the remnant mass $M_{\rm {rem}}$ using the relationship described in Gao Zhang, & Lü (Reference Gao, Zhang and Lü2016). We then calculate whether that merger would produce a gamma-ray burst and afterglow marginalising over piecewise polytropic equations of state with $M_{\rm {TOV}}$ between $2.2{-}2.4\,{\rm M}_{\odot}$ motivated by current estimates (e.g., Ai, Gao, & Zhang Reference Ai, Gao and Zhang2020).

Beyond the prompt gamma-ray emission and broadband afterglow, r-process nucleosynthesis from the neutron-rich ejecta of a binary neutron star merger is expected to power a week-long thermal transient, that is, a kilonova. This was likely first seen in GRB130603B (Tanvir et al. Reference Tanvir, Levan, Fruchter, Hjorth, Hounsell, Wiersema and Tunnicliffe2013; Berger, Fong, & Chornock Reference Berger, Fong and Chornock2013), but has since been studied in considerably more detail through multi-wavelength observations of AT2017gfo, the kilonova counterpart to GW170817 (e.g., Evans et al. Reference Evans2017; Smartt et al. Reference Smartt2017; Villar et al. Reference Villar2017). There are more kilonova models in the literature than there are observations. Many models do a reasonable job of explaining the observations of AT2017gfo (see Breschi et al. Reference Breschi, Perego, Bernuzzi, Del Pozzo, Nedora, Radice and Vescovi2021, for a quantitative model comparison). Some models are predominantly phenomenological while others attempt to connect the progenitor’s parameters to the resulting lightcurve. We use the latter class of models. In particular, we use the kilonova model from Metzger (Reference Metzger2017) relating it to the intrinsic binary parameters through fits of the ejecta properties derived from numerical-relativity simulations (Dietrich & Ujevic Reference Dietrich and Ujevic2017) with a piecewise polytropic equation of state (using the software presented in Hernandez Vivanco et al. Reference Hernandez Vivanco, Smith, Thrane and Lasky2020) and a baryonic mass relation (e.g., Coughlin et al. Reference Coughlin, Dietrich, Kawaguchi, Smartt, Stubbs and Ujevic2017; Gao et al. Reference Gao, Ai, Cao, Zhang, Zhu, Li, Zhang and Bauswein2020).

The afterglow emission is broadband and will be detectable with various electromagnetic telescopes and surveys at different wavelengths. However, for the purposes of detectability here, we focus specifically on the optical afterglow with the Vera Rubin Observatory. The
$5\sigma$
magnitude threshold for single images in Vera Rubin in *r*-band is
$24.5$
with an equivalent flux density of
$4.79\times10^{-4}\mathrm{mJy}$
, which requires a
$30\mathrm{s}$
integration time (Ivezić et al. Reference Ivezić2019). Therefore, we consider any optical afterglow signal with an *r*-band peak flux greater than the single-image threshold to be detectable. We use the same threshold for a kilonova signal. For gamma-ray prompt emission, we consider a gamma-ray burst to be detectable if the peak gamma-ray flux in the 30–150 keV regime is above
$3 \times 10^{-8}\mathrm{erg\ cm^{-2}\ s^{-1}}$
. This is comparable to the expected sensitivity of Theseus (Amati et al. Reference Amati2018). We could also use Theseus to study the detectability of the x-ray afterglow, however the x-ray instrument proposed for Theseus does not have the wide-field survey capability of Vera Rubin and will not serendipitously detect as many afterglows.

We simulate the expected electromagnetic emission using the models described above using the population of binary neutron star mergers described in Section 2.1. The prior ranges of the various parameters in the relevant models are motivated by analyses of GRB170817A (Ryan et al. Reference Ryan, van Eerten, Piro and Troja2019) and AT2017gfo (Coughlin et al. Reference Coughlin, Dietrich, Kawaguchi, Smartt, Stubbs and Ujevic2017) and summarised in Table B.1 in Appendix B. We calculate the number of these events that produce a detectable electromagnetic counterpart at a given distance, which is shown in the top left panel of Figure 1 for kilonovae, bottom left panel for prompt emission, and bottom right panel for the broadband afterglows. The bands correspond to the uncertainties in detectable fraction due to the uncertain model parameters, nuclear equation of state, and binary neutron star merger rate. The absolute merger rate is shown in each panel as the dashed curves.

## 3. Associating gravitational-wave and electromagnetic signals

In the following, we calculate coincident detection rates under different scenarios. We generalise the Bayesian calculation of coincidences between any two data sets first presented in Ashton et al. (Reference Ashton2018) to multiple data sets. We provide details of the formalism in Appendix A.

### 3.1. A GW170817-like event

We first walk through a fiducial example of a coincident multimessenger event akin to GW170817. We include a single NEMO-like gravitational-wave detector, with electromagnetic observations that provide prompt gamma-ray burst identification, broadband afterglow and kilonovae observations, and host-galaxy identification. This is truly a jackpot event; in subsequent sections we consider scenarios with only a subset of these observations.

We simulate an event coined GW301116, which is a face-on (inclination angle $\iota=0$ ) binary neutron star merger at luminosity distance $D_L=40\,\mathrm{Mpc}$ with component masses $m_1=1.5\mathrm{M_{\odot}}$ and $m_2=1.4\mathrm{M_{\odot}}$ , dimensionless tidal deformabilities $\Lambda_1=400$ and $\Lambda_2=450$ , dimensionless in-plane spins $\chi_{1,2}=0.02$ , at ${\rm RA}=1.375$ and ${\rm Dec}=0.55$ . We inject this signal into a single NEMO detector using the IMRPhenomPv2_NRTidal (Dietrich et al. Reference Dietrich2019) waveform, and recover the system parameters using the same waveform. We utilise Parallel Bilby (Ashton et al. Reference Ashton2019; Romero-Shaw et al. Reference Romero-Shaw2020; Smith et al. Reference Smith, Ashton, Vajpeyi and Talbot2020) with the dynesty nested sampler (Speagle Reference Speagle2020). The sky map showing the 50% and 90% posterior credible intervals is shown in Figure 2, and the one-dimensional marginalised posteriors for the time of coalescence and luminosity distance are shown in green in Figures 3 and 4, respectively.

The event GW301116 being on axis and at
$40\,\mathrm{Mpc}$
implies the prompt gamma-ray signal would be observed by any telescope not occulted by the Earth (e.g., Howell et al. Reference Howell, Ackley, Rowlinson and Coward2019). We consider capabilities of future gamma-ray telescopes SVOM (Wei et al. Reference Wei2016), and Theseus (Amati et al. Reference Amati2018), which are all capable of localising the gamma-ray burst to
${\lesssim}15$
arcmin precision (Amati et al. Reference Amati2018). We thus model the sky-localisation posterior as a two-dimensional distribution mimicking current gamma-ray detectors such as *Swift* using the GWCelery package (Singer et al. Reference Singer2020). This sky map is shown as the inset in Figure 2.

A prompt gamma-ray detection also provides a strong constraint on the merger time, even accounting for the uncertainty in the delay between the merger and prompt emission signal. By the late 2020 s/early 2030 s, we estimate that the uncertainty on the merger time from the identification of a prompt gamma-ray signal will be $\mathcal{O}(2 \mathrm{s})$ . We therefore draw samples for the time of coalescence from a uniform distribution spanning 2 s to derive a posterior distribution for $t_c$ , which is shown as the purple curve in Figure 3.

An on-axis binary neutron star merger at $40\,\mathrm{Mpc}$ will also be detectable as a kilonova with observatories such as Vera Rubin (Ivezić et al. Reference Ivezić2019) and the Zwicky Transient Facility (Bellm et al. Reference Bellm2019), and as an afterglow with various broadband electromagnetic telescopes. These observations will enable the identification of the host galaxy and redshift, as well as an indirect measurement of the merger time. The precision for sky localisation of such an event will therefore be sub-arcsecond (Wei et al. Reference Wei2016; Ivezić et al. Reference Ivezić2019; Amati et al. Reference Amati2018). This is such a small localisation area by comparison to gravitational waves that it is not visible in Figure 2. We estimate that by the late 2020 s/early 2030 s, the error on the redshift will be $\Delta z\sim\mathcal{O}({10^{-4}})$ . We note that the error on the redshift for the host galaxy of GW170817, NGC4993, is $\mathcal{O}({10^{-3}})$ (Hjorth et al. Reference Hjorth2017; Abbott et al. Reference Abbott2017b), and improvements in electromagnetic surveys and better understanding of systematic uncertainties like peculiar velocities will improve this limit (Howlett & Davis Reference Howlett and Davis2020). We further note that for on-axis events it may be difficult to disentangle the optical afterglow from the kilonova but we ignore this complication. We assume a Gaussian posterior of width $\sigma_z=10^{-4}$ for the redshift, which we convert to luminosity distance assuming Planck 2018 cosmology (Aghanim et al. Reference Aghanim2020). The luminosity distance posterior is shown in Figure 4 as the red curve.

Using these various posterior distributions, we calculate the odds using Equations (A1) and (A6) as derived in Appendix A. For a jackpot event like GW301116, the Bayes factor is $4.3\times10^{11}$ . To estimate the prior odds, we consider how many events with a short gamma-ray burst, afterglow, and kilonova we expect in the four-dimensional combined uncertainty region (i.e., $t_c$ , $D_L$ , RA, and Dec.) given the binary neutron star merger rate. For a merger rate of $320^{+490}_{-240}\mathrm{Gpc^{-3} yr^{-1}}$ , and uncertainties on the coalescence time, luminosity distance, and sky position, we expect one merger every $\approx 12$ yr in the four-dimensional uncertainty region, implying a prior odds of $9.8\times 10^{11}$ . This implicitly assumes that every binary neutron star merger at $D_L = 40\,\mathrm{Mpc}$ will produce a detectable prompt gamma-ray emission, afterglow or kilonova, which is true for an on-axis merger such as GW301116, but may not be for other mergers. We explore the effect of the inclination in Section 3.3. For prior odds of $9.8\times 10^{11}$ , the odds is therefore $4.2\times 10^{23}$ , indicating it is $4.2\times 10^{23}$ times more likely that the observations come from the same astrophysical event than they do not. We note that in more realistic scenarios, the prior odds will be dictated by the rate of false positives for a given observatory, the true astrophysical rate, and the potential for event misidentification. We elaborate on these more realistic prior odds calculations in the next section. We note that the above scenario of observing an on-axis event at $D_L = 40\,\mathrm{Mpc}$ is going to be exceedingly rare. In the next sections we consider more common types ofevents.

### 3.2. Coincident detections as a function of distance

The expected detection rate of GW301116-like events that include coincident gravitational-wave, prompt emission, afterglow, and kilonovae observations as well as host-galaxy identification is far less than one per year. More realistically, we expect most multimessenger gravitational-wave events to be observed only with one of the aforementioned emissions and at varying distances. To wit, we repeat the calculation of Section 3.1 but with a number of events at different distances from $D_L=40$ to $400\,\mathrm{Mpc}$ . All other parameters and the calculation remain the same as detailed in the previous section.

The afterglow and kilonova provide only an indirect measure on the time of coalescence $t_{c}$ . The peak timescales for an on-axis merger are dictated by the deceleration of the relativistic jet for the afterglow, and the timescale for the expansion of the ejecta to balance the diffusion timescale for the kilonova. These peak timescales are $\mathcal{O}(\mathrm{1-2}\mathrm{s})$ and $\mathcal{O}(1\mathrm{d})$ , respectively. By the late 2020 s/early 2030 s, we predict that advancements in cadence of optical telescopes and theoretical modelling of these phenomena will enable identification of afterglows and kilonovae within approximately 3 and $6\mathrm{h}$ from the merger time, respectively. We emphasize that this is speculation, and realistic values will depend on detector design and survey choices too difficult to predict. Depending on the scenario used, we take this time uncertainty as the width of the posterior on the time of coalescence $t_{c}$ .

The prior odds depend on the scenario being considered. For a jackpot event like GW170817, that is, an event with an observation of a short gamma-ray burst, multi-wavelength afterglow and kilonova, the prior odds is the inverse of the number of events expected in the given four-dimensional combined uncertainty region given the binary neutron star merger rate itself. Given the low merger rate, this implies *a priori* that if we see a confident short gamma-ray burst, kilonova, gravitational-wave signal, and afterglow in the same four-dimensional combined uncertainty region, it is overwhelmingly likely they are coincident (see previous section for the calculation). The same is not true for other scenarios. For example, the prior odds on a gamma-ray burst in coincidence with a gravitational-wave signal is dictated by the rate of triggers (astrophysical or otherwise) from gamma-ray telescopes. For Theseus, this is predicted to be
${1{-}8\,\mathrm{week}^{-1}}$
(Amati et al. Reference Amati2018), which dwarfs the astrophysical rate of detectable short gamma-ray bursts themselves.

Optical telescopes will find many potential candidates for both kilonovae and afterglows. Even if these candidates are astrophysical, they could easily be mistaken as supernovae or afterglows of long gamma-ray bursts for example. A realistic rate of triggers in an observatory such as Vera Rubin is not known and will depend on the survey strategy (Setzer et al. Reference Setzer, Biswas, Peiris, Rosswog, Korobkin and Wollaeger2019; Cowperthwaite et al. Reference Cowperthwaite, Villar, Scolnic and Berger2019). We therefore take estimates from current wide-field optical telescopes, using particularly the Gravitational-Wave Optical Transient Observer (GOTO; Dyer et al. Reference Dyer2020), which has recently announced results from a survey that searched for the optical afterglows of gamma-ray bursts, and had a trigger rate of ${\approx} 1{-}20\, {\mathrm{week}^{-1}}$ (Mong et al. Reference Mong2021). This is significantly larger than the astrophysical rate of kilonova-like or afterglow-like optical transients, and can be used to determine a conservative estimate of the prior odds for scenarios involving kilonovae and afterglows. We use these trigger estimates, alongside the rate of binary neutron star mergers and potential contaminants such as long gamma-ray bursts and their afterglows (Ghirlanda et al. Reference Ghirlanda2021) and supernovae and fast-blue optical transients (Ho et al. Reference Ho2020), to calculate the prior odds.

We note that alongside binary neutron star mergers, some fraction of neutron star-black hole mergers are expected to produce a kilonova, short gamma-ray burst and afterglow. In fact, there may already be evidence of this contamination in the observed sample of short gamma-ray bursts (Siellez et al. Reference Siellez, Boer, Gendre and Regimbau2016; Gompertz, Levan, & Tanvir Reference Gompertz, Levan and Tanvir2020; Hamburg et al. Reference Hamburg2020). However, given the astrophysical rate of neutron star-black hole mergers is small (Abbott et al. Reference Abbott2021e), ignoring contamination from neutron star-black hole mergers does not affect our result.

The green ‘GW170817-like’ curve in Figure 5 shows the Bayes factor (top panel) and odds (bottom panel) as a function of distance for events where we detect the gravitational-wave signal as well as the prompt emission, broadband afterglow and kilonovae, as well as identify the host galaxy. At all distances, the odds is overwhelmingly in favour of coincidence, going from $8.3\times10^{27}$ at $D_L=40\,\mathrm{Mpc}$ , to $1.7\times10^{24}$ at $D_L=400\,\mathrm{Mpc}$ , largely due to the strong a priori odds for a coincident event in that scenario. The purple ‘GRB + Afterglow’ curve shows the same, however assuming that the kilonova was not observed. The red ‘Afterglow + kilonova’ curve assumes only the prompt emission is missed, while the turquoise ‘Kilonova only’ curve assumes only the gravitational-wave and kilonova signal are detected (i.e., the prompt and afterglow are missed).

The horizontal grey band in the bottom panel shows an odds of $\log_{10}\mathcal{O}_{C/R}=6$ , corresponding approximately to a five-sigma threshold for the association of the gravitational-wave signal with the electromagnetic counterpart. Using this relatively arbitrary threshold, we see that all face-on binary neutron star mergers will have confident coincidences between the gravitational-wave and electromagnetic counterpart, with the exception of signals that are only detected through their kilonova at distance $D_L\gtrsim250\,\mathrm{Mpc}$ . We note that this threshold is somewhat arbitrary, and deriving a realistic threshold requires understanding the distribution of $\log_{10}\mathcal{O}_{C/R}$ .

### 3.3. The effect of inclination angle

The previous two sections assume the gravitational-wave and electromagnetic signal(s) are observed, and quantifies the confidence we could get that those signals originate from the same source. Not all electromagnetic signals will be detectable, in part due to the inclination angle of the source implying, for example, the gamma-ray burst jet is pointed away from Earth. To understand this, we look at the number of detectable gamma-ray burst, afterglow, kilonovae, and gravitational-wave sources as a function of both distance and inclination angle.

At fixed values of distance, we create a population of neutron star mergers that are isotropically distributed in inclination angle (i.e., uniform in $\cos\iota$ ) with the same intrinsic parameters as described in Section 2. We assume the same jet structure and kilonova models as described in the previous sections. For each event, we calculate whether the signal is detectable with each messenger and, if detectable, whether any two or more messengers can be confidently associated. We marginalise over systematic uncertainties due to the unknown equation of state, jet structure and prior odds discussed in Sections 2.2 and 3, which directly affect the expected kilonova, gamma-ray burst signature and odds calculation, respectively.

In Figure 6, we show the total coincident detection rate per year as a function of luminosity distance. As an example, the top left panel shows in the turquoise-coloured band the total detection rate for which the gravitational-wave signal and kilonova can be confidently identified (above a $5{-}\sigma$ threshold) as coming from the same source, marginalising over all previously mentioned uncertainties, as well as the uncertainty in the overall merger rate. This coloured band can be contrasted with the total binary neutron star merger rate, where the 90% confidence intervals are shown as the dashed curves in each panel.

One can inspect the curves in Figure 6 to derive an overall merger fraction $\mathcal{F}$ of mergers that will be detected with gravitational waves and various electromagnetic messengers, as well as the overall rate $\mathcal{R}$ of detections. These numbers are provided in Table 1. These numbers show $14^{+25}_{-11}\mathrm{yr^{-1}}{}$ (or ${\approx}40\%$ ) of binary neutron star mergers per year will have confident associations between their independently identified electromagnetic signatures and gravitational-wave signal with only a NEMO detector.

We note that the above numbers, in particular the fraction of binary neutron star mergers with any confident association, should be taken with some caution. These numbers are sensitive to three assumptions: (1) The true distribution of $\log_{10}\mathcal{O}_{C/R}$ is distributed in a way such that our threshold for a confident association is too weak. (2) The modelling assumptions and priors for the electromagnetic models are too optimistic. (3) The fraction $\mathcal{F}$ represents the fraction of mergers detectable with both an electromagnetic telescope and with NEMO that can also be confidently associated with one another. However, whether a given electromagnetic telescope can actually detect a counterpart (even if it is above the threshold for detection) will depend on the survey strategy and specifics of the detector. These strategies and detector specifics are too difficult to predict for observatories that do not yet exist, and a calculation of the true distribution of $\log_{10}\mathcal{O}_{C/R}$ is not computationally feasible as it requires simulating the complete population multiple times to create a distribution of $\log_{10}\mathcal{O}_{C/R}$ for each binary neutron star merger. However, both of these factors will likely be better understood in the near future, accompanied by a better understanding of electromagnetic models.

## 4. Conclusions

The future of gravitational-wave astronomy relies on sequential technological and infrastructure upgrades, as well as the development of objectively expensive new facilities (albeit far less than some space-based infrared telescopes). One plausible scenario sees a period of time for which a designated kilohertz gravitational-wave detector like the proposed NEMO instrument (Ackley et al. Reference Ackley2020) is the only operational ground-based observatory. In this work, we show that such an instrument has significant scientific utility on its own due to the multimessenger capabilities expected in the same era. For example, in operation with the Vera Rubin observatory (Ivezić et al. Reference Ivezić2019) and Theseus (Amati et al. Reference Amati2018), we show that approximately 40% of binary neutron star mergers out to $300\,\mathrm{Mpc}$ can be confidently identified as multimessenger coincident observations. Given current estimates for the local merger rate of neutron stars (Abbott et al. Reference Abbott2021b), this implies a total number of $14^{+25}_{-11}\mathrm{yr^{-1}}{}$ multimessenger detections per year.

Coincident multimessenger observations of binary neutron star mergers will enable several precise tests into the nature of nuclear matter. The ${\gtrsim}500\mathrm{Hz}$ gravitational-wave sensitivity enables simultaneous measurements of the neutron star’s masses and tidal deformability, which can be further constrained when coupled with kilonovae observations that inform the ejecta mass, and hence the progenitor’s mass ratio and complementary information about the equation of state (e.g., Metzger Reference Metzger2017). For nearby mergers (see Ackley et al. Reference Ackley2020) this will be combined with measurements of gravitational waves from the hot post-merger remnant, which can further inform the equation of state (e.g., Bauswein et al. Reference Bauswein2019), including the potential to discern temperature-dependent phase transitions (Bauswein, Stergioulas, & Janka Reference Bauswein, Stergioulas and Janka2016). Combining gravitational-wave information from the inspiral and post-merger remnant with constraints on the jet-launching timescale from timing of the prompt emission (Ren et al. Reference Ren, Lin, Zhang, Wang, Li, Wang and Liang2020; Beniamini et al. Reference Beniamini, Duran, Petropoulou and Giannios2020b) can constrain the maximum mass of hot and cold neutron stars (e.g., Chatziioannou Reference Chatziioannou2020; Sarin & Lasky Reference Sarin and Lasky2021); also a key indicator of the nuclear equation of state. This information can also be determined from complementary observations of, for example, the colour and energetics of the kilonova (Margalit & Metzger Reference Margalit and Metzger2019) and the duration of x-ray plateaus from surviving supramassive neutron stars (Rowlinson et al. Reference Rowlinson, O’Brien, Metzger, Tanvir and Levan2013; Ravi & Lasky Reference Ravi and Lasky2014; Sarin, Lasky, & Ashton Reference Sarin, Lasky and Ashton2020).

Associating gamma-ray bursts and kilonovae with gravitational-wave signals will also enrich our understanding of the physics driving these enigmatic transients. In particular, constraints on the delay time between merger and prompt emission will deepen our understanding of ultra-relativistic jet launching and propagation (Ren et al. Reference Ren, Lin, Zhang, Wang, Li, Wang and Liang2020). This will also shed insight into the prompt emission mechanism. Detection of post-merger gravitational waves combined with prompt emission observations will provide a smoking-gun observation into the contentious issue of whether a black hole central engine is required to launch an ultra-relativistic jet capable of generating short gamma-ray burst emission. This will also provide an opportunity to understand the impact of the remnant on r-process enrichment (e.g., Perego et al. Reference Perego, Rosswog, Cabezón, Korobkin, Käppeli, Arcones and Liebendörfer2014; Metzger Reference Metzger2017; Bernuzzi Reference Bernuzzi2020). Coincidentmultimessenger observations will also improve our understanding of the jet structure and beaming of gamma-ray burst jets (e.g., Beniamini et al. Reference Beniamini, Duque, Daigne and Mochkovitch2020a; Biscoveanu et al. Reference Biscoveanu, Thrane and Vitale2020).

The capabilities of survey telescopes in the 2030s will enable host-galaxy identification of neutron star mergers, and hence precise redshift measurements (e.g., Howlett, Staveley-Smith, & Blake Reference Howlett, Staveley-Smith and Blake2017). Required peculiar-velocity measurements will be enabled by next-generation spectroscopic surveys (e.g., Palmese et al. Reference Palmese2019; Howlett & Davis Reference Howlett and Davis2020), and possibly percent-level measurements of the Hubble constant with single events (Coughlin et al. Reference Coughlin2020; Calderón Bustillo et al. Reference Calderón Bustillo, Leong, Dietrich and Lasky2021).

Dedicated kilohertz gravitational-wave observatories have the power to produce valuable science *without* the requirement of living in a global, heterogeneous network of other gravitational-wave observatories, provided the capabilities of electromagnetic telescopes do not get worse. In this paper, we show this is true only considering the guaranteed science case of binary neutron star mergers. Neutron star-black hole binaries, millisecond pulsars, and supernovae are all targets for NEMO-like observatories that would provide rich multimessenger information about the physics and astrophysics of these extreme objects. However, the unknown neutron star-black hole tidal disruption rate, the unknown ellipticity of rapidly rotating isolated neutron stars, and the unknown amplitude of gravitational waves from supernovae prohibits us from making concrete statements on the likelihood of coincident detections. Regardless, we know that the science outcomes from such events would likely be bountiful.

## Acknowledgements

We are grateful to David Ottaway and Eric Howell for their thoughtful comments on the manuscript. This work was supported through Australian Research Council (ARC) Centre of Excellence CE170100004, ARC Future Fellowship FT160100112, and ARC Discovery Project DP180103155. This work was performed on the OzSTAR national facility at Swinburne University of Technology. The OzSTAR programme receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government.

## A. Nitty-gritty formalism

We wish to determine the probability that observed electromagnetic and gravitational-wave signals are coming from the same source. This involves two questions: (1) are the two events coincident in space and time, and (2) is it physically possible/probable that the signals were generated in the same event. In a Bayesian sense, the former question is quantified by the Bayes factor, while the latter is quantified by the prior odds. Together, they enable us to calculate the odds that any two or more pieces of data were generated by the same source.

To determine whether two or more independent observations are coincident (i.e., come from the same physical source), we generalise the Bayesian formalism of Ashton et al. (Reference Ashton2018) from two to *N* coincident data sets. We calculate the *odds*
$\mathcal{O}_{C/R}$
(Ashton et al. Reference Ashton2018, Reference Ashton, Ackley, MagaÑa Hernandez and Piotrzkowski2020), where *C* refers to the hypothesis that the two independent detections share a common origin and *R* refers to the hypothesis that they are random. This is given by

where

is the Bayes factor,
${P\left(C \mid d_{1}, d_{2}, \ldots,d_{N}, I\right)}$
is the probability of a common hypothesis conditional on data sets
$d_{1}, d_{2}, \ldots,d_{N}$
, while the denominator refers to the random hypothesis conditioned on the same data. The second term in Equation (A1),
$\pi_{C / R}$
is the *a priori* probability of a coincident event, which conservatively can be estimated as one over the number of expected events in the given overlapping region.

For a given gravitational-wave and electromagnetic transient to share a common origin, they must have occurred at the same location in the sky, at the same distance, and within some coincident time-span as specified by theoretical models. This implies they must have the same sky-location parameters $\Omega$ , same luminosity distance $D_L$ , and same inferred time of coalescence $t_{c}$ . In the case of observations with a common source, some other parameters will also be the same (e.g., source inclination). However, given theoretical limitations in current electromagnetic models, particularly in connecting to the progenitor parameters measured by the gravitational-wave signal, we ignore these common parameters. The confidence gained by including these additional parameters is proportional to the information gained on these parameters relative to the prior.

Following the framework from Ashton et al. (Reference Ashton2018), the numerator of the Bayes factor is,

Here, $\theta$ is the set of parameters common to all data sets, $\pi(\theta \mid C)$ is the common hypothesis prior on $\theta$ , and

Substituting Equation (A4) into Equation (A3) gives,

where
$\mathcal{I}_\theta$
is the posterior overlap integral for the parameter
$\theta$
generalised for *N* different data sets

We now move from the generalised abstract formalism described above to our specific study. For example, we consider the scenario of three different data sets: a gravitational-wave signal, kilonova, and afterglow. The overlap integral can be written as

where $d_{\mathrm{gw}}$ , $d_{\mathrm{kn}}$ , and $d_{\mathrm{ag}}$ are the gravitational-wave, kilonova and afterglow data, respectively. The odds given our common parameters $\theta=\{D_{L},\,\Omega,\,t_{c}\}$ is

In main body of the paper, we discuss what constraints a detection of kilonovae, prompt gamma-ray and afterglow emission and a gravitational-wave signal detected by a kHz-band gravitational-wave observatory can provide on these joint parameters. We use the framework described above to determine how confidently we can associate a given electromagnetic counterpart with the gravitational-wave signal. We note that since these transients are independently identified, we do not rely on low-latency estimates to calculate the odds and associate events.

## B. Priors

Here we list the priors used to model the prompt emission, afterglow, and kilonova signal from a binary neutron star merger. These priors are motivated by analysis of the afterglow of GRB170817A (Ryan et al. Reference Ryan, van Eerten, Piro and Troja2019) and the kilonova AT2017gfo (Coughlin et al. Reference Coughlin, Dietrich, Kawaguchi, Smartt, Stubbs and Ujevic2017). The priors are displayed inTable B.1.