Observing cosmological binary mergers with next generation neutrino and gravitational wave detectors

TL;DR

This paper explores the potential for detecting thermal neutrinos from binary mergers using next-generation neutrino and gravitational wave detectors, enabling multi-messenger observations up to redshift 2 and providing insights into merger dynamics.

Contribution

It proposes a multi-messenger observational strategy combining neutrino and gravitational wave data to study binary mergers and estimates expected neutrino event rates for future detectors.

Findings

Expected neutrino events range from 0.1 to 10 for typical merger rates.

Extreme merger scenarios could produce up to 100 neutrino events.

Detection up to redshift z~0.5 is feasible with current planned detector sensitivities.

Abstract

We discuss the potential of detecting thermal neutrinos from matter-rich binary mergers, via a decades-long multi-messenger campaign involving a Mt-scale water Cherenkov neutrino detector and one or more next generation gravitational wave detectors, capable of observing mergers up to redshift $z\sim 2$. The search of neutrinos in time-coincidence with gravitational wave detections will allow to identify single neutrinos from individual mergers above the background, and to study their distributions in energy, redshift and type (double neutron-star or neutron-star-black hole merger) of the candidate sources. We find that, for merger rates consistent with current LIGO-Virgo constraints, and for a $100~{\rm Mt\cdot yr}$ exposure, between ${\mathcal O(10^{-1})}$ and ${\mathcal O(10)}$ neutrino events are expected. For extreme cases of mergers with more than $10^{52}$ ergs emitted in $\mathrel{{\bar \nu}_e}$, the number of events can be as large as $\sim 100$, with sensitivity to mergers up to redshift $z\sim 0.5$ or so. Such scenarios can already be tested with a $10~{\rm Mt\cdot yr}$ exposure, resulting in constraints on the post-merger evolution of the systems being considered.

Tables

Table 1. Table 1: Parameters and neutrino emission properties for the merger models considered in this study. The columns from left to right contain: (1) time duration, T 99 , subscript 𝑇 99 T_{99}, within which 99 % percent 99 99\% of the total energy in ν ¯ e subscript ¯ 𝜈 𝑒 \mathrel{{\bar{\nu}}_{e}} is emitted; (2) total energy emitted in ν ¯ e subscript ¯ 𝜈 𝑒 \mathrel{{\bar{\nu}}_{e}} ; (3) average energy of the ν ¯ e subscript ¯ 𝜈 𝑒 \mathrel{{\bar{\nu}}_{e}} (time-integrated) flux; (4) type of merger (BNS or NSBH); (5) type of central remnant (CR) which could be a hypermassive neutron star (HMNS) or a black-hole (BH); (6) the central remnant mass (CRM) in units of solar mass, M ⊙ subscript 𝑀 direct-product M_{\odot} ; (7) the disk mass (DM; note that for model M of Sekiguchi et al., this quantity is not available); (8) the model name and (9) corresponding reference (first author only; see text for full reference).

$T_{99}\left(s\right)$	${ℰ}_{{\overline{\nu }}_e}\left({10}^{51}erg\right)$	$⟨E_{{\overline{\nu }}_e}⟩\left(MeV\right)$	type	CR	CRM ($M_{\odot }$)	DM ($M_{\odot }$)	Model	Ref
0.58	4.4	18	BNS	HMNS	3	0.03	Hinf	J.Lippuner(2017)
0.40	2.0	16.5	BNS	BH	3	0.03	B090
0.30	1.8	15.4	NSBH	BH	8.1	0.1	BF15
0.10	1.0	17.8	BNS	BH	3	0.03	M3A8m03a5	O.Just(2015)
0.27	11.2	16	NSBH	BH	6	0.3	M6A8m3a5
0.99	14	10	BNS	HMNS	2.7	0.2	DD2-1351350-On-H	S.Fujibayashi(2017)
0.58	40	20	BNS	HMNS	3		M	Y.Sekiguchi(2011)
0.08	19.8	24	NSBH	BH	6.6	0.49	A5	H.T.Janka (1999)
0.16	23.2	28	NSBH	BH	11.6	0.47	B10

Equations

$$\frac{dV}{dz}$$

$$D_c$$

$$N(z_{max})={\int }_0^{z_{max}}\frac{R(z)dV}{(1+z)dz}dz.$$

$$F_{{\bar{\nu }}_e}(E)=\frac{{\mathcal{E}}_{{\bar{\nu }}_e}}{⟨E_{{\bar{\nu }}_e}⟩}\frac{2}{3T^3\zeta (3)}\frac{E^2}{e^{E/T}+1},$$

$${\Phi }_{{\bar{\nu }}_e}(E)=f_{osc}\frac{c}{H_0}{\int }_0^{z_{max}}R(z)F_{{\bar{\nu }}_e}(E^{\prime })\frac{dz}{\sqrt{{\Omega }_m(1+z{)}^3+{\Omega }_{\Lambda }}},$$

$$N_s=\sum _{i=1}^N{n}_s(t_i)\Delta t_i$$

$$N_s\simeq N_t\Delta T{\int }_{E_{th}}^{E_{max}}\eta (E){\Phi }_{\bar{e}}(E)\sigma (E)dE.$$

$$F(T)=1-\frac{(RT{)}^0{\exp }^{-RT}}{0!}=1-e^{-RT},$$

$$f(T)=\frac{dF(T)}{dT}=R{e}^{-RT}.$$

$$\Lambda (T)=\frac{f(T;R_1)}{f(T;R_0)}=\frac{R_1}{R_0}{e}^{-(R_1-R_0)T}.$$

$$P_D={\int }_{T:\Lambda (T)>k}f(T;R_1)dT,$$

$$P_F={\int }_{T:\Lambda (T)>k}f(T;R_0)dT,$$

Full text

I Introduction

Gravitational wave (GW) astronomy is rapidly developing. Since the first discovery of GW from a binary black hole merger in 2016 Abbott et al. (2016a), several binary mergers of black holes and neutron stars have been observed by LIGO-Virgo (see, e.g. Abbott et al. (2018)). One event in particular, a binary neutron star merger, was heralded as the first multi-messenger observation of a merger, with detections in GW and at electromagnetic frequencies Abbott et al. (2017a); Abbott et al. (2017b, c); Soares-Santos et al. (2017); Hallinan et al. (2017). It delivered a wealth of scientific information, largely confirming theoretical predictions of the merger dynamics and GW emission, and establishing neutron star mergers as sources of short gamma ray bursts Abbott et al. (2017c) and sites of r-process nucleosynthesis Drout et al. (2017); Chornock et al. (2017); Pian et al. (2017). Several open questions remain, especially about the nature of the compact object which is born in the merger and its short- and long-term evolution.

