Introduction

We have entered a new phase of astronomical observations, the so-called multimessenger astronomy era. Experiments and observatories are more than ever now able to observe the sky in different energy regions (from eV to EeV) with different messengers (photons, cosmic rays, neutrinos, or gravitational waves).

Since 2019 April, the LIGO/Virgo collaboration (LVC) has been publicly releasing their alerts for gravitational waves (GW) directly through their own GraceDB service and through the GCN system. 58 Within a few minutes after the first detection, the first alert is sent with a precise time stamp and a rough sky localization allowing quick follow-ups from other observatories. More precise information on localization is provided in the following days.

The detected GW emitters are categorized by LVC into several types, for which high-energy neutrino (HE-ν) emission is also expected from relativistic outflows and hadronic interactions within these sources: binary neutron star mergers (BNS; Kimura et al. 2018), neutron star-black hole mergers (NSBH; Kimura et al. 2017), or binary black hole mergers (BBH; Kotera & Silk 2016). Such astrophysical objects may also emit MeV neutrinos (LE-ν; for BNS, see Foucart et al. 2016).

However, such a joint observation of GWs and neutrinos is yet to be observed. Even a single event of this type would provide useful information on the underlying mechanisms. Furthermore, high-energy neutrinos would allow improving the localization in the sky of a single GW event, increasing the chance for a pointing observatory (e.g., follow-up telescopes) to observe a third correlated signal if the alert is provided promptly.

The IceCube (Countryman et al. 2019) and ANTARES/ KM3NeT (Dornic et al. 2019) experiments are already taking part in such a follow-up program for every single reported GW event. Nevertheless, these neutrino telescopes are mainly 55 Also at BMCC/CUNY, Science Department, New York, NY 10007, USA. 56 Now at IFIC (CSIC-U. Valencia), Spain. 57 Deceased.

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

covering HE-ν above 100 GeV. Super-Kamiokande (SK) can complement such follow-up studies, as, despite its much smaller size, it is sensitive to lower energies (MeV-TeV). In the past, SK has performed such studies, but only for a few of the detected GW events: GW150914/GW151226 in Abe et al. (2016a) and GW170817 in Abe et al. (2018). For the MeV region, searches have also been carried out in KamLAND (Abe et al. 2021b) and Borexino (Agostini et al. 2017).

This paper is focused on the follow-ups of GW triggers detected during the first half of the third observing run (O3a) of LVC, from 2019 April to 2019 September, and described in the GWTC-2 catalog (Abbott et al. 2021). Each GW was classified as a BNS, BBH, or NSBH based on the measured masses of the two objects (m < 3 M e = neutron star, m > 3 M e = black hole, where M e is the solar mass).

This article is organized as follows. Section 2 describes SK and the used data sets. In Section 3, the search method and basic results will be described. Section 4 focuses on the extraction of flux limits and signal significance out of each of individual followup, while Section 5 describes how the results can be combined to constrain different source populations. Section 6 summarizes and concludes the discussion. The data release accompanying this article (doi:10.5281/zenodo.4724822) includes all the figures, the tables of observations and calculated flux limits, and the SK effective area.

Super-Kamiokande and Event Samples

SK (Fukuda et al. 2003) is a water Cerenkov detector located in the Mozumi mine in Gifu Prefecture, Japan. It lies under Mount Ikeno (Ikenoyama) with a total of 2700 m.w.e (meters water equivalent) mean overburden, reducing the cosmic-ray muon rate at the detector by a factor of ∼10 5 with respect to the surface. The detector consists of a cylindrical stainless-steel tank of 39 m diameter and 42 m height, filled with 50 kt of water. It is optically separated into an inner detector (ID) and an outer detector (OD) by a structure at ∼2 m from the wall. The ID is instrumented with 11,129 photomultiplier tubes (PMTs) to observe the Cerenkov light emitted by charged particles produced in neutrino interactions. The OD, instrumented with 1885 PMTs, is primarily used as a veto for external background. SK is sensitive to neutrinos with energies ranging from several MeV to TeV.

