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

K. Abe  Kamioka Observatory, Institute for Cosmic Ray Research, University of
Tokyo, Kamioka, Gifu 506-1205, Japan  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), The University of Tokyo Institutes for
Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583, Japan  C.
Bronner  Kamioka Observatory, Institute for Cosmic Ray Research, University of
Tokyo, Kamioka, Gifu 506-1205, Japan  Y. Hayato  Kamioka Observatory,
Institute for Cosmic Ray Research, University of Tokyo, Kamioka, Gifu
506-1205, Japan  Kavli Institute for the Physics and Mathematics of the
Universe (WPI), The University of Tokyo Institutes for Advanced Study,
University of Tokyo, Kashiwa, Chiba 277-8583, Japan  M. Ikeda  S. Imaizumi
Kamioka Observatory, Institute for Cosmic Ray Research, University of Tokyo,
Kamioka, Gifu 506-1205, Japan  J. Kameda  Kamioka Observatory, Institute for
Cosmic Ray Research, University of Tokyo, Kamioka, Gifu 506-1205, Japan  Kavli
Institute for the Physics and Mathematics of the Universe (WPI), The
University of Tokyo Institutes for Advanced Study, University of Tokyo,
Kashiwa, Chiba 277-8583, Japan  Y. Kanemura  Y. Kataoka  S. Miki  Kamioka
Observatory, Institute for Cosmic Ray Research, University of Tokyo, Kamioka,
Gifu 506-1205, Japan  M. Miura  S. Moriyama Kamioka Observatory, Institute for
Cosmic Ray Research, University of Tokyo, Kamioka, Gifu 506-1205, Japan  Kavli
Institute for the Physics and Mathematics of the Universe (WPI), The
University of Tokyo Institutes for Advanced Study, University of Tokyo,
Kashiwa, Chiba 277-8583, Japan  Y. Nagao  Kamioka Observatory, Institute for
Cosmic Ray Research, University of Tokyo, Kamioka, Gifu 506-1205, Japan  M.
Nakahata  Kamioka Observatory, Institute for Cosmic Ray Research, University
of Tokyo, Kamioka, Gifu 506-1205, Japan  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), The University of Tokyo Institutes for
Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583, Japan  S.
Nakayama  Kamioka Observatory, Institute for Cosmic Ray Research, University
of Tokyo, Kamioka, Gifu 506-1205, Japan  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), The University of Tokyo Institutes for
Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583, Japan  T. Okada
K. Okamoto  A. Orii  G. Pronost Kamioka Observatory, Institute for Cosmic Ray
Research, University of Tokyo, Kamioka, Gifu 506-1205, Japan  H. Sekiya  M.
Shiozawa  Kamioka Observatory, Institute for Cosmic Ray Research, University
of Tokyo, Kamioka, Gifu 506-1205, Japan  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), The University of Tokyo Institutes for
Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583, Japan  Y. Sonoda
Y. Suzuki Kamioka Observatory, Institute for Cosmic Ray Research, University
of Tokyo, Kamioka, Gifu 506-1205, Japan  A. Takeda  Kamioka Observatory,
Institute for Cosmic Ray Research, University of Tokyo, Kamioka, Gifu
506-1205, Japan  Kavli Institute for the Physics and Mathematics of the
Universe (WPI), The University of Tokyo Institutes for Advanced Study,
University of Tokyo, Kashiwa, Chiba 277-8583, Japan  Y. Takemoto  A. Takenaka
Kamioka Observatory, Institute for Cosmic Ray Research, University of Tokyo,
Kamioka, Gifu 506-1205, Japan  H. Tanaka  Kamioka Observatory, Institute for
Cosmic Ray Research, University of Tokyo, Kamioka, Gifu 506-1205, Japan  S.
Watanabe  T. Yano  Kamioka Observatory, Institute for Cosmic Ray Research,
University of Tokyo, Kamioka, Gifu 506-1205, Japan  S. Han  Research Center
for Cosmic Neutrinos, Institute for Cosmic Ray Research, University of Tokyo,
Kashiwa, Chiba 277-8582, Japan  T. Kajita  Research Center for Cosmic
Neutrinos, Institute for Cosmic Ray Research, University of Tokyo, Kashiwa,
Chiba 277-8582, Japan  Kavli Institute for the Physics and Mathematics of the
Universe (WPI), The University of Tokyo Institutes for Advanced Study,
University of Tokyo, Kashiwa, Chiba 277-8583, Japan  K. Okumura  Research
Center for Cosmic Neutrinos, Institute for Cosmic Ray Research, University of
Tokyo, Kashiwa, Chiba 277-8582, Japan  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), The University of Tokyo Institutes for
Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583, Japan  T.
Tashiro  R. Wang  J. Xia  Research Center for Cosmic Neutrinos, Institute for
Cosmic Ray Research, University of Tokyo, Kashiwa, Chiba 277-8582, Japan  G.
D. Megias  Institute for Cosmic Ray Research, University of Tokyo, Kashiwa,
Chiba 277-8582, Japan  D. Bravo-Berguño  L. Labarga  Ll. Marti  B. Zaldivar
Department of Theoretical Physics, University Autonoma Madrid, 28049 Madrid,
Spain  B. W. Pointon  Department of Physics, British Columbia Institute of
Technology, Burnaby, BC, V5G 3H2, Canada  F. d. M. Blaszczyk  Department of
Physics, Boston University, Boston, MA 02215, USA  E. Kearns  Department of
Physics, Boston University, Boston, MA 02215, USA  Kavli Institute for the
Physics and Mathematics of the Universe (WPI), The University of Tokyo
Institutes for Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583,
Japan  J. L. Raaf  Department of Physics, Boston University, Boston, MA 02215,
USA  J. L. Stone  Department of Physics, Boston University, Boston, MA 02215,
USA  Kavli Institute for the Physics and Mathematics of the Universe (WPI),
The University of Tokyo Institutes for Advanced Study, University of Tokyo,
Kashiwa, Chiba 277-8583, Japan  L. Wan  Department of Physics, Boston
University, Boston, MA 02215, USA  T. Wester  Department of Physics, Boston
University, Boston, MA 02215, USA  J. Bian  N. J. Griskevich W. R. Kropp  S.
Locke  S. Mine  Department of Physics and Astronomy, University of California,
Irvine, Irvine, CA 92697-4575, USA  M. B. Smy  H. W. Sobel  Department of
Physics and Astronomy, University of California, Irvine, Irvine, CA
92697-4575, USA  Kavli Institute for the Physics and Mathematics of the
Universe (WPI), The University of Tokyo Institutes for Advanced Study,
University of Tokyo, Kashiwa, Chiba 277-8583, Japan  V. Takhistov  Department
of Physics and Astronomy, University of California, Irvine, Irvine, CA
92697-4575, USA  Kavli Institute for the Physics and Mathematics of the
Universe (WPI), The University of Tokyo Institutes for Advanced Study,
University of Tokyo, Kashiwa, Chiba 277-8583, Japan  P. Weatherly  Department
of Physics and Astronomy, University of California, Irvine, Irvine, CA
92697-4575, USA  J. Hill  Department of Physics, California State University,
Dominguez Hills, Carson, CA 90747, USA  J. Y. Kim  I. T. Lim  R. G. Park
Institute for Universe and Elementary Particles, Chonnam National University,
Gwangju 61186, Korea  B. Bodur  Department of Physics, Duke University, Durham
NC 27708, USA  K. Scholberg C. W. Walter Department of Physics, Duke
University, Durham NC 27708, USA  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), The University of Tokyo Institutes for
Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583, Japan  L.
Bernard  A. Coffani  O. Drapier  S. El Hedri  A. Giampaolo  M. Gonin  Th. A.
Mueller  P. Paganini  B. Quilain  Ecole Polytechnique, IN2P3-CNRS, Laboratoire
Leprince-Ringuet, F-91120 Palaiseau, France  T. Ishizuka  Junior College,
Fukuoka Institute of Technology, Fukuoka, Fukuoka 811-0295, Japan  T. Nakamura
Department of Physics, Gifu University, Gifu, Gifu 501-1193, Japan  J. S. Jang
GIST College, Gwangju Institute of Science and Technology, Gwangju 500-712,
Korea  J. G. Learned  Department of Physics and Astronomy, University of
Hawaii, Honolulu, HI 96822, USA  L. H. V. Anthony  D. G. R. Martin  A. A.
Sztuc  Y. Uchida  Department of Physics, Imperial College London , London, SW7
2AZ, United Kingdom  V. Berardi  M. G. Catanesi E. Radicioni  Dipartimento
Interuniversitario di Fisica, INFN Sezione di Bari and Università e
Politecnico di Bari, I-70125, Bari, Italy  N. F. Calabria  L. N. Machado  G.
De Rosa  Dipartimento di Fisica, INFN Sezione di Napoli and Università di
Napoli, I-80126, Napoli, Italy  G. Collazuol  F. Iacob  M. Lamoureux N. Ospina
Dipartimento di Fisica, INFN Sezione di Padova and Università di Padova,
I-35131, Padova, Italy  L. Ludovici  INFN Sezione di Roma and Università di
Roma “La Sapienza”, I-00185, Roma, Italy  Y. Maekawa  Y. Nishimura  Department
of Physics, Keio University, Yokohama, Kanagawa, 223-8522, Japan  S. Cao  M.
Friend  T. Hasegawa  T. Ishida  M. Jakkapu  T. Kobayashi  T. Matsubara  T.
Nakadaira  High Energy Accelerator Research Organization (KEK), Tsukuba,
Ibaraki 305-0801, Japan  K. Nakamura  High Energy Accelerator Research
Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan  Kavli Institute for the
Physics and Mathematics of the Universe (WPI), The University of Tokyo
Institutes for Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583,
Japan  Y. Oyama  K. Sakashita  T. Sekiguchi  T. Tsukamoto  High Energy
Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan  Y.
Kotsar  Y. Nakano H. Ozaki  T. Shiozawa  Department of Physics, Kobe
University, Kobe, Hyogo 657-8501, Japan  A. T. Suzuki  Department of Physics,
Kobe University, Kobe, Hyogo 657-8501, Japan  Y. Takeuchi Department of
Physics, Kobe University, Kobe, Hyogo 657-8501, Japan  Kavli Institute for the
Physics and Mathematics of the Universe (WPI), The University of Tokyo
Institutes for Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583,
Japan  S. Yamamoto  Department of Physics, Kobe University, Kobe, Hyogo
657-8501, Japan  A. Ali  Y. Ashida  J. Feng  S. Hirota  T. Kikawa  M. Mori
Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan  T.
Nakaya  Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The
University of Tokyo Institutes for Advanced Study, University of Tokyo,
Kashiwa, Chiba 277-8583, Japan  R. A. Wendell Department of Physics, Kyoto
University, Kyoto, Kyoto 606-8502, Japan  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), The University of Tokyo Institutes for
Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583, Japan  K.
Yasutome  Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502,
Japan  P. Fernandez  N. McCauley  P. Mehta  A. Pritchard  K. M. Tsui
Department of Physics, University of Liverpool, Liverpool, L69 7ZE, United
Kingdom  Y. Fukuda  Department of Physics, Miyagi University of Education,
Sendai, Miyagi 980-0845, Japan  Y. Itow  Institute for Space-Earth
Environmental Research, Nagoya University, Nagoya, Aichi 464-8602, Japan
Kobayashi-Maskawa Institute for the Origin of Particles and the Universe,
Nagoya University, Nagoya, Aichi 464-8602, Japan  H. Menjo  T. Niwa  K. Sato
Institute for Space-Earth Environmental Research, Nagoya University, Nagoya,
Aichi 464-8602, Japan  M. Tsukada  Institute for Space-Earth Environmental
Research, Nagoya University, Nagoya, Aichi 464-8602, Japan  P. Mijakowski
National Centre For Nuclear Research, 02-093 Warsaw, Poland  J. Jiang  C. K.
Jung  C. Vilela  M. J. Wilking  C. Yanagisawa  also at BMCC/CUNY, Science
Department, New York, New York, 1007, USA.  Department of Physics and
Astronomy, State University of New York at Stony Brook, NY 11794-3800, USA  K.
Hagiwara  M. Harada  T. Horai  H. Ishino S. Ito  Department of Physics,
Okayama University, Okayama, Okayama 700-8530, Japan  Y. Koshio  Department of
Physics, Okayama University, Okayama, Okayama 700-8530, Japan  Kavli Institute
for the Physics and Mathematics of the Universe (WPI), The University of Tokyo
Institutes for Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583,
Japan  H. Kitagawa  W. Ma  N. Piplani  S. Sakai  Department of Physics,
Okayama University, Okayama, Okayama 700-8530, Japan  Y. Kuno  Department of
Physics, Osaka University, Toyonaka, Osaka 560-0043, Japan  G. Barr  D. Barrow
Department of Physics, Oxford University, Oxford, OX1 3PU, United Kingdom  L.
Cook  Department of Physics, Oxford University, Oxford, OX1 3PU, United
Kingdom  Kavli Institute for the Physics and Mathematics of the Universe
(WPI), The University of Tokyo Institutes for Advanced Study, University of
Tokyo, Kashiwa, Chiba 277-8583, Japan  A. Goldsack  Department of Physics,
Oxford University, Oxford, OX1 3PU, United Kingdom  Kavli Institute for the
Physics and Mathematics of the Universe (WPI), The University of Tokyo
Institutes for Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583,
Japan  S. Samani  Department of Physics, Oxford University, Oxford, OX1 3PU,
United Kingdom  C. Simpson  D. Wark  Department of Physics, Oxford University,
Oxford, OX1 3PU, United Kingdom  STFC, Rutherford Appleton Laboratory, Harwell
Oxford, and Daresbury Laboratory, Warrington, OX11 0QX, United Kingdom  F.
Nova Rutherford Appleton Laboratory, Harwell, Oxford, OX11 0QX, UK  T. Boschi
F. Di Lodovico J. Migenda  S. Molina Sedgwick  now at IFIC (CSIC - U.
Valencia), Spain.  M. Taani  S. Zsoldos  Department of Physics, King’s College
London, London, WC2R 2LS, UK  J. Y. Yang  Department of Physics, Seoul
National University, Seoul 151-742, Korea  S. J. Jenkins  M. Malek  J. M.
McElwee O. Stone  M. D. Thiesse  L. F. Thompson  Department of Physics and
Astronomy, University of Sheffield, S3 7RH, Sheffield, United Kingdom  H.
Okazawa  Department of Informatics in Social Welfare, Shizuoka University of
Welfare, Yaizu, Shizuoka, 425-8611, Japan  S. B. Kim  I. Yu  Department of
Physics, Sungkyunkwan University, Suwon 440-746, Korea  K. Nishijima
Department of Physics, Tokai University, Hiratsuka, Kanagawa 259-1292, Japan
M. Koshiba  Deceased.  The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan
K. Iwamoto  Department of Physics, University of Tokyo, Bunkyo, Tokyo
113-0033, Japan  Y. Nakajima  Department of Physics, University of Tokyo,
Bunkyo, Tokyo 113-0033, Japan  Kavli Institute for the Physics and Mathematics
of the Universe (WPI), The University of Tokyo Institutes for Advanced Study,
University of Tokyo, Kashiwa, Chiba 277-8583, Japan  N. Ogawa  Department of
Physics, University of Tokyo, Bunkyo, Tokyo 113-0033, Japan  M. Yokoyama
Department of Physics, University of Tokyo, Bunkyo, Tokyo 113-0033, Japan
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The
University of Tokyo Institutes for Advanced Study, University of Tokyo,
Kashiwa, Chiba 277-8583, Japan  K. Martens  Kavli Institute for the Physics
and Mathematics of the Universe (WPI), The University of Tokyo Institutes for
Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583, Japan  M. R.
Vagins  Kavli Institute for the Physics and Mathematics of the Universe (WPI),
The University of Tokyo Institutes for Advanced Study, University of Tokyo,
Kashiwa, Chiba 277-8583, Japan  Department of Physics and Astronomy,
University of California, Irvine, Irvine, CA 92697-4575, USA  S. Izumiyama  M.
Kuze M. Tanaka  T. Yoshida  Department of Physics,Tokyo Institute of
Technology, Meguro, Tokyo 152-8551, Japan  M. Inomoto  M. Ishitsuka  H. Ito R.
Matsumoto  K. Ohta  M. Shinoki  Department of Physics, Faculty of Science and
Technology, Tokyo University of Science, Noda, Chiba 278-8510, Japan  J. F.
Martin  H. A. Tanaka  T. Towstego  Department of Physics, University of
Toronto, ON, M5S 1A7, Canada  R. Akutsu  M. Hartz  A. Konaka  P. de Perio  N.
W. Prouse  TRIUMF, 4004 Wesbrook Mall, Vancouver, BC, V6T2A3, Canada  S. Chen
B. D. Xu  Department of Engineering Physics, Tsinghua University, Beijing,
100084, China  M. Posiadala-Zezula  Faculty of Physics, University of Warsaw,
Warsaw, 02-093, Poland  D. Hadley  B. Richards  Department of Physics,
University of Warwick, Coventry, CV4 7AL, UK  B. Jamieson  J. Walker
Department of Physics, University of Winnipeg, MB R3J 3L8, Canada  A. Minamino
K. Okamoto  G. Pintaudi  S. Sano  R. Sasaki  Department of Physics, Yokohama
National University, Yokohama, Kanagawa, 240-8501, Japan  A. K. Ichikawa  K.
Nakamura  Department of Physics, Faculty of Science, Tohoku University,
Sendai, Miyagi, 980-8578, Japan