At this time, neutrinos are still missing from the multi-messenger picture of mergers Albert et al. (2017); Hayato et al. (2018); Albert et al. (2019); Avrorin et al. (2018). Neutrino production is expected in compact-object mergers, i.e., binary neutron star (BNS) and neutron-star-black hole (NSBH) mergers. One site of neutrino production is the remnant of mergers and its accretion disk, where the dense and hot conditions allow efficient cooling via neutrinos similarly to a core collapse supernova. Therefore, a burst of $\sim$10 MeV, thermal neutrinos is a generic prediction (see, e.g. Janka and Ruffert (1996); Ruffert et al. (1997); Rosswog and Liebendoerfer (2003)). Furthermore, another possible channel of neutrino emission is due to the pion and muon decays generated through the photomeson production. The latter takes place when accelerated cosmic ray particles propagate through background photons and produce high energy (multi-TeV) neutrinos Gao et al. (2013); Fang and Metzger (2017); Kimura et al. (2017); Murase et al. (2018); Biehl et al. (2018). This scenario predicts a high energy neutrino flare accompanied by a gamma ray burst. Strong neutrino emission is possible if the remnant is a (quasi-)stable magnetar, whereas little or no high energy neutrino production is expected if a black hole is directly produced (e.g., Gao et al. (2013); Fang and Metzger (2017); Kimura et al. (2017)).

The theme of this paper is the thermal neutrinos from BNS and NSBH mergers and their detectability. A detection of these neutrinos would be very important to advance our understanding of mergers. Indeed, while GW data are mostly sensitive to the physics before and during the merger, neutrinos can give unique information on the post-merger phase. Models show that the duration of the neutrino burst, its luminosity and its energy spectrum depend strongly on the type of remnant produced (stable neutron star, transient neutron star or black hole), on the mass and accretion rate of the accretion disk surrounding the remnant, etc. (see e.g., Janka et al. (1999); Sekiguchi et al. (2011); Just et al. (2015); Fujibayashi et al. (2017); Lippuner et al. (2017) for examples, more details in Sec. II.1). The features of the neutrino burst in turn have an important impact, through their effect on the electron fraction, on the rates of formation of heavy elements via the r-process, see, e.g., Wanajo et al. (2014); Foucart et al. (2016); Roberts et al. (2017). Moreover, studying the neutrino emission could be important to understand the process of neutrino-antineutrino annihilation into $e^{+}\ e^{-}$ pairs. This mechanism transfers part of the gravitational energy of the accretion disk to electromagnetic radiation, thus possibly powering a short gamma ray burst Dessart et al. (2009); Janka et al. (1999); Fujibayashi et al. (2017). Finally, the flavor energy spectra of the merger neutrino burst at Earth will carry an imprint of the complex pattern of neutrino flavor oscillations inside the accretion disk, with the potential to probe its density and temperature profiles (e.g., Malkus et al. (2016); Zhu et al. (2016); Frensel et al. (2017); Wu and Tamborra (2017); Tian et al. (2017)).

The thermal neutrinos from mergers are very challenging to detect. Their diffuse flux is overwhelmed by the diffuse supernova neutrino background – which has a similar spectrum – and by experimental backgrounds; only in extreme cases, its observation might be possible Caballero et al. (2009); Schilbach et al. (2018). Recently, a long-term strategy has been proposed by Kyutoku and Kashiyama Kyutoku and Kashiyama (2018), which is based on exploiting the time-coincidence with GW detections to strongly reduce backgrounds, and identify individual neutrino events from mergers. Assuming the projected performance of the future Advanced LIGO Abadie et al. (2010), with its sensitivity up to $\sim 200$ Mpc distance (i.e., redshift $z\simeq 0.05$), it was found that a Mt-scale water Cherenkov detector like the planned HyperKamiokande Abe et al. (2011, 2018) might obtain a high significance detection of ${\mathcal{O}}(1)$ neutrinos from mergers, over about a century of operation. The operating time needed would be significantly reduced if a multi-Mt detector like DeepTITAND Kistler et al. (2011) and MICA Böser et al. (2015) is constructed in the next few decades.

In parallel to the development of neutrino observatories, next generation GW detectors, like the proposed Voyager B. P. Abbott et al. (2016), Einstein Telescope Punturo et al. (2010), and Cosmic Explorer Abbott et al. (2017d), will become a reality. They will allow a high-efficiency sensitivity to binary mergers as far as $z\sim 2$ Mills et al. (2018), thus opening up a much larger cosmic volume for searches of neutrinos in coincidence with GW detections. Inspired by this prospect, here we elaborate on the time-coincidence method for a future scenario where both a Mt-scale neutrino detector and at least one GW detector with cosmological ($z\gtrsim 1$) sensitivity are operating in synergy for several decades. One main element of novelty of this work is that we examine the detection of neutrinos produced at cosmological distances, studying the contributions of mergers at different redshifts to the total signal. We also discuss the significance and physics potential of a very low statistics signal, at the level of a single neutrino per decade or so. Another new element is that results are shown for different assumptions on the merger rate, and for several different models of neutrino emission from mergers.

The paper is structured as follows: in Sec. II, scenarios of neutrino production in BNS and NSBH mergers are illustrated, and the calculation of the neutrino flux at Earth is presented. In Sec. III, the method of detection is described, and results are shown for the number of events of signal and background, and for the probability of detection as a function of the detector exposure. In the discussion section, Sec. IV, the results are examined in their broader implications, and conclusions are drawn.

II Neutrinos from binary mergers: fluxes and spectra

II.1 Neutrino emission in mergers

We consider mergers of binary systems made of two neutron stars or a neutron star and a stellar-mass black hole. We assume that the black holes are of stellar origin, i.e., the product of core collapse supernovae.