The experiment has been operating since 1996, and data taking can be separated into six distinct periods, from SK-I to SK-VI, with the latter starting in 2020 July, being the first phase where gadolinium sulfate has been dissolved into the otherwise pure water. In this paper focused on O3a GW events, only data from SK-V (2019 January-2020 July) were used for analysis, as this is covering the full O3 period.

HE-ν Samples

The high-energy samples correspond to neutrinos with E ν > 100 MeV (which is linked to an electron equivalent energy / visible energy greater than 30 MeV). The neutrino is detected thanks to the outgoing lepton produced in the neutrino charged-current interaction. Data are further divided into three subsamples based on event topology.

The fully contained (FC) and partially contained (PC) samples have a reconstructed neutrino interaction vertex inside the fiducial volume of the ID. 59 The separation between FC and PC is based on the number of effective PMT hits in OD (<16 hits for FC, 16 hits for PC).

The muons entering the detector can originate from muon neutrino interactions in the rock surrounding the detector. As such events would be indistinguishable from downward-going cosmic-ray muons, only events with upward-going direction are considered, hence the name UPMU (for "Upward-going muons"). Events are either through-going (with a requirement on track length >7 m) or stopping in the detector (with a requirement on reconstructed muon momentum >1.6 GeV). Further details are documented in Ashie et al. (2005).

Typical neutrino energy for FC (PC) will be between 0.1 GeV and 10 GeV (1 and 100 GeV), and these samples are sensitive to ν μ , ν e , nm , and ne . The UPMU sample is only sensitive to muon neutrinos and muon antineutrinos, but it covers energies from  (GeV) to  (TeV). The contribution of tau neutrinos is subdominant and was therefore neglected in the following, even though it may improve the final limits in a next iteration of the analysis.

LE-ν Sample

The low-energy sample corresponds to events with energy between 7 and 100 MeV. The largest cross section in this range is the inverse beta decay (IBD) of electron antineutrinos (n e + p → e + + n), and the second most dominant is neutrino elastic scattering (ν + e -→ ν + e -), which is sensitive not only to electron neutrinos but also to other flavors. Interactions on 16 O are neglected in the following analysis.

There are two existing data samples in the SK low-energy analysis. In the 7-15.5 MeV range, the selection tuned for the solar neutrino analysis (Abe et al. 2016b) is applied, while the supernova relic neutrino (SRN) search (Bays et al. 2012) is focused on the 15.5-100 MeV range. The main background below 20 MeV is spallation products from cosmic-ray muons (Zhang et al. 2016); above 20 MeV, it is dominated by interactions from atmospheric neutrinos and electrons from muon decays.

The solar neutrino analysis is in principle sensitive down to 3.5 MeV (Abe et al. 2016b). However, to reduce the background originating from radioactive decays (especially Rn Nakano et al. 2020), only events with reconstructed energy above 7 MeV are considered in this paper.

The expected background is higher than for the HE-ν samples (see the next section), and, as opposed to the latter, the reconstructed neutrino direction cannot be reliably used to identify spatial coincidence with the GW localization, IBD being mostly insensitive to the original direction.

Search Method and Results

The information related to O3a GW triggers is extracted from the FITSfoot_3 files in the data release accompanying Abbott et al. (2021). The main input for the SK analysis is the trigger time t GW : it is used to define a ±500 s time window centered on t GW . The choice of this window is based on the conservative considerations proposed by Baret et al. (2011). The SK data in this window were collected and divided into the four samples (three HE-ν, one LE-ν) described in Section 2. Downtime periods, due to detector calibration or other maintenance (e.g., preparation for Gd-loading in early 2020), prevent the follow-up of some GW triggers. Out of the 39 confirmed events from O3a, SK was able to perform the analysis for 36 of them, with a live time within the 1000 s window of ∼99.5% for each (not 100% because of cosmic muon vetoes and other trigger dead times). Additionally, one of those (GW190512_180714) was not suited for low-energy analysis since there were large noise fluctuations in SK near the GW time due to high-voltage issues.

