Search for neutrinos in coincidence with gravitational wave events from the LIGO-Virgo O3a Observing Run with the Super-Kamiokande detector

The Super-Kamiokande detector can be used to search for neutrinos in time coincidence with gravitational waves detected by the LIGO-Virgo Collaboration (LVC). Both low-energy ($7-100$ MeV) and high-energy ($0.1-10^5$ GeV) samples were analyzed in order to cover a very wide neutrino spectrum. Follow-ups of 36 (out of 39) gravitational waves reported in the GWTC-2 catalog were examined; no significant excess above the background was observed, with 10 (24) observed neutrinos compared with 4.8 (25.0) expected events in the high-energy (low-energy) samples. A statistical approach was used to compute the significance of potential coincidences. For each observation, p-values were estimated using neutrino direction and LVC sky map ; the most significant event (GW190602_175927) is associated with a post-trial p-value of $7.8\%$ ($1.4\sigma$). Additionally, flux limits were computed independently for each sample and by combining the samples. The energy emitted as neutrinos by the identified gravitational wave sources was constrained, both for given flavors and for all-flavors assuming equipartition between the different flavors, independently for each trigger and by combining sources of the same nature.


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 (Barthelmy 2021). 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.
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. 2020) experiments are already taking part in such follow-up program for every single reported GW event. Nevertheless, these neutrino telescopes are mainly 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 gravitational wave 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. 2020b). Each GW was classified as a BNS, BBH or NSBH based on the measured masses of the two objects (m < 3 M = Neutron Star, m > 3 M = Black Hole, where M is the solar mass).
This article is organized as follows. Section 2 describes SK and the used datasets. 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 follow-up, 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 (Abe et al. 2021a) 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 Mt. 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 sub-samples 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 1 . 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 GeV and 100 GeV) and these samples are sensitive to ν µ , ν e ,ν µ andν e . The UPMU sample is only sensitive to muon neutrinos and muon antineutrinos, but it covers energies from O(GeV) to O(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 (ν 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 are extracted from the FITS files (NASA/GSFC 2021) in the data release accompanying (Abbott et al. 2020b). 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 vetos 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 2 .
For 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, right ascension (RA) and declination (Dec), 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 (E ν O(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 dataset from 2019 February to 2020 March. The expected number of background events in a 1000-second 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 ), the 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 dataset, 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.

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: p time = p(N ≥ 1) = 1 − e −n B . 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 log-likelihood 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); Hussain et al. (2020) has been adapted to SK.
For each sample (k = FC, PC or UPMU), the point-source likelihood L S , γ; Ω S ) is defined as: where n (k) S is the number of signal events in the sample (to be fitted), γ is the spectral index of the signal neutrino spectrum (dn/dE ν ∝ E −γ , to be fitted as well), Ω S is the probed source direction, n (k) B is the expected number of Table 1. Summary of all GW triggers from the O3a observation run. The first four columns summarize GW information (name, time, event type, and mean distance), the fifth column corresponds to SK live time in the 1000 s time window, and the last columns present the observed number of events in the four SK samples. b The detector was not taking data due to calibrations or maintenance. showing the distribution of SK events, superimposed with the GW probability distribution, for the ten GW triggers with one observed event in the O3a observation run. The region representing the 1σ angular resolution is indicated in blue for FC and in green for UPMU. The dark red contour shows the 90% containment of the GW probability. The shaded area shows the sky region that is below SK horizon (where the UPMU sample is sensitive). background events in the time window, and N (k) is the observed number of events.
is the background pdf, which depends solely on event information. The S (k) and B (k) functions were computed and tuned for k = FC, PC, UPMU independently, using atmospheric neutrino Monte Carlo simulation datasets. They are both written as the product of an angular and an energy component: and S , γ; Ω S ) were obtained using iminuit (Dembinski et al. 2020) 3 . The log-likelihood ratio Λ(Ω S ) was then computed: where k sums over the three considered samples and P GW (Ω S ) is the spatial prior given directly by the GW sky map. The test statistic T S 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 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 T S 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, T S was computed and the probability distribution function P bkg (T S) 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 best space = 1.72% (p best = 0.22%), corresponding to 2.1σ (respectively 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 pseudo-experiments and checking how often one gets min{p i } i=1...N < p best . This gives post-trial values P best space = 15.9% (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, that 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 φ 0 = E 2 ν dn/dE ν [in GeV cm −2 ], for the different samples and neutrino flavors (ν µ ,ν µ , ν e ,ν e ).

Sample-by-sample approach
For a given sample s, flavor f and source position Ω, the neutrino flux E 2 dn/dE is related to the number of events: where A (s,f ) eff (E ν , Ω) is the SK detector effective area for the selected sample and neutrino flavor, the integration range is 0.1-10 5 GeV. The quantity 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: where n (k) B 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 φ 0 and finding the 90% percentile, for each sample and flavor (this effectively corresponds to the Bayesian limit with flat prior on φ 0 ): The effective areas have been computed explicitly as a function of neutrino energy and zenith angle, using 500 years 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. 2020c) 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 GeV 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 P 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 P n S (T S) 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 c(Ω) = s c (s,f ) (Ω) 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 φ 0 = E 2 ν dn/dE ν is then simply obtained as in Equation 10. The procedure can be repeated independently for each neutrino flavor or also combining flavors, e.g., ν µ +ν µ (in the latter case, both P n S (T S) 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 T S data and thus impact P k (T S data ) used in the 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. Table 3. Obtained 90% C.L. upper limits on E 2 dn/dE for GW190425 and GW190521. For HE-ν, limits on E 2 dn/dE [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 Eiso [in erg] (assuming isotropic emission) are also presented: one limit per flavor and limits for νe +νe, νµ +νµ 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,νe, νx = νµ + ντ ,ν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 Eiso [in erg] for the Fermi-Dirac scenario.

NEUTRINO EMISSION LIMITS AND POPULATION CONSTRAINTS
None of the joint observations has 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, E −2 spectrum is assumed, the 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): The quantity c (r, Ω) 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 P GW (r, Ω) is the three-dimensional 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 GW 3D localization).
One can derive E iso limits independently for each GW trigger as it 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ν µ ; the limit on the total energy emitted in neutrinos assuming equipartition is then dominated by the latter contributions: E all−flavors iso,90% 3 × E νµ+νµ 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 ∼ 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.  . The limits are following two lines E 90% iso ∝ distance 2 based on geometrical considerations, one of the lines shows events dominated by UPMU νν /νµ contributions (giving more stringent limits) while the other line contains GW triggers that are less constrained. The two GW used in Table 3 are labelled in the plots. The complete figure set (5 images, one per considered neutrino flavor + all-flavors) is available in the online journal.
A more realistic toy scenario would be that the neutrino emission scales with the total mass M tot of the binary system: E ν iso = f ν × M tot . One can then use the following likelihood to constrain f ν : where f ν , in erg/M 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 = 62%), and p GW (M (i) 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 4 , 2 NSBH (GW190426 152155 and GW190814), 33 BBH (all other events in O3a). The all-flavor 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 all objects have similar emission. This turns to E iso < 9.73 × 10 56 erg for 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 of 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ν e limit. To cover the distance uncertainty, the following likelihood was defined:  Figure 5. 90% C.L. upper limits on the isotropic energy emitted in neutrinos by combining GW triggers with the same nature, for νµ, νµ +νµ, νe +νe and all-flavor emission (assuming equipartition). The left panel shows the results assuming that all selected sources are emitting the same Eiso while the right panel is assuming that neutrino emission is scaling with the total mass of the binary system.
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, p GW (r) is the p.d.f. 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 gravitational waves 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-ν, 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 the 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) (predicted luminosity L model iso ∼ 4-7 × 10 53 erg s −1 ), even though the exact shape of the neutrino spectrum (beyond the assumed simple Fermi-Dirac distribution with E ν = 20 MeV) may modify slightly the obtained upper limits.
For high-energy neutrino emissions, the obtained limits on E iso assuming 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 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.