As a result of a merger, a massive remnant is formed, surrounded by an accretion disk. For NSBH mergers, the outcome is a black hole-torus system. Matter in the torus is accreted on the black-hole, so that part of its gravitational binding energy is converted into internal energy. For a mass $M_{accr}\sim 0.1M_{\odot}$ of material being accreted, the internal energy budget is $E_{tot}\gtrsim 10^{53}$ erg, of which about half is thermal energy and is released via neutrino emission Dessart et al. (2009); Perego et al. (2014). For BNS mergers, the post merger phase could be different depending on the binary parameters and the properties of the nuclear equation of state. It could lead to the direct formation of a black hole, in which case neutrino production would proceed similarly to the NSBH case. Other possibilities are either the formation of a stable NS in the central remnant or the birth of a transient hypermassive neutron star (HMNS). The latter might eventually collapse into a black hole. In the presence of a newly formed neutron star, the neutrino luminosity receives contributions from both the accretion process and the cooling of the neutron star itself. The neutrino luminosity from a HMNS or a stable neutron star is powered by an internal energy reservoir, which is approximately the difference between the internal energy of a hot and of a cold HMNS/NS. For a HMNS of $\sim 2.5M_{\odot}$ and $\sim 15$ MeV , the internal energy is E_{th}\sim 3\times\ 10^{52}\ ergs Perego et al. (2014).

For both BNS and NSBH mergers, typically the neutrino emission is dominated by $\mathrel{{\nu_{e}}}$ and $\mathrel{{\bar{\nu}}_{e}}$, with a certain preponderance of $\mathrel{{\bar{\nu}}_{e}}$ due to the matter of the disk – and of the compact object, if a neutron star is formed – being neutron-rich. In scenarios where a HMNS collapses into a black hole, the $\mathrel{{\nu_{e}}}$ and $\mathrel{{\bar{\nu}}_{e}}$ luminosities only decrease gradually, around the time of collapse, because these species continue to be produced in the accretion disk via charged current processes. On the other hand, the muon and tau neutrino luminosity decreases sharply after the black hole formation because their inner, hotter surface of last scattering is swallowed by the black hole. The duration of a merger neutrino  burst may be of the order of a fraction of a second, up to several seconds, depending on the time-scale of the evolution of the central object Sekiguchi et al. (2011, 2015); Janka et al. (1999); Just et al. (2015); Fujibayashi et al. (2017).

A detailed description of the physics of mergers, and their neutrino emission requires large scale numerical modeling. Here we adopt a simplified description of a burst of $\mathrel{{\bar{\nu}}_{e}}$ with thermal (Fermi-Dirac) spectrum with duration of up to $\sim 1$s. The effect of flavor oscillations on the $\mathrel{{\bar{\nu}}_{e}}$ spectrum and luminosity at Earth will be accounted for by a energy-independent suppression factor (see Sec. III). The total energy emitted, as well as the burst duration and average energies are drawn from numerical simulations; as summarized in Table 1. These simulations suffer from various approximations that are made necessary by the limited computational resources. Still, they are state-of-the-art, and give reasonable estimations on ejected matter as well as r-process nucleosynthesis properties consistent with previous works Sekiguchi et al. (2016); Martin et al. (2015); Kyutoku et al. (2013); Kyutoku and Kashiyama (2018). Therefore, the model results in Table 1 are credible, and can be used to define a space of the parameters that can be considered plausible for neutrinos from mergers. Here, we will focus mainly on the two-dimensional space of the total energy in $\mathrel{{\bar{\nu}}_{e}}$ and the average energy of the $\mathrel{{\bar{\nu}}_{e}}$ spectrum (see Fig. 7).

A brief description of each model follows below.

• •

Lippuner et al. Lippuner et al. (2017), computed a long term ($\sim 10$ s simulation time) evolution of the accretion disk surrounding a HMNS of variable lifetime, including a stable one (infinite lifetime, model Hinf), which we use here. In this two-dimensional simulation, the HMNS (of mass $\sim$3 $M_{\odot}$) is approximated by using a reflecting inner boundary condition Metzger and Fernández (2014), and a parameterized isotropic outward neutrino flux is imposed on the boundary. The neutrino luminosity has a time dependence of $\approx t^{-\frac{1}{2}}$, which is general of the cooling of a proto-neutron star (see, e.g., Pons et al. (1999)). When the HMNS collapses into a black hole, the boundary is changed from reflecting to absorbing, and the neutrino emission is set to zero. Lippuner et al. also produced additional models with spinning black holes at the center; among these, we use models B090 and BF15.

• •

Sekiguchi, Kiuchi, Kyutoku and Shibata Sekiguchi et al. (2011), simulated BNS mergers, using the Shen et al. Shen et al. (1998) equation of state, and including effects of general relativity and neutrino cooling. A long-lived HMNS is formed after the merger. Among their results, we use the one for a system of equal-mass neutron stars with an intermediate mass value, $M_{NS}=1.5~{}M_{\odot}$ (“middle” model, or M). Because of the continuous collision between HMNS and spiral arms, the outer region of the HMNS is heated up, and, as a result, neutrino emission is strong Rosswog and Liebendoerfer (2003); Dessart et al. (2009). Sekiguchi et al. found that a HMNS with mass $\lesssim 3M_{\odot}$ has a life time much longer than its dynamical time scale ($t_{dyn}\sim$10 ms) and determined by the neutrino cooling process, with time-scale of the order of seconds. The simulations by Sekiguchi et al. end at about 20 ms after the merger, not long enough to estimate the total energy emitted in $\mathrel{{\bar{\nu}}_{e}}$, which is needed here. To circumvent this difficulty, we used the similarity between model M by Sekiguchi et al. and model Hinf by Lippuner et al., which share similar HMNS mass and neutrino emission time scale. Considering the generality of the luminosity time dependence found in Lippuner et al. (see above), we used it to extend the results of Sekiguchi et al beyond 20 ms (with the appropriate normalization of the luminosity to match the numerical result at 20 ms, $L_{\mathrel{{\bar{\nu}}_{e}}}\approx 2.3\times 10^{53}\,erg\,s^{-1}$).

• •