HE-ν Samples

The events passing the selection described in Section 2.1 were stored. For FC and PC events, the total visible energy E vis is a good estimator of the incoming neutrino energy. For UPMU events, the reconstructed muon momentum p μ is not an accurate estimator because the original neutrino energy cannot be inferred as it interacted in the surrounding rock; it is, however, a lower bound for the neutrino energy. 61For FC and PC, the direction of each event was estimated by reconstructing the Cerenkov rings in the ID, while the direction of UPMU event is determined using the OD hit information. This local direction was converted to equatorial coordinates, R. A. and decl., so that it can easily be compared with GW localization. The associated angular uncertainty was obtained by comparing the reconstructed direction with the true neutrino direction in atmospheric Monte Carlo samples of similar energies. For the lower energies ( n   E (GeV)), the angular resolution is limited by the scattering angle between the neutrino and the lepton (e.g., σ ∼ 20°for FC, E ν = 2 GeV). Resolution of the order of the degree can be achieved with the UPMU sample, as detailed in Hagiwara (2020).

The expected background rate in the high-energy samples was stable over the full data period and therefore can be extracted from data using the full data set from 2019 February to 2020 March. The expected number of background events in a 1000 s time window is 0.112, 0.007, and 0.016 events for FC, PC, and UPMU, respectively (with negligible statistical uncertainties).

The numbers of observed events in the different samples for each individual follow-up are presented in Table 1. Out of the 36 performed follow-ups, 10 of them have associated SK HE-ν events in time coincidence (8 FC, 0 PC, 2 UPMU); this can be compared to the expected background over the 36 GWs: 4.0, 0.3, and 0.6 events, respectively, for FC, PC, and UPMU. For each selected neutrino event, the timing (in particular Δt = t νt GW ), energy, direction, and its related angular uncertainty are provided. The latter information is presented in Table 2, and the angular distributions are shown in Figure 1.

LE-ν Sample

The events within the 1000 s time window passing the selection described in Section 2.2 were extracted. As in the case of HE-ν samples, the total number of observed events is compared to the background expectation. The latter is based on the average event rate computed using the total SK-V data set, which corresponds to 0.729 expected events in 1000 s; this background level has been found to be stable over the whole period. The results for all follow-ups are summarized in Table 1. No significant excess was observed with respect to the expected Poisson statistics, with 24 observed events and 25.0 expected.

Event-by-event Statistical Analysis

Observation Significance

The significance of a given observation in HE-ν samples can be quantified in terms of p-value. The latter can be divided into a temporal component p time that is evaluating the probability to observe at least one SK event in time coincidence with the GW and a spatial component p space comparing the direction of reconstructed neutrinos with the GW localization: p = p time × p space . This discrimination allows separating the discrete time component (due to the low expected background) from the continuous spatial component.

The term p time is simply the Poisson probability to observe at least one event in the selected time window:

We have p time ; 12.6% for n B ; 0.13 (total number of expected events in 1000 s). The term p space is obtained using a maximum likelihood method with the GW localization used as a spatial prior. The best-fit sky position of the potential joint source is obtained by maximizing the loglikelihood ratio, and the obtained test statistic value is compared to its expected distribution from background events to extract p space as the probability for the observation to be compatible with the background-only hypothesis given that at least one SK event in time coincidence has been observed. The method presented in Aartsen et al. (2020) and Hussain et al. (2019) has been adapted to SK.

For each sample (k = FC, PC, or UPMU), the point-source likelihood ( ) k is the number of signal events in the sample (to be fitted), γ is the spectral index of the signal neutrino spectrum ( µ n g dn dE E , to be fitted as well), Ω S is the probed source direction, ( ) n B k is the expected number of background events in the time window, and ( ) N k is the observed number of events. ( )