(Received April 20, 2021; Revised June 15, 2021; Accepted June 16, 2021)

###### Abstract

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 ($7100\text{\,}\mathrm{MeV}$) and high-energy
($0.1{10}^{5}\text{\,}\mathrm{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.

††journal: ApJ

September 13, 2021

##  1 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.

The detected GW emitters are categorized by LVC into several types, for which
high-energy neutrino (HE-$\nu$) 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 $\mathrm{MeV}$ neutrinos
(LE-$\nu$) (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., 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-$\nu$ above $100\text{\,}\mathrm{GeV}$. Super-Kamiokande (SK) can
complement such follow-up studies, as, despite its much smaller size, it is
sensitive to lower energies ($\mathrm{MeV}$-$\mathrm{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 $\mathrm{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\text{\,}\mathrm{M_{\odot}}$$ =
Neutron Star, $m>$3\text{\,}\mathrm{M_{\odot}}$$ = Black Hole, where
$\mathrm{M_{\odot}}$ 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.

##  2 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\text{\,}\mathrm{m.w.e.}$ (meters water equivalent) mean
overburden, reducing the cosmic-ray muon rate at the detector by a factor of
$\sim 10^{5}$ with respect to the surface. The detector consists of a
cylindrical stainless-steel tank of $39\text{\,}\mathrm{m}$ diameter and
$42\text{\,}\mathrm{m}$ height, filled with $50\text{\,}\mathrm{kt}$ of water.
It is optically separated into an inner detector (ID) and an outer detector
(OD) by a structure at $\sim$2\text{\,}\mathrm{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
$\mathrm{MeV}$ to $\mathrm{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.

###  2.1 HE-$\nu$ samples

The high-energy samples correspond to neutrinos with
$E_{\nu}>$100\text{\,}\mathrm{MeV}$$ (which is linked to an electron
equivalent energy / visible energy greater than $30\text{\,}\mathrm{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
ID111The fiducial volume for this analysis is defined as the region in the ID
more than $1\text{\,}\mathrm{m}$ ($2\text{\,}\mathrm{m}$ for PC) from any
wall.. The separation between FC and PC is based on the number of effective
PMT hits in OD ($<16$ hits for FC, $\geq 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\text{\,}\mathrm{m}$$) or stopping in the detector (with a requirement on
reconstructed muon momentum $>$1.6\text{\,}\mathrm{GeV}$$). Further details
are documented in Ashie et al. (2005).

Typical neutrino energy for FC (PC) will be between
$0.1\text{\,}\mathrm{GeV}10\text{\,}\mathrm{GeV}$
($1\text{\,}\mathrm{GeV}100\text{\,}\mathrm{GeV}$) and these samples are
sensitive to $\nu_{\mu}$, $\nu_{e}$, $\bar{\nu}_{\mu}$ and $\bar{\nu}_{e}$.
The UPMU sample is only sensitive to muon neutrinos and muon antineutrinos,
but it covers energies from $\mathcal{O}($\mathrm{GeV}$)$ to
$\mathcal{O}($\mathrm{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.

###  2.2 LE-$\nu$ sample

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

There are two existing data samples in the SK low-energy analysis. In the
$715.5\text{\,}\mathrm{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.5100\text{\,}\mathrm{MeV}$ range. The main background below
$20\text{\,}\mathrm{MeV}$ is spallation products from cosmic-ray muons (Zhang
et al., 2016) ; above $20\text{\,}\mathrm{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\text{\,}\mathrm{MeV}$ (Abe et al., 2016b). However, to reduce the
background originating from radioactive decays (especially $\mathrm{Rn}$
Nakano et al. (2020)), only events with reconstructed energy above
$7\text{\,}\mathrm{MeV}$ are considered in this paper.

The expected background is higher than for the HE-$\nu$ 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.

##  3 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_{\rm GW}$: it is used to
define a $\pm$500\text{\,}\mathrm{s}$$ time window centered on $t_{\rm 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-$\nu$, one LE-$\nu$) 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\text{\,}\mathrm{s}$
window of $\sim 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.

###  3.1 HE-$\nu$ samples

The events passing the selection described in section 2.1 were stored. For FC
and PC events, the total visible energy $E_{\rm vis}$ is a good estimator of
the incoming neutrino energy. For UPMU events, the reconstructed muon momentum
$p_{\mu}$ 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 energy222In the following, we will use the
notation $E_{\rm reco}$ to refer to $E_{\rm vis}$ and $p_{\mu}$ for FC/PC and
UPMU respectively, with the related caveats..

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_{\nu}\lesssim\mathcal{O}($\mathrm{GeV}$)$), the angular
resolution is limited by the scattering angle between the neutrino and the
lepton (e.g. $\sigma\sim$20\text{\,}\mathrm{\SIUnitSymbolDegree}$$ for FC,
$E_{\nu}=$2\text{\,}\mathrm{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-$\nu$ 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 $\Delta t=t_{\nu}-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.

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\text{\,}\mathrm{s}$ time window, and the last columns present the observed number of events in the four SK samples. Trigger name | Alert time | GW Type | Distance | Live time | HE-$\nu$ | LE-$\nu$
---|---|---|---|---|---|---
| (UTC) |  | $\mathrm{Mpc}$ | seconds | FC | PC | UPMU |
GW190408_181802 | 2019-04-08 18:18:02 | BBH | 1547.5 | 993 | 0 | 0 | 0 | 3
GW190412 | 2019-04-12 05:30:44 | BBH | 734.1 | 993 | 0 | 0 | 0 | 0
GW190413_052954 | 2019-04-13 05:29:54 | BBH | 4189.6 | 993 | 0 | 0 | 0 | 0
GW190413_134308 | 2019-04-13 13:43:08 | BBH | 5181.6 | 993 | 0 | 0 | 0 | 0
GW190421_213856 | 2019-04-21 21:38:56 | BBH | 3165.5 | 993 | 0 | 0 | 0 | 3
GW190424_180648 | 2019-04-24 18:06:48 | BBH | 2568.4 | 993 | 1 | 0 | 0 | 1
GW190425 | 2019-04-25 08:18:05 | BNS | 156.8 | 993 | 0 | 0 | 0 | 1
GW190426_152155 | 2019-04-26 15:21:55 | NSBH | 377.2 | 993 | 0 | 0 | 1 | 0
GW190503_185404 | 2019-05-03 18:54:04 | BBH | 1527.3 | 993 | 0 | 0 | 0 | 0
GW190512_180714aaThe low-energy sample is not used because of high-voltage issues. | 2019-05-12 18:07:14 | BBH | 1462.5 | 994 | 0 | 0 | 0 | -
GW190513_205428 | 2019-05-13 20:54:28 | BBH | 2189.7 | 994 | 1 | 0 | 0 | 0
GW190514_065416 | 2019-05-14 06:54:16 | BBH | 4987.6 | 994 | 0 | 0 | 0 | 1
GW190517_055101bbThe detector was not taking data due to calibrations or maintenance. | 2019-05-17 05:51:01 | BBH | 2270.5 | 0 | - | - | - | -
GW190519_153544 | 2019-05-19 15:35:44 | BBH | 3023.5 | 994 | 0 | 0 | 0 | 1
GW190521 | 2019-05-21 03:02:29 | BBH | 4566.9 | 994 | 0 | 0 | 0 | 3
GW190521_074359 | 2019-05-21 07:43:59 | BBH | 1244.2 | 994 | 0 | 0 | 0 | 0
GW190527_092055 | 2019-05-27 09:20:55 | BBH | 3562.9 | 994 | 1 | 0 | 0 | 1
GW190602_175927 | 2019-06-02 17:59:27 | BBH | 3138.1 | 994 | 1 | 0 | 0 | 0
GW190620_030421 | 2019-06-20 03:04:21 | BBH | 3210.9 | 995 | 0 | 0 | 1 | 1
GW190630_185205 | 2019-06-30 18:52:05 | BBH | 956.2 | 992 | 0 | 0 | 0 | 2
GW190701_203306 | 2019-07-01 20:33:06 | BBH | 2152.4 | 992 | 0 | 0 | 0 | 0
GW190706_222641 | 2019-07-06 22:26:41 | BBH | 5184.0 | 992 | 0 | 0 | 0 | 2
GW190707_093326 | 2019-07-07 09:33:26 | BBH | 790.8 | 992 | 0 | 0 | 0 | 0
GW190708_232457 | 2019-07-08 23:24:57 | BBH | 887.9 | 993 | 0 | 0 | 0 | 0
GW190719_215514 | 2019-07-19 21:55:14 | BBH | 4786.3 | 993 | 0 | 0 | 0 | 1
GW190720_000836 | 2019-07-20 00:08:36 | BBH | 906.0 | 993 | 0 | 0 | 0 | 0
GW190727_060333 | 2019-07-27 06:03:34 | BBH | 3608.9 | 992 | 0 | 0 | 0 | 0
GW190728_064510 | 2019-07-28 06:45:10 | BBH | 857.6 | 993 | 1 | 0 | 0 | 2
GW190731_140936 | 2019-07-31 14:09:36 | BBH | 4033.7 | 993 | 0 | 0 | 0 | 1
GW190803_022701 | 2019-08-03 02:27:01 | BBH | 3749.6 | 993 | 0 | 0 | 0 | 0
GW190814 | 2019-08-14 21:10:38 | NSBH | 240.7 | 994 | 1 | 0 | 0 | 0
GW190828_063405 | 2019-08-28 06:34:05 | BBH | 2160.3 | 542 | 0 | 0 | 0 | 1
GW190828_065509bbThe detector was not taking data due to calibrations or maintenance. | 2019-08-28 06:55:09 | BBH | 1657.8 | 0 | - | - | - | -
GW190909_114149 | 2019-09-09 11:41:49 | BBH | 4923.7 | 994 | 0 | 0 | 0 | 0
GW190910_112807 | 2019-09-10 11:28:07 | BBH | 1670.1 | 994 | 1 | 0 | 0 | 2
GW190915_235702bbThe detector was not taking data due to calibrations or maintenance. | 2019-09-15 23:57:02 | BBH | 1714.6 | 0 | - | - | - | -
GW190924_021846 | 2019-09-24 02:18:46 | BBH | 572.4 | 994 | 1 | 0 | 0 | 0
GW190929_012149 | 2019-09-29 01:21:49 | BBH | 3901.5 | 994 | 0 | 0 | 0 | 0
GW190930_133541 | 2019-09-30 13:35:41 | BBH | 785.9 | 994 | 0 | 0 | 0 | 0
Table 2: List of selected SK HE-$\nu$ events in time coincidence with GW triggers from O3a observation run. The p-values obtained with the statistical method presented in section 4.1 are also listed.  Trigger name | SK sample | $\Delta T$ | $E_{\rm reco}$ | RA | Dec | $\sigma_{\rm ang}$ | $p_{\rm space}$ | $p$
---|---|---|---|---|---|---|---|---
|  | seconds | GeV | degree | degree | degree | % | %
GW190424_180648 | FC | $104.03$ | $0.57$ | $210.82$ | $-58.74$ | $52.08$ | $48.55$ | $6.12$
GW190426_152155 | UPMU | $278.99$ | $9.52$ | $352.37$ | $-8.46$ | $2.15$ | $100.00$ | $12.59$
GW190513_205428 | FC | $-183.27$ | $0.68$ | $279.34$ | $-37.27$ | $41.19$ | $8.59$ | $1.08$
GW190527_092055 | FC | $248.41$ | $0.48$ | $54.09$ | $18.80$ | $52.08$ | $58.93$ | $7.43$
GW190602_175927 | FC | $-286.52$ | $2.75$ | $93.67$ | $-38.90$ | $16.22$ | $1.72$ | $0.22$
GW190620_030421 | UPMU | $-327.70$ | $2.33$ | $177.69$ | $-35.59$ | $8.04$ | $100.00$ | $12.62$
GW190728_064510 | FC | $102.99$ | $0.19$ | $300.45$ | $29.71$ | $92.51$ | $21.02$ | $2.65$
GW190814 | FC | $250.36$ | $1.21$ | $157.59$ | $-9.47$ | $28.26$ | $100.00$ | $12.61$
GW190910_112807 | FC | $301.42$ | $1.08$ | $160.13$ | $-22.70$ | $32.09$ | $57.11$ | $7.20$
GW190924_021846 | FC | $411.87$ | $0.30$ | $281.38$ | $-54.52$ | $73.58$ | $50.49$ | $6.37$
Figure 1: Sky maps (in equatorial coordinates) 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\sigma$ 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).

###  3.2 LE-$\nu$ sample

The events within the $1000\text{\,}\mathrm{s}$ time window passing the
selection described in section 2.2 were extracted. As in the case of HE-$\nu$
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\text{\,}\mathrm{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.

##  4 Event-by-event statistical analysis

###  4.1 Observation significance

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

The term $p_{\rm time}$ is simply the Poisson probability to observe at least
one event in the selected time window: $p_{\rm time}=p(N\geq 1)=1-e^{-n_{B}}$.
We have $p_{\rm time}\simeq 12.6\%$ for $n_{B}\simeq 0.13$ (total number of
expected events in $1000\text{\,}\mathrm{s}$). The term $p_{\rm 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_{\rm 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
$\mathcal{L}^{(k)}_{\nu}(n_{S}^{(k)},\gamma;\Omega_{S})$ is defined as:

$\mathcal{L}^{(k)}_{\nu}(n_{S}^{(k)},\gamma;\Omega_{S})=\frac{e^{-(n_{S}^{(k)}+n_{B}^{(k)})}(n_{S}^{(k)}+n_{B}^{(k)})^{N^{(k)}}}{N^{(k)}!}\prod_{i=1}^{N^{(k)}}\frac{n_{S}^{(k)}\mathcal{S}^{(k)}(\boldsymbol{x_{i}},E_{i};\Omega_{S},\gamma)+n_{B}^{(k)}\mathcal{B}^{(k)}(\boldsymbol{x_{i}},E_{i})}{n_{S}^{(k)}+n_{B}^{(k)}},$
(1)

where $n_{S}^{(k)}$ is the number of signal events in the sample (to be
fitted), $\gamma$ is the spectral index of the signal neutrino spectrum (${\rm
d}n/{\rm d}E_{\nu}\propto E^{-\gamma}$, to be fitted as well), $\Omega_{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. $\mathcal{S}^{(k)}(\boldsymbol{x_{i}},E_{i};\Omega_{S},\gamma)$ is the
signal probability density function (pdf), which depends on reconstructed
event direction $\boldsymbol{x_{i}}$, reconstructed event energy $E_{i}$,
source spectral index, and direction.
$\mathcal{B}^{(k)}(\boldsymbol{x_{i}},E_{i})$ is the background pdf, which
depends solely on event information.

The $\mathcal{S}^{(k)}$ and $\mathcal{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:

$\displaystyle\mathcal{S}^{(k)}(\boldsymbol{x_{i}},E_{i};\Omega_{S},\gamma)$
$\displaystyle=\mathcal{A}_{S}(\boldsymbol{x_{i}};E_{i},\Omega_{S})\mathcal{E}_{S}(E_{i};\gamma)$
(2) $\displaystyle\text{and }\mathcal{B}^{(k)}(\boldsymbol{x_{i}},E_{i})$
$\displaystyle=\mathcal{A}_{B}(\boldsymbol{x_{i}})\mathcal{E}_{B}(E_{i}),$ (3)

where the point-spread function
$\mathcal{A}_{S}(\boldsymbol{x_{i}};E_{i},\Omega_{S})$ is characterizing the
angular resolution of the detector at the considered energy $E_{i}$ (it is
maximal for $\boldsymbol{x_{i}}=\Omega_{S}$, i.e. neutrino in the direction of
the probed point source), and the energy part $\mathcal{E}_{S}(E_{i};\gamma)$
is the convolution of the energy spectrum ${\rm d}n/{\rm d}E_{\nu}\propto
E^{-\gamma}$ with the energy response function $f(E_{i};E_{\nu})$. The
functions $\mathcal{A}_{B}(\boldsymbol{x_{i}})$ and $\mathcal{E}_{B}(E_{i})$
characterize the expected background distribution with direction and energy.

For each sample $k$ and direction $\Omega_{S}$, the values
$\widehat{n_{S}^{(k)}}$ and $\widehat{\gamma^{(k)}}$ maximizing the likelihood
$\mathcal{L}^{(k)}_{\nu}(n_{S}^{(k)},\gamma;\Omega_{S})$ were obtained using
iminuit (Dembinski et al., 2020)333In the implementation, as the maximization
is performed independently for each sample k, the $\widehat{\gamma^{(k)}}$ may
differ, even though from the physical point of view, there should be only one
common value. In practice, this has almost no impact as, in most of the cases,
only zero or one event are observed in SK in the time window, so that Equation
4 is much simpler and only one sample contributes.. The log-likelihood ratio
$\Lambda(\Omega_{S})$ was then computed:

$\displaystyle\Lambda(\Omega_{S})$
$\displaystyle=2\sum_{k}\ln\left[\dfrac{\mathcal{L}_{\nu}(\widehat{n_{S}^{(k)}},\widehat{\gamma^{(k)}};\Omega_{S})}{\mathcal{L}_{\nu}(n_{S}^{(k)}=0;\Omega_{S})}\right]+2\ln\mathcal{P}_{GW}(\Omega_{S}),$
(4)

where $\sum_{k}$ sums over the three considered samples and
$\mathcal{P}_{GW}(\Omega_{S})$ 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 $\Lambda(\Omega)$ 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):

$TS=\underset{\Omega}{\max}\left[\Lambda(\Omega)\right].$ (5)

Finally, this number can be used to compute $p_{\rm space}$. First, the
observation in SK was used to compute $TS_{\rm 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
$\mathcal{P}_{\rm bkg}(TS)$ was obtained and compared to the data value:

$p_{\rm space}=\int_{TS_{\rm data}}^{\infty}\mathcal{P}_{\rm bkg}(TS)\,{\rm
d}TS.$ (6)

The value $p_{\rm 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_{\rm space}$ for the GW triggers with at
least one SK event (for the other triggers, we trivially have $p=p_{\rm
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^{\rm best}_{\rm
space}=1.72\%$ ($p^{\rm best}=0.22\%$), corresponding to $2.1\sigma$
(respectively $2.9\sigma$). However, one needs to take into account the total
number of trials involved in the catalog search ($N=10$ for $p_{\rm 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\ldots N}<p^{\rm best}$. This gives post-trial
values $P^{\rm best}_{\rm space}=15.9\%$ ($1.0\sigma$) and $P^{\rm
best}=7.8\%$ ($1.4\sigma$), which are fully consistent with the background-
only hypothesis.

###  4.2 Flux limits using high-energy samples

Because no statistically significant event excess was observed within the
$1000\text{\,}\mathrm{s}$ time window in the HE-$\nu$ 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 $\gamma=2$, that is commonly used in such
astrophysical searches (e.g. Abe et al. (2018)). The neutrino flux can then be
written as:

$\dfrac{{\rm d}n}{{\rm d}E_{\nu}}=\phi_{0}E_{\nu}^{-2}.$ (7)

In the following, we will report the upper limits on
$\phi_{0}=E_{\nu}^{2}\,{\rm d}n/{\rm d}E_{\nu}$ [in
$\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$], for the different samples and
neutrino flavors ($\nu_{\mu}$, $\bar{\nu}_{\mu}$, $\nu_{e}$, $\bar{\nu}_{e}$).

####  4.2.1 Sample-by-sample approach

For a given sample $s$, flavor $f$ and source position $\Omega$, the neutrino
flux $E^{2}{\rm d}n/{\rm d}E$ is related to the number of events:

$\displaystyle N_{\rm sig}$ $\displaystyle=\phi_{0}^{(s,f)}\times\int
A^{(s,f)}_{\rm eff}(E_{\nu},\Omega)\times E_{\nu}^{-2}\,{\rm d}E_{\nu}$
$\displaystyle=\phi_{0}^{(s,f)}\times c^{(s,f)}(\Omega),$ (8)

where $A^{(s,f)}_{\rm eff}(E_{\nu},\Omega)$ is the SK detector effective area
for the selected sample and neutrino flavor, the integration range is
$0.1{10}^{5}\text{\,}\mathrm{GeV}$. The quantity $c^{(s,f)}(\Omega)$ is the
detector acceptance, which takes into account the source direction. To
marginalize over the source localization, the following Poisson likelihood is
then defined:

$\mathcal{L}^{(s,f)}(\phi_{0};n_{B}^{(k)},N^{(k)})=\displaystyle\int\dfrac{(c^{(s,f)}(\Omega)\times\phi_{0}+n_{B}^{(k)})^{N^{(k)}}}{N^{(k)}!}e^{-(c^{(s,f)}(\Omega)\times\phi_{0}+n_{B}^{(k)})}\times\mathcal{P}_{\rm
GW}(\Omega)\,{\rm d}\Omega,$ (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 $\phi_{0}$
and finding the 90% percentile, for each sample and flavor (this effectively
corresponds to the Bayesian limit with flat prior on $\phi_{0}$):

$\dfrac{\int_{0}^{\phi_{0}^{\rm
90\%}}\mathcal{L}^{(s,f)}(\phi_{0};n_{B}^{(k)},N^{(k)})\,{\rm
d}\phi_{0}}{\int_{0}^{\infty}\mathcal{L}^{(s,f)}(\phi_{0};n_{B}^{(k)},N^{(k)})\,{\rm
d}\phi_{0}}=0.90.$ (10)

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
$\theta$ for FC and PC, while UPMU covers only efficiently below the horizon
($\theta>$90\text{\,}\mathrm{\SIUnitSymbolDegree}$$), with a nonnegligible
variation with $\theta$, as shown in Figure 2. The UPMU sample has very
limited sensitivity at and above the horizon
($0<\theta<$90\text{\,}\mathrm{\SIUnitSymbolDegree}$$), 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.

Figure 2: Effective area of the SK event selection, for the different HE-$\nu$
samples and neutrino flavors, as a function of neutrino energy: dark (light)
green for muon (anti-)neutrinos and dark (light) orange for electron
(anti-)neutrinos. The different line styles show the variation of effective
area for different ranges in zenith angle.

The full results for $\nu_{\mu}$ 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.

Figure 3: Obtained 90% C.L. upper limits on $E^{2}{\rm d}n/{\rm d}E$ for $\nu_{\mu}$ and for the different GW triggers, using the methods presented in section 4.2.1 (sample-by-sample) and section 4.2.2 for the combined analysis. The complete figure set (4 images, one per considered neutrino flavor) is available in the online journal. Table 3: Obtained 90% C.L. upper limits on $E^{2}{\rm d}n/{\rm d}E$ for GW190425 and GW190521. For HE-$\nu$, limits on $E^{2}{\rm d}n/{\rm d}E$ [in $\mathrm{GeV}\text{\,}{\mathrm{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_{\rm iso}$ [in erg] (assuming isotropic emission) are also presented: one limit per flavor and limits for $\nu_{e}$ \+ $\bar{\nu}_{e}$, $\nu_{\mu}$ \+ $\bar{\nu}_{\mu}$ and on the total energy in all flavors assuming equipartition (including unseen tau neutrinos). For LE-$\nu$, limits on the total neutrino fluence $\Phi$ [in ${\mathrm{cm}}^{-2}$] are given for $\nu_{e}$, $\bar{\nu}_{e}$, $\nu_{x}=\nu_{\mu}+\nu_{\tau}$, $\bar{\nu}_{x}=\bar{\nu}_{\mu}+\bar{\nu}_{\tau}$ assuming Fermi-Dirac spectrum (with average energy of $20\text{\,}\mathrm{MeV}$) and flat spectrum (within the range $7100\text{\,}\mathrm{MeV}$), as well as upper limits on $E_{\rm iso}$ [in erg] for the Fermi-Dirac scenario.  Trigger name | Sample | $\nu_{e}$ | $\bar{\nu}_{e}$ |  $\nu_{\mu}$ ($\nu_{x}$) |  $\bar{\nu}_{\mu}$ ($\bar{\nu}_{x}$)
---|---|---|---|---|---
| HE $E^{2}\dfrac{{\rm d}n}{{\rm d}E}$ | FC | $2.22\cdot 10^{3}$ | $4.32\cdot 10^{3}$ | $3.91\cdot 10^{3}$ | $9.42\cdot 10^{3}$
|  | PC | $3.32\cdot 10^{4}$ | $1.12\cdot 10^{5}$ | $4.81\cdot 10^{3}$ | $8.74\cdot 10^{3}$
|  | UPMU | $-$ | $-$ | $-$ | $-$
|  | Combined | $2.09\cdot 10^{3}$ | $4.28\cdot 10^{3}$ | $2.16\cdot 10^{3}$ | $4.20\cdot 10^{3}$
| HE $E_{\rm iso}$ | Per-flavour | $1.98\cdot 10^{56}$ | $3.85\cdot 10^{56}$ | $1.96\cdot 10^{56}$ | $3.69\cdot 10^{56}$
GW190425 |  | $\nu+\bar{\nu}$ | $2.62\cdot 10^{56}$ | $2.52\cdot 10^{56}$
|  | All | $3.47\cdot 10^{56}$
| LE $\Phi$ | Flat | $1.49\cdot 10^{9}$ | $1.83\cdot 10^{7}$ | $9.35\cdot 10^{9}$ | $1.11\cdot 10^{10}$
|  | Fermi-Dirac | $3.92\cdot 10^{9}$ | $9.57\cdot 10^{7}$ | $2.43\cdot 10^{10}$ | $2.87\cdot 10^{10}$
| LE $E_{\rm iso}$ | Per-flavour | $3.92\cdot 10^{59}$ | $9.59\cdot 10^{57}$ | $2.43\cdot 10^{60}$ | $2.87\cdot 10^{60}$
|  | All | $5.54\cdot 10^{58}$
| HE $E^{2}\dfrac{{\rm d}n}{{\rm d}E}$ | FC | $2.27\cdot 10^{3}$ | $4.71\cdot 10^{3}$ | $3.76\cdot 10^{3}$ | $9.60\cdot 10^{3}$
|  | PC | $3.66\cdot 10^{4}$ | $3.68\cdot 10^{4}$ | $4.89\cdot 10^{3}$ | $8.35\cdot 10^{3}$
|  | UPMU | $-$ | $-$ | $4.48\cdot 10^{1}$ | $5.04\cdot 10^{1}$
|  | Combined | $2.21\cdot 10^{3}$ | $4.60\cdot 10^{3}$ | $3.75\cdot 10^{1}$ | $4.82\cdot 10^{1}$
| HE $E_{\rm iso}$ | Per-flavour | $1.69\cdot 10^{59}$ | $3.46\cdot 10^{59}$ | $2.58\cdot 10^{57}$ | $3.72\cdot 10^{57}$
GW190521 |  | $\nu+\bar{\nu}$ | $2.26\cdot 10^{59}$ | $3.00\cdot 10^{57}$
|  | All | $8.94\cdot 10^{57}$
| LE $\Phi$ | Flat | $2.63\cdot 10^{9}$ | $3.22\cdot 10^{7}$ | $1.65\cdot 10^{10}$ | $1.95\cdot 10^{10}$
|  | Fermi-Dirac | $6.89\cdot 10^{9}$ | $1.68\cdot 10^{8}$ | $4.27\cdot 10^{10}$ | $5.04\cdot 10^{10}$
| LE $E_{\rm iso}$ | Per-flavour | $5.85\cdot 10^{62}$ | $1.43\cdot 10^{61}$ | $3.63\cdot 10^{63}$ | $4.28\cdot 10^{63}$
|  | All | $8.26\cdot 10^{61}$

####  4.2.2 Combination of the samples (using test statistic)

As the neutrino spectrum is expected to span the full range from
$0.1\text{\,}\mathrm{GeV}$ to ${10}^{5}\text{\,}\mathrm{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 $\mathcal{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 $\mathcal{P}_{n_{S}}(TS)$ 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:

$\mathcal{L}(\phi_{0};TS_{\rm
data},\mathcal{P}_{GW})=\int\sum_{k=0}^{2}\Big{[}\dfrac{\left(c(\Omega)\phi_{0}\right)^{k}}{k!}e^{-c(\Omega)\phi_{0}}\times\mathcal{P}_{k}(TS_{\rm
data})\Big{]}\times\mathcal{P}_{GW}(\Omega)\,{\rm d}\Omega,$ (11)

where $c(\Omega)=\sum_{s}c^{(s,f)}(\Omega)$ 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 $\phi_{0}=E_{\nu}^{2}{\rm d}n/{\rm d}E_{\nu}$ is
then simply obtained as in Equation 10. The procedure can be repeated
independently for each neutrino flavor or also combining flavors, e.g.,
$\nu_{\mu}+\bar{\nu}_{\mu}$ (in the latter case, both
$\mathcal{P}_{n_{S}}(TS)$ terms and $c(\Omega)$ 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_{\rm data}$ and thus impact $\mathcal{P}_{k}(TS_{\rm data})$ used
in the Equation 11.

###  4.3 Flux limits using Low-energy sample

The flux limit calculation for the low-energy sample is similar to HE-$\nu$,
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-$\nu$
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-$\nu$ case.

The upper limit on the total fluence ($\Phi=\int{\rm d}n/{\rm d}E\;{\rm d}E$)
is then simply computed as:

$\Phi_{90}=\frac{N_{90}}{N_{T}\int\lambda(E_{\nu})\sigma(E_{\nu})R(E_{e},E_{\rm
vis})\epsilon(E_{\rm vis})\,{\rm d}E_{\nu}},$ (12)

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, $\sigma$ is the combined cross section for all
interactions, $\epsilon$ is the detection efficiency, $\lambda$ is the energy
density assuming a given spectrum (${\rm d}n/{\rm
d}E_{\nu}=\Phi\times\lambda(E)$) and $R$ is the response function to convert
electron or positron energy ($E_{e}$) to visible energy in SK ($E_{\rm vis}$).
The response function and the detection efficiency ($\epsilon$) are calculated
using SK detector Monte Carlo simulations, and related systematic
uncertainties are neglected as in the HE-$\nu$ case.

In this analysis, two types of spectra were considered: flat spectrum
($\lambda=\textrm{constant}$) and Fermi-Dirac spectrum with an average energy
of $20\text{\,}\mathrm{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 $\bar{\nu}_{e}$ case, given that the main interaction
channel in the detector is inverse beta decay of $\bar{\nu}_{e}$, as described
in section 2.

##  5 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_{\rm iso}$ emitted by neutrinos from a source at a distance $d$ is
directly related to the detected flux at Earth:

$E_{\rm iso}=4\pi d^{2}\int\frac{{\rm d}n}{{\rm d}E}\times E\,{\rm d}E.$ (13)

Knowing $d$, one can then constrain $E_{\rm iso}$, as described in the
following.

###  5.1 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:

$\dfrac{E_{\rm iso}}{4\pi d^{2}}=\int\phi_{0}E^{-2}\times E\,{\rm
d}E=\phi_{0}\times\ln\left(\dfrac{E_{\rm max}}{E_{\rm min}}\right).$ (14)

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

$\mathcal{L}(E_{\rm iso};TS_{\rm
data}^{(i)},\mathcal{V}_{GW}^{(i)})=\int\sum_{k=0}^{2}\Big{[}\dfrac{\left(c^{\prime}(r,\Omega)E_{\rm
iso}\right)^{k}}{k!}e^{-c^{\prime}(r,\Omega)E_{\rm
iso}}\times\mathcal{P}^{(i)}_{k}(TS^{(i)}_{\rm
data})\Big{]}\times\mathcal{V}^{(i)}_{GW}(r,\Omega){\rm d}r{\rm d}\Omega.$
(15)

The quantity $c^{\prime}(r,\Omega)$ is the conversion factor from $E_{\rm
iso}$ to the expected number of signal events in SK for known source distance
$r$ and direction $\Omega$. The test statistic distributions
$\mathcal{P}^{(i)}_{k}(TS)$ and the measured test statistic $TS^{(i)}_{\rm
data}$ for trigger $i$ are the same as defined in section 4.2.2.
$\mathcal{V}_{GW}^{(i)}(r,\Omega)$ 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_{\rm iso}$ limits independently for each GW trigger as it
has been done for the flux limits. For a given flavor (e.g., $\nu_{\mu}$), 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 $(\nu_{e}:\nu_{\mu}:\nu_{\tau})$ equal to $(1:2:0)$,
which would become $\sim(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 $\nu_{\mu}$ and $\bar{\nu}_{\mu}$ ; the
limit on the total energy emitted in neutrinos assuming equipartition is then
dominated by the latter contributions: $E_{\rm iso,90\%}^{\rm all-
flavors}\simeq 3\times E_{\rm iso,90\%}^{\nu_{\mu}+\bar{\nu}_{\mu}}$. 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 $\sim
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_{\rm iso}$
limits.

Figure 4: 90% C.L. upper limits on the isotropic energy emitted in neutrinos
for the 36 GW triggers followed up by SK, as a function of source distance.
The distance and its error, as well as the source type (indicated by the
different colors and markers), are provided using the data from Abbott et al.
(2020b) ($m<$3\text{\,}\mathrm{M_{\odot}}$$ = NS,
$m>$3\text{\,}\mathrm{M_{\odot}}$$ = BH). The limits are following two lines
$E_{\rm iso}^{90\%}\propto\text{distance}^{2}$ based on geometrical
considerations, one of the lines shows events dominated by UPMU
$\nu_{\nu}/\bar{\nu}_{\mu}$ 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.

The combination of a meaningful set of GW events to constrain further $E_{\rm
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:

$\mathcal{L}^{\rm Pop}(E_{\rm iso};\\{TS_{\rm
data})^{(i)}\\},\\{\mathcal{V}_{GW}^{(i)}\\})=\prod_{i=1}^{N}\mathcal{L}(E_{\rm
iso};TS_{\rm data})^{(i)},\mathcal{V}_{GW}^{(i)}),$ (16)

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 $\mathcal{M}_{\rm tot}$ of the binary system: $E_{\rm
iso}^{\nu}=f_{\nu}\times\mathcal{M}_{\rm tot}$. One can then use the following
likelihood to constrain $f_{\nu}$:

$\mathcal{L}^{\rm Pop}(f_{\nu};\\{TS_{\rm
data}^{(i)}\\},\\{\mathcal{V}_{GW}^{(i)}\\},\\{\mathcal{M}_{\rm
tot}^{(i)}\\})=\prod_{i=1}^{N}\int\mathcal{M}_{\rm
tot}^{(i)}\times\mathcal{L}(f_{\nu}\mathcal{M}_{\rm tot}^{(i)};TS_{\rm
data}^{(i)},\mathcal{V}_{GW}^{(i)})\times p_{\rm GW}(\mathcal{M}_{\rm
tot}^{(i)})\times{\rm d}\mathcal{M}_{\rm tot}^{(i)},$ (17)

where $f_{\nu}$, in $\mathrm{erg}\text{/}\mathrm{M_{\odot}}$ is to be
constrained (simplifying the units, $f_{\nu}$ can be expressed as the
proportion of the total mass converted in neutrinos: e.g.
$f_{\nu}=${10}^{54}\text{\,}\mathrm{erg}\text{/}\mathrm{M_{\odot}}$=62\%$),
and $p_{\rm GW}(\mathcal{M}_{\rm tot}^{(i)})$ 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
candidate444In this case, the result is the same as using directly GW190425
event, as it is the only identified BNS in O3a, 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_{\rm
iso}<$4.16\text{\times}{10}^{55}\text{\,}\mathrm{erg}$$ assuming all objects
have similar emission. This turns to $E_{\rm
iso}<$9.73\text{\times}{10}^{56}\text{\,}\mathrm{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.

Figure 5: 90% C.L. upper limits on the isotropic energy emitted in neutrinos
by combining GW triggers with the same nature, for $\nu_{\mu}$,
$\nu_{\mu}+\bar{\nu}_{\mu}$, $\nu_{e}+\bar{\nu}_{e}$ and all-flavor emission
(assuming equipartition). The left panel shows the results assuming that all
selected sources are emitting the same $E_{\rm iso}$ while the right panel is
assuming that neutrino emission is scaling with the total mass of the binary
system.

###  5.2 Low-energy neutrino emission

As for the flux limits, the low-energy case is much simpler. $E_{\rm 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 $\bar{\nu}_{e}$ limit. To cover the distance uncertainty, the
following likelihood was defined:

$\mathcal{L}(E_{\rm iso};N_{\rm obs},N_{\rm
bkg})=\displaystyle\int_{0}^{\infty}\dfrac{\left(N_{\rm bkg}+c^{\rm
LE}(r)\times E_{\rm iso}\right)^{N_{\rm obs}}}{N_{\rm obs}!}e^{-\left(N_{\rm
bkg}+c^{\rm LE}(r)\times E_{\rm iso}\right)}\times p_{\rm GW}(r){\rm d}r,$
(18)

where $N_{\rm obs}$ and $N_{\rm bkg}$ are the observed and expected number of
LE-$\nu$ events, $c^{\rm LE}(r)$ is the conversion factor from $E_{\rm iso}$
to number of signal events assuming Fermi-Dirac spectrum and source at
distance $r$, $p_{\rm 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.

##  6 Discussion and conclusions

The results of the follow-up of LVC O3a gravitational waves with the SK
detector have been presented. In the $\pm$500\text{\,}\mathrm{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-$\nu$, one for LE-$\nu$). Upper limits on the incoming neutrino flux were
computed for the different neutrino flavors. For HE-$\nu$, $E^{-2}$ spectrum
was assumed, while for LE-$\nu$ limits, Fermi-Dirac emission with average
energy of $20\text{\,}\mathrm{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_{\rm 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_{\rm iso}^{\rm
model}\sim$47\text{\times}{10}^{53}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$),
even though the exact shape of the neutrino spectrum (beyond the assumed
simple Fermi-Dirac distribution with $\langle
E_{\nu}\rangle=$20\text{\,}\mathrm{MeV}$$) may modify slightly the obtained
upper limits.

For high-energy neutrino emissions, the obtained limits on $E_{\rm iso}$
assuming $E^{-2}$ spectrum are barely covering the nonphysical region where
the total mass of the binary system is converted to neutrinos
($f_{\nu}\lesssim${10}^{54}{10}^{56}\text{\,}\mathrm{erg}\text{/}\mathrm{M_{\odot}}$\simeq
60-6000\%$), while the region currently probed by IceCube is $f_{\nu}\lesssim
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 $\mathrm{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+$\nu$ source population and better understand the underlying mechanisms.

We gratefully acknowledge the cooperation of the Kamioka Mining and Smelting
Company. The Super‐Kamiokande experiment has been built and operated from
funding by the Japanese Ministry of Education, Culture, Sports, Science and
Technology, the U.S. Department of Energy, and the U.S. National Science
Foundation. Some of us have been supported by funds from the National Research
Foundation of Korea NRF‐2009‐0083526 (KNRC) funded by the Ministry of Science,
ICT, and Future Planning and the Ministry of Education (2018R1D1A3B07050696,
2018R1D1A1B07049158), the Japan Society for the Promotion of Science, the
National Natural Science Foundation of China under Grants No.11620101004, the
Spanish Ministry of Science, Universities and Innovation (grant
PGC2018-099388-B-I00), the Natural Sciences and Engineering Research Council
(NSERC) of Canada, the Scinet and Westgrid consortia of Compute Canada, the
National Science Centre, Poland (2015/18/E/ST2/00758), the Science and
Technology Facilities Council (STFC) and GridPPP, UK, the European Union’s
Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-
Curie grant agreement no.754496, H2020-MSCA-RISE-2018 JENNIFER2 grant
agreement no.822070, and H2020-MSCA-RISE-2019 SK2HK grant agreement no.
872549.

##  Appendix A Additional material and data release

This section details all the results not presented in the main text of the
paper. The Table 4 contains the number of observed and expected events in the
different samples and for each follow-up, as well as computed flux and $E_{\rm
iso}$ limits.

Additionally, the data release Abe et al. (2021a) contains the effective areas
that have actually been involved in the computation of flux upper limits, as
presented in section 4.2.1. These can be used to derive again the upper limits
with a specific source position or a different spectrum.

Table 4: Content of the detailed data release table.

Col.# | Label | Unit | Description
---|---|---|---
1 | GW_NAME | - | Name of the GW trigger
2 | GW_UTC | - | UTC time of the trigger
3 | GW_SKYAREA90 | deg2 | Surface of the 90% containment of GW localization
4 | GW_DISTANCE | $\mathrm{Mpc}$ | Mean estimate of the distance to GW source
5 | SK_LIVETIME | $\mathrm{s}$ | Live time of Super-Kamiokande over the selected $1000\text{\,}\mathrm{s}$ time window
6 | SK_FC_OBSERVED | - | Number of observed events in the HE-$\nu$/FC sample in the time window
7 | SK_FC_EXPECTED | - | Number of expected background events in the HE-$\nu$/FC sample in the time window
8 | SK_PC_OBSERVED | - | Same for HE-$\nu$/PC
9 | SK_PC_EXPECTED | - | Same for HE-$\nu$/PC
10 | SK_UPMU_OBSERVED | - | Same for HE-$\nu$/UPMU
11 | SK_UPMU_EXPECTED | - | Same for HE-$\nu$/UPMU
12 | SK_LOWE_OBSERVED | - | Same for LE-$\nu$
13 | SK_LOWE_EXPECTED | - | Same for LE-$\nu$
14 | E2PHI90_NUE_FC | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | 90% U.L. on $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{e}}$ using the HE-$\nu$/FC sample, assuming $E^{-2}$ spectrum
15 | E2PHI90_NUEB_FC | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\bar{\nu}_{e}}$
16 | E2PHI90_NUENUEB_FC | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{e}+\bar{\nu}_{e}}$
17 | E2PHI90_NUMU_FC | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{\mu}}$
18 | E2PHI90_NUMUB_FC | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\bar{\nu}_{\mu}}$
19 | E2PHI90_NUMUNUMUB_FC | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{\mu}+\bar{\nu}_{\mu}}$
20 | E2PHI90_NUE_PC | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | 90% U.L. on $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{e}}$ using the HE-$\nu$/PC sample, assuming $E^{-2}$ spectrum
21 | E2PHI90_NUEB_PC | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\bar{\nu}_{e}}$
22 | E2PHI90_NUENUEB_PC | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{e}+\bar{\nu}_{e}}$
23 | E2PHI90_NUMU_PC | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{\mu}}$
24 | E2PHI90_NUMUB_PC | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\bar{\nu}_{\mu}}$
25 | E2PHI90_NUMUNUMUB_PC | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{\mu}+\bar{\nu}_{\mu}}$
26 | E2PHI90_NUMU_UPMU | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | 90% U.L. on $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{\mu}}$ using the HE-$\nu$/UPMU sample, assuming $E^{-2}$ spectrum
27 | E2PHI90_NUMUB_UPMU | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\bar{\nu}_{\mu}}$
28 | E2PHI90_NUMUNUMUB_UPMU | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{\mu}+\bar{\nu}_{\mu}}$
29 | E2PHI90_NUE_COMBINED | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | 90% U.L. on $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{e}}$ using all HE-$\nu$ samples, assuming $E^{-2}$ spectrum
30 | E2PHI90_NUEB_COMBINED | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\bar{\nu}_{e}}$
31 | E2PHI90_NUENUEB_COMBINED | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{e}+\bar{\nu}_{e}}$
32 | E2PHI90_NUMU_COMBINED | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{\mu}}$
33 | E2PHI90_NUMUB_COMBINED | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\bar{\nu}_{\mu}}$
34 | E2PHI90_NUMUNUMUB_COMBINED | $\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}$ | Same for $E^{2}\left.{\rm d}n/{\rm d}E\right|_{\nu_{\mu}+\bar{\nu}_{\mu}}$
35 | EISO90_NUE_COMBINED | $\mathrm{erg}$ | 90% U.L. on $\left.E_{\rm iso}\right|_{\nu_{e}}$ using all HE-$\nu$ samples, assuming $E^{-2}$ spectrum
36 | EISO90_NUEB_COMBINED | $\mathrm{erg}$ | Same for $\left.E_{\rm iso}\right|_{\bar{\nu}_{e}}$
37 | EISO90_NUENUEB_COMBINED | $\mathrm{erg}$ | Same for $\left.E_{\rm iso}\right|_{\nu_{e}+\bar{\nu}_{e}}$
38 | EISO90_NUMU_COMBINED | $\mathrm{erg}$ | Same for $\left.E_{\rm iso}\right|_{\nu_{\mu}}$
39 | EISO90_NUMUB_COMBINED | $\mathrm{erg}$ | Same for $\left.E_{\rm iso}\right|_{\bar{\nu}_{\mu}}$
40 | EISO90_NUMUNUMUB_COMBINED | $\mathrm{erg}$ | Same for $\left.E_{\rm iso}\right|_{\nu_{\mu}+\bar{\nu}_{\mu}}$
41 | EISO90_ALL_COMBINED | $\mathrm{erg}$ | Same for all-flavors $E_{\rm iso}$ (assuming equipartition between flavors)
42 | EISO90_ALL_COMBINED_GAMMA3 | $\mathrm{erg}$ | Same but assuming $E^{-3}$ spectrum
43 | FLUENCE90_LOWE_NUE_FERMIDIRAC | ${\mathrm{cm}}^{-2}$ | 90% U.L. on $\Phi_{\nu_{e}}$ using LE-$\nu$ sample, assuming Fermi-Dirac spectrum
44 | FLUENCE90_LOWE_NUEB_FERMIDIRAC | ${\mathrm{cm}}^{-2}$ | Same for $\Phi_{\bar{\nu}_{e}}$
45 | FLUENCE90_LOWE_NUX_FERMIDIRAC | ${\mathrm{cm}}^{-2}$ | Same for $\Phi_{\nu_{\mu}+\nu_{\tau}}$
46 | FLUENCE90_LOWE_NUXB_FERMIDIRAC | ${\mathrm{cm}}^{-2}$ | Same for $\Phi_{\bar{\nu}_{\mu}+\bar{\nu}_{\tau}}$
47 | FLUENCE90_LOWE_NUE_FLAT | ${\mathrm{cm}}^{-2}$ | 90% U.L. on $\Phi_{\nu_{e}}$ using LE-$\nu$ sample, assuming flat spectrum
48 | FLUENCE90_LOWE_NUEB_FLAT | ${\mathrm{cm}}^{-2}$ | Same for $\Phi_{\bar{\nu}_{e}}$
49 | FLUENCE90_LOWE_NUX_FLAT | ${\mathrm{cm}}^{-2}$ | Same for $\Phi_{\nu_{\mu}+\nu_{\tau}}$
50 | FLUENCE90_LOWE_NUXB_FLAT | ${\mathrm{cm}}^{-2}$ | Same for $\Phi_{\bar{\nu}_{\mu}+\bar{\nu}_{\tau}}$
51 | EISO90_LOWE_NUE_FERMIDIRAC | $\mathrm{erg}$ | 90% U.L. on $\left.E_{\rm iso}\right|_{\nu_{e}}$ using LE-$\nu$ sample, assuming Fermi-Dirac spectrum
52 | EISO90_LOWE_NUEB_FERMIDIRAC | $\mathrm{erg}$ | Same for $\left.E_{\rm iso}\right|_{\bar{\nu}_{e}}$
53 | EISO90_LOWE_NUX_FERMIDIRAC | $\mathrm{erg}$ | Same for $\left.E_{\rm iso}\right|_{\nu_{\mu}+\nu_{\tau}}$
54 | EISO90_LOWE_NUXB_FERMIDIRAC | $\mathrm{erg}$ | Same for $\left.E_{\rm iso}\right|_{\bar{\nu}_{\mu}+\bar{\nu}_{\tau}}$
55 | EISO90_LOWE_ALL_FERMIDIRAC | $\mathrm{erg}$ | Same for all-flavors $E_{\rm iso}$ (assuming equipartition between flavors)

Note. — Only the description of the table is shown here. The complete table is
available as MRT format in the online journal and as CSV format in Abe et al.
(2021a).

## References

  * Aartsen et al. (2020) Aartsen, M., et al. 2020, Astrophys. J. Lett., 898, L10, doi: 10.3847/2041-8213/ab9d24
  * Abbasi et al. (2021) Abbasi, R., et al. 2021, Probing neutrino emission at GeV energies from compact binary mergers with the IceCube Neutrino Observatory.  https://arxiv.org/abs/2105.13160
  * Abbott et al. (2020a) Abbott, B., et al. 2020a, Astrophys. J. Lett., 892, L3, doi: 10.3847/2041-8213/ab75f5
  * Abbott et al. (2020b) Abbott, R., et al. 2020b, GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run.  https://arxiv.org/abs/2010.14527
  * Abbott et al. (2020c) —. 2020c, Phys. Rev. Lett., 125, 101102, doi: 10.1103/PhysRevLett.125.101102
  * Abe et al. (2016a) Abe, K., et al. 2016a, Astrophys. J. Lett., 830, L11, doi: 10.3847/2041-8205/830/1/L11
  * Abe et al. (2016b) —. 2016b, Phys. Rev. D, 94, 052010, doi: 10.1103/PhysRevD.94.052010
  * Abe et al. (2018) —. 2018, Astrophys. J. Lett., 857, L4, doi: 10.3847/2041-8213/aabaca
  * Abe et al. (2021a) —. 2021a, Data release, Zenodo, doi: 10.5281/zenodo.4724822
  * Abe et al. (2021b) Abe, S., et al. 2021b, Astrophys. J., 909, 116, doi: 10.3847/1538-4357/abd5bc
  * Agostini et al. (2017) Agostini, M., et al. 2017, Astrophys. J., 850, 21, doi: 10.3847/1538-4357/aa9521
  * Ashie et al. (2005) Ashie, Y., et al. 2005, Phys. Rev. D, 71, 112005, doi: 10.1103/PhysRevD.71.112005
  * Baret et al. (2011) Baret, B., et al. 2011, Astropart. Phys., 35, 1, doi: 10.1016/j.astropartphys.2011.04.001
  * Barthelmy (2021) Barthelmy, S. 2021, GCN: The Gamma-ray Coordinates Network, https://gcn.gsfc.nasa.gov/
  * Bays et al. (2012) Bays, K., et al. 2012, Phys. Rev. D, 85, 052007, doi: 10.1103/PhysRevD.85.052007
  * Countryman et al. (2019) Countryman, S., Keivani, A., Bartos, I., et al. 2019, Low-Latency Algorithm for Multi-messenger Astrophysics (LLAMA) with Gravitational-Wave and High-Energy Neutrino Candidates.  https://arxiv.org/abs/1901.05486
  * Dembinski et al. (2020) Dembinski, H., et al. 2020, scikit-hep/iminuit: v1.5.2, v1.5.2, Zenodo, doi: 10.5281/zenodo.4047970
  * Dornic et al. (2020) Dornic, D., Coleiro, A., Colomer Molla, M., Kouchner, A., & Pradier, T. 2020, PoS, ICRC2019, 872
  * Foucart et al. (2016) Foucart, F., Haas, R., Duez, M. D., et al. 2016, Phys. Rev. D, 93, 044019, doi: 10.1103/PhysRevD.93.044019
  * Fukuda et al. (2003) Fukuda, Y., et al. 2003, Nucl. Instrum. Meth. A, 501, 418, doi: 10.1016/S0168-9002(03)00425-X
  * Gorski et al. (2005) Gorski, K., Hivon, E., Banday, A., et al. 2005, Astrophys. J., 622, 759, doi: 10.1086/427976
  * Hagiwara (2020) Hagiwara, K. 2020, PhD thesis, Graduate School of Natural Science and Technology, Okayama University
  * Hussain et al. (2020) Hussain, R., Vandenbroucke, J., & Wood, J. 2020, PoS, ICRC2019, 918, doi: 10.22323/1.358.0918
  * Kimura et al. (2018) Kimura, S. S., Murase, K., Bartos, I., et al. 2018, Phys. Rev. D, 98, 043020, doi: 10.1103/PhysRevD.98.043020
  * Kimura et al. (2017) Kimura, S. S., Murase, K., Mészáros, P., & Kiuchi, K. 2017, Astrophys. J. Lett., 848, L4, doi: 10.3847/2041-8213/aa8d14
  * Kotera & Silk (2016) Kotera, K., & Silk, J. 2016, Astrophys. J. Lett., 823, L29, doi: 10.3847/2041-8205/823/2/L29
  * Nakano et al. (2020) Nakano, Y., Hokama, T., Matsubara, M., et al. 2020, Nucl. Instrum. Meth. A, 977, 164297, doi: 10.1016/j.nima.2020.164297
  * NASA/GSFC (2021) NASA/GSFC. 2021, Flexible Image Transport System (FITS), https://fits.gsfc.nasa.gov/
  * Singer et al. (2016) Singer, L. P., et al. 2016, Astrophys. J. Lett., 829, L15, doi: 10.3847/2041-8205/829/1/L15
  * Veske et al. (2020) Veske, D., Márka, Z., Bartos, I., & Márka, S. 2020, JCAP, 05, 016, doi: 10.1088/1475-7516/2020/05/016
  * Zhang et al. (2016) Zhang, Y., et al. 2016, Phys. Rev. D, 93, 012004, doi: 10.1103/PhysRevD.93.012004