Fujibayashi, Sekiguchi, Kiuchi and Shibata Fujibayashi et al. (2017), performed long term, two dimensional simulations of the evolution of the remnant of a BNS mergers – with general relativity effects – using initial conditions taken from a three-dimensional, numerical relativity simulationSekiguchi et al. (2015). In their fiducial model (called DD2-135135-On-H, see Fujibayashi et al. (2017) for its detailed meaning) the $\mathrel{{\bar{\nu}}_{e}}$ luminosity is $L_{\mathrel{{\bar{\nu}}_{e}}}\sim 10^{53}\,erg\,s^{-1}$ at the beginning of simulation, and decreases rapidly over $\sim$100 ms, setting to a nearly constant value of $L_{\mathrel{{\bar{\nu}}_{e}}}\sim 10^{52}\,erg\,s^{-1}$ after 300 ms. The simulation ends at 400 ms. To estimate the total energy in $\mathrel{{\bar{\nu}}_{e}}$, we extrapolate the neutrino luminosity up to $t\sim 1$ s, similarly to the approach used by Kyutoku and Kashiyama Kyutoku and Kashiyama (2018).

• •

Just, Bauswein, Pulpillo, Goriely and Janka Just et al. (2015), studied NSBH and BNS merger resulting in BH-torus systems, using a relativistic smooth-particle-hydrodynamic code Oechslin et al. (2002), with Newtonian hydrodynamics, and including viscosity effects. They found that the masses of black hole and of the torus are important parameters, strongly influencing both the neutrino-driven wind and the neutrino luminosities. The neutrino energy loss rates can reach $L_{\nu}\sim few~{}10^{52}\,erg\,s^{-1}$, and even exceed $L_{\nu}\sim 10^{53}\,erg\,s^{-1}$ for each of $\nu_{e}$ and $\mathrel{{\bar{\nu}}_{e}}$ for a few hundreds of ms. For more detail on these models, we refer to Table 1 (see also Table 2 of Just et al. (2015)).

• •

Janka, Eberl, Ruffert and Fryer Janka et al. (1999), performed Newtonian hydrodynamic simulations of NSBH and BNS mergers using a relativistic EoS Lattimer and Swesty (1991). In their simulations, the maximum temperature $T_{max}$ of the accretion flow can reach several tens of MeV, resulting in efficient neutrino emission at densities $\rho\sim 10^{10}-10^{13}\,{\rm g\,cm^{-3}}$. Total energies ${\mathcal{E}}_{\nu}\approx 3\times 10^{52}\,erg$ are found for the neutrino emission. The parameters of their models are given in details in Table 1 (see also Table 2 in Janka et al. (1999)). Note that the time duration of the $\mathrel{{\bar{\nu}}_{e}}$ luminosity in the Janka et al. model is approximated by the authors to be the same as torus life, which is around $50-150$ ms. We assume this to be the time $T_{99}$ shown in the Table.

II.2 Cosmological merger neutrino flux

For a given class of mergers (BNS or NSBH), the cosmological neutrino flux depends on the volumetric merger rate (number of mergers per unit of comoving volume, per unit time), as a function of the redshift, $R(z)$. This quantity depends on the core collapse supernova rate, which in turn tracks the star formation rate with its characteristic redshift evolution like $\sim(1+z)^{3}$ (for $z\lesssim 1$). It also depends on the post-collapse evolution of the remnants, which is only partially understood; therefore uncertainties on the merger rate are large (see, e.g. Dominik et al. (2013); Belczynski et al. (2016); Chruslinska (2018); Mapelli and Giacobbo (2018) for overviews).

For each merger type, we use two representative models for $R(z)$, whose predictions are consistent with the LIGO-Virgo  observational bounds. The first, more optimistic (“optimistic” from here on), is from the work of Eldridge, Stanway and Tang Eldridge et al. (2019). Specifically, we use the results for the kick model there (see also Bray and Eldridge (2018) for a description of the model itself) 111In Eldridge et al. (2019), the redshift dependence of the rate for the kick model is not explicitly given. So, we constructed an approximate analytical form by taking an analytical expression used for other models in the same paper, with parameters adjusted to fit the tabulated rates (Table 2 in Eldridge et al. (2019)) for the kick model. . The second model, more conservative (“moderate” from here on), is from Mapelli and Giacobbo Mapelli and Giacobbo (2018) (case labeled as $\alpha=5,~{}$ low $\sigma$ in Mapelli and Giacobbo (2018)). Both models are illustrated in Fig. 1 for BNS and NSBH mergers.

To estimate the neutrino flux from mergers in a detector, one needs to use the comoving volume enclosed between redshift $z$ and $z+dz$, which can be written in terms of the comoving distance, $D_{c}$, as follows:

[TABLE]

The rate of mergers with $z<z_{max}$, in the frame on an observer on Earth, is then given by:

[TABLE]

Here, $\Omega_{m}=0.3$ and $\Omega_{\Lambda}=0.7$ are the fractions of the cosmic energy density in matter and dark energy respectively; $c$ is the speed of light and $H_{0}$ is the Hubble constant. We assume the $F_{\mathrel{{\bar{\nu}}_{e}}}(E)$ follows Firmi-Dirac distribution with zero chemical potential. Therefore, the (time-integrated) number of $\mathrel{{\bar{\nu}}_{e}}$ emitted per unit energy by an individual merger, $F_{\mathrel{{\bar{\nu}}_{e}}}(E)$, is:

[TABLE]

where ${\mathcal{E}}_{\mathrel{{\bar{\nu}}_{e}}}$ is the total neutrino energy emission, $\langle E_{\mathrel{{\bar{\nu}}_{e}}}\rangle$ is the mean neutrino energy, $T$ is the neutrino temperature, and $2/3T^{3}\zeta(3)$ serves as a normalization factor of the Fermi-Dirac distribution. Given $F_{\mathrel{{\bar{\nu}}_{e}}}$,the total, time averaged flux (differential in energy, surface and time) at Earth due to all mergers with $z<z_{max}$ can be expressed as:

[TABLE]

where $E^{\prime}=E(1+z)$, describes the redshift of the neutrino energy during propagation. The constant $f_{osc}$ is a phenomenological factor accounting for the effect of flavor oscillations. Detailed studies of neutrino oscillations in the post-merger environment for BNS and NSBH mergers have shown a complex pattern that depends strongly on the relative intensity of the different flavor fluxes and other parameters describing the compact object and the accretion disk. The $\mathrel{{\bar{\nu}}_{e}}$ survival probability varies, roughly, between $0.5$ and 1 (corresponding to no oscillations), see, e.g., Malkus et al. (2012, 2016); Zhu et al. (2016); Tian et al. (2017); Vlasenko and McLaughlin (2018). Therefore, here we set $f_{osc}=1$.