is the signal probability density function (pdf), which depends on reconstructed event direction x i , reconstructed event energy E i , source spectral index, and direction.

( )

i is the background pdf, which depends solely on event information.

The ( )  k and ( )  k functions were computed and tuned for k = FC, PC, and UPMU independently, using atmospheric neutrino Monte Carlo simulation data sets. They are both written as the product of an angular and an energy component:

where the point-spread function (

S is characterizing the angular resolution of the detector at the considered energy E i (it is maximal for x i = Ω S , i.e., neutrino in the direction of the probed point source), and the energy part

is the convolution of the energy spectrum µ

were obtained using iminuit (Dembinski et al. 2020). 62 The loglikelihood ratio Λ(Ω S ) was then computed:

where ∑ k sums over the three considered samples and

is the spatial prior given directly by the GW sky map.

The test statistic TS was defined by finding the direction in the sky maximizing Λ(Ω) while scanning the full sky, after it has been divided into equal-area pixels using the HEALPix method (Gorski et al. 2005; same pixelization method as used by LVC for GW sky maps):

Finally, this number can be used to compute p space . First, the observation in SK was used to compute TS data . Over 10 5 background toys were generated with neutrino events distributed according to the background distribution (empty toys with zero events are not considered). For each toy, TS was computed and the probability distribution function ( )  TS bkg was obtained and compared to the data value:

The value p space is the probability for the observation to be compatible with the background-only hypothesis given that at least one SK event in time coincidence has been observed. Table 2 presents the obtained p space for the GW triggers with at least one SK event (for the other triggers, we trivially have p = p space = 100%). No significant deviations from the background hypothesis (uniform distribution) were observed. The most significant coincidence is associated with GW190602_175927, with a p-value = p 1.72% space best (p best = 0.22%), corresponding to 2.1σ (2.9σ). However, one needs to take into account the total number of trials involved in the catalog search (N = 10 for p space , as the analysis is restricted to GW with at least one SK event in time coincidence, or N = 36 for p). The trial factor correction is computed by performing 10 5 background-only pseudoexperiments and checking how often one gets

best . This gives post-trial values = P 15.9% space best

(1.0σ) and P best = 7.8% (1.4σ), which are fully consistent with the background-only hypothesis.

Flux Limits Using High-energy Samples

Because no statistically significant event excess was observed within the 1000 s time window in the HE-ν samples, the observation can be converted to an upper limit on the incoming neutrino flux. In the first approach, this was done separately for FC, PC, and UPMU samples, using a similar procedure to that in Abe et al. (2018), while a second approach used the test statistic defined in Section 4.1 to combine the samples.

In the following, the neutrino energy spectrum is assumed to follow a simple power law with a spectral index γ = 2, which is commonly used in such astrophysical searches (e.g., Abe et al. 2018). The neutrino flux can then be written as

In the following, we will report the upper limits on f =

), for the different samples and neutrino flavors (ν μ , nm , ν e , ne ).

Sample-by-sample Approach

For a given sample s, flavor f, and source position Ω, the neutrino flux E dn dE 2 is related to the number of events:

is the SK detector effective area for the selected sample and neutrino flavor, and the integration range is 0.110 5 GeV. The quantity ( ) ( ) W c s f , is the detector acceptance, which takes into account the source direction. To marginalize over the source localization, the following Poisson likelihood is then defined:

, 9

where ( ) n B k and ( ) N k are, respectively, the expected and observed number of events in sample s. One can then derive 90% confidence level (C.L.) upper limits by computing the likelihood as a function of f 0 and finding the 90th percentile, for each sample and flavor (this effectively corresponds to the Bayesian limit with flat prior on f 0 ):