III Detectability

III.1 The method; expected backgrounds

We consider a near-future scenario where: (i) a neutrino detector of mass ${\mathcal{O}(100)}$ kt or higher exists, with good timing resolution and low energy threshold, like the upcoming HyperKamiokande; and (ii) one or more next generation GW detectors are available, with the capability to observe mergers with $\sim 100\%$ efficiency, up to a redshift $z_{GW}\gtrsim 1$. We also assume that the GW data will allow to establish the time of the merger and its distance/redshift with good precision (e.g., tens of per cent)Safarzadeh et al. (2019).

In this scenario, a realistic detection method would be the one proposed by Kyutoku and Kashiyama Kyutoku and Kashiyama (2018), which is generalized here to account for non-zero redshift. The method consists of considering the $N$ mergers that are observed in GW over a long period of time, $\Delta T$, with their merger times $t_{i}$ ($i=1,2,...,N$), and redshifts $z_{i}$, and restricting the neutrino data analysis to time windows of width $\Delta t_{i}\sim T_{99}(1+z_{i})$ after each merger, where $T_{99}$ is the time duration of neutrino burst within which $99\%$ of $\nu_{\bar{e}}$ energy has emitted. Here the factor $(1+z_{i})$ ensures that each time windows is appropriately adjusted for the cosmological time dilation. Thus, the number of neutrino events collected in the detector over $\Delta T$ is:

[TABLE]

where, $n_{s}$ is the (average) neutrino event rate per merger. For the same observation time, the expected number of background events in the detector has an expression analogous to Eq. (6), where the rate of background events, $n_{b}(t)$, can be approximated as constant in time: $N_{b}\simeq n_{b}\sum^{N}_{i=1}\Delta t_{i}$. The effect of restricting to the effective observation time, $\Delta T_{eff}=\sum^{N}_{i=1}\Delta t_{i}$ results in a strong enhancement of the signal-to-background ratio Kyutoku and Kashiyama (2018).

For simplicity, we neglect the possible corrections due to the neutrino mass and to gravitational lensing on the expected time window $\Delta t_{i}$ Stodolsky (2000); Baker and Trodden (2017) (see also discussion in Kyutoku and Kashiyama (2018)), as well as effects of possible physics beyond the Standard Model of particle physics. We also neglect differences in the neutrino emission between individual mergers, thus making the somewhat simplistic assumption that all mergers of the same class are identical. We estimate that the corrections due to these neglected effects would be sub-dominant compared to the very large uncertainties on the neutrino emission models and merger rates.

Let us consider a water Cherenkov detector of the next generation, and its main channel of detection, inverse beta decay: $\mathrel{{\bar{\nu}}_{e}}+p\rightarrow n+e^{+}$. One can approximate the sum for $N_{s}$ as an integral of the diffuse flux over the observation time $\Delta T$:

[TABLE]

In this approximation, it is assumed that the choice of the time windows, $\Delta t_{i}$, allows to capture the entire post-merger neutrino flux, i.e., $\Delta t_{i}/(1+z_{i})>T_{99}$. In this work, a rather conservative time window, $\Delta t_{i}/(1+z_{i})=1$ s (to be compared with $T_{99}$ for various models, see Table 1), will be used in all calculations. Here $\sigma(E)$ is the detection cross section and $N_{t}$ is the number of protons in the detector; $\eta(E)$ is the detector efficiency. The interval $E=E_{th}-E_{max}$ is a suitable energy window, determined by the need to minimize the background, see Sec. III.2.1.

Let us now discuss the experimental background. We consider two different detector configurations, the first with pure water, and the second with water with the addition of Gadolinium (Gd) Beacom and Vagins (2004), for better background reduction. The option with Gd is being realized, for the first time, in the upgraded SuperKamiokande (SuperK-Gd; see, e.g., Marti-Magro (2018)), and is envisioned as a possibility for a second phase of HyperKamiokande Abe et al. (2018). For pure water, the region of visible energy $E_{vis}\lesssim 18$ MeV is unaccessbile due to the overwhelming spallation background. At higher energy, the background is dominated by invisible muons (Fig. 2, bottom panel): these are sub-Cherekov atmospheric muons that become visible in the detector only when they decay, thus mimicking inverse beta decay. We take $\eta\simeq 0.9$ as a realistic efficiency Hirata et al. (1988); Kyutoku and Kashiyama (2018); Abe et al. (2011).

For water with Gd, the background due to spallation can be effectively reduced, with only a small residual contribution due to the production of a short-lived Lithium isotope (${\rm{}^{9}Li}$). Therefore, the window of sensitivity extends down to $E_{vis}\sim 11$ MeV, below which events due to reactor $\mathrel{{\bar{\nu}}_{e}}$ dominate. With Gd, the invisible muon background is suppressed by a factor of $\sim 5$ Beacom and Vagins (2004), and the total background at $E_{vis}\gtrsim 11$ MeV has comparable contributions from different sources. Specifically, besides invisible muons, one should include inverse beta decay events due to the diffuse supernova neutrino background and to atmospheric $\mathrel{{\bar{\nu}}_{e}}$, as well as events due to neutral current scattering of atmospheric neutrinos. For illustration of these backgrounds, we refer to the specific set of background cuts that was discussed for the diffuse supernova neutrino background at SuperK-Gd and HyperKamiokande Kunxian (2015); Abe et al. (2018). These result in a signal efficiency $\eta\simeq 0.67$, and produce the residual background spectrum shown in Fig. 2, top panel.

III.2 Results

III.2.1 Event rates and physics potential

In this section, results will be shown for the signal and background event rates. Unless stated otherwise, they refer to a detector exposure $T_{exp}=100~{}{\rm Mt~{}yr}$, corresponding, e.g. to a 1 Mt (10 Mt) mass detector running for a century (decade). While such exposure does not seem entirely realistic for HyperKamiokande as currently planned (total fiducial mass of 374 kt, see Abe et al. (2018)), it is technlogically feasible, and within the realm of possibility for the next-to-next generation water Cherenkov detectors. For illustration, here several results will be shown for neutrino parameters (${\mathcal{E}}_{\mathrel{{\bar{\nu}}_{e}}}$ and $\langle E_{\mathrel{{\bar{\nu}}_{e}}}\rangle$) of the specific models in Table 1. Considering the uncertainties that affect numerical simulations (see Sec. II.1), results that are intermediate between models should also be considered possible.

To assess the detectability prospects, let us first examine the energy spectra of signal (in terms of observed positron energy, $E_{e}$) and background events, and establish the interval of energy where the signal exceeds the background. The spectra are shown in Fig. 3 for mergers of both types (BNS and NSBH), with redshift in the interval $z=0-0.05$. For water with Gd, it appears that the energy window of signal dominance is $E_{e}\sim 14-34$ MeV for the most conservative models. It extends down to $E_{e}\sim 11$ MeV and beyond 50 MeV for the most luminous mergers of either type. For pure water, the energy window is bound from below by the threshold due to spallation (see Sec. 2); its upper edge is as low as $\sim 25$ MeV or as high as 50 MeV or beyond depending on the signal flux model.

In principle, the energy window depends on the redshift of the merger of interest, because of the effect of the redshift of energy on the signal spectrum. We explored the possibility of using a $z-$dependent window to further reduce backgrounds, and found that it provides some advantage for the models with lower neutrino emission, but is slightly detrimental for the most optimistic models, where the energy window extends beyond 50 MeV for a wide range of $z$.

For definiteness, we settle on a moderately conservative and $z-$independent energy window, $E_{e}=18-50$ MeV ($E_{e}=11-50$ MeV) for pure water (water with Gd)222The decision to limit the analysis to $E_{e}<50$ MeV is in part motivated by the fact that the background is poorly known at higher energy, where contributions from the inelastic scattering of high energy atmospheric neutrinos becomes increasingly important..

When searching for neutrinos in time coincidence with a GW detection, a relevant question is the distance (or, equivalently, redshift) range of sensitivity of the neutrino detector. Fig. 4 shows an example of distribution in redshift bins (of width $\Delta z=0.05$) of the number of mergers expected to be observed at Earth (assuming 100% efficiency in GW detection), and of the neutrino signal and background event rates. The most optimistic models of neutrino emission (BNS merger with long lived HMNS from model M, see Table 1) and of merger rates were chosen here (see figure caption). Additionally, a detector with Gd was assumed here.

We see that the number of mergers grows rapidly with $z$ for $z<0.5$, reflecting the growth of the merger rate and of the cosmic volume with redshift. For $z<0.05$ ($z<0.5$), $N\sim 10^{4}$ ($N\sim 10^{6}$) mergers per century are expected, in agreement with existing predictions for GW  detectors of the next generation, see, e.g., Mills et al. (2018).

The number of background events in a neutrino detector roughly follows the growth of the number of mergers, with a faster increase with $z$ due to our accounting for time dilation when choosing the time window of observation for coincidence (see Sec. III.1). The number of signal events can be as high as $N_{s}\sim 10$ for $z<0.05$. It decreases with $z$, because the increase of the merger rate with $z$ is overcompensated by effects of the redshift of energy (as the redshift increases, a larger fraction of the neutrino flux falls below the energy window), and of the volume-dilution of the neutrino flux from a single source. The rate of signal events exceeds the background in the first 3-4 bins, corresponding to $z<z_{S/B}\simeq 0.2$, where $z_{S/B}$ is the redshift at which the number of background events becomes comparable to the number of signal events.

Figs. 5 and 6 show the number of signal and background events for the interval $z<z_{max}$, as a function of $z_{max}$, for the different models considered in this work. Results are shown for different combinations of merger types, merger rates and detector configuration. For the case of water with Gd, the figures confirm the results shown in Fig. 4 for the most optimistic scenario. In addition, they show that for other models of BNS mergers, and for all the models of NSBH merger, the numbers of signal events are at the level of $N_{s}\sim 0.1-2$ for $z_{max}\sim 0.1$. Backgrounds overwhelm the signal already at $z_{S/B}\sim 0.1-0.2$. For pure water, the figures reveal that the situation is noticeably worse due to the larger background rate. $z_{S/B}<0.1$ is expected for the more conservative flux models, and in certain cases the background exceeds the signal at all $z$.

Let us briefly discuss the statistical significance of a detection. As a criterion for statistical significance we establish that the total predicted number of events, $N_{tot}=N_{S}+N_{B}$, should be higher than a threshold $N_{T}$ beyond which the hypothesis of background-only is unlikely to be realized: $\sum_{n=N_{T}}^{\infty}P(n|N_{B})<0.0015$. Here $P(n|N_{B})$ is a Poisson distribution of $n$, with the mean being $N_{B}$. The criterion we apply here naturally becomes the 3-sigma-limit when the statistics is high.

The statistical significance of a detection is dependent on both the effective observation range, $z_{max}$, and the operation time of a network of neutrino/GW detectors. Because the signal-to-background ratio decreases with $z_{max}$, for a given flux model and exposure, one can find the maximum $z_{max}$ within which the signal could be distinguished from noise with high statistical significance. For example, for the scenario in the left upper (left bottom) panel of Fig. 5, where BNS mergers predicted by model M with optimistic (moderate) merger rate, a detector with Gd, and a $100$ Mt year exposure are used, we find that this extreme value is $z_{max}=0.74$ ($z_{max}=0.29$), for which the number of signal events is $N_{S}=71$ ($N_{S}=11$). The corresponding number of background events is $N_{B}=529$ ($N_{B}=10$). Including mergers with higher redshift in the analysis will worsen the significance of the signal, due to the overwhelming background.