The effective areas have been computed explicitly as a function of neutrino energy and zenith angle, using 500 yr of atmospheric Monte Carlo simulations. There is a very small dependency on the local zenith angle θ for FC and PC, while UPMU covers only efficiently below the horizon (θ > 90°), with a nonnegligible variation with θ, as shown in Figure 2.

The UPMU sample has very limited sensitivity at and above the horizon (0 < θ < 90°); only near-horizontal neutrinos with slightly upgoing muons can be identified. No systematic uncertainties are applied to the detector effective area estimate, as the detector response is relatively stable and well understood and the analysis is strongly dominated by limited statistics.

The full results for ν μ are presented in Figure 3. They show a wide variety of limits: in particular, UPMU upper limits are only reported if the GW localization is mainly below the horizon, where this sample has sensitivity. Detailed numbers for GW190425 (Abbott et al. 2020a) and GW190521 (Abbott et al. 2020b) are presented in Table 3. These two events are illustrating the two scenarios and are the only BNS candidate in GWTC-2 and the heaviest BBH, respectively. Results for all the triggers are given in Table 4.

Combination of the Samples (Using Test Statistic)

As the neutrino spectrum is expected to span the full range from 0.1 to 10 5 GeV, it is worth combining the different samples that have varying sensitivities (in energy, flavor, and direction). The method initially presented in Veske et al. (2020) was implemented using the test statistic defined previously.

Signal simulations were performed, assuming E -2 spectrum and that at most two signal neutrinos are detected in SK; the source direction is chosen randomly based on GW sky map  GW ; and the distribution of signal toy events between the samples is done according to the relative effective areas. As with the background toys in Section 4.1, this allows computing the pdf ( )

for a given number of signal events n S = 0, 1, 2. Assuming that at most two signal neutrinos will be observed for a given GW trigger, the following flux likelihood is defined, based on the observed test statistic and GW sky map:

where ( ) ( )

is the total detector acceptance (summing all samples) assuming E -2 spectrum and the other quantities have already been defined above. The likelihood is composed of a sum of Poisson terms that quantify the relation between number of events and the flux, weighted by the probability to observe the measured test statistic given the different signal hypotheses.

The 90% C.L. upper limit on

is then simply obtained as in Equation (10). The procedure can be repeated independently for each neutrino flavor or also combining flavors, e.g., n n + m m (in the latter case, both ( )

terms and c(Ω) are computed assuming equally distributed flux between the different flavors). The results are presented in Figure 3 and Table 3 for the two examples mentioned above and in Table 4 for all the events.

The combined limits are usually close to the limits obtained by the most constraining individual sample. If the UPMU sample is used (GW localized mainly below the horizon), the combined limit is similar to the UPMU limit. Otherwise, it is consistent with the result of FC+PC. In the case of GW190602_175927, the combined limit is slightly worse than the individual UPMU because of the observed FC event in the same direction as the GW, which gives higher TS data and thus impact ( )  TS k data used in Equation (11).

Flux Limits Using Low-energy Sample

The flux limit calculation for the low-energy sample is similar to HE-ν, except that the effective area is parameterized as in Abe et al. (2018). As there is no direction dependence of the latter and there is only one LE-ν sample, there is no need to define a likelihood in order to perform a combination or to marginalize over the sky like in the HE-ν case.

The upper limit on the total fluence ( ò F = dn dE dE) is then simply computed as

T e 90 90 vis vis

where N 90 is the 90% C.L. upper limit on the number of signal events (calculated from a Poisson distribution), N T is the number of target nuclei in SK fiducial volume, σ is the combined cross section for all interactions, ò is the detection efficiency, λ is the energy density assuming a given spectrum ( ( ) l = F ń dn dE E ), and R is the response function to convert electron or positron energy (E e ) to visible energy in SK (E vis ). The response function and the detection efficiency (ò) are calculated using SK detector Monte Carlo simulations, and related systematic uncertainties are neglected as in the HEν case.