Let us briefly discuss the physics potential of a search for neutrinos from mergers. Given a measured event rate, it will be possible to place constraints on the main parameters governing the merger neutrino flux, namely the total energy $\mathcal{E}_{\bar{\nu}_{e}}$ and the average energy of the $\mathrel{{\bar{\nu}}_{e}}$ spectrum,$\langle{E_{\mathrel{{\bar{\nu}}_{e}}}}\rangle$. As an illustration of possible constraints on these parameters, in Fig. 7 we show two sets of iso-contours of the total number of events (signal plus background) in the space of $\mathcal{E}_{\mathrel{{\bar{\nu}}_{e}}}$ and $\langle{E_{\mathrel{{\bar{\nu}}_{e}}}}\rangle$, for the two redshift intervals of reference, $z\leq 0.05$ and $z\leq 0.74$.The optimistic BNS merger rate was adopted here. For comparison, the predictions from each of the models in Table 1 for BNS mergers are shown. We assume that uncertainties on the merger rates and on the background rates are negligible, so errors on the number of events are only statistical. In each of the panels of the figure, the contours are chosen so that they form two bands of $90\%$ confidence level statistical uncertainty centered around two hypothetical measured event rates, one pessimistic and one optimistic ($N_{tot}=1,8$ for $z\leq 0.05$ and $N_{tot}=537,600$ for $z\leq 0.74$, see figure caption for more information). These bands represent the region of the parameter space that can be considered consistent with the measurement. We see that a measurement of a relatively large event rate would disfavor a wide portion of the parameter space, and result in a lower bound on ${\mathcal{E}}_{\mathrel{{\bar{\nu}}_{e}}}$, under minimal assumptions on the average energy (e.g., if one assumes ${\mathcal{E}}_{\mathrel{{\bar{\nu}}_{e}}}<24$ MeV as a theoretical prior). In the case of a low measured number of events, most of the parameter space would be allowed (including ${\mathcal{E}}_{\mathrel{{\bar{\nu}}_{e}}}=0$, the background-only case), and only the most extreme corner of the parameter space would be disfavored. In this latter case, most or even all the models would be consistent with the data.

In addition to bounds from the event rate, it will be possible to further constrain the parameter space using the energy distribution of the events. Roughly, these will correspond to vertical lines in Fig. 7, which we do not show for simplicity. More realistically, a fit of the data energy spectrum and overall normalization will give characteristic boomerang-shaped confidence level contours similar to those obtained for neutrinos from the supernova SN1987A (see, e.g. Lunardini and Smirnov (2001) and references therein). A detailed illustration of the potential of data fits is beyond the scope of the present work.

III.2.2 Probability of first detection

As was shown in the previous section (see Fig. 7), even the detection of a single candidate event, or even a null result (no detection) can be significant, and provide useful constraints. This can be especially relevant in the short term, when the detector exposure is still low. As a further illustration, we follow Kyutoku and Kashiyama Kyutoku and Kashiyama (2018), and discuss how the waiting time – i.e., the time that elapses from the beginning of the experimental search until the first candidate event is observed – is a useful statistical observable to distinguish between different merger models.

We consider the total number $N_{tot}$ of events that are observed in time-coincidence with GW detections, and suppose that it follows a Poisson process with event rate $R$. Then the length of time until the first arrival, denoted by $T$, is a continuous random variable. In a simple scenario where $z_{max}$ is fixed (rather than adjusted depending on the model, to optimize the statistical significance), the rate is a constant, $R=N_{tot}/T$, which is determined by the merger flux model. The background rate can be assumed as known here.

The cumulative distribution function of $T$, which is the probability that the first arrival is observed within the time interval $[0,T]$, is:

[TABLE]

and the probability density that the first arrival is observed within the time interval $[T,T+dT]$ is333Note that $f(T)dT$ is the probability that no events are recorded in the interval $[0,T]$, and one event is recorded in the interval $[T,T+dT]$.:

[TABLE]

With a single measurement of $T$ (waiting time of the first event arrival) with probability density function given by $f(T)$, one can use a likelihood ratio test to distinguish between two hypotheses, $H_{1}$, and $H_{0}$, according to the Neyman-Pearson lemma 444When we want to distinguish an alternative theory $H_{1}$ from a default (“null”) theory $H_{0}$, a type I error occurs if we reject $H_{0}$ when it is true. The power of a hypothesis test is the probability of making a correct decision if $H_{1}$ is true. The Neyman-Pearson lemma Kay (1993) states that a likelihood ratio test is the test that has maximum power for a given type I error. .

To fix the ideas, let us consider an example drawn from Fig. 7 (upper panel): a search for neutrinos from BNS mergers, for the optimistic merger rate, a 1 Mt detector (with Gd) and $z_{max}=0.05$. In this example, hypothesis $H_{1}$ predicts a rate $R_{1}=8\cdot 10^{-2}~{}{\rm yr^{-1}}$ (black solid curve in the upper panel of Fig. 7), while $H_{0}$ predicts $R_{0}=10^{-2}~{}{\rm yr^{-1}}$ (blue solid curve in the same panel).

Let us introduce the likelihood ratio:

[TABLE]

When a measurement gives $\Lambda\gg 1$ ($\Lambda\ll 1$), hypothesis $H_{1}$ ($H_{0}$) is strongly favored over the other one. Focusing on the case $\Lambda>1$ for definiteness, one can choose a threshold value, $k>1$, such that if $\Lambda>k$ hypothesis $H_{1}$ is accepted as true, and $H_{0}$ is discarded (as unlikely to be true).

Let us also define the probabilities:

[TABLE]

where the integration is performed in the interval $[0,T_{k}]$, where $T_{k}$ is the value of $T$ such that $\Lambda=k$ is satisfied.

Here $P_{D}$ (“probability of detection”, in common jargon) is the probability that hypothesis $H_{1}$ is correctly chosen. In other words, it is the probability that the analysis will correctly identify the true hypothesis, and represents the power of the hypothesis test. $P_{F}$ is a “false alarm” probability; it is the probability that $H_{1}$ is (wrongly) chosen while $H_{0}$ is true. It represents the probability that the analysis will give an incorrect answer, and is a type I error (see, e.g. Kay (1993)). Both $P_{D}$ and $P_{F}$ increase with decreasing likelihood ratio threshold $k$, or, equivalently, with increasing waiting time threshold $T_{k}$. A detector with good performance to test one model against another is expected to have high $P_{D}$ while low $P_{F}$. This behavior can be seen in the Receiver Operating Characteristic (ROC) curve, which is defined as the path in the space of $P_{F}$ and $P_{D}$ which is obtained by varying $k$.