In this analysis, two types of spectra were considered: flat spectrum (l = constant) and Fermi-Dirac spectrum with an average energy of 20 MeV. The results for a selection of triggers are shown in Table 3, while the rest are detailed in Table 4. The limits are more stringent for the ne case, given that the main interaction channel in the detector is IBD of ne , as described in Section 2.

Neutrino Emission Limits and Population Constraints

None of the joint observations have a significance high enough in order to classify them as detection (as presented in Table 2), and the flux limits provided in the previous sections do not directly constrain the physical quantities related to the astrophysical objects. In this section, the neutrino emission at the source is assumed to be isotropic, so that the intrinsic energy E iso emitted by neutrinos from a source at a distance d is directly related to the detected flux at Earth:

Knowing d, one can then constrain E iso , as described in the following.

High-energy Neutrino Emission

If, as in Section 4.2, the E -foot_6 spectrum is assumed, Equation (13) can then be integrated under this particular assumption:

To use the GW sky map as an input, the following likelihood is defined (Veske et al. 2020):

is the conversion factor from E iso to the expected number of signal events in SK for known source distance r and direction Ω. The test statistic distributions ( )

and the measured test statistic ( ) TS i data for trigger i are the same as defined in Section 4.2.2.

( )

is the threedimensional LVC sky map provided for trigger i, taking into account both the direction localization and the distance to the source (see Singer et al. 2016 for details on GW3D localization).

One can derive E iso limits independently for each GW trigger as has been done for the flux limits. For a given flavor (e.g., ν μ ), the obtained limit is on the isotropic energy emitted from the source and that would be detected with this given flavor on Earth (with no assumptions on the flavor distribution). Limits on the total energy emitted by neutrinos of all flavors can be obtained by considering all detectable flavors in SK and assuming equal proportions of them at Earth. This is a reasonable assumption in the most common source scenario, where neutrinos are produced in pion decays in a flavor ratio (ν e : ν μ : ν τ ) equal to (1:2:0), which would become ∼(1: 1: 1) at Earth, after oscillations.

The results are detailed for a selection of triggers in Table 3 and are plotted in Figure 4; the full results are shown in Table 4. In the example of GW190521, the UPMU sample contributed to the observation so that the most constraining limits are obtained for ν μ and nm ; the limit on the total energy emitted in neutrinos assuming equipartition is then dominated by the latter contributions:

all flavors iso,90% . Instead, for GW190425, the limit has similar contributions from all neutrino flavors, as the UPMU sample is not contributing to the limit.

It is worth mentioning that the only BNS in the catalog, GW190425, is located in a sky region for which the observation with the UPMU sample is not possible, as already mentioned in Section 4.2.2. If it had been located in a more favorable region, the upper limit would have improved by a factor of ∼30. This is promising for future observations. Furthermore, if the spectrum is assumed to follow a E -3 spectrum, all the limits presented above are getting less constraining, due to this less favorable spectrum (shifted to lower energies where associated effective areas are smaller), as detailed in Table 4 for the combined all-flavor E iso limits.

The combination of a meaningful set of GW events to constrain further E iso is also worthwhile to infer information about the common physical processes involved in a given source population. This can be performed for different sets of triggers, based on the classification provided in the GW catalog. The relevant categories are BBH, BNS, and NSBH. If emission from all objects of the same nature is assumed to be similar (independently of their individual characteristics), one can define the likelihood:

where the sum runs over the selected GW triggers to be combined.

A more realistic toy scenario would be that the neutrino emission scales with the total mass  tot of the binary system:

One can then use the following likelihood to (in Gev cm -2 ) are presented for the different neutrino flavors, assuming E -2 spectrum. Upper limits on the total energy emitted by the source as neutrinos E iso (in erg) (assuming isotropic emission) are also presented: one limit per flavor and limits for ν e + ne , ν μ + nm , and on the total energy in all flavors assuming equipartition (including unseen tau neutrinos). For LE-ν, limits on the total neutrino fluence Φ (in cm -2 ) are given for ν e , ne , ν x = ν μ + ν τ , ¯¯n n n = + m t x assuming Fermi-Dirac spectrum (with average energy of 20 MeV) and flat spectrum (within the range 7-100 MeV), as well as upper limits on E iso (in erg) for the Fermi-Dirac scenario. constrain f ν :

where f ν , in erg/M e is to be constrained (simplifying the units, f ν can be expressed as the proportion of the total mass converted in neutrinos, e.g., f ν = 10 54 erg/M e = 62%), and

 p i GW tot is the posterior distribution of the total mass of the binary system, as obtained from the LVC data release.

Figure 5 presents the results for the three categories defined above: 1 BNS candidate,foot_7 2 NSBH (GW190426_152155 and GW190814), and 33 BBH (all other events in O3a). The allflavor limit values are indicated on the figures, with the most constraining results obtained for the BBH population: E iso < 4.16 × 10 55 erg assuming that all objects have similar emission. This turns to E iso < 9.73 × 10 56 erg for the E -3 spectrum.

Despite the objects being closer, the BNS and NSBH limits are worse than the ones for BBH because of the limited statistics for these two samples and the fact that the three corresponding GW events have localization above the SK horizon.

Low-energy Neutrino Emission

As for the flux limits, the low-energy case is much simpler. E iso limits are directly obtained by scaling the flux limits using the source distance estimate. In case per-flavor limits are combined, the limit on the total energy emitted in all flavors, assuming equipartition, is, however, dominated by the ne limit.

To cover the distance uncertainty, the following likelihood was defined: where N obs and N bkg are the observed and expected number of LE-ν events, c LE (r) is the conversion factor from E iso to number of signal events assuming Fermi-Dirac spectrum and source at distance r, and p GW (r) is the pdf of distance estimation provided by LIGO-Virgo (Singer et al. 2016). Detailed results for selected triggers are shown in Table 3.

Discussion and Conclusions

The results of the follow-up of LVC O3a GWs with the SK detector have been presented. In the ±500 s time windows centered on the triggers, no excess with respect to the background hypothesis was observed in any of the four considered samples (three for HE-ν, one for LE-ν). Upper limits on the incoming neutrino flux were computed for the different neutrino flavors. For HE-ν, the E -2 spectrum was assumed, while for LE-ν limits, Fermi-Dirac emission with average energy of 20 MeV was considered. In both cases, detailed results are presented in Table 4. Assuming isotropic emissions and equipartition between the different flavors, upper limits on the total energy as neutrinos E iso were derived, both individually for each trigger and by combining the different triggers of the same type, assuming the same emission or that the neutrino emission is scaling with the total mass of the binary system.

For low-energy neutrino emissions, the upper limits on the isotropic energy are not yet constraining enough to probe existing models such as Foucart et al. (2016) ), even though the exact shape of the neutrino spectrum (beyond the assumed For high-energy neutrino emissions, the obtained limits on E iso assuming the E -2 spectrum are barely covering the nonphysical region where the total mass of the binary system is converted to neutrinos ( f ν  10 54 -10 56 erg/M e ; 60%-6000%), while the region currently probed by IceCube is f ν  1% (Aartsen et al. 2020). However, this depends greatly on the assumed spectrum; if the latter happens to be different from the E -2 standard scenario or features a cutoff, the limits would be changed as illustrated in Section 5.1 for the E -3 spectrum. A larger GeV component would favor detection and precise reconstruction of such neutrinos at SK as compared to larger neutrino telescopes like IceCube (Abbasi et al. 2021).

Even though the present paper has focused on the O3a catalog and the analysis was performed offline, the selections and techniques could also be used for real-time follow-up in the O4 observation period and beyond. With these constantly increasing statistics, it may finally be possible to probe the GW+ν source population and better understand the underlying mechanisms.