For our example, the ROC curve is shown in the bottom panel of Fig. 8. From this figure we learn that, if a relatively short waiting time threshold is set, e.g., $T_{k}=10$ yr, the power of the test is moderate, $P_{D}\sim 0.55$, although the type I error associated to it is rather low, $P_{F}\simeq 0.1$. Qualitatively, this can be understood by thinking that, if a short arrival time is observed, $T<T_{k}=10$ yr, then $H_{1}$ will be favored with a good level of confidence, but the probability of this situation to be realized is relatively low if $H_{0}$ is true.

If a longer waiting time threshold is established, e.g. $T_{k}=20$ yr, the power of the test increases to $P_{D}\simeq 0.82$, however the probability to falsely choose $H_{1}$ when instead $H_{0}$ is true also increase, $P_{F}\simeq 0.24$. Intuitively, this means that, by setting a weaker criterion on the likelihood ($\Lambda\geq 1.52$), we obtain a test which is more likely to favor $H_{1}$ (because the condition on the likelihood is rather likely to be met), however this comes at the price of having a larger probability to come to the incorrect conclusion.

If a very long waiting time is observed, $T\sim 100$ yr, then $H_{0}$ will be strongly favored over model $H_{1}$, because $\Lambda(T)\ll 1$, and for even larger $T$ both models might be poor interpretations of the data, and a third model (e.g., background only) might have to be considered as an alternative. Indeed, we must emphasize that the likelihood ratio method discussed here is in the spirit of discriminating between two models only, and therefore might be suitable to use if the theoretical panorama evolves towards a somewhat bipolar situation (two fairly established classes of models only).

IV Summary and discussion

We found that the detection of thermal neutrinos from mergers is possible, over the next several decades, if the synergy between Mt-scale neutrino detectors and next generation GW detectors, capable of observing mergers up to redshift $z\sim 2$, is fully exploited.

A very important element is the time coincidence between GW and neutrino detection, which will have two important benefits. The first is a strong reduction of the background at the neutrino detector, without which the identification of the neutrino flux from mergers would be almost impossible. The second is that each neutrino candidate event can be analyzed as a single, stand-alone signal, with measured energy and arrival time, and with an identified (potential) parent merger, of known merger time, redshift, type (BNS or NSBH) and morphology.

About the detectability of neutrinos from mergers, we find that, among the neutrino emission models we have considered, the strongest neutrino signal is expected for BNS mergers (see Table 1), if their rate is in the higher part of the interval currently allowed by the LIGO-Virgo data, and if the compact object emits more than $\sim 10^{52}$ ergs in $\mathrel{{\bar{\nu}}_{e}}$s, due to its being a long-lived or hypermassive neutron star. In this scenario, up to ${\mathcal{O}}(10^{2})$ neutrino events are predicted for an exposure of $100~{}{\rm Mt\cdot yr}$, for an energy window $E_{e}=10-50$ MeV and a redshfit window $z\simeq 0-0.5$. Several tens of percent of the candidate events (those with $z\lesssim 0.1$ or so) will be attributed to parent mergers of known redshift with high statistical significance (see Fig. 4). And part of them might even be individually statistically significant over the background.

For mergers with more conservative parameters, i.e. lower merger rates, and lower total neutrino luminosity – $E_{\mathrel{{\bar{\nu}}_{e}}}\sim few\cdot 10^{51}$ ergs – the event rates can be ${\mathcal{O}}(1)$ event or less for the same exposure ($100~{}{\rm Mt\cdot yr}$), and the minimum exposure needed for a first detection would be several tens of ${\rm Mt\cdot yr}$. Because the parameters and models considered here are only a subset of all the possibilities, our results should be interpreted as a broad (and tentative) interval of possible neutrino detection rates: $\dot{N}\sim(10^{-3}-1)~{}{\rm Mt^{-1}yr^{-1}}$.

On the front of what can be learned about mergers from a neutrino detection, data from a neutrino detector – after the background reduction allowed by the time-coincidence method – will probably be analyzed statistically to extract basic, low precision, information on the merger neutrino flux parameters, mainly the energy spectrum and the total energy emitted. These will then be used to constrain numerical models of neutrino emission, thus indirectly testing the features (e.g., mass of the compact object and/or of the accretion disk) of the underlying merger models. The methods used in such analyses will draw from the experience of other high-background searches, for example those for the diffuse supernova neutrino background, where the energy distribution and total count of the data are fit with a mixed model (signal plus background) with both signal and background parameters to be found by fitting (see, e.g. Malek et al. (2003)). Advanced analyses where the time window $\Delta t_{i}$ is also optimized by fitting may further constrain the underlying model, and also enhance the signal significance by further reducing the background contamination.

In the fortunate case (see above) where a few single neutrino detections are individually significant, there might be potential to search for correlations with the physical features of the mergers provided by the GW observation, using methods similar to those used at IceCube in analyses of high energy neutrino data in a multi-messenger context (see, e.g. Aartsen et al. (2018)). A full exploration of the potential of such event-by-event analyses is left for future work.

We emphasize that even a single neutrino detected in coincidence with a merger, or even a non-detection, can place important constraints on the neutrino luminosity in a merger, and therefore on the type and physics of the newly born compact object, if the detector exposure exceeds $\sim$10 Mt yr for a configuration with water plus Gadolinium. In this context, the waiting time before the first detection is a useful statistical observable. Considering HyperKamiokande (in its futuristic upgrade to water plus Gd, with the currently planned fiducial mass of $\sim 370$ kt) a waiting time of $\sim$30 years is needed to start constraining theoretical flux models. Such waiting time is realistic for a long-term, broad-scope project like HyperKamiokande.

In closing, we believe that the detection of thermal neutrinos in time-coincidence with gravitational waves from mergers is a new, realistic goal to be added to the agenda of Mt-scale neutrino detectors with $\sim$ MeV energy thresholds. It may be the first observation of $\sim 10$ MeV-scale neutrinos from individually resolved cosmological sources, thus drastically expanding the horizon of possibilities of low energy neutrino astronomy in the context of multi-messenger studies.

Acknowledgments

The authors acknowledge funding from the National Science Foundation grant number PHY-1613708. They thank the authors of ref. Lippuner et al. (2017) for sharing their numerical results and for discussions. They are also grateful to Liliana Caballero, Laura Cadonati and Mohammadtaher Safarzadeh for useful comments.