Search for Neutrinos Associated with
Gravitational Waves inSuper-Kamiokande

2021, March
Wenjie Ma

Graduate School of Natural Science and Technology(Master’s Course)

Okayama University

Abstract
This thesis focuses on neutrino searching associated with gravitational waves in Super-Kamiokande(SK)which is the largest Water-Cherenkov detector with the sensitivity of MeV to TeV. We get gravitational information from gravitational detector LIGO/VIRGO which was restarted in April 2019 called O3. O3 has improved
sensitivity, and 56 gravitational waves quick alerts have already been observed and ended in April 2020.And at
the catalogs of LIGO/VIRGO,they update data from Apr.2019 to Sep.2019 called O3a, reported 39 gravitational
waves events. Gravitational waves can be generated by various processes such as the combination merger of black
holes and black holes, the merger of neutron stars. Among them, the merger of neutron stars is considered to
produce neutrino signals. If this kind of neutrino signal can be observed by us, it will have extremely profoundsignificance for the study of the specific physical process of neutron star assembly and synthesis products. Immediate observation of neutrinos is considered to be important for transmitting the neutrino search results tothe optical telescope in conjunction with the gravitational wave signal. For this, an automatic analysis system isneeded. The purpose of this research is to create this automatic analysis system and search for neutrinos derived
from gravitational waves.
In order to cover a very wide range of neutrino spectra, both low-energy (7–100 MeV) and high-energy(0.1–105GeV) samples were analyzed. O3a reported the tracking of 36 gravitational waves (39 in total) andchecked 2 catalogs.From the result in this study,in high energy region (from 100 M eV to 10 T eV ), 10 neutrinos
were observed, and the expected event was 4.8 ( 2σ) in high-energy .In low-energy(from 3.5 M eV to 100 M eV
) for each events, the max observed signals are 3 events/1000s as the expected event was 0.729 /1000s ( 2σ).
In the first chapter, this thesis will discuss Einstein’s theory of relativity and the generation of gravitational
waves. In addition, the binary star system that can generate gravitational waves is explained. The neutron starmerger is considered to be capable of generating neutrinos. According to different theoretical predictions, theenergy of the generated neutrinos is also different. The main goal of this thesis is to detect neutrinos so that wecan better understand the process of binary star system. The second chapter of this paper mainly introduces thesuper-Kamiokade(SK) detector, the main detector for detecting neutrinos. Introduce the history of the detector,
the principle of detection (via Cerenkov light), and the operating system in SK. Chapter 3 mainly introduces the
LIGO/VIRGO detectors of gravitational waves. The detection of gravitational waves is mainly carried out by
the Michelson interferometer, and LIGO/VIRGO has the most sensitive Michelson interferometer (4km arm) in
the world. In this chapter, I will introduce the history of LIGO/VIRGO and the detection principle of Michelsoninterferometer, as well as how to obtain the mass, distance and other data of the satellite system through thedetected frequency. Chapter 4 mainly introduces the work of this paper: how to build an automated system that
can automatically obtain gravitational wave data from LIGO/VIRGO and analyze it from the two region of highenergy and low energy. Among them, in the low energy range (3.5MeV to 100MeV), because the backgroundis not stable, after analyzing the energy higher than 3.5MeV, 5MeV, 6MeV, and 7MeV, it is decided to cut thedata at 7MeV as the result of this detection. Chapter 5 mainly introduces the result processing and analysis in
the two ranges of high energy and low energy respectively, and analyzes the sensitivity of the detected events.
First, the data in the two energy ranges are analyzed by simple Poisson distribution. The results show that
there is no neutrino signal. Then, in the high-energy range, considering the direction, energy, and location ofgravitational waves, the Lambda and chi-square methods were used to analyze them, and the results showedthat the neutrino signal was still not found. Although the neutrino signal could not be found, the upper limit ofneutrino was still calculated, and the isotropic total neutrino energy was calculated at the same time. Chapter
6 mainly introduces the output results of the automation system. Chapter 7 and Chapter 8 respectively show
the results of O3 and O3a of this research. The results show that there is no significant signal found.

1

Contents
1 Introduction
1.1 Gravitational waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.1 Newtonian gravity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.2 Special relativity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.3 Relativistic gravity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Binary star system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Neutrino . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4
4
4
4
5
5
6

2 Super-Kamiokande(SK)
2.1 Detector Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 SK history . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Detector Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Photomultiplier Tubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 Water Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6 Radon Free Air Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.7 Front-end Electronics and Data Acquisition System . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.8 Monitoring system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
7
7
8
8
10
10
12
12

3 LIGO/VIRGO Detector
14
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2 Detector Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3 History . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.4 Signal Display . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4 Automatic system and methods
21
4.1 Alert receiver and alert database . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.1.1 GCN Notice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.1.2 GCN Notice receiver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.2 High Energy analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2.1 Event selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.2.2 Background estimation and characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.2.3 Signal characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.2.4 Outputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.3 Low Energy analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.3.1 Event selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.3.2 Background estimation and characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.3.3 Outputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
5 Statistical analysis
36
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
5.1.1 Super-Kamiokande approach in previous analyses . . . . . . . . . . . . . . . . . . . . . . . . . 36
5.2 Signal significance:first approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.3 Flux limits:sample-by-sample, position dependent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.4 Flux limits:sample-by-sample, marginalised . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.5 Combined statistical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.5.1 Lambda method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.5.2 Chi-square method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.6 Flux limits: combined . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.6.1 Signal simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.6.2 Marginalised limits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.7 Isotropic total neutrino energy limits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.7.1 High Energy sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.7.2 Low Energy sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

2

6 Result merger and result database53
6.1 Inventory of different results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.1.1 Alert receiver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.1.2 High Energy analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.1.3 Low Energy analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
6.2 Output processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
6.3 Output formats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
6.3.1 Results database . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
6.3.2 Web page . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
6.3.3 PDF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
7 Follow-up of O3 real-time events58
7.1 All GW triggers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
8 Follow-up of O3a catalogued events61
8.1 All GW triggers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
8.2 Flux limits and Fluence limits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
8.3 Total energy limits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
9 Summary and Outlook
68
9.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
9.2 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

3

1

Introduction

1.1

Gravitational waves

Gravitational waves were proposed by Henry Poincaré in 1905. They are disturbances in the curvature of space andtime produced by accelerating masses and propagate from their source at the speed of light. Then in 1916, Albert
Einstein made a prediction based on his theory of relativity. Gravitational wave is a wave similar to electromagnetic
radiation, which transmits energy as gravitational radiation. In classical mechanics, according to Newton’s law of
universal gravitation, their existence is not stipulated, because Newtonian mechanics is based on the assumptionthat physical interaction is instantaneous (infinite speed) propagation.
1.1.1

Newtonian gravity

First, we conduct a qualitative analysis of gravitational waves. Gravitational radiation is a natural result ofgravitational interaction. Then, compare the characteristics of gravitational waves and electromagnetic waves
predicted by general relativity
Among the four basic interactions known in nature, the first to be discovered, described and modeled is universal
gravitation. The law of universal gravitation was first published in 1687. The law states that two large point-shaped
objects are attracted to each other by a force F~ , and the force paradigm is |F~ |. = Gm1 m2 /r2 is proportional totheir masses m1 and m2 , and inversely proportional to the square of their interval r, where G is a universal constant. Recall that this force originates from the local potential Φ. A common form of Newton’s law is the Poisson
equation:∇2 Φ = 4πGρ

(1)with ρ the mass density of matter, acting as the source of the gravitational potential Φ . Hence, in Newtonian
gravity, the gravitational interaction acts instantaneously. This was already of some concern to Newton himself,
but it clearly became a significant problem with the advent of Einstein’s theory of special relativity.
1.1.2

Special relativity

In 1887, Abraham Michelson and Edward Morley conducted an experiment aimed at detecting the relative movementof substances relative to luminescent ether. Using what is now called the Michelson interferometer, Michelson and
Molly measured the speed of light in two orthogonal directions from a common light source.
The result of this experiment is negative, because it obtains the same speed of light, regardless of the position
and motion of the earth around the sun. This posed a major problem in physics, and its solution led to the proposalof the theory of relativity in 1905. Einstein’s theory is based on the following two assumptions:(1)Principle of relativity: the equations describing the laws of physics have the same form in all inertial reference
frames;(2)Invariant light speed: in a vacuum, light propagates at a constant speed c, irrespective of the state of motionof the source.
Although the principle of relativity has been realized in Galileo and Newtonian mechanics, the second hypothesis
is a drastic amendment to our understanding of time and space.
A central concept of special relativity is the space-time interval between two events. Let ∆t, ∆x, ∆y and ∆z
denote the coordinate difference between the two events p and q relative to the global inertial reference system.
Then, the space-time interval between these events is∆s2 = −c2 (∆t)2 + (∆x)2 + (∆y)2 + (∆z)2(2)

The form of the interval (2) is quadratic in the differences of the coordinates, and invariant under the Poincare
group(translations, rotations, boosts [20]), thus ensuring that the speed of light is indeed the same in all inertial
frames. This observation suggests that, in full analogy with the Euclidean geometry of three-dimensional space,special relativity can be formulated as a theory of the Lorentzian geometry of four-dimensional spacetime.
In addition, the space-time interval can be used to explore the causal structure of space-time. See Figure 1.
Given the event p, the light cone Cp is the set of all events q, such that ∆s2 = 0. These events are called light
related to p, because all these events can be reached by light passing through p. All events in Cp make ∆s2 < 0. It
is said that these events are time-dependent with p, because the massive particles passing through p can reach anyof them at least in principle. The remaining events, that is, events other than Cp , make ∆s2 > 0. Those events are4

Figure 1: In special relativity , the causal structure of space-time defines a notion of light-cone Cp at any event p.
All events on Cp are lightlike related to p, while all events within (respectively, outside) Cp are timelike (respectively,
spacelike) related to p.[24]called p-related spatial events because they cannot be reached without mass particles or any light passing throughp. Two events related to space cannot have a causal relationship with each other.
1.1.3

Relativistic gravity

Next, let us explore whether gravitational interaction can adapt to the model of special relativity. For Poisson’sequation (1), we extend it to replace the elliptic Laplace operator ∇2 with the hyperbolic hyperbolic d’Alembert
operator, and replace the mass density ρ with the Lorentz covariate source . Therefore, we can get the gravitational
field equation as follows:2Φ = −

4πG
T
c2

(3)

1 ∂2
Where 2 = − 2 2 + ∇2 is the usual plane wave operator, and T is Trace of the energy momentum tensor ofc ∂tmatter. This scalar gravity theory follows the principle of special relativity and reproduces Poisson’s equation (1)in the non-relativistic limit, where c−1 → 0.
This shows that gravitational propagation needs to be combined with a limited speed, just like the form ofelectromagnetic wave propagation, gravitational waves should also be in the form of traveling waves and propagate
at the speed of light. Of course, the generation of gravitational waves is due to the vibration of space. There aremany ways to generate gravitational waves. The most common is the merger of binary star systems. Of course,
supernova explosions will also generate gravitational waves. This thesis mainly discusses the gravitational waves
generated in the merger of binary star systems, especially the gravitational waves generated by the merger of neutronstars.

1.2

Binary star system

Binary neutron stars(BNS) Neutron stars are one of the strangest objects in the universe. The mass of theneutron star is about 1.5 of that of the sun, and the radius is about 10 km. It is the most compact stable structure.
The degenerate pressure of the material can still balance gravity. Further compression will cause gravity to collapse
and form a black hole. Since gravity is extreme, its rotation is extreme: A neutron star is the fastest known rotatingstar with a period as short as milliseconds. When a neutron star and a normal star are in a binary star system,
its huge gravitational field can capture the material of the companion star accumulated on the neutron star. The
huge release of the gravitational force of the accretion material makes the system’s luminosity in the X-ray bandas high as 100,000 suns. Neutron stars also have very strong surface magnetic fields, ranging from 107 Gauss to1015 Gauss.The magnetic field can transfer the magnetic poles to the material near the surface. Under the actionof such a high magnetic field of neutron stars, considering that most neutron stars are born with a certain angular5momentum, this periodic magnetic field changes are modulated during rotation. Make the neutron star a pulsar.
Because neutron stars have a huge moment of inertia( 1045 g/cm2 ), neutron stars have become the most stable clockin the universe and beacons in cosmic navigation.
In the old binary star system, when the magnetic field of the neutron star decays to 108 ∼ 109 Gauss, theaccumulation of matter and angular momentum can accelerate the neutron star to a very fast spin cycle, which
can reach several hundred revolutions per second. The maximum rotation speed of a neutron star depends on its
compactness, that is, its mass radius ratio. In other words, it depends on the equation of state of the ultra-dense
matter in the neutron star. The search for the fastest spinning neutron star may impose important constraints on
the equation of state for the densest matter configuration in the universe.
In this research, the main research goal is the neutrino signal that can be received in the binary neutron starmerger (BNS), but because it is an automated system, when a gravitational wave signal is detected, we will also
proceed accordingly Neutrino detection.
Binary black hole(BBH) A black hole is a place in space where the gravitational pull is so great that evenlight cannot be emitted. The force of gravity is so great because the matter is squeezed into a small space. Thiscan happen when a star is about to die.
Because there is no light emission, people cannot see black holes. They are invisible. Black holes can only beobserved indirectly by observing the shadow of the black hole and the nearby accretion disk. This black hole photois the first photo of the black hole through the accretion disk around the black hole. Observing a black hole requires
a telescope with an aperture equivalent to that of the Earth to effectively observe it. The telescope used for thisobservation is the Event Horizon Telescope, which can be observed simultaneously by combining radio telescopes
around the world. The black hole photo shooting also won the Nobel Prize.
Because of the big mass, we can easier get the signal which from the black hole merger than BNS. But as weknow there are nothing can escape from the black hole, also neutrinos. So in our research, the expected of neutrinosignals which emitted by BBH is 0.

1.3

Neutrino

Pauli predicted in 1930 that neutrinos would retain the momentum of beta decay. Reines and Cowan conducted
experimental observations in 1956. Initially, neutrinos were considered massless. Now, since neutrino oscillations
are observed, they have mass. They have 3 flavors corresponding to light protein. Because their quality intrinsic
state and flavor intrinsic state are different, the flavor oscillates or changes during driving. In the Standard Model,
neutrinos are elementary particles. Since the interaction of gravity is very small, it can be ignored. It is electricallyneutral and will not interact with strong forces. Weak nuclear force is the main interaction channel. Neutrinos
pass through large objects with little interaction. Its experimental research requires large detectors and advanced
technology.
Neutrinos are produced through nuclear reactions, such as the fusion of the sun, the fission of the crust andmantle of the earth, and reactors. When cosmic rays interact with atomic nuclei in the atmosphere, they also produce
neutrinos, which produces unstable particles that decay. There are also artificial neutrino beams. The propertiesof neutrinos have been studied from these various sources. Neutrinos can also be used to study astronomical
events. Many astronomical events produce or are expected to produce large numbers of neutrinos. Due to the lowinteraction rate of neutrinos, it can bring different information to the earth without scattering or shielding. The
only astronomical neutrinos observed so far are from the sun and supernova 1987A. This thesis is to detect theneutrinos produced during the merger of celestial bodies when gravitational waves are generated, especially theneutrinos produced during the merger of neutron stars. According to different theoretical predictions and equationof state, the energy of neutrinos produced during the merger of neutron stars will be different to a certain extent (thefigure Figure 2 shows the prediction of neutrino production with different energies according to different equationof state). If the generated neutrinos can be detected The energy of neutrinos helps us to speculate on the reactions
that take place during the merger of neutron stars, and provides the corresponding data basis for the merger ofneutron stars.

6

Figure 2: Observed at infinity neutrino average energies for electron neutrinos and electron antineutrinos for threedifferent EoS[14]Figure 3: Super-Kamiokande Detector

2Super-Kamiokande(SK)2.1

Detector Overview

Super-Kamiokande (SK) is a 50kilo-ton Cherenkov water detector located in the Kamioka mining area. The cosmic
ray muon reduced the 1,000m rock shield (2700 m.w.e.). The detector is cylindrical with a diameter of 39.3 m and aheight of 41.4 m. The stainless steel frame structure divides the detector into two parts, the inner detector (ID) andthe outer detector (OD). The structure forms an optical barrier and supports the photo-multiplier tube (PMT). In
the OD, there are 1885 8-inch PMTs that monitor particles entering from the outside, and 2.5 m of water shields
gamma rays from the rocks surrounding the detector.
The inner side of the frame structure is surrounded by a polyethylene terephthalate black sheet to reduce opticalreflection. The ID is formed by 11,129 20-inch PMT, providing a photo coverage of ∼ 40%. The volume within 2
m from the ID wall is not used for SK analysis. The remaining volume of 22.5 tons is called the fiducial volume(FV). The space above the detector is a research dome with 4 electronic huts and a central control room.

2.2

SK history

SK’s predecessor, the Kamioka Nuclear Decay Experiment (KamiokaNDE) detector, has a water tank with adiameter of 15.6 m and a height of 16 m. It contains about 3 tons of pure water and 948 tons of PMT, 123 tonsof PMT in the OD. Observations were made from July 1983 to April 1996, during which time the neutrinos ofthe supernova (SN1987A) were observed. As the successor to the KamiokaNDE experiment, the construction of
the SK detector began in 1991. For SK storage tanks, holes in the mine were dug until mid-1994. Thereafter, theconstruction of the water tank began and the PMT was installed until the end of 1995. Water injection started in7

January 1996. The experimental data was collected after a month of test operation, starting from April 1, 1996.
From the beginning of operation to the present, SK is divided into six stages (ie SK-I, SK-II, SK-III, SK-IV, SK-V
and SK-Gd (SK-VI)).
1. SK-I started observation on April, 1996. This period ended on July, 2001.
2. After regular maintenance in July, 2001, a shock wave from the implosion of a PMT caused the chain implosionof almost half of the PMTs. Then, the remaining PMTs were re-distributed and SK-II started on October, 2002.
From SK-II, PMTs are protected from shock waves by acrylic covers and fiberglass reinforced plastic (FRP). SK-II
was stopped to mount new PMTs on October, 2005.
3.The data cycle from July 2006 to August 2008 is SK-III. In October 2005, SK re-installed the PMT andcompleted it in July 2006. SK-III has 11129 ID PMTs, and the coverage rate of photos is 40%.
4.Since September 2008, as all front-end electronics have been upgraded, the detector has been shut down. Thencontinue with SK-IV from September 2008 to June 2018. During this period, the T2K experiment using SK as afar detector was started, and the first event of the T2K experiment was successfully observed in SK on February
24, 2010.
5. In preparation for the SK-Gd phase, SK repaired the leak of the water tank, installed new pipes, replaced
the failed ID and OD PMT, and cleaned the water tank. This stage is SK-V, which starts in January 2019 andends in July 2020.
6.The SK-Gd phase has stared in July 2020, which dissolve gadolinium in ultra-pure water at SK. Gadolinium
has the largest cross-section of thermal neutron capture among natural elements, and a total of about 8 MeV
gamma rays are emitted. Therefore, it is possible to distinguish between neutrino reactions with neutrons andthose without.

2.3

Detector Principle

SK detects Cherenkov light, which is emitted when charged particles travel faster than the speed of light in themedium. Cherenkov light radiates in a cone shape along the path of charged particles. The angle between the
direction of the charged particle and the direction of Cherenkov light is called Cherenkov angle θc , which can becalculated by the following formula1(4)
cos θc =βn(λ)
Where β is the speed of the charged particle relative to the speed of light in vacuum, and n(λ) is the refractive
index at the wavelength λ. In pure water, the refractive index at 589 nm is 1.33. For relativistic particles likeelectrons and positrons, θc is 42° also a muon with 1 GeV. The energy threshold of Cerenkov light emitted isdetermined by the following formula:
Ethr = pm
1 − (1/n)2(5)

where m is the rest mass of the charged particle. The number of Cherenkov photons along the trajectory per
wavelength is calculated byd2 N
2πz 2 α2πz 2 α1=)=sin2 θc
(1
−
dλdx
λ2
β 2 n2 (λ)λ(6)

where z is the particle charge in units of e and α is the fine structure constant.

2.4

Photomultiplier Tubes

The 20-inch PMT (R3600) [29]for the ID was developed by Hamamatsu Photonics. The schematic is shown in
Figure 4.
Its photocathode is made of bailkali (Sb-K-Cs) and is sensitive to photons with wavelengths of 300 to 600
nanometers. Its quantum efficiency is about 20% at the maximum(Figure 5). Photoelectrons (pe) generated by thephotoelectric effect are amplified 106 ∼ 107 by the dynode chain. The PMT dark noise is 4 kHz. The Helmholtz
coil around the detector reduces the geomagnetic field that affects the PMT response from 450 mG to 50 mG.

8Figure 4: A schematic of the Hamamatsu R3500 PMT.[19]Figure 5: The quantum efficiency of the Hamamatsu R3600 as a function of wavelength.[19]9Figure 6: The SK water purification system.[19]2.5

Water Systems

The water purification system (Figure 6) was developed to maintain the quality of water transparency and radioactivity. The higher transparency gives less attenuation of the Cherenkov light. Radioactivity can produce
backgrounds in the low energy analysis.
The water system consists of the following components:•1µm filter: Removes large particles.
• Heat exchanger: Maintains the water temperature at 13◦ C. Water temperature variations cause convection,
while higher temperatures increase PMT dark noise and the growth of bacteria.
• Cartridge polisher: Removes heavy ions.
• UV sterilizer: Kills bacteria.
• Vacuum de-gasifier: Removes oxygen and radon in water.
• Ultra filter: Removes particles larger than 10 nm.
• Membrane de-gasifier: Removes dissolved gas.
• Reverse osmosis: Removes large particles with molecular weight> 100.
The water circulates at a rate of 60 tons per hour. It fills from the bottom of the detector and drains from thetop (Figure 7).
The water temperature is monitored in the ID and OD. It is uniform from the bottom to 11m below the centerand slightly increases with height (Figure 8). The difference between top and bottom is 0.2◦ C .

2.6

Radon Free Air Systems

The air in the mine includes radon in the surrounding rocks. Radon will affect the low energy background rate.
SK’s system can draw in air from outside the mine and reduce radon in the air. The schematic diagram of the air
system without the radon is shown in Figure 9 . It consists of the following components:• Compressor: Provides air pressure up to 7-8.5 atm.
• Air filter: Removes dust in the air. The air goes through 3 types of filters, 0.3 µm, 0.1 µm, and 0.01 µm ,from large to small.
• Air drier: Removes water and CO2 in the air.
• Carbon column: Absorbs the radon in the air.8 m3 of charcoal is used in total.
10Figure 7: The SK water flow.[28]Figure 8: The z-dependence of the water temperature.[28]11Figure 9: The radon free air system.[19]• Cooled charcoal: Has better random removal efficiency. The volume is 50 L.
The radon level increases to ∼ 2000 Bq/m3 in summer and decreases to ∼ 260 Bq/m3 in winter. The air system
without radon is maintained in the experimental area ∼ 40 Bq/m3 . In order to prevent the air with radon in themine from entering the tank, SK will make the pressure of the radon free air in the tank higher than the outsidepressure (0.3kP a

2.7

Front-end Electronics and Data Acquisition System

The front-end electronics and data acquisition system (DAQ) was upgraded in August 2008. This upgrade canachieve more stable data collection and greater data throughput.
Each PMT hit is recorded without loss and selected by a trigger. Because SK has multiple physical purposes,
multiple types of triggers are installed. Low energy (LE) trigger, high energy (HE) trigger and super high energy(SHE) trigger have a timing window of 40 µs. Before the trigger, 5 µs is stored for pre-event γ research, and thesubsequent 35 µs is used for post-event research. At the beginning of a new DAQ period, the trigger threshold is
47, 50, and 70 hits in 200 nsec for LE, HE, and SHE, respectively. The SHE threshold was changed to 58 hits inthe summer of 2011.
Figure 10 shows the schematic of the DAQ system. The signals from PMTs go to the new front-end electronics,
QTC (charge-to-time converter) Based Electronics with Ethernet (QBEE, Figure 11). The QBEE has a 24 analog
input channel with a dynamic range of 0.2 to 2500 pC. It works 5-times faster than the previous electronics[28]. 8
QTC chips mounted in the QBEE integrate the input charge and output a pulse of proportional width. The pulse
width and timing is digitized by a time-to-digital converter (TDC). The digital information is processed by a Field
Programmable Gate Array (FPGA)[27].
TCP(Transmission Control Protocol) which is one of the main protocols of the Internet protocol suite packets
from 500 QBEE are sent to 20 front-end PCs. The front-end PC will sort the data in time, and then send the data
to 10 combined PCs. The merged PC also sorts the data by time, and then applies the event construction process
to the data in the same time zone. During the event building process, the trigger will scan data from 20 front-end
PCs. The triggered event will be sent to the organizer’s PC and recorded on the disk. Gigabit Ethernet is used fordata transmission between PCs.

2.8

Monitoring system

The status of data acquisition is monitored by various tools. The ”Slow Control” monitor checks the high voltage
status and the temperature of the electronic equipment. The data flow, PMT conditions and trigger rate aremonitored by the distributed process. It also monitors data transfer and offline processes. The online program
displays visual images of events in real time. This display helps to monitor DAQ system and PMT.

12

Figure 10: The schematic of the DAQ system.[32]Figure 11: QBEE.[28]13Figure 12: LIGO Livingston

3

LIGO/VIRGO Detector

3.1

Introduction

LIGO is the abbreviation of Laser Interferometer Gravity Wave Observatory. It is currently the most sensitiveand largest gravitational wave detector in operation in the world. LIGO first observed gravitational waves in 2015
and won the Nobel Prize for it. LIGO has two observation points, which are located in Hanford, Washington,
and Livingston, Louisiana, USA, which are 3000km apart. Simultaneous observation of the two observation points
can greatly reduce the possibility of misjudgment of gravitational waves, and pass The time difference between
the arrival of the gravitational wave signal and the distance between the two can roughly determine the sourceand direction of the signal. LIGO mainly uses the Michelson Observer to observe the gravitational wave. The
observation principle will be described in the next section.
LIGO has super large interference arms (each interference arm has 4km), which can improve the accuracy ofobservation. Since the observation of gravitational waves is different from other astronomical observations, and there
is no electromagnetic interference in traditional observations, the components of LIGO’s detectors are completely
isolated. The detector consists of 1.2m wide steel vacuum tubes arranged in an L-shape and covered by a concrete
bunker 10 feet wide and 12 feet high.

3.2

Detector Principle

The LIGO detector is designed to detect the differential strain of the incident GW based on a simple Michelsoninterferometer (see Figure 12). This design starts from the input laser incident on the beam splitter located at theorigin, and its face forms an angle of 45 degrees between the x-axis and the y-axis. The beam splitter directs 50%of the light in two vertical directions x and y. The light travels to the mirrors located at the end of the arm (Lx , 0)and (0, Ly ), which redirect the light back to the beam splitter. The cumulative phase of the light from the beam
splitter to the end mirror is given by:
Z L
Φ=
0

2πL
2π
dx =
λ
λ

(7)

Where L is the distance traveled and λ is the wavelength of light. The total phase of such a round trip willbe the sum of the integral drop and the integral return. When the light returns to the beam splitter, the two
beams interfere with each other under such conditions.In this way, depending on the length of the differential arm,
the light will either return to the laser (symmetric port), or toward the photodetector (antisymmetric port), or a14

combination of both. Then, by looking at the power of the light from the antisymmetric port PAS , the differential
arm length Lx − Ly can be measured in a way that the mode is a wavelength. This is caused by
4π(Lx − Ly ))(8)
λWhere P0 is the input laser power, φx is the cumulative phase of the beam in the x arm, and φy is the cumulative
phase of the beam in the y arm. This shows that PAS depends on the length of the differential arm. Let us assume
that the length of the arm is the same L and the amplitude of the incident GW is h. It stretches the detector along
the x-axis and squeezes it along the x-axis. This effect causes The distance the light needs to travel is given by
L0 ≈ L(1 ± h/2), and φx and φy are
PAS = P0 sin2 (φx − φy ) = P0 sin2 (

Z L0

Z 0
2π
2π
dx −
dx
λ
0
L λ
0
4πL
h
≈
(1 + )
λ
2
4πL
h
φy ≈
(1 − )
λ2With these phase shifts, we find that PAS is then given by
φx =

(9)

(10)

hL 2
)
(11)
λ
Thus, for the same change of PAS , the strain h is inversely proportional to the length L of the arm. However,
although the antisymmetric port of this interferometer has a signal, it is a The method of detecting GW is poor.
First, since the above power on the asymmetric port is proportional to sin2 (4φ), Taylor can be expanded toapproximately 4φ2 about zero, so the signal we are measuring It is square and proportional to the small GW signal
we want to sense. In addition, there is no information that we need to move the mirror to bring the signal back to
the dark fringe.
In order to solve these problems, LIGO adopted the ”heterodyne detection” technology, which combines thephase modulation input laser with Schnupp asymmetric technology. In this setting, before the laser enters theinterferometer, the light is phase-modulated at the frequency Ω. In this way, instead of the three superimposedlaser beams incident at three frequencies, they all enter the detection at one frequency. The electric field amplitude:
the original beam with frequency ω and amplitude E0 , also called carrier, and two sidebands with frequency ω ± Ω
and amplitude E1 .
When the arm lengths are exactly equal, all these beams will exit at the symmetrical port, which will produce
the same answer as above. However, there will be a different response when Schnupp asymmetry is introduced
between the two arms. This happens when the unequal amount of the arm of the interferometer is an integer
multiple of the wavelength of the carrier, rather than an integer multiple of the wavelength of the sideband. In thiscase, in the absence of a GW signal, the carrier light will be emitted from the symmetric port, and the sideband
will be transmitted to the antisymmetric port. When a GW signal is present, a phase shift will be introduced on
all three beams so that a mixture of all three beams will exit the antisymmetric port. Performing demodulation
processing on the signal from the asymmetric port as needed will cause the error signal to be linear in the GW
amplitude.
In order to improve the sensitivity of the detector, LIGO modified the standard Michelson interferometer twice.
The first is to use Fabry-Perot cavities in the two arms to increase the storage time of light in the arms, effectively
increasing the length of the arms, thereby increasing the phase shift of the light by the same factor. The second is
to increase the power recovery mirror to enhance the light circulating in the interferometer under a given laser input
power. The design schematic diagram of LIGO, including all the additional functions of the standard Michelsoninterferometer, is shown in Figure 13.
PAS = P0 sin (4π3.3History

Till now, there are three phases in LIGO history called O1,O2,O3.
O1 It started on September 18, 2015 and ended on January 12, 2016. The data from the surrounding engineering
period is of sufficient quality to be included in the analysis, which means that the observation data was collectedfrom September 12, 2015 to January 19, 2016. The operation involves Hanford (H) and Livingston (L) detectors.
15

Figure 13: A simple schematic of the LIGO detectors. The mirror labels are as follows: PRM (power recycling
mirror), BS (beam splitter), ITMX/Y (x-/y- arm input test mass), ETMX/Y (x-/y- arm end testmass), AS (antisymmetric port). The width of the laser beam denotes power in different portions of the detector. A simple
Michelson interferometer includes the mirrors BS and ETMX/Y. The addition of the ITMX/Y mirrors makes useof Fabry-Perot cavities for the arms. The use of the PRM boosts the power circulating in the whole interferometer16

Figure 14: LIGO phases
O2 The project started on November 30, 2016 and will end on August 25, 2017. Prior to this, the project began on
October 31, 2016 in Livingston and on November 14, 2016 in Hanford. Hanford’s delay was to facilitate additional
commissioning activities. The sensitivity achieved during the whole operation is usually in the BNS range of 80-100
Mpc. AdV interferometer joined O2 on August 1, 2017, forming a network of three detectors in the last month ofoperation. The target is a BNS range of 40 Mpc. Since the vacuum contamination problem has been solved, AdV
uses steel wire instead of fused silica fiber to suspend the test substance. This limits the maximum BNS range of
AdV; in O2, the BNS range reaches 30 Mpc.
O3 The third observing run started on April 1, 2019 and was expected to end on April 30, 2020, with a commissioning break from October 1, 2019 to November 1, 2019. While this article was in review the COVID-19 Pandemic
led to suspension of the observing run on March 27, 2020. The increase in sensitivity of the LIGO detectors (whosetarget sensitivity was expected to be 120 Mpc) comes from a variety of changes, chiefly from increasing the input
laser power, adding a squeezed vacuum source at the interferometer output and mitigating noise arising from scattered light. Additionally, end test-mass optics with lower-loss coatings, along with new reaction masses, have been
installed in each interferometer. The Livingston instrument began the run with an average BNS range of 130 Mpc
and the Hanford instrument typically operates with an average range of 110 Mpc.Fused silica fibers were installed on
the AdV test mass suspensions in preparation for O3. Other improvements included reduction of technical noises,
increasing the input laser power and installation of a squeezed vacuum source. The result was a BNS range of 50
Mpc at the start of O3.

3.4

Signal Display

In the gravitational waves detection, the most important indicator is frequency.By using Numerical Relativity,weµm1 m2
can define a Chirp-Mass[12] as M = M η 3/5 which η =,µ =, M = m1 + m2 .And we can also computeMM
this chirp-Mass from frequency likeM=

c3 5 −8/3 −11/3 ˙ 3/5
[ π
f
f]
G 96

(12)

where f and f˙ are the observed frequency and its time derivative and G and c are the gravitational constantand speed of light.By using this equation, the chirp-Mass can be computed by frequency. For GW150914, LIGO
estimates f and f˙ based on the data in Figure 16[4]. We get the chi mass of M ' 30M , which shows that the
total mass M = m1 + m2 in the detector frame is > 70M . This limits the sum of the Schwarzschild radii of the17

binary components to 2GM/c2 > 210km. In order to reach the 75 Hz orbital frequency (half the frequency of thegravity wave), objects must be very close and very compact; orbits equal to the mass of the Newtonian point are
only 350 km apart at this frequency. Although a pair of neutron stars is compact but will not have the requiredmass, a black hole neutron star binary with inferential chi mass will have a very large total mass and therefore will
merge at a much lower frequency. This leaves a black hole because the only known object is compact enough toreach an orbital frequency of 75 Hz without contact.
Next, we can use the max amplitude to get the distance of the source byhmax =

4(GM)5/3 (πf )2/3
4G2 M µ
=
4
c Dr
c4 D

(13)

Where D is the distance , f is frequency, G is gravitational constant,M is chirp-Mass.
For the emitted energy, we can compute by the wave function thatdE
2c3 dA 2=() , h = A(u) sin2 θ/rdt
15G duAnd we can also compute merger time by using those parameters liket0 =5256(Gc−3 M5/3 (πf0 )8/3 )

(14)

(15)

From Figure 15 ,comparing the t0 of each signal ,So we can easily distinguish Binary Neutron Star or Black
Hole Merger.
For each GW event , we get information from frequency and merger times and we can compute distances ,
mass,and so on.

18

Figure 15: merger time for Binary Neutron star(up) and Black hole Merger(down)

19

Figure 16: The gravitational-wave event GW150914 observed by the LIGO Hanford (H1, left column panels) and
Livingston (L1, right column panels) detectors. Times are shown relative to September 14, 2015 at 09:50:45 UTC.
For visualization, all time series are filtered with a 35–350 Hz bandpass filter to suppress large fluctuations outside
the detectors’ most sensitive frequency band, and band-reject filters to remove the strong instrumental spectral
lines seen in the Fig. 15 spectra. Top row, left: H1 strain. Top row, right: L1 strain. GW150914 arrived first+0.5
at L1 and 6.9−0.4
ms later at H1; for a visual comparison, the H1 data are also shown, shifted in time by thisamount and inverted (to account for the detectors’ relative orientations). Second row: Gravitational-wave strain
projected onto each detector in the 35–350 Hz band. Solid lines show a numerical relativity waveform for a systemwith parameters consistent with those recovered from GW150914 [26][17] confirmed to 99.9% by an independentcalculation based on[15]. Shaded areas show 90% credible regions for two independent waveform reconstructions.
One (dark gray) models the signal using binary black hole template waveforms[6]. The other (light gray) does
not use an astrophysical model, but instead calculates the strain signal as a linear combination of sine-Gaussian
wavelets[17][5]. These reconstructions have a 94% overlap, as shown in[6]. Third row: Residuals after subtracting
the filtered numerical relativity waveform from the filtered detector time series. Bottom row:A time -frequency
representation [16]of the strain data, showing the signal frequency increasing over time.[3]20Figure 17: (a): Location and orientation of the LIGO detectors at Hanford, WA (H1) and Livingston, LA (L1).(b):The instrument noise for each detector near the time of the signal detection; this is an amplitude spectral density,
expressed in terms of equivalent gravitational-wave strain amplitude. The sensitivity is limited by photon shot noiseat frequencies above 150 Hz, and by a superposition of other noise sources at lower frequencies[3]. Narrow-bandfeatures include calibration lines (33–38, 330, and 1080 Hz), vibrational modes of suspension fibers (500 Hz andharmonics), and 60 Hz electric power grid harmonics.[3]4Automatic system and methods

4.1

Alert receiver and alert database

4.1.1

GCN Notice

Alerts of the gravitational wave detection by LIGO and Virgo are public since the beginning of their third operation
period, so-called O3. They are distributed to the follow-up observers of these events via the Gamma-ray Coordinates
Network (GCN). GCN has been originally developed by NASA for quick distribution of gamma-ray burst alertsfrom satellite experiments to the follow-up telescopes on the ground, and now it extends its activity into many typesof astronomical transient events, high-energy neutrino events and supernova neutrino events, as well as GW events.
The information and results of the follow-up observations are shared with the registered individuals and institutes
via GCN also. The former function, alert distribution, is called as Notice, and the later one is called as Circular.
GCN Notice messages are distributed via a packet connection between a GCN server in NASA and a clientmachine of a follow-up observer. Several protocols, 160 bytes binary and XML text, are available for Notice, andthe format is well defined for each type of the alerts. The messages of GCN Circulars are distributed via E-mail,and they are human-readable descriptions about follow-up results. All Notice and Circular messages are archived
in the official web page of GCN.
4.1.2

GCN Notice receiver

A receiver of GCN Notice is implemented by using the python library, gcn. The library provides functions for
initialization of a packet connection to the NASA server and listening of notices, and they enable us to implementa GW alert receiver by a very simple code. A GW notice message is received as the VOEvent format written in
XML, and we can extract several information from the message such as event ID, trigger time, classification of
the event, a type of the alert, and links to a web page and a fit file of the event in GraceDB (describe later).
These extracted information of the event is filled out into a data base written with a CSV format, and the original
VOEvent message(Figure 18) is also saved into a single file.
As an example of the notice, Fig. 18 shows the parameters of a first GW event in O3, which is publishedon April 8, 2019. Each event is given an unique identification number, such as ”S190408an” (see the parameter,
TRIGGER NUM). Usually more than one notice are published for a single GW event based on a LIGO/Virgo alert
sequence. The first notice is called as Preliminary. It is sent fully automatically within some minutes from theGW event. The others are Initial and Update. They include additions or updates of the information of the event
considering results of other analysis pipelines, and they are published in a few hour and a few days for Initial andUpdate, respectively. The type is found in the NOTICE TYPE row in Fig. 18. Figure 19 describe the alert time line21

of LIGO/VIRGO. Sometimes, they find that the event already published by a Preliminary notice is background. In
this case, Retraction notice is sent in order to inform it to the follow-up observers. In case we receive a Retraction
notice, we remove the event from our database. A GW notice also includes the information of GW event type, suchas BBH, BNS and NSBH(neutron star black hole merger), which is provided as event probability of each type (seePROB XXX). The most important parameter for our analysis is the time of the event, given in TRIGGER TIME.
Another important parameter, especially for the High Energy analysis, is the location of the event in the sky. This
is not provided as a simple parameter but a probability-density map of the event location on the sky. It can beaddressed via a Fits-formatted file downloaded from the link provided in the notice, SKYMAP FITS URL. As an
example, the probability density distribution of the GW event S190408an is shown in Fig. 20. A header of the fits
file contains some additional information of the event such as the estimated distance between the GW object andthe earth, as well as the basic information of the sky map, e.g. the pixel size.
The all public information of the GW candidate events is achieved into a dedicated database for GW, GravitationalWave Candidate Event Database (GraceDB [25]). In the web page, all Notice messages, sky map, and reports offollow-up observations for each event are available. If a notice is missed due to a maintenance of the software, theSK server and etc., we can download the notice message from the GraceDB web page (or the NASA GCN webpage). An automatic script to do it have already been prepared.

Figure 18: The message of GCN Notice for the GW event, S190408an [11]. This was a first notice for the event in
O3, and the information was updated by the second notice labelled as Initial, which was published 1.5 hour later
than the first notice.

Figure 19: The types of GW alerts and related timelines22Figure 20: The probability density map of the GW event, S190408an[25]4.2High Energy analysis

The High Energy(ATMPD which notates ATMospheric neutrino and Proton Decay analysis)framework is built
independently of the rest of the framework. It consists in a watcher script waiting for new gravitational wave eventsin the alert database. When it is received, the script is then looking for the processing of the Super-Kamiokande
data of the few minutes around the time of the external trigger, there are three possible scenarios:• data already processed: the High Energy analysis can be launched directly• data not processed: the script is waiting for the data to be ready. More precisely, we wait till the last processedrun (according to the run summary file on the analysis machine) has an end time later than (talert + 1000 s).
When this condition is fulfilled:– if there is no run covering the ±500 s time window centered on the external trigger, the High Energy
analysis is anyway launched, in order to produce the minimal output.
– if there are some runs covering the time window, we wait for the end of the processing of FC, PC
and UPMU files

SK data

Expectedbackgroundwaiting for
online reduction

New alert

Watcher

Compute> 0s
livetime= 0s

END

Find data
eventsFlux limits in any caseSignal significance if observed events

Distributions,
event displays. . .

Figure 21: Overview of the High Energy analysis

23

4.2.1

Event selection

In High Energy region , we have three samples:• Fully-Contained (FC): neutrino interacts in the inner detector and the outgoing lepton has a short range,leading to a clear ring in the inner detector and nothing in the outer detector. The typical neutrino energy
range is 0.1 GeV to 10 GeV.
• Partially-Contained (PC): neutrino interacts in the inner detector and the outgoing lepton exits the detector,
leading to a disk in the inner detector and one cluster in the outer detector. The typical neutrino energy
range is 0.1 GeV to 100 GeV.
• Up-going muons (UPMU): neutrino interacts in the rock or in the water in the outer detector. One or two
clusters are detected in the outer detector as well as light in the inner detector. Most of the events arethrough-going, which means it can be fitted with a straight line crossing the full detector. Only events comingfrom above horizon are selected, in order to reject most of the down-going atmospheric muon background.
The typical neutrino energy range is 1.6 GeV to 10 000 GeV.
When an external trigger comes, only the files containing the runs covering the 1000 s time window are processed.
As runs are typically 24 h-long, we need to further pre-select the events within the time window. This is done byusing jd and fjd variables (where jd and fjd are respectively the integer and the fractional part of the Julian dateof event).
The final selection cuts applied in the analysis are following the standard cuts used in the official ATMPD
software:• FC:– wall>100 (cut on the distance to the wall in cm)– evis>30 (cut on visible energy in MeV)
– nhitac<16 (cut on number of OD hits)• PC:– wall>200 (cut on the distance to the wall in cm)– evis>350 (cut on visible energy in MeV)
– nhitac>=16 (cut on number of OD hits)• UPMU:– Fit directional z>0 (up-going events)– Stopping muons: Fit pid==1 and Fit mom>1.6 (cut on fitted momentum in GeV)
– Or through-going muons: Fit pid==2 and Fit len>700 (cut on fitted track length in cm)
We store all events in the time window of each sample in ROOT files, and make official sample reductions. This
paper will use these samples for analysis, regardless of whether they passed the official sample reduction or failed(because the latter can be used for further inspection or adjustment). In addition to the official variables, three
new variables have been added:• selected is a boolean variable indicating if the event passed the official sample cuts or not. The events withselected==True are the ones to be considered for final analysis.
• delta t is the time difference [in seconds] between the neutrino event and the external trigger time.
• sigma ang is the expected 1 − σ angular resolution of the event. It is estimated as a function of reconstructedenergy.

24

4.2.2

Background estimation and characteristics

The three High-Energy samples are dominated by atmospheric neutrinos, that will therefore be the dominantbackground for the analysis. Fortunately the rate of such events is quite stable with time and their angular
distribution only depends on zenith angle θ. However, it will have to be carefully treated in the subsequentstatistical analysis, as the Earth rotation makes the background at a given position (δ, α) is varying with time,
because the corresponding θ changes with time .

sin δ = cos φSK cos a cos A + sin φSK sin a
sin a − sin δ sin φSKcos h =cos δ cos φSK− sin A cos asin h =cos δ
α = θL − h(16)

For the integrated expected number of background event, due to its stability, it can easily estimated by multiplying Super-Kamiokande livetime in the selected time window by the average background rate.
Background: For each sample, the background rate is estimated using the full dataset of the current SK period,
where background is assumed to be stable. This is done by counting the total number of events passing the selection
cut and dividing by the total livetime.
For the real-time processing (and therefore during SK data acquisition), the full SK dataset is not available. We
then use a sufficiently large sub-sample of data that have already been processed. When the dataset is officialised(after e.g. eye-scans and use of official calibration constants), we can use it in order to estimate the background
event rate. Results are shown inTable 1.
Table 1: Background event rate using realtime processing (runs 80027 ∼ 82389) and official SK5 dataset of May
2020 (runs 80031 ∼ 82585) with the associated statistical errors.

Realtime

Official

SampleEvent rate [s−1 ]Error on event rate [s−1 ]

FC
PC
UPMU
FC
PC
UPMU

1.109 × 10−4
7.262 × 10−6
2.102 × 10−5
1.120 × 10−4
7.279 × 10−6
1.637 × 10−5

0.021 × 10−4
0.540 × 10−6
0.092 × 10−50.020 × 10−40.490 × 10−6
0.073 × 10−5

Livetime: In SK analysis,we usually take 24 hours data as a normal run, and divide into subruns. The livetime isestimated using the specific run(s) overlapping with the selected time window. For each run, contains the computed
livetime Ls and deadtime Ds for each subrun, but the subrun start time ts is only precise to the minutes andtherefore it is not enough to know which fraction of a given subrun is overlapping with the time window.
Taking into account the fact that each subrun is approximately 2 min-long, we select the subruns that are locatedwithin few minutes of the selected time window [t1 , t2 ]:t1 − 3 min < ts < t2 + 1 min(17)

For these subruns, the reformatted root files are used to extract the precise start and end time of the subrun(using GPS time, precise to the nanoseconds). We can then compute the fractional overlapping Os of the subrun
with the time window. We have Os = 0 (resp. Os = 1) if the subrun is completely outside (resp. inside) the time
window.
The total livetime of Super-Kamiokande over the selected time window [t1 , t2 ] is:X
L=Os × Ls(18)
s25

Figure 22: Angular and energy dependence of atmospheric neutrino background, using SK4 Monte Carlo. The
distributions are normalised to unity.
Background characteristics: For the statistical analysis, it will be needed to quantify if a given event in SuperKamiokande is compatible with background or not. This is done by using two discriminatory variables: directionand energy. Background will have one behaviour (atmospheric-like spectrum, isotropic angular distribution) while
signal will behave differently (different spectrum, pointing to the source).
The Fig.21 shows the angular and energy distribution for background using SK4 Monte Carlo and the selection
cuts described in Sec.4.2.1.
4.2.3

Signal characteristics

In order to extract some limits, sensitivity and significance from the observations, it is needed to characterise thebehaviour of the expected astrophysical neutrino signal.
Neutrino effective area:
events in the detector.

Neutrino effective area allows converting from an incoming flux to a number of expectedZd2 nnsig = Aeff (E, Ω)× dEdΩ(19)dEdΩ

For an instantaneous fluxdn
[cm−2 GeV−1 ] = φ0 [cm−2 GeV−1+γ ]E −γ from a source at position Ωsrc :dEZnsig = Aeff (E, Ωsrc ) × φ0 E −γ × dE

(20)

The effective area(Aeff) depends on the considered sample, on the considered neutrino flavour (νµ , νe , ν̄µ , ν̄e ),on neutrino energy and on direction (more precisely the zenith angle θ, for simplicity we neglect the dependencywith azimuth)
The effective area can be estimated using Monte Carlo simulations. For Super-Kamiokande, the biggest available
sample is the standard atmospheric MC, corresponding to 500-years-equivalent. In the analysis, we use the SK4
Monte-Carlo.
For a given neutrino flavour, we know the atmospheric neutrino flux Φatm (E, θ) [GeV−1 cm−2 s−1 sr−1 ]. For agiven bin [E, E + dE], [θ, θ + dθ] (integrating over azimuth angle) in terms of truth information, for a given livetimebin
T (e.g. 500 years), we then define Mtrue[in cm−2 ]:
bin
Mtrue
= Φatm (E, θ) × dE × sin θdθ × 2π × T(21)

For a given sample and neutrino flavour, we can look for the number of events after simulation + reconstructionbin+ selection in the same bin (truth information) Nselected. The effective area for this particular sample is estimatedto be:
Af,seff (E, θ) =26

binNselected
bin
Mtrue

(22)

By averaging over the zenith angles[?] [?], we can get the effective area as a function of energy only, as it isoften pictured in astrophysical publications(e.g[1]).
FC, νµ
−3

10

2.5

3.0

10−4

2.5

2.0

10−5

2.0

1.5

−6

10

1.5

1.0

10−7

1.0

0.5

10−8

0.5

0.0
−1.0

10−9
−0.5

0.0

0.5

PC, νµFC, νµ

3.0

1.0

1.5

10−3

3.0
10−3

2.5
10−4
2.0

10−4
10−5

10−6

1.5
10−5

1.0
0.5

10−6

10−7
0.0
−1.0

−0.5

0.0

(a) FC νµ

0.5

1.0

0.0

1.5

0.2

0.4

0.6

(b) FC ν̄µFC, νe

0.8

−2

10−3

1.6

1.8

102
1

3.0
−3

10
2.5

2.0

1.4

UPMU, νµ

3.0

10
2.5

1.2

(c) PC νµPC, νe

3.0

1.0

10−5
2.5

10−10

−4

10
2.0

10−15

2.0

10−20
1.5

−4

−5

1.5

10

1.5

1.0

10−6

1.0

10−25

10
1.0

10−5

0.5

0.5

10−30
10−35
10−40

0.5
10−7

0.0
−1.0

−0.5

0.0

0.5

1.0

(d) FC νe

1.5

0.0

0.2

0.4

0.6

0.8

1.0

1.2

(e) PC νe

1.4

1.6

1.8

10−45
0.0
0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

(f) UPMU νµFigure 23: Effective area as a function of neutrino energy (log(E[GeV] on x-axis) and zenith angles (rad on y-axis)
for a selection of samples and neutrino flavours.
We see a small dip in the FC plot, around log(E[GeV])=-0.6 corresponding to E ∼ 250 MeV.
Angular resolution: For a given neutrino detection in Super-Kamiokande, the reconstructed direction is notpointing exactly to the neutrino direction. There are two reasons for that:• During a neutrino-nucleus scattering, the outgoing muon is not aligned with the incoming neutrino direction
because there is a momentum exchange with the nucleus. This effect is particularly important for the lower
momentum, as shown in Fig.25: 30°at 1 GeV and 0.01°at 3 TeV.
• The Super-Kamiokande reconstruction is not perfect and therefore the angle between the reconstructed direction of the muon and its true direction is not zero. This effect depends on the considered sample but is more
independent of the initial neutrino energy.
As the goal in this analysis is to reconstruct the position of the source (a.k.a. the neutrino direction) based on the
only information we have, i.e. the direction of reconstructed event. Therefore we only care about the compositionof the two effects: resolution for the angle between true neutrino and reconstructed direct of the event.
In order to estimate the resolution, we use SK4 Monte Carlo. For each sample (FC, PC, UPMU), we apply thenormal event selection and plot ∆Φ = Angle(dtrue; drecoνevt ) vs reconstructed energy (visible energy for FC/PC, fitted
momentum for UPMU).
Then, for different bins in reconstructed energy (sufficiently large to contain enough statistics, we project theprevious plot in one-dimension:Fig.26. This allows to obtain the p.d.f. describing the distribution of ∆Φ for a given
reconstructed energy.
This will be of particular importance in two occasions:• quantify the agreement between the reconstructed direction of a given event and an assumed source direction• simulate fake signal event from a given point source, taking into account resolution effect27

effective area [cm2]

10−2
10−3
10−4
10−5
10−6

νµ
νµ
νe
νe

10−7
10−8
10−9
10−10
−1.0

−0.5

0.0

0.5

1.0

1.5

2.0

log (Eν[GeV])
10effective area [cm2](a) FC

10−2
10−3
10−4
10−5
10−6

νµ
νµ
νe
νe

10−7
10−8
10−9
10−10
0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

log (Eν[GeV])
10effective area [cm2](b) PC

103
102
10
1
−1

10

10−2

νµ
νµ
νe
νe

10−3
10

−4

10−5
10−6
10−7
0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

log (Eν[GeV])
10(c) UPMU

Figure 24: Effective area as a function of neutrino energy for the three ATMPD samples. The different lines
correspond to the different neutrino flavours.

28

Figure 25: Distribution of true neutrino energy for true νµ events, separated between events reconstructed as µ-likeor e-like.

Figure 26: Root mean square of the angle between the true direction of the neutrino and the true direction of themuon, for different values of energy.

(a) FC(b) PC(c) UPMU

Figure 27: Probability density function of cos ∆Φ for different bins in energy.

29

(a) FC(b) PC(c) UPMU

Figure 28: Distribution of true neutrino energy (quantity of interest) vs reconstructed event energy (measured
quantity) for the different samples.
For these purposes, it is needed to smoothen the distribution: this is done by fitting it.
For a given energy, we expect:∆Φ2P (∆Φ; E) = A sin ∆Φ exp −2σ(E)2(23)

where σ(E) is the angular resolution. Therefore:(arccos x)2
P (x = cos ∆Φ; E) = A exp −2σ(E)20(24)

As the plots are done for a range of energy instead of one, the p.d.f. are better fitted with the sum of two
distributions Eq.24 with different sigma values (to cover the difference in angular resolutions over the energy
range).
Energy reconstruction: The true neutrino energy spectrum of the signal is often thought as being of the form:dn= φ0 E −γdE(25)

where γ is the spectral index. While atmospheric neutrinos correspond to γatm. = 3, astrophysical publications areoften referred to γastro. = 2 when giving upper limits on flux or significance, even though there are some hints that
this is not fitting perfectly the current data.
Even if the true neutrino spectrum is fixed, it is still needed to convert to reconstructed energy in order topredict the expected energy distribution of observed events from the source. The Fig.27 shows how true neutrinoenergy relates to the reconstructed energy.
As we have done for the angular resolution, we can project in 1D and plot R = Ereco /Etrue for different trueenergy bins. Some example plots are shown in Fig.28.
To convert from true energy Etrue to reconstructed energy Ereco , we select the correct 1D plot and select one
initialvalue of R using this plot and use Ereco = Rthrown × Etrue.
4.2.4

Outputs

We can separate the output in three categories:• Basic information about SK data-taking: Was the detector running? What is the coverage of the sky if
we consider only upgoing muons (relevant for the UPMU sample)? . . .
• Event information: How many events are observed in SK? what is the expected background? What are
the energy/direction of these events?
• SK-GW information: What is the distribution of the time difference between neutrino detection andtrigger? Where are the events localised with respect to the gravitational wave skymap?
In the following, I want to make the example of S190513bm.
30

(a) FC(b) PC(c) UPMU

Figure 29: Distribution of Ereco /Etrue for different true neutrino energies and different samples.
Detector status The list of runs/subruns that are covering the 1000 s time window around a given gravitational
event is given, with their corresponding livetime and fractional overlapping with the time window (according to themethod presented in Sec.4.2.2), as shown in Tab.2. The main output is then the total livetime of Super-Kamiokandeover the 1000 seconds.
Table 2: Status of Super-Kamiokande on the 1000 seconds time window around S190513bm.
Run

Subrun

81163
81163
81163
81163
81163
81163
81163
81163
81163
81163

334
335
336
337
338
339
340
341
342
343

Start time (JST)Livetime (seconds)

Overlapping (%)

190514/05:44
190514/05:46
190514/05:48
190514/05:50
190514/05:52
190514/05:54
190514/05:56
190514/05:58
190514/06:00
190514/06:02

116.44116.11114.80114.74113.10114.68117.77117.79116.60115.6030.80
100.00100.00100.00100.00100.00100.00100.00100.0027.78

TOTAL (over 1000s)

993.57The Super-Kamiokande status can then be classified into three main categories:• Livetime > 950 s: Super-Kamiokande is taking good quality data over the full time window (sign of normal
operation of the detector).
• 0 < Livetime < 950 s: Super-Kamiokande is taking good quality data only for a fraction of the time window(sign that there may have been a change of run in the middle or some calibration ongoing).
• Livetime = 0 s: Super-Kamiokande was down during the time window (because of an ongoing calibration
period or other issues)
The upgoing muon sample can only be used to detect neutrinos coming from a negative altitude. The coverageskymap is showing which part of the sky is accessible at the time of the alert, e.g. Fig.29. This plot is obtained bysimply computing the altitude a corresponding to each point of the sky at the time of the alert and checking if a ispositive (above horizon, not seen by UPMU, coverage=0 in the plot) or negative (below horizon, seen by UPMU,
coverage=1 in the plot). The boundary between the two regions are intermediary values, because they correspondto sky positions that are below the horizon only for a fraction of the time window ; as the time window is shortwith respect to one day (1000 s over 86 400 s), this effect is relatively small.
Events in Super-Kamiokande All the events passing the selection of Sec.4.2.1 are saved in a list. It is thensaved in a text file with all relevant information.
31

Figure 30: Coverage of the sky with UPMU sample, for the S190513bm event.
Table 3: List of Super-Kamiokande events for S190513bm
Event Time

Delta T

Event Type2019-05-13 20:51:28-183.273

FC

Energy [GeV]

RA [deg]

Dec [deg]

0.678

278.722-37.455

Plots showing their time distribution and the spatial distribution are also obtained, as shown in Figures 31a
and 31b. The spatial distribution is superimposed with the skymap provided by the LIGO/Virgo collaboration,
allowing direct comparison of sky positions.
Based on the reconstructed energy, the angular uncertainty for each SK events can also be provided.

(a) Time distribution(b) Skymap

Figure 31: Plots of Super-Kamiokande events for S190513bm.

Summary of number of observed events and expected background The number of observed events inSuper-Kamiokande can be compared to the expected background (coming mainly from standard neutrino interactions).

4.3

Low Energy analysis

The LOWE framework is also built independently of the rest of the framework and also consists in a watcher script
waiting for new gravitational wave events in the alert database. When it is received, the script is then looking
for the processing of the Super-Kamiokande data of the few minutes around the time of the external trigger, the
LOWE energy range is 3.5 Mev to 100 MeV.
32• data already processed: the LOWE analysis can be launched directly• data not processed: the script is waiting for the data to be ready. More precisely, we wait till the last processedrun (according to the run summary file on sukap) has an end time later than (talert + 1000 s). When thiscondition is fulfilled:– if there is no run covering the ±500 s time window centered on the external trigger, the LOWE
analysis is anyway launched, in order to produce the minimal output– if there are some runs covering the time window, we wait for the end of the processing of makelowe sample

Figure 32: LOWE framework.

4.3.1

Event selection

Assuming that the flux is approximately equally distrubuted among flavors ,as for a core-collapse supernova, thedominant channei in the 3.5 MeV to 100 MeV range is the inverse beta decay of electron antineutrinos ν e + p =e+ + n.The second-most dominant one is neutrino elastic scatteringν + e− = ν + e− ,which is sensitive to all neutrino
flavors but dominated by electron neutrinos.Positrons or electrons from these interactions can produce observablesignals in the SK detector. There are other charged-current and neutral-current interactions with 16 O nuclei thatare subdominant.

Figure 33: LOWE background event rate.
But in our analysis we choose to use 7 MeV to 100 MeVV because we compared with the background of lowe
event(with over 3.5 MeV, 5 MeV ,6 MeV,7 MeV).
By compared with background event rate in lowe cases,we found that event rate which is over 7 MeV is moststable in this analysis.The reason is when in the beginning of SK5 state, the water system is not stable and morenoise is appeared like decay of Radon. It is necessary to cut with 7 MeV that we can get a more stable data.

33

4.3.2

Background estimation and characteristics

As we discussed in lowe events selection ,we cut both background and signal event in 7 MeV and
For the integrated expected number of background event, due to its stability, it can easily estimated by multiplying Super-Kamiokande livetime in the selected time window by the average background rate.
Background rate For each sample, the background rate is estimated using the full dataset of the current SK
period, where background is assumed to be stable. This is done by counting the total number of events passing theselection cuts and dividing by the total livetime.
Background selection In Lowe cases, we try to compared with the background only in one SK normal run
and a week before GW event.The data shows that one week background is more stable than one SK normal run
data.We calculate the average of background event rate is 0.729(events/1000s) by using all of the data from SK5
beginning,and compared average of all data with recent 7 days data.

(a) background for 1 SK run(b) background for 7days

Figure 34: Estimation of the background rate
If the Recent background data is 3σ greater than average background ,we choose to use the Recent backgroundand if not we use the average background to calculate.

(a) Time plot(b) List of selected events

Figure 35: Example of S200225q.

34

Table 4: Fluence limit sample for S200316bj.

4.3.3

GW alert
S200316bj

FL(Fer-Dir)νe
1.1E+10FL(Fer-Dir)ν e
6.8E+10FL(Fer-Dir)νx
9.7E+08

FL(Fer-Dir)ν x
8.0E+10GW alert
S200316bj

FL(Flat)νe
1.2E+10FL(Flat)ν e
7.2E+10FL(Flat)νx
1.1E+09

FL(Flat)ν x (%)
8.4E+10Outputs

Like High Energy work we can separate the output in two categories:• Event information: How many events are observed in SK(Energy from 3.5MeV to 100MeV)? what is the
expected background? What are the energy/direction of these events?
• Calculate information: 90 percent Confidence level and fluence limit for Fermi-Dirac spectrum(Eave =20 MeV), Flat spectrum, and some energies(5,7,10,14,20,30,50,80,100MeV)
Event information We make plot for each event in ±500s in lowe region(Energy from 3.5MeV to 100MeV).And
we select event over 7MeV in N90 calculation and fluence limit calculation.

Figure 36: Fermi-Dirac Spectrum (E = 20MeV).

Calculate information We calculate N90 and fluence limit event by event by using low energy event (E¿7MeV).For
fluence limit, we calculate in 2 types of energy spectrum(Fermi-Dirac,flat) and in mono-energetic(5,7,10,14,20,30,50,80,100MeV).

Figure 37: Fluence limit as a function of energy.

35

5

Statistical analysis

5.1

Introduction

The main goal of the Gravitational Wave follow analysis is to conclude whether or not we have detected a joint
source of gravitational waves and neutrinos. If there is no conclusive observations, we want to be able to put upper
limits on the incoming neutrino flux from the GW source.
There have been numerous publications in Super-Kamiokande and in other collaborations where different methods have been used in order to obtain such results, going from the simple flux limits to some extensive frequentist
or bayesian studies.
5.1.1

Super-Kamiokande approach in previous analyses

In the recent years, the Super-Kamiokande collaboration has published several papers related to the search forneutrinos in coincidence with detected gravitational waves or other transient events: Gamma-Ray Bursts from 1996
to 2000 [18], GRB 080319B [30], GW150914 & GW151226 [8], GW170817 [9], Blazar TXS 0505+056 [21].
In all the papers mentioned above, there have been no significant enough observation to claim for discovery fromthe various tests performed in the analyses (searching for time coincidence only, searching for spatial coincidence
in all data, . . . ). However, all of them use a very similar approach to extract limits on the neutrino flux/fluence.
Fluence limits Assume an observation of Nobs events in a time window Tw centered on the alert time. The
expected number of events from background-only can be simply expressed as: Nb = δt(Tw ) × r, where r is the
background event rate. These two information is converted to an upper limit on the number of signal events Ns
by using Poisson statistics where Nobs is a realisation of the Poisson process with λ = Ns + Nb . There are different
methods to obtain an upper limit NsCL at a given confidence level CL directly from the Poisson probability (we getfor instance Ns90% = 2.3 for Nobs = Nb = 0).
As we have seen in Eq.19, the number of signal events can be related to the flux by integrating over energies withan additional factor Aeff that is taking into account all the effects related to the detector and to neutrino interaction
properties. In all previous Super-Kamiokande publications, the notion of effective area is not used explicitely but
its different components are explicited.
It is worth noting that, depending on the considered sample (FC, PC, UPMU, LOWE), the formula to be used
is not the same because these different samples do not have the same energy ranges and efficiencies.
In the following, we write the flux from a point-source for a given neutrino flavour f :dn[GeV−1 cm−2 ] = Φf [cm−2 ] × λ(E)[GeV−1 ]
(26)
dE
R
where λ(E) is a shape function such that λ(E)dE = 1 and Φf is the fluence (in cm−2 ). For an assumed
neutrino energy spectrum with spectral index γ, we have λ(E) ∝ E −γ and we can also rewrite the flux as:
dn= φ0 × E −γ
dE
where φ0 is the flux constant and can easily be related to the fluence:
Φf =

(27)


φ0  1−γ
1−γEmax − Emin
1−γ

(28)

for γ 6= 1, using simple integration of the spectrum between its boundaries Emin and Emax .
In high energy,we have:x +ν̄x
ΦνFC
=

NT

R

FC
N90
dEν (σ νx (Eν )ενx (Eν ) + σ ν̄x (Eν )εν̄x (Eν )) λ(Eν )(x = e, µ),

PC
N90
,
ν
ν
µ
µ
NT dEν (σ (Eν )ε (Eν ) + σ ν̄µ (Eν )εν̄µ (Eν )) λ(Eν )UPMU
N90νµ +ν̄µ
R
.
ΦUPMU
=
ν
ν
Aeff (θ) dEν (P µ (Eν )S µ (θ, Eν ) + P ν̄µ (Eν )S ν̄µ (θ, Eν )) λ(Eν )ν +ν̄µ

µ
ΦPC

=

R

36

(29)
(30)
(31)

where NT is the number of nucleons (targets) in the fiducial volume, σ is the neutrino-nucleus cross section, εrefers to detection efficiency, P is the probability for a neutrino to produce a muon in the surrounding rock thatwill reach the detector, S covers the neutrino flux attenuation from the Earth and Aeff (θ) is the muon geometricaleffective area (function only of zenith angle).
In this particular example, neutrino and anti-neutrinos of the same flavour are considered together but they canalso be considered separately.
In [21], the integration boundaries are 5.1 GeV to 10 GeV for FC, 1.8 GeV to 100 GeV for PC, 1.6 GeV-10 TeV
for UPMU.In [9],they are 0.1 GeV to 10 GeV for FC-PC, 1.6 GeV-100 PeV for UPMU.
And also in Low energy , we have:ΦLowE =

NT

R

LowE
N90
dEv λ(Ev )σ(Ev )R(Ev , Evis )(Evis )(32)

Where R is the response function to convert electron or positron energy (Ee ) to kinetic energy in SK (Evis ).
The response function and the detection efficiency () are calculated using SK detector Monte Carlo simulation
λ is the density of Eµ assuming an energy spectrum.Here we assume the spectrum is fermi-dirac distribution ,flat and mono-energetic.(See in Fig.35 Fig.36)
We have seen we can rewrite all the previous equations in a simpler form using neutrino effective area:
s
N90
s,f
dEν Aeff (Eν , θ)λ(Eν )Φs,f90 = R(33)

for all samples s (∈ FC, PC, UPMU, . . . ) and all flavours f (∈ νµ , νe , ν̄µ , ν̄e , νµ + ν̄µ , νe + ν̄e ). The effective area
computation was already presented in Sec.4. It is worth noting that effective area depends on the zenith angle θ, sothat fluence limits will depend on the position of the source. Previous analyses have either considered a fixed source
position (because the position was well known from external observations) or scan over the sky (and therefore overdifferent zenith angles).
Some of the previous analyses have also reported mono-energetic fluence limits, i.e. assuming λ(E) = δ(E − E0 )rather than a continuous spectrum.In [18],such limits Φmono (E) are computed and can therefore be used to obtain
fluence limits assuming a given continuous spectrum λ by doing a convolution:
Z
Φλ =

λ(E)
dE
Φmono (E)

−1
(34)

Observation significance In order to quantify if an observation done with Super-Kamiokande is consistent withthe existence of a correlation between detected neutrino flux and external trigger, there have been different methods
in the different publications.
The Table 5 summarises the analyses from these publications.
As can be seen, beyond simple Poisson tests and different methods for fluence limits, there have been no special
statistical treatment to quantify if data indicate the discovery of a GW+ν or not. We will therefore suggest new
methods in the following sections.

5.2

Signal significance:first approach

If we want to quickly say if data are consistent or not with background-only hypothesis without any complicatedstatistical treatment, we can proceed with simple counting.
In Sec.4.2.2, we saw the expected number of events in a 1000 s time window is around nB = 0.111, 0.007 and0.021, respectively for FC, PC and UPMU. In the background-only hypothesis, for a given sample, the number ofobserved events N should be distributed according to a Poisson distribution Poisson(N, nB ).
For a given GW event, we get specific values of N and nB (nB can vary from above values if SK livetime is not100%). We can therefore compute the p-value of the observation for each sample s:
ps =

∞
X

Poisson(k, nsB )

(35)

k=Ns

If we have N = 0, we get p = 100% (observation fully consistent with background). If we have NFC = 1 (resp.
NUPMU = 1), we get p ∼ 10.5% (resp. p ∼ 2.1%).
1 special sample of events with > 1.75 × 106 p.e in the ID

37

Table 5: Summary of the different strategies to evaluate neutrino+external trigger coincidence.
Ref.

[18]

[30]

[8]

Evt.type

SK samples

Time/spatial window

Tests

Spectrum

GRBs 1996 ∼ 2000

LE
HE-e
HE-µ
UPMU

±10 − 100 − 1000s±10 − 100 − 1000s±10 − 100 − 1000s± 1000 s, < 15°

Poisson p-value

E −2 , mono-E

GRB 080319B

LowE
FC
PC
UPMU
UHE1

±1000s-E −2

GW150914
GW151226

LowE
FC+PC
UPMU

±500s

-

flat, mono-E
E −2
E −2

GW170817

LowE
FC+PC
UPMU

±500s
±500s, +14d
±500s, +14d < 5°-flat, Fermi-Dirac
E −2
E −2

Blazar TXS
0506+056

FC,PC
UPMU

1996 ∼ 2018, < 10°±500s, < 5°event rate evolution(dev. from const.rate)

[9]

[21]

E

−2

E −2
, mono-E

We can also do the test with the sum of the contribution of the different channels:p=∞
X

Poisson(k,

X

nsB )

(36)

s

P
k= s Ns

If we have N = 0 (resp. N = 1), we get p = 100% (resp. p ∼ 12.9%).
Also in low energy,we use Poisson-distribution to caclulate the p-value.
These tests are quite limited because they only consider the time correlation of the SK events with the GW.
Moreover, they only probe the agreement of data with background hypothesis and not with the signal+backgroundhypothesis.
However, it can be quite interesting to do the test with the combination of all GW events. In the following
section, we will consider more advanced methods which take into account both time and spatial correlations inorder to quantify signal significance.

5.3

Flux limits:sample-by-sample, position dependent

In a first time, we will reproduce the limits with a similar method with respect to what has been done in previous
Super-Kamiokande publications. The inputs to the analysis will be the same:s• number of observed events in a given sample Nobss
• number of expected background events in the same sample Bexp

We produce a flux limit map following the current procedure:
1. Given an external trigger at time t0 , for each sky position (δ, α), we compute the corresponding zenith angle
θ. We obtain a map as the one presented in Figure 38.
ss
2. For each SK sample, we compute the upper limit on the number of signal events from Nobs
and Bexp
using
the Kraft, Burrows, and Nousek method [23]. We will consider four different samples: FC, PC, UPMU and
FC+PC (where we sum the contribution of the two samples). In the case of each flavour, we will consider
this upper limit to be an upper limit on the number of signal neutrinos of this given flavour, even though
we expect signal neutrinos to be distributed among the different flavours (with mixture 1:1:1 after standard
oscillations).

38

3. For each neutrino flavour, SK sample and for each pixel of the map, we compute limit on E 2 dn/dE, assuming
E −2 spectrum:
s
N90
φs,f
(37)
90 = R
s,fdEν Aeff (Eν , θ)Eν−2
where the integration boundaries used in the analysis are covering the full range E = 0.1 GeVto105 GeV.
4. This allows to obtain maps of flux limits. For UPMU, limits are computed only for zenith angles above π/2i.e. positions below the horizon.
5. We can extract important numbers from these maps: the minimum and maximum values of the flux limits
within the 90% contour. An example is illustrated in Table 6.

Figure 38: Example of zenith angle map for LST = 0 h.

Table 6: Output of the analysis in terms of flux 90% upper limits (in GeV.cm− 2) for the different samples andflavours, for the GW event S190412m. We assume E −2 spectrum.
SampleFlavour

Min 90%GW

Max 90%GW

Marginalised

FC

νµ
νe
ν̄µ
ν̄e
νµ + ν̄µνe + ν̄eall

3.89 · 103
2.16 · 103
9.22 · 103
4.76 · 103
5.47 · 103
2.97 · 103
3.88 · 103

4.03 · 103
2.36 · 103
9.58 · 103
4.99 · 103
5.67 · 103
3.18 · 103
4.06 · 103

3.97 · 103
2.28 · 103
9.45 · 103
4.88 · 103
5.60 · 103
3.10 · 103
3.99 · 103

PC

νµ
νe
ν̄µ
ν̄e
νµ + ν̄µνe + ν̄eall

4.83 · 103
3.35 · 104
7.23 · 103
8.93 · 104
5.82 · 103
5.27 · 104
1.05 · 104

4.93 · 103
3.92 · 104
8.76 · 103
1.27 · 105
6.31 · 103
5.76 · 104
1.14 · 104

4.86 · 103
3.78 · 104
8.44 · 103
1.10 · 105
6.16 · 103
5.59 · 104
1.11 · 104

UPMU

νµ
ν̄µ
νµ + ν̄µ3.16 · 101
3.65 · 101
3.67 · 101

3.73 · 101
4.51 · 101
3.81 · 101

3.35 · 101
4.27 · 101
3.73 · 101

The second computation that can be performed in the same fashion is by assuming mono-energetic emission, as
described in Sec.5.2. The method is the same as above except that λ(E) = δ(E − E0 ) and we scan different valuesof E0 between 0.1 GeV and 100 TeV.

39

(a) FC(b) PC(c) UPMU(d) FC+PC

Figure 39: Energy-dependent fluence 90% upper limits for νµ and ν̄µ for FC, PC, UPMU and FC+PC.

N90
dEν Aeff (Eν , θ)δ(Eν − E)
N90=
Aeff (E, θ)Φmono (E) = R(38)

As this procedure has to be repeated independently for all pixels of the sky map, it is not easy to visualise2 . In
the Figure 39, we show the result of such computation at the best-fit position of S190412m.
Several remarks on the results can be made:• For FC and PC, the result does not depend very much on the sky position: the skymap is more or less uniform
and the results from the Table 6 are similar for average/min/max... This is expected because the efficiencyof the selection and the acceptance of Super-Kamiokande does not depend much on event direction3 .
• For UPMU, as the sample is only covering the sky below the horizon, we get an “incomplete” flux limits
map. It can then happen that some of the numerical applications from Table 6 cannot be done for UPMU, as
2 In [18], it was more straightforward as the external trigger was well localised an therefore only one possible position was considered
3 This was indeed assumed implicitly in previous analyses where the flux formula was not depending on θ at all40

Figure 40: Example of flux likelihood using Equation 39, for the different samples (FC/PC/UPMU) and neutrino
flavours.
e.g. the best-fit of the GW contour is above the horizon. In the worst case, if the full GW 90% is above thehorizon, the UPMU sample can only be used to constrain the sky-averaged (below the horizon) flux limits,
but this quantity is of little relevance as the astrophysical event is most probably located in the “unseen”
region.
• The energy-dependent fluence limits are a simple reflection of the effective area as expected from Equation 38.

5.4

Flux limits:sample-by-sample, marginalised

The approach presented in previous section is ideal if the position of the source is known precisely, so that the limit
at this position can be computed directly. However, the GW sky localisation are generally quite large.
The first strategy is to simply quote a range of limits covering the full sky localisation by using the values
obtained above, while a second strategy consists in marginalising over the localisation to get only one value at theend.
The N90 limit in the previous section is obtained by using the Poisson likelihood with n = Nobs and λ = Nexp .
Similarly, we define the following likelihood based on Poisson statistics and weighted by GW localisation probability:Z(c(Ω)φ0 + nB )N −(c(Ω)φ0 +nB )L(φ0 ; nB , N ) =e
PGW (Ω)dΩ(39)
N!
where c(Ω) is the conversion factor from E 2 dn/dE to number of events (this is the denominator in Equation 37),
nB is the expected background, N is the observed number of events and PGW (Ω) is the GW localisation.
The 90% upper limit is simply obtained by solving:Z φ0,upL(φ0 )dφ0 = 0.90

(40)

0

assuming L is correctly normalised. This effectively corresponds to a Bayesian upper limit, with flat prior on φ0 .
This has to be repeated independently for each SK sample and neutrino flavour.
The Figure 40 shows an example of likelihood as a function of φ0 and the results for S190412m are presented inTable 6, in comparison with the min-max range computed in the previous section.
It is also possible to compute upper limits on νµ + ν̄µ , νe + ν̄e and νµ + ν̄µ +νe + ν̄e (=ALL), assuming equipartitionof the incoming flux between the considered flavours.

41

Conversion from flux to source energy release For each GW event, we have provided upper limits on thedetected neutrino flux. However, the true physical quantity is rather the total energy emitted in neutrinos from thesource. We often use Eiso to refer to this energy, assuming spherical symmetry.
If the position of the source is perfectly known (distance r, sky position Ω), we have:
Eiso
=
4πr2

Z Emax
Emin

dntot
EdE
dE

Z Emax= φ0,tot(41)

E −2 EdE

Emin

Using this formula, it is easy to convert from previously computed flux limits to limits on Eiso . However; thedistance to the source has some uncertainty and there are some correlations between distance and angles. It is thenneeded to use a more sophisticated approach, as it will be described later.

5.5

Combined statistical analysis

We have discussed the observation significant that we compute it by using Poisson distribution.This method is only
depended on observed number, But in high-energy region, we have angular-distribution and energy-distribution.To
get a more accurate data,we use the follow two method to compute the significance and flux limits.
Here we discuss this following parameter:• Angular distribution: The neutrino events should cluster around a given direction that is compatible with
the GW, with a spread that is related to the detector angular resolution (both due to experimental resolution
and to the scattering angle between the neutrino and the outgoing muon).
• Energy distribution: While the background from atmospheric neutrinos have an energy spectrum following
approximately a power law E −2.7 , the signal should have a distinct shape. In the case of simple Fermi
acceleration mechanisms, we can expect a power law E −2 . In the following, we will assume a generic powerlaw E −γ .
We want to use our data D in order to distinguish between two hypotheses:• H0 : data can be completely described by background atmospheric neutrinos• Hs : data consists in a superposition of atmospheric neutrinos and astrophysical neutrinos
To do so, we will use as test statistic the following likelihood ratio:
P (D|Hs )Λ = 2 ln
P (D|H0 )(42)

where P (D|Hs ) and P (D|H0 ) are the likelihood associated to signal+bkg and bkg hypotheses respectively, using
all our present knowledge on angular and energy distributions for both atmospheric background and astrophysicalevents. This definition of Λ is convenient because large values would correspond to the case where data favours Hs
while small values correspond to the case where data favours H0 .
5.5.1

Lambda method

Likelihood definition In a point-source search, assuming a direction ~xS , we can build the following likelihood,as
it has been done in [22, 2]:NLν (nS , γ; x~S ) =e−(nS +nB ) (nS + nB )N Y nS S(x~i , Ei ; x~S , γ) + nB B(x~i , Ei )N!
nS + nBi=1(43)

where nS is the number of signal events (variable of interest), x~S is the test source direction, γ is the spectral
index of the signal (dn/dE ∝ E −γ ) the nB is the expected number of background events, N is the number of observed
events in the selected time window centered on GW trigger, the product is over these events (characterised by theirdirection x~i and their energy Ei ) and S (resp. B) are signal (resp. background) p.d.f. that encodes expectedbehaviour as a function of energy and direction.
42

It is possible to rewrite P (D|Hs ) and P (D|H0 ) using this likelihood and assuming that the signal hypothesis is
the one corresponding to nS = n
cS and γ = γb, where ncS and γ
b are the values maximising the likelihood:
P (D|Hs ) = Lν (cnS , γb; x~S )P (D|H0 ) = Lν (nS = 0; x~S )


Lν (c
nS , γb; x~S )Λν (x~S ) = 2 ln
Lν (nS = 0; x~S )(44)

In the case of Gravitational Wave follow-up analysis, the direction of the test source is not known perfectly, but
we rather have a probability density distribution represented by a skymap wL (x~S ). We can inserted in the previous
formula as an additional weight.
Lν (cnS , γb; x~S )+ 2 ln wL (x~S )Λ(x~S ) = 2 ln
Lν (nS = 0; x~S )
"
#

N
X
n
cS Si
= 2 −cnS ++ 2 ln wL (x~S )ln 1 +
nB Bi
i=1


(45)
(46)

This procedure can be repeated for every pixels of the GW skymap, effectively scanning all the sky in search
for a point-source. We will then use as a test statistic the best fit position in the sky:
T S = max [Λ(x~S )]
x~S

(47)

Signal and background p.d.f. We have used Si = S(x~i , Ei ; x~S , γ) and Bi = B(x~i , Ei ) to quantify how likely agiven event i is compatible with signal and background. It is trivial that Si depends on the characteristics of thesource (nS , γ) while Bi is independent from them.
To ensure proper mathematical and statistical properties, we should have: ZS(E, Ω; x~S , γ) dE dΩ = 1 for all values of x~S , γ
Z
(48)


B(E, Ω) dE dΩ = 1
where we integrate over the full energy-angle phase-space.
We have already characterised in details how they are both distributed with respect to energy and angle andcan therefore define them as it has been done in [13].
For background, assuming angular and energy dependencies are uncorrelated, we can then write:
B(Ei , x~i = [θ, φ]) = Bene (Ei ) × Bang (θi )(49)

For signal, the same decomposition can be performed:
S(Ei , x~i ; x~S , γ) = Sene × Sang

(50)

The angular component is directly obtained from the angular resolution,~ , x~S ); Ei )Sang (x~i ; Ei , x~S ) = P (∆Φ = Angle(xastroi(51)

It is worth noting that the angular component depends on detected neutrino energy and on the test source~ is used to highlight the fact that we should compare the neutrino and the test source
direction. The notation xastroidirection in the same coordinate system4 . As the SK standard ROOT files after reduction already contain the
direction of each events in equatorial coordinates (αi , δi ), we will use the system by convenience5 .
Finally the signal energy component is the result of the correlation between the detector energy resolution (seeEreco
Figure 29) and the assumed signal energy spectrum. More precisely, if we note F (r = E; Etrue ) the distributions
trueof Figure 29 (we have one such distribution for each value of Etrue ), we can then write:Z
Ei; Eν ) × Eν−γ dEν
(52)Sene (Ei ; γ) = C × F (r =Eνwhere C ensures normalisation to unity.
4 Reminder: we can convert from (θ, φ) (or rather altitude-azimuth (a, A)) to (δ, α) (and reverse). If we know the current local
sidereal time θL
5 The choice does not matter as the angular distance between two directions is an invariant in the change of coordinate system.

43

Combination of different samples In the case of High-Energy analyses, we have defined three different andindependent samples: FC, PC and UPMU. Of course, these samples will have different angular resolution, energy
reconstruction, . . . , we may wonder then how to combine them efficiently. The likelihood from Equation 43 caneasily be expressed independently for each sample s:
(s)

(s)

L(s)
ν (nS , γ; x~S ) =

(s)

(s)

(s)

e−(nS +nB ) (nS + nB )N
N (s) !

(s)

(s)
N
Y

nS S (s) (x~i , Ei ; x~S , γ) + nB B (s) (x~i , Ei )i=1

nS + nB

(s)

(s)

(s)

(s)

(53)

where we only consider events in the specified sample s.
(s)
We can then simply define a combined Λ where nS and γ are maximised independently for each sample6 :



d
N
(s) (s)
d
(s) )
Xd
X
n
S
(.
.
.
,
γ
(s)
 + 2 ln wL
Λ(x~S ) = 2 −
nS +
(54)
ln 1 + S i (s) (s)
n
B
i=1
i
B
where the sum is over all the events.
The test statistic T S is defined as before.
Obtaining p-value In an ideal scenario where both background and signal have a very predictive behaviour, we
could simply compute the test statistic for data and compare it to the expected distributions for background andfor signal+background.
However, both the strength (nS ) and the shape (γ) are unconstrained parameters7 . What is fixed, nonetheless,is the background behaviour: we know the event rate, its energy distribution and its angular distribution.
We can then compute a p-value, which quantifies how likely the observation is compatible with the backgroundonly hypothesis. For a given gravitational wave event, this can be obtained in two steps:
1. Compute T S background distribution using randomised neutrinos (random energy, direction and time) andthe fixed GW sky map ⇒ Pbkg (T S)
2. Compute the observed T Sdata and compare to the background distribution:
Z ∞p=
Pbkg (T S) dT S

(55)

T SdataThis procedure is repeated independently for all the GW events. In the case where data are completely compatible with background-only hypothesis, the distribution of the p-values between 0 and 1 should be uniform. If
there is some signal evidence, there should be an excess toward p = 0.
To handle with the combination of the different samples, the distribution Pbkg (T S) is simply obtained using
randomised neutrinos with representative proportions of the different samples.
The zero event toys One particular issue has to be taken care of: the background event rate in SK is relatively
small, such that in any given 1000 s time window, only 0.14 events are expected. This rises a P0 ∼ 87% probability
to observe no event. In such case, it is irrelevant to compute signal significance, as we would always obtain
T S = max 2 ln [wL (x~S )] ≡ κ.
x~S

Nevertheless, these scenarii should in principle affect Pbkg (T S): when doing background toys to estimate
Pbkg (T S), 87% of the toys would have no event and therefore Pbkg (T S) would have a big peak containing 87% of
the probability at T S = κ. Then, when computing the p-value for a given GW follow-up, there are two possibilities:• if there are no SK events in the 1000 s time window, we have T Sdata = κ and p = 100%
• if there are some SK events in the 1000 s time window, we have T Sdata > κ and p < 1 − P0 % (as we are surely
located after the peak containing P0 of Pbkg (T S) integral)
6 In principle, there is no need to maximise γ independently if the spectrum follows an unique power-law E −γ , but doing so allow tohandle the possibility e.g. that FC sample sees a different spectrum than UPMU sample between they have different energy ranges.
7 Even if we had a perfect modelling of all astrophysical sources, they would depend on the distance to the source and the type ofthe source.

44This leaves a big unfilled gap of values between 1 − P0 and 1, that has no physical meaning.
However, we can get rid of this situation by simply changing point of view: as stated before, it only makes sense
to compute the p-value in the case where we have detected some events in Super-Kamiokande. Then, what wereally need is the probability for the observation to be compatible with background given that we have at least
one event in time coincidence. Effectively, it is equivalent to the suppression of zero event toys in Pbkg (T S)
distribution.
Then, the p-value will be distributed between 0 and 1 continuously, as expected and if an observation is fullyinconsistent (resp. consistent) with the signal direction, we will get p close to 1 (resp. 0).
The additional interesting feature is that, if one wants to obtain the p-value taking into account the zerobackground toys (i.e. the fact that observing at least one event in time coincidence is already quite “rare”), one
can simply multiply the above p-value by 1 − P0 (for a simple order of magnitude, let’s say divide by a factor 10),that would effectively reintroduce the zero event toys in the game.
Example of application Before presenting the general implementation of the analysis, let’s discuss one specificexample: S190513bm.
Results for this particular event have already been presented in section 4.2.4. The most important piece of
information is presented in Figure 31: there was one SK event (in FC sample), with a reconstructed energy of0.68 GeV and with a reconstructed direction close to the GW contour.
For every pixel of the sky, we compute Λ using Equation 54. We obtain Figure 41.

Figure 41: Lambda map for S190513bm.
The final value of the test statistic correspond to the maximum value of above map: T Sdata ∼ −2.1415.
Using the rates from Table 1, we throw randomly the number of observed events in FC/PC/UPMU (with theconstraint that NFC + NPC + NUPMU > 0). For each event, its energy and direction are chosen randomly from theexpected background distribution obtained from SK4 MC8 .
After throwing ∼ 50k toys, we can plot the distribution presented in Figure 42. From this, we can easily extract
the p-value by comparing it with T Sdata .
We obtain p = 4.8%. This should only be interpreted in frequentist fashion: if there is truly no signal and if wewere able to reproduce the same conditions (same detector status, same GW localisation) N times, 4.8% of such
hypothetical experiments would have a test statistic at least as high as the one from real data.
Final results For reference, the Table 7 shows the results for all gravitational waves from April 2019 to March
2020 in real-time, with a SK event detected in the 1000 s time window.
We see some minor disagreements between the two methods, which are simply due to statistical fluctuations
and precision loss in the operations.
5.5.2

Chi-square method

As a cross-check, it can be interesting to consider a second approach, with a different method to obtain a p-value.
We will follow the method presented in [10].
8 This could be performed using data directly.

45

Figure 42: Background TS distribution for S190513bm.

Table 7: Results of the likelihood statistical analysis for all gravitational waves from April 2019 to March 2020 witha SK event detected in the 1000 s time window. We remind the date/time of the event, the number of detected
events in the different samples. Two set of p-values are provided: pold corresponds to the one without using thepre-computation trick (where we take the time to do the full computation from scratch); ppre. is obtained using thepre-computation trick.
TriggerDate/time

FC
Nobs

PC
Nobs

UPMU
Nobs

pold (%)

ppre. (%)

S190426c
S190513bm
S190602aq
S190728q
S190814bv
S190924h
S191213g
S200316bj

2019-04-26 15:21:55
2019-05-13 20:54:28
2019-06-02 17:59:27
2019-07-28 06:45:10
2019-08-14 21:10:39
2019-09-24 02:18:46
2019-12-13 04:34:08
2020-03-16 21:57:56

0
1
1
1
1
1
1
0

0
0
0
0
0
0
0
0

1
0
0
0
0
0
0
1

100.0
5.2
2.5
17.0
100.0
53.2
16.0
57.4

100.0
4.9
2.6
17.9
100.0
51.5
15.4
56.1

46

First, we define the likelihood:
Z
L=

L(x~S )dx~S

with L(x~S ) =(56)
j~j , x~S ; Ei )j Sν (xQ j
BGW (x~S ) j Bν (x~j , x~S )SGW (x~S )

Q

(57)

where the products is over all SK events in the 1000 s time window and the SX / BX are normalised to unity(when integrating over the full sky) and are respectively characterising signal and background angular distribution
for the component X.
More precisely, Sνj is defined in the same way as in Equation 51 and Bνj is defined as being the angular component
from Equation 49. SGW (x~s ) is the GW probability distribution from LIGO-Virgo and BGW (x~s ) = 1/(4π).
From this likelihood, we can directly define a p-value that is covering only the possible angular correlations
between neutrinos and the gravitational wave:
Z ∞Pbkg (L0 ) dL0(58)
psky =
Ldata

where Pbkg (L0 ) is the distribution of likelihood obtained with a sample of O(100)k background toys. This is obtainedin the same way as for the Λ method, with the exact same consideration of removing the toys with zero event in
the time window.
We then define several p-values that will cover different points:• pGW = 1 − Poisson(0, T · FAR) is the probability for such a GW event to raise from background. FAR is thefalse alarm rate of the given event and T is the considered coincidence time window.
R∞(j)• pν = Ej Pbkg (E 0 ) dE 0 is the p-value obtained by comparing reconstructed neutrino event energy with thebackground distribution.
PN −2• pcluster = 1 − i=0 Poisson(i, nB ), where nB is the expected number of background events. This will favour
the case where more than 1 event are detected in the time window (pcluster = 1 if N = 1).
The different p-values are then combined in a joint χ2 :
χ2data = −2 ln psky × pGW ×


Y


p(j)
ν × pcluster

(59)

j

From this new test statistic, we can compute the final p-value:
Z ∞p=Pbkg (χ2 ) dχ2(60)

χ2dataIn [10], the background distribution Pbkg (χ2 ) is obtained using time-scrambled neutrino data and time-scrambled
GW data. As the goal here is to do a cross-check with the previously described Λ method, with a given GW event,in the following, we will consider the GW event fixed, even for the background estimation. Therefore, FAR will
always have the same issue, so that we can remove pGW from the Equation 59.
Combination of different samples As for the previous method, we may wonder how to combine FC, PC and
UPMU efficiently. Here, it is quite straightforward that we can consider all the events from the different samples
P (s)(with their different p.d.f.) together in Equation 57 and nB = s nB for pcluster .
Example of application Let’s consider again S190513bm. The Figure 43 shows the background distribution
and the value for data of the likelihood from Equation 57, while the Figure 44 shows the background distribution
and the value for data of χ2 .
We have Ldata = 2.17 and χ2data = 62.66. From the background distributions, we directly obtain the final
p-value: 12.9%

47

Figure 43: Background likelihood distribution compared to the value obtained for data for S190513bm.

Figure 44: Background χ2 distribution compared to the value obtained for data for S190513bm.

48

5.6

Flux limits: combined

5.6.1

Signal simulation

In order to compute flux limits, we will need to be able to simulate neutrino signal events from a given source
direction. Fortunately, we already have most of the needed ingredients.
We have defined how signal and background are distributed with energy and angles. We can use this to do
simulations.
For a given sample and neutrino flavour, we can generate a signal event by applying the following procedure:
1. generate Etrue following Aeff (E, θ) × λ(E)Erecofollowing F (r; Eν ) (already used in Equation 52)
2. obtain Ereco by generating r = E
true

3. generate the angular distance ∆Φ between true neutrino direction (α, δ) and reconstructed direction, following
the angular resolution p.d.f. P (∆Φ; Ei ) (already used in Equation 51)
4. generate the actual reconstructed direction randomly on the circle of center (α, δ) and of radius ∆ΦFor the combined analysis (combining the samples), every event is assigned randomly to a SK sample based on
the relative weights:
Z
Aeff (s, f )(E, θ) × λ(E)dE for s ∈ {F C, P C, U P M U }(61)

In the case of combining flavours, e.g. νµ + ν̄µ , the flavours is also assigned randomly to the event by using
similar weights.
5.6.2

Marginalised limits

A second approach using the test statistic (Λ method) has been implemented. It follows the procedure describedin [31]. Instead of computing flux limits for each direction of the sky, the direction is marginalised by defining theflux likelihood:
Z X
2 h
k
i
(c(Ω)φ0 )exp (−c(Ω)φ0 ) × Pk (T Sdata ) × PGW (Ω)dΩ(62)
L(φ0 ; T Sdata , PGW ) =k
k=0

where c(Ω) is the conversion factor from φ0 to the expected number of signal events in SK, Pk (T S) is the test
statistic distribution for k injected signal events, simulated using the strategy described in section 5.6.1), T Sdata is
the measured test statistic and PGW (Ω) is the GW skymap provided by LIGO-Virgo.
The Figure 45 shows the distributions Pi (T S) (i = 0, 1, 2) for the simulation of νµ events. From these normalised
distribution, we can extract Pi (T Sdata ).
One example of L(φ0 ; T Sdata , PGW ) for S190408an is shown in Figure 46.
φ0 = E 2 dn/dE limits are simply obtained by solving:Z ∞Z [φ0 ]90
L(φ0 ) × dφ0 = 0.90
0

L(φ0 ) × dφ0

(63)

0

The example of S190408an and S190412m is presented in Table 8.
dn
Table 8: Computation of the limits on φ0 = E 2 dE
for S190408an and S190412m [in GeV.cm−2 ], obtained using the
combined method compared with the sample-by-sample results.

Trigger

Flavour

FC

PC
3

3

UPMU

Combined

5

S190408an

νµν̄µ4.00 · 10
9.19 · 103

5.16 · 10
8.94 · 103

2.47 · 10
7.48 · 105

2.27 · 103
4.53 · 103

S190412m

νµν̄µ3.97 · 103
9.45 · 103

4.86 · 103
8.44 · 103

3.35 · 101
4.27 · 101

3.38 · 101
4.31 · 101

We can see that the marginalised results is consistent with the min-max range, but it has the advantage ofconsisting in one value instead of two
49

Figure 45: Test statistic distribution for nS = 0 (only background), nS = 1 (background + 1 νµ signal events),
nS = 2 (background + 2 νµ signal events).

Figure 46: Flux likelihood after marginalising over the sky, for S190408an.

50

5.7

Isotropic total neutrino energy limits

5.7.1

High Energy sample

Individual upper limits The physical quantity of interest is the total energy emitted in neutrinos by the source(correcting by distance), assuming isotropy. We can use the same method than for flux limits presented in section 5.6
but using the 3D skymap (r in Mpc) VGW (r, Ω) as in [31]:
#
"
Z X
2
k
(c0 (r, Ω)Eiso ) −c0 (r,Ω)Eiso
× Pk (T Sm ) VGW (r, Ω)r2 drdΩ
(64)
L(Eiso ) =
e
k!
k=0

0

where c (r, Ω) is the factor to convert Eiso to number of signal events (combining all the samples) for known source
position Ω and distance r.
upUpper limit Eiso
is obtained by using:up
Z EisoL(Eiso )dEiso = 0.90
(65)
0assuming the likelihood is correctly normalised.
For a given neutrino flavour, we can then get an upper limit on the energy in neutrinos that would be detectedν
as this particular flavour at Earth. For instance, using νµ , this allows to put a limit on Eisoµ .
In order to extrapolate to the total energy emitted in neutrinos by the source, we can assume equipartition on
Earth, i.e. νe : νµ : ντ : ν̄e : ν̄µ : ν̄τ = 1 : 1 : 1 : 1 : 1 : 1. This assumption is reasonable is the main production
channel for neutrinos is pion decay in the astrophysical source vicinity.
νall
In the case of the example above, we can extrapolate as Eiso
= 6 × Eisoµ .
If the likelihood of Equation 64 is computed for νµ + ν̄µ , one can constrain jointly the two flavours and can
ν +ν̄
all
extrapolate to total energy: Eiso
= 3 × Eisoµ µ . The situation is similar for νe + ν̄e . Finally, if all flavour areν +ν̄ +ν +ν̄allconsidered together, i.e. νµ + ν̄µ + νe + ν̄e , we get Eiso
= 3/2 × Eisoµ µ e e .
The Table 9 summarises the results for GW190521. As the flux limits for νe and ν̄e are dominated by FC and areν +ν̄ ,up
all,upallmuch less stringent than νµ and ν̄µ limits, the final Eiso
are dominated by muon limits and Eiso
∼ 3 × Eisoµ µ .
Table 9: Example of summary of flux and Eiso limits, for the different samples and flavours.
Name

SampleHE E 2dn
dE

GW190521
HE Eiso5.7.2

νµ (νx )ν̄µ (ν̄x )

νe

ν̄e

FC
PC
UPMU
Combined

3.96 · 103
4.89 · 103
4.24 · 101
3.73 · 101

9.60 · 103
8.35 · 103
4.96 · 101
4.93 · 101

2.28 · 103
3.99 · 104
2.21 · 103

4.71 · 103
1.02 · 105
4.60 · 103

Per-flavourν + ν̄
All

2.61 · 1057 3.70 · 1057 1.70 · 1059 3.43 · 1059
3.22 · 1057
2.26 · 1059
57
9.30 · 10

Low Energy sample

As High-Energy sample, the Eiso can be defined as
Eiso /(4πd2 ) =Z EmaxZ EmaxΦlowe × λ(E) × EdE = Φlowe

Emin

λ(E) × EdE

(66)

Emin

In Low-energy cases , we assume the energy spectrum as Fermi-Dirac spectrum.That we can get λ(E) = λ(F D).
And compute Eiso like:
Eiso /(4πd2 ) =Z 100

Φtotallowe × λ(F D) × EdE

(67)

0

cE 2

We know that the Fermi-Dirac distribution iseE
Eave (20)/3.15)and we can get :51(here we assume the average energy is 20 M eV+1

Eiso = (4πd2 ) ×

Z 100

Φtotallowe ×

0

CE 2

× EdE

E
Eave (20)
3.15

(68)

e
+1
where C is Integral constant. As we discussed in High energy , the Fluence-limit of all flavor can be 6 times one
flavor(In low energy cases , the most of neutrino is anti-election neutrino,So here we use 6 × Φ(ν e ))Limit with likelihood Here we compute the isotropic with a simple method.Like the High-Energy, here we must
consider the distance uncertainty.So we must consider it with a likelihood. Here we use Poisson likelihood that:
N

L(Nsig ; Nobs , Nbkg ) =(Nbkg + Nsig )obs −(Nbkg +Nsig )eNobs !

(69)

And we can get the likelihood with limit is:N(Nbkg + CΦ)obs −(Nbkg +CΦ)e
Nobs !
R
where C is the conversion from fluence to Nsig : C = Aeff (E)λ(E)dE
We can get Likelihood with Eiso for known distance d is :L(Φ; Nobs , Nbkg ) =

(70)

N
Nbkg + C 0 /d2 × Eiso obs −(Nbkg +C 0 /d2 ×Eiso )eL(Eiso ; Nobs , Nbkg ) =Nobs !
00

where C is the conversion from Eiso at 1 Mpc to Nsig : C =

4πd2

R

(71)

R
Aeff (E)λ(E)dE
1R×C =Eλ(E)dE
4π(1M pc)2 Eλ(E)dE

And likelihood with distance error can be :Z ∞L(Eiso ; Nobs , Nbkg ) =0

Nbkg + C 0 /r2 × Eiso
Nobs !

Nobs

0

2

e−(Nbkg +C /r ×Eiso ) × P (r)dr

(72)



(r − dmean )2
where P (r) ∝ exp −for r ≥ 0
2σd2
RP
So we can get the final Function: F (Eiso , r)×P (r)dr = i F (Eiso , ri ), where ri are randomly thrown following
the gaussian distribution.
The Table 10 summarises the results for GW 190408 181802,GW 190424 180648.
Table 10: Example of Low energy summary Fluence ,Eiso limits, for the different samples and flavours.
Name

GW 190408 181802

GW 190412 180648

Obs

Type

Fluence limit(Fermi-Dirac)

Eiso

Eiso (DistanceError)

3

νe
νe
νx
νx
Combine

6.9 · 109
1.7 · 1084.3 · 1010
5.0 · 1010
9.7 · 1086.7 · 1061
1.6 · 1060
4.2 · 1062
4.9 · 1062
9.5 · 1062

9.5 · 1061
2.3 · 1060
5.9 · 1062
7.0 · 1062
1.3 · 1061

1

νe
νe
νx
νx
Combine

3.9 · 1099.6 · 107
2.4 · 1010
2.9 · 1085.5 · 1081.1 · 1062
2.6 · 1060
6.5 · 1062
7.7 · 1062
1.5 · 1061

2.0 · 1062
4.7 · 1060
1.3 · 1063
1.4 · 1063
2.9 · 106152

6Result merger and result database

6.1

Inventory of different results

Let’s detail all the outputs from the different part of the analysis. In each case, the basic information are stored ina csv file, in the form of a dictionary. Some additional information are stored as plot in png format or as list in
dat files.
6.1.1

Alert receiver

The basic information from the GW event is obtained from the GCN notice file (xml format):• ALERT FILE: path to the alert file
• ALERT NATURE: nature of the alert (“GW”)• NOTICE TYPE: type of the GCN notice (“Preliminary”, “Initial”, “Update”)• NOTICE UTC: date and time of the GCN notice (datetime object)• TRIGGER NAME: name of the trigger (e.g. S190513bm)• TRIGGER UTC: date and time of the trigger (datetime object)• TRIGGER DATE: date of the trigger (string object)• TRIGGER TIME: time of the trigger (string object)• TRIGGER JD: Julian date of the trigger• GW DISTANCE: best-fit distance to the Gravitational Wave (in Mpc)
6.1.2

High Energy analysis

The basic information are coming from atmpd results.csv file (one created for each GW follow-up). For simplicity,
there is some redundancy and complementarity with the data from Sec.6.1.1.
• ALERT FILE: path to the alert file
• GW TYPE: type of GW event (BBH, NSBH, BNS...)• GW SKYMAP: path to the skymap file
• GW SKYAREA: area of the 90% contour (in %)• TIMEWINDOW: size of the used time window (in seconds)• SK LIVETIME: livetime of the SK detector (in seconds)• SK RUNS: list of runs and subruns covering the time window (with start time + livetime + overlap with thetime window, see Sec.4.2.2• SK SAMPLE OBS9 : observed number of events in the sample• SK SAMPLE EXP: expected number of background events in the sample• SK UPMU SKYCOV: which fraction of the GW probability is covered by upgoing sky (altitude < 0)
• SK UPMU SKYCOV 90: which fraction of the 90% GW contour is covered by upgoing sky (altitude < 0)
• SK PVALUE: what is the p-value obtained by simple Poisson counting as described in Sec.5.2
9 In the following, SAMPLE will already refer to either FC or PC or UPMU or FCPC, variables with this notation are defined forthe 4 samples.

53

All the flux limits are also provided in terms of float values, for the different samples (FC, PC, UPMU, combined)
and different neutrino flavours (νµ , ν̄µ , νe , ν̄e and combining flavours νµ + ν̄µ , νe + ν̄e , νµ + ν̄µ + νe + ν̄e ).
Additional data are provided in the form of additional files:• atmpd events.csv: list of events in the time window (with their respective sample, time, energy, position...
• atmpd time.png: plot of the time distribution of events in the time window• atmpd skymap.png: skymap distribution of events in the time window, on top of the GW probability• skymap SK coverage.png: skympa showing the part of the sky covered by UPMU sample (altitude < 0)• event displays: one for each selected SK events• flux limit plots: one for each sample and neutrino flavour
Finally all the outputs from statistical analysis (only relevant if at least one event is observed) are provided
through a set of variables and plots:• SK PVALUE Lambda: p-value from Lambda method• SK PVALUE Chi2: p-value from chi-square method• Lambda lambdaData.png and Lambda lambdaBkg.png: plots like in and• Chi2 chi2Bkg.png: plot like in
6.1.3

Low Energy analysis

The basic information are coming from one csv file (one created for each GW follow-up). When imported, thenames are modified in order to match above naming conventions• NUMBER OF OBSERVED EVENTS → SK LOWE OBS: number of observed events• BACKGROUND RECENT and BACKGROUND CONTS: expected background from LOWE (respectively estimated using
recent data and using average value)• BACKGROUND USED → SK LOWE EXP: expected number of background events, as explained in lowe section• N90 → SK LOWE UP: 90% upper limit on the number of signal events
In addition:• event *.dat: list of events with time and energy• plot *.pdf: plot of ∆t vs energy for selected events

6.2

Output processing

When a new gravitational wave alert is received, the result merger code is triggered, waiting for High energy and
Low energy results, as presented in Figure 47.
There are two possible outcomes:• while we are waiting for High Energy/LOWE results, the GW event is retracted by LIGO-Virgo: we can giveup waiting for High Energy/LOWE (in parallel, both analyses are stopped as well)• if there is no retraction, where High Energy and LOWE results are available, we process them
The diagram in Figure 48 is showing the general flow of the result merger script, with its main tasks: collect
the results from the different branches (GCN receiver, High Energy, LOWE) ; merge them in the results database; broadcast the results in different formats (more explanations in the next section).

54

GCN notice

GW alert

SK datasocket connection

GCN receiverappend

ATMPD(HighEnergy) frameworkatchResults databaseerge

LOWE frameworkmm

Alert database

watch

automatic email + web

Expert check
manual email

SK GCN circular
Figure 47: Overview of the general framework

GW Notice
New or Update

GW Notice
Retraction

Remove fromresults database

if retractionin the meantime

Waiting forHigh Energy/LOWE processing

GCN
HighEnergy

Getting
results

LOWE
Writting toresults database

Update webpage
Write summary PDF

Send internal email
Figure 48: Flow chart of the result merger process

55

6.3

Output formats

6.3.1

Results database

The different variables described in section 6.1 are stored in the results database, created and updated using pandas.
The csv format is used for storing the database. The plots and additional data files are stored in a dedicated folder(one for each GW trigger).
6.3.2

Web page

The panda database is also exported as HTML format.
The Figure 49 shows a screenshot of the page.

Figure 49: Screenshot of the web page on SK intranet summarising the results of GW follow-ups.

6.3.3

PDF

Using PyLaTeX, a PDF is created for each GW trigger, that is containing all relevant information and plots.

56

3Super-Kamiokande events in the 1000 seconds time window
Table 1: Number of observed events (compared to background expectation)

Follow-up of S190513bm
The Super-Kamiokande collaboration
April 30, 2020

FC

PC

UPMU

LOWE

Observed
Expected

1
0.110

0
0.007

0
0.021

0
0.729

Summary of the alert

100

Energy (MeV)

1

Sample• Status of the notice: Initial• Date of the notice: 2019-05-13 21:22:17• Unique name of the trigger: S190513bm• Date of the trigger: 2019-05-13 20:54:28.747089

90
80
70
60
50

• JD of the trigger: 2458617.37116

4030• Type of event: BBH (94.3%)

20• Distance: 1987.0 +/- 501.2 Mpc

10
0

2

Status of Super-Kamiokande

05/14
05:30

05/14
05:45

05/14
06:00

05/14
06:15

Figure 2: Time distribution of ATMPD (left) and LOWE (right) events in the time window

During the 1000 seconds time window:
Run

Subrun

81163
81163
81163
81163
81163
81163
81163
81163
81163
81163

334
335
336
337
338
339
340
341
342
343

Start time (JST)Livetime (seconds)

Overlapping (%)

190514/05:44
190514/05:46
190514/05:48
190514/05:50
190514/05:52
190514/05:54
190514/05:56
190514/05:58
190514/06:00
190514/06:02

116.44116.11114.80114.74113.10114.68117.77117.79116.60115.6030.16
100.00100.00100.00100.00100.00100.00100.00100.0028.43

TOTAL (over time window)

993.57 seconds

Figure 3: GW skymap and distribution of ATMPD events

Figure 1: Sky coverage of Super-Kamiokande UPMU
Figure 4: SK ATMPD events display (in the order of appearance in the above table)• Fraction of the GW skymap probability covered by Super-Kamiokande UPMU: 0.00%

No observed LOWE events• Fraction of the GW 90% area covered by Super-Kamiokande UPMU: 0.00%

1

2

Table 2: List of Super-Kamiokande ATMPD events
Time

Delta T

2019-05-13 20:51:28

4-183.272552

Type

RA [deg]

Dec [deg]

Energy [GeV]

Ang.Unc [deg]FC

278.917694-36.721024

0.707847

28.276404

Fluence limits
Table 3: Constraints on neutrino fluence using Super-Kamiokande samples

SampleObs.

min-max 90%UL in GW contourExp.

Φ(νe ) [cm−2 ]Φ(νµ ) [cm−2 ]Φ(νe ) [cm−2 ]

FC
PC
UPMU
LOWE

1
0
0
0

0.110
0.007
0.021
0.729

3.79
2.30
2.30
-

8.79 · 104 − 8.92 · 104
1.07 · 105 − 1.14 · 105-

6.46 · 104 − 6.62 · 104
9.06 · 105 − 1.50 · 106-

2.49 · 105 − 2.52 · 105
1.83 · 105 − 2.02 · 105-

1.69 · 105 − 1.76 · 105
1.63 · 106 − 3.24 · 106-

Λ (> T Sbkg,50% )Λ (> T Sbkg,90% )Λ (> T Sbkg,99% )

-

-

-

1.93 · 104 − 2.98 · 104
2.25 · 104 − 3.69 · 104
2.80 · 104 − 4.79 · 104

-

-

-

(b) FC+PC νe(c) UPMU νµ(a) FC+PC νµS

Φ(νµ ) [cm−2 ](a) FC νµ(b) PC νµ(c) UPMU νµFigure 6: Energy-dependent fluence limits [in cm−2 ](d) UPMU ν̄µFigure 5: Fluence limits maps [in cm−2 ]5

ATMPD statistical analysis

Figure 7: Fluence limits for νµ using the Lambda method (T Sthr > T Sbkg,99% ) [in cm−2 ]The probability to observe data at least as extreme as the real data, given the expected background, using simple Poisson statistics
is: 12.92%
Table 4: p-values from statistical analyses
MethodLambdaChi-squarep-value

5.67%

12.78%(a) Lambda method(b) Chi-square method

Figure 8: TS background distribution compared to the data value for the two methods

7

Follow-up of O3 real-time events

From April 2019 to March 2020, the LIGO-Virgo collaboration has sent 80 public real-time alerts, 24 of which havebeen retracted, leaving 56 GW triggers that are interesting for follow-up observations.

7.1

All GW triggers

The Table 11 summarises the follow-ups of O3 events and the Table 13 shows the full description of all followups performed with GW triggers from O3 observation run, including the ones where SK was down. For more
information, you can go to the PDF that is created for each trigger, as described in section 6.3.3, on the web page.
Name

Value

Total number of GW triggers
Number of non-retracted GW triggers
Number of triggers with SK livetime
Number of triggers with SK events

80
56
46
8

Table 11: Summary of the number of GW triggers in O3 observation run.
We can sum the number of observed and expected events for each sample to search for any deviation frombackground expectation.
SampleObserved

Expected

p-value (%)

FC
PC
UPMU
LOWE

6
0
2
38

5.06
0.33
0.74
31.98

39.48
100.00
16.98
25.67

Table 12: Total number of observed and expected events for all O3 observation run, along with the associatedPoisson p-value.

58

Table 13: Summary of all GW triggers for the O3 observation run. The different columns are: Trigger name ; Alert time in UTC ; Most probable GW
event type (with its probability, classifcation based on the masses) ; area of the 90% GW contour ; best-fit distance to the astrophysical object ; SK detector
livetime ; number of observed events and number of expected background events for the four samples.

59

Name

Alert UTC

Type

2
A90%
GW [deg ]d [Mpc]

L [s]S190408an
S190412m
S190421ar
S190425z
S190426c
S190503bf
S190510g
S190512at
S190513bm
S190517h
S190519bj
S190521g
S190521r
S190602aq
S190630ag
S190701ah
S190706ai
S190707q
S190718y
S190720a
S190727h
S190728q
S190814bv
S190828j
S190828l
S190901ap
S190910d
S190910h
S190915ak
S190923y
S190924h
S190930s
S190930t
S191105e

2019-04-08 18:18:02
2019-04-12 05:30:44
2019-04-21 21:38:56
2019-04-25 08:18:05
2019-04-26 15:21:55
2019-05-03 18:54:04
2019-05-10 02:59:39
2019-05-12 18:07:14
2019-05-13 20:54:28
2019-05-17 05:51:01
2019-05-19 15:35:44
2019-05-21 03:02:29
2019-05-21 07:43:59
2019-06-02 17:59:27
2019-06-30 18:52:05
2019-07-01 20:33:06
2019-07-06 22:26:41
2019-07-07 09:33:26
2019-07-18 14:35:12
2019-07-20 00:08:36
2019-07-27 06:03:33
2019-07-28 06:45:10
2019-08-14 21:10:39
2019-08-28 06:34:05
2019-08-28 06:55:09
2019-09-01 23:31:01
2019-09-10 01:26:19
2019-09-10 08:29:58
2019-09-15 23:57:02
2019-09-23 12:55:59
2019-09-24 02:18:46
2019-09-30 13:35:41
2019-09-30 14:34:07
2019-11-05 14:35:21

BBH (100.0%)
BBH (100.0%)
BBH (96.7%)
BNS (99.9%)
BNS (49.3%)
BBH (96.3%)
Terrestrial (58.0%)
BBH (99.0%)
BBH (94.3%)
BBH (98.3%)
BBH (95.6%)
BBH (96.6%)
BBH (99.9%)
BBH (99.0%)
BBH (94.3%)
BBH (93.4%)
BBH (99.0%)
BBH (100.0%)
Terrestrial (97.9%)
BBH (98.9%)
BBH (92.2%)
BBH (95.4%)
NSBH (99.8%)
BBH (100.0%)
BBH (100.0%)
BNS (86.1%)
NSBH (97.6%)
BNS (61.2%)
BBH (99.5%)
NSBH (67.8%)
MassGap (100.0%)
MassGap (95.1%)
NSBH (74.3%)
BBH (95.3%)

386
157
1443
7461
1131
448
1166
252
691
939
967
765
488
1172
1483
49
825
921
7246
443
152
105
24
228
358
14753
2482
24265
318
2107
303
1747
24221
643

1472.9812.2
1628.4
156.1
376.7
421.1
227.2
1387.7
1987.02950.03153.53931.41136.1
797.3
925.7
1848.95262.9
781.4
226.6
868.5
2838.6
873.8
267.4
1946.31528.3
240.9
631.5
229.5
1583.5
438.1
547.9
708.9
108.5
1182.8

993.1
993.1
993.1
993.4
992.6
993.2
0.0
993.7
993.6
0.0
993.7
994.0
994.0
993.6
992.4
992.5
992.4
992.4
993.0
992.8
992.5
992.5
993.9
542.5
0.0
993.9
994.1
994.0
0.0
994.2
994.1
994.1
994.3
991.6

Nobs

FC
Nexp

Nobs

PC
Nexp

0
0
0
0
0
0
0
1
0
0
0
1
0
0
0
0
0
0
0
1
1
0
0
0
0
0
1
0
0
0

0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.061
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.004
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007

UPMU
Nobs Nexp
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.009
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016

LOWE
Nobs Nexp
3
1
3
1
0
0
1
0
1
3
0
0
2
0
2
0
1
0
0
2
0
1
1
1
0
1
0
0
2
0

0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.955
0.729
0.962
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729

Table 13: Summary of all GW triggers for the O3 observation run. The different columns are: Trigger name ; Alert time in UTC ; Most probable GW
event type (with its probability, classifcation based on the masses) ; area of the 90% GW contour ; best-fit distance to the astrophysical object ; SK detector
livetime ; number of observed events and number of expected background events for the four samples.

60

Name

Alert UTC

Type

2
A90%
GW [deg ]d [Mpc]L [s]S191109d
S191129u
S191204r
S191205ah
S191213g
S191215w
S191216ap
S191222n
S200105ae
S200112r
S200114f
S200115j
S200128d
S200129m
S200208q
S200213t
S200219ac
S200224ca
S200225q
S200302c
S200311bg
S200316bj

2019-11-09 01:07:17
2019-11-29 13:40:29
2019-12-04 17:15:26
2019-12-05 21:52:08
2019-12-13 04:34:08
2019-12-15 22:30:52
2019-12-16 21:33:38
2019-12-22 03:35:37
2020-01-05 16:24:26
2020-01-12 15:58:38
2020-01-14 02:08:18
2020-01-15 04:23:09
2020-01-28 02:20:11
2020-01-29 06:54:58
2020-02-08 13:01:17
2020-02-13 04:10:40
2020-02-19 09:44:15
2020-02-24 22:22:34
2020-02-25 06:04:21
2020-03-02 01:58:11
2020-03-11 11:58:53
2020-03-16 21:57:56

BBH (100.0%)
BBH (100.0%)
BBH (100.0%)
NSBH (93.2%)
BNS (76.8%)
BBH (99.7%)
BBH (99.1%)
BBH (100.0%)
Terrestrial (97.3%)
BBH (100.0%)
Unmodeled
MassGap (100.0%)
BBH (96.9%)
BBH (100.0%)
BBH (99.3%)
BNS (62.9%)
BBH (96.4%)
BBH (100.0%)
BBH (95.8%)
BBH (89.0%)
BBH (100.0%)
MassGap (99.6%)

1487
852
104
6378
4480
361
253
1850
7373
6199
403
765
2520
54
27
2325
781
73
22
5656
35
507

1810.1
742.4
677.5
385.0
200.9
1769.6375.92518.1282.8
1136.3340.2
4031.1908.4
2142.0200.93533.11575.0994.9
1820.11114.61178.0

0.0
992.1
992.1
991.9
978.6
992.0
992.0
992.0
991.8
992.1
0.0
991.9
991.6
991.6
991.9
0.0
0.0
0.0
991.4
0.0
992.5
992.6

Nobs

FC
Nexp

Nobs

PC
Nexp

0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0

0.111
0.111
0.111
0.110
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007

UPMU
Nobs Nexp
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1

0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016

LOWE
Nobs Nexp
2
0
2
2
2
0
0
0
1
0
3
1
0
1
0
1

0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.998
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729

0

0
GW 8_18
G
1
1
GW 904 W19 802
13 04
1
9
GW 04 _05 12
1 2
GW1904 3_13 954
19 21_ 430
04 21 82 3
GW 4_18 856
GW 06
1
9
GW 04 19 48
2 0
GW1905 6_15 425
0 2
1
9GW 05 3_18 155
1 5
1
GW 905 2_18 404
1 0
1
9GW 05 3_20 714
1 5
1
GW 905 4_06 428
19 17_ 541
05 05 6
1 5
GW 9_15 101
GW 35
1
9
GW 05 19 44
2 0
GW1905 1_07 521
2 41
9GW 06 7_09 359
0 2
1
GW 906 2_17 055
2 5
1
9GW 06 0_03 927
3 01
GW 907 0_18 421
0 51
9GW 07 1_20 205
0 3
1
GW 907 6_22 306
0 2
1
9GW 07 7_09 641
0 3
1
GW 907 8_23 326
1 2
19GW 07 9_21 457
2 5
1
GW 907 0_00 514
2 0
1
9
GW 07 7_06 836
2 0
1
GW 907 8_06 333
19 31_ 451
08 14 00 0
GW 3_02 936
GW 27
1
9
GW 08 19 01
2 0
GW1908 8_06 814
2 3
1
9GW 09 8_06 405
0 5
1
GW 909 9_11 509
1 4
19
GW 09 0_11 149
1 2
1GW 909 5_23 807
2 5
19
GW 09 4_02 702
19 29_ 184
09 01 630 21
_13 49
35
41

04

19

GW

Super-Kamiokande livetime [s]Follow-up of O3a catalogued events

8.1
All GW triggers

400

FC

200

PC
UPMU
SampleObserved
Expected
p-value (%)

FC
PC
UPMU
8
0
2
3.95
0.26
0.58
4.84
100.00
11.47

1000

600
expected

0.1

Figure 50: Summary of observations in O3a catalog events.

61

Number of SK events

8

In the new GWTC-2 catalog [7], 39 gravitational waves were reported. Out of these, 36 have SK detector runningand 10 have SK events in time coincidence.
The Figure 50 and the Table 15 are showing the final results for all the events. We can sum the number ofobserved and expected events for each sample to search for any deviation from background expectation.

Table 14: Total number of observed and expected events for O3a observation run, along with the associated Poisson
p-value.

1

800

0.01

Table 15: Summary of all GW triggers for the O3a observation run. The different columns are: Trigger name ; Alert time in UTC ; Most probable GW
event type (with its probability, classifcation based on the masses) ; area of the 90% GW contour ; best-fit distance to the astrophysical object ; SK detector
livetime ; number of observed events and number of expected background events for the four samples.

62

Name

Alert UTC

Type

2
A90%
GW [deg ]d [Mpc]L [s]GW190408 181802
GW190412
GW190413 052954
GW190413 134308
GW190421 213856
GW190424 180648
GW190425
GW190426 152155
GW190503 185404
GW190512 180714
GW190513 205428
GW190514 065416
GW190517 055101
GW190519 153544
GW190521
GW190521 074359
GW190527 092055
GW190602 175927
GW190620 030421
GW190630 185205
GW190701 203306
GW190706 222641
GW190707 093326
GW190708 232457
GW190719 215514
GW190720 000836
GW190727 060333
GW190728 064510
GW190731 140936
GW190803 022701
GW190814
GW190828 063405
GW190828 065509
GW190909 114149

2019-04-08 18:18:02
2019-04-12 05:30:44
2019-04-13 05:29:54
2019-04-13 13:43:08
2019-04-21 21:38:56
2019-04-24 18:06:48
2019-04-25 08:18:05
2019-04-26 15:21:55
2019-05-03 18:54:04
2019-05-12 18:07:14
2019-05-13 20:54:28
2019-05-14 06:54:16
2019-05-17 05:51:01
2019-05-19 15:35:44
2019-05-21 03:02:29
2019-05-21 07:43:59
2019-05-27 09:20:55
2019-06-02 17:59:27
2019-06-20 03:04:21
2019-06-30 18:52:05
2019-07-01 20:33:06
2019-07-06 22:26:41
2019-07-07 09:33:26
2019-07-08 23:24:57
2019-07-19 21:55:14
2019-07-20 00:08:36
2019-07-27 06:03:34
2019-07-28 06:45:10
2019-07-31 14:09:36
2019-08-03 02:27:01
2019-08-14 21:10:38
2019-08-28 06:34:05
2019-08-28 06:55:09
2019-09-09 11:41:49

BBH
BBH
BBH
BBH
BBH
BBH
BNS
NSBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
BBH
NSBH
BBH
BBH
BBH

41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253
41253

1547.5734.1
4189.65181.63165.52568.4156.8
377.2
1527.31462.52189.74987.62270.53023.54566.91244.23562.93138.13210.9956.2
2152.45184.0790.8
887.9
4786.3906.0
3608.9857.64033.73749.6240.7
2160.31657.84923.7

993.1
993.1
993.0
993.1
993.1
993.4
993.4
992.6
993.2
993.7
993.6
993.8
0.0
993.7
994.0
994.0
993.7
993.6
994.5
992.4
992.5
992.4
992.4
992.5
993.0
992.8
992.5
992.5
993.0
993.4
993.9
542.5
0.0
994.0

Nobs

FC
Nexp

Nobs

PC
Nexp

0
0
0
0
0
1
0
0
0
0
1
0
0
0
0
1
1
0
0
0
0
0
0
0
0
0
1
0
0
1
0
0

0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.111
0.061
0.111

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.007
0.004
0.007

UPMU
Nobs Nexp
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0

0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.016
0.009
0.016

LOWE
Nobs Nexp
3
0
0
0
3
1
1
0
0
1
0
1
1
0
0
1
0
1
2
0
2
0
0
1
0
0
2
1
0
0
1
0

0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.729
0.955
0.729
0.729
0.729
0.962
0.729

Table 15: Summary of all GW triggers for the O3a observation run. The different columns are: Trigger name ; Alert time in UTC ; Most probable GW
event type (with its probability, classifcation based on the masses) ; area of the 90% GW contour ; best-fit distance to the astrophysical object ; SK detector
livetime ; number of observed events and number of expected background events for the four samples.
Name

Alert UTC

Type

2
A90%
GW [deg ]d [Mpc]L [s]GW190910 112807
GW190915 235702
GW190924 021846
GW190929 012149
GW190930 133541

2019-09-10 11:28:07
2019-09-15 23:57:02
2019-09-24 02:18:46
2019-09-29 01:21:49
2019-09-30 13:35:41

BBH
BBH
BBH
BBH
BBH

41253
41253
41253
41253
41253

1670.1
1714.6572.4
3901.5
785.9

994.1
0.0
994.1
994.2
994.1

Nobs

FC
Nexp

Nobs

PC
Nexp

1
1
0
0

0.111
0.111
0.111
0.111

0
0
0
0

0.007
0.007
0.007
0.007

UPMU
Nobs Nexp
0
0
0
0

0.016
0.016
0.016
0.016

LOWE
Nobs Nexp
2
0
0
0

0.729
0.729
0.729
0.729

63

8.2

Flux limits and Fluence limits

we can compute flux limits for all GW triggers taking benefit of the three ATMPD samples and marginalising over
GW localisation. The Table 17 shows a summary of the results for FC+PC and UPMU, νµ + ν̄µ and νe + ν̄e .
Similarly, the Table 16 shows the results with the combined (Λ) method.
Table 16: E 2 dn/dE limits [in GeV cm−2 ] for all GW triggers with the combined method.
Name
GW190408 181802
GW190412
GW190413 052954
GW190413 134308
GW190421 213856
GW190424 180648
GW190425
GW190426 152155
GW190503 185404
GW190512 180714
GW190513 205428
GW190514 065416
GW190519 153544
GW190521
GW190521 074359
GW190527 092055
GW190602 175927
GW190620 030421
GW190630 185205
GW190701 203306
GW190706 222641
GW190707 093326
GW190708 232457
GW190719 215514
GW190720 000836
GW190727 060333
GW190728 064510
GW190731 140936
GW190803 022701
GW190814
GW190828 063405
GW190909 114149
GW190910 112807
GW190924 021846
GW190929 012149
GW190930 133541

νµ

ν̄µ
3

2.31 · 10
3.14 · 101
2.23 · 103
2.20 · 103
3.42 · 101
2.67 · 103
2.16 · 103
2.35 · 103
4.29 · 101
2.32 · 103
3.66 · 103
2.27 · 103
2.06 · 103
3.75 · 101
3.52 · 101
3.47 · 101
5.61 · 101
2.11 · 103
2.28 · 103
2.34 · 103
1.29 · 103
1.87 · 103
2.15 · 103
1.92 · 103
3.36 · 101
3.58 · 101
3.59 · 101
3.42 · 101
2.33 · 103
2.16 · 103
1.31 · 103
2.10 · 103
2.36 · 103
2.46 · 103
1.68 · 103
2.34 · 103

νe3

4.56 · 10
4.39 · 101
4.41 · 103
4.55 · 103
4.51 · 101
4.92 · 103
4.20 · 103
4.08 · 103
9.02 · 101
4.53 · 103
6.68 · 103
4.45 · 103
4.31 · 103
4.82 · 101
4.04 · 101
4.49 · 101
7.58 · 101
4.14 · 103
4.34 · 103
3.75 · 103
2.25 · 103
3.85 · 103
4.18 · 103
3.78 · 103
4.76 · 101
5.02 · 101
4.00 · 101
4.46 · 101
4.22 · 103
4.58 · 103
2.03 · 103
4.29 · 103
4.46 · 103
4.09 · 103
3.32 · 103
5.47 · 103

64

ν̄e
3

1.82 · 10
2.11 · 103
1.83 · 103
1.84 · 103
2.24 · 103
2.99 · 103
2.09 · 103
1.97 · 103
2.10 · 103
1.82 · 103
2.93 · 103
1.86 · 103
1.90 · 103
2.21 · 103
2.26 · 103
2.77 · 103
3.89 · 103
2.16 · 103
2.08 · 103
2.13 · 103
2.23 · 103
2.11 · 103
2.09 · 103
2.09 · 103
2.22 · 103
2.17 · 103
3.30 · 103
2.19 · 103
2.11 · 103
1.82 · 103
2.25 · 103
1.97 · 103
3.09 · 103
2.26 · 103
2.25 · 103
1.96 · 103

3.13 · 103
4.64 · 103
3.31 · 103
3.53 · 103
4.87 · 103
5.81 · 103
4.28 · 103
3.89 · 103
4.20 · 103
3.14 · 103
5.03 · 103
3.29 · 103
3.39 · 103
4.60 · 103
4.86 · 103
6.21 · 103
8.28 · 103
4.39 · 103
4.46 · 103
4.74 · 103
3.33 · 103
3.28 · 103
4.13 · 103
3.32 · 103
4.54 · 103
4.53 · 103
6.95 · 103
4.86 · 103
4.82 · 103
3.62 · 103
3.48 · 103
3.40 · 103
6.08 · 103
5.05 · 103
3.65 · 103
4.52 · 103

Table 17: Summary of E 2 dn/dE limits [in GeV cm−2 ] for all GW triggers. In the UPMU columns, empty values corresponds to cases where the GW
localisation is mainly above the horizon.
FC

Name
νµ

65

GW190408 181802
GW190412
GW190413 052954
GW190413 134308
GW190421 213856
GW190424 180648
GW190425
GW190426 152155
GW190503 185404
GW190512 180714
GW190513 205428
GW190514 065416
GW190519 153544
GW190521
GW190521 074359
GW190527 092055
GW190602 175927
GW190620 030421
GW190630 185205
GW190701 203306
GW190706 222641
GW190707 093326
GW190708 232457
GW190719 215514
GW190720 000836
GW190727 060333
GW190728 064510
GW190731 140936
GW190803 022701
GW190814
GW190828 063405
GW190909 114149
GW190910 112807
GW190924 021846
GW190929 012149
GW190930 133541

ν̄µ34.02 · 10
3.95 · 103
4.00 · 103
3.97 · 103
3.93 · 103
6.50 · 103
3.95 · 103
4.01 · 103
4.03 · 103
4.02 · 103
6.61 · 103
4.01 · 103
3.96 · 103
3.96 · 103
3.91 · 103
6.39 · 103
6.29 · 103
3.95 · 103
3.94 · 103
3.95 · 103
4.01 · 103
4.00 · 103
3.96 · 103
3.98 · 103
3.94 · 103
3.96 · 103
6.46 · 103
3.92 · 103
3.95 · 103
6.48 · 103
4.01 · 103
3.99 · 103
6.49 · 103
6.49 · 103
4.00 · 103
3.81 · 103

PC
νe

3

9.00 · 10
9.40 · 103
9.23 · 103
9.62 · 103
9.05 · 103
1.57 · 104
9.47 · 103
8.92 · 103
1.00 · 104
8.96 · 103
1.47 · 104
9.15 · 103
9.84 · 103
9.59 · 103
9.20 · 103
1.49 · 104
1.54 · 104
9.50 · 103
9.29 · 103
9.28 · 103
1.01 · 104
1.01 · 104
9.56 · 103
1.00 · 104
9.67 · 103
9.63 · 103
1.53 · 104
9.05 · 103
9.08 · 103
1.61 · 104
1.01 · 104
9.75 · 103
1.57 · 104
1.51 · 104
1.01 · 104
9.84 · 103ν̄e3

2.10 · 10
2.17 · 103
2.13 · 103
2.11 · 103
2.32 · 103
3.70 · 103
2.25 · 103
2.25 · 103
2.14 · 103
2.11 · 103
3.52 · 103
2.13 · 103
2.15 · 103
2.28 · 103
2.34 · 103
3.59 · 103
3.53 · 103
2.25 · 103
2.24 · 103
2.40 · 103
2.33 · 103
2.28 · 103
2.25 · 103
2.24 · 103
2.27 · 103
2.24 · 103
3.85 · 103
2.28 · 103
2.29 · 103
3.44 · 103
2.37 · 103
2.19 · 103
3.70 · 103
3.85 · 103
2.34 · 103
2.01 · 103νµ33.26 · 10
4.76 · 103
3.49 · 103
3.72 · 103
5.06 · 103
7.32 · 103
4.45 · 103
3.95 · 103
4.19 · 103
3.27 · 103
5.61 · 103
3.43 · 103
3.59 · 103
4.72 · 103
4.94 · 103
8.21 · 103
7.90 · 103
4.53 · 103
4.62 · 103
4.75 · 103
3.42 · 103
3.43 · 103
4.32 · 103
3.46 · 103
4.63 · 103
4.65 · 103
8.07 · 103
5.05 · 103
4.87 · 103
6.30 · 103
3.63 · 103
3.57 · 103
7.40 · 103
7.96 · 103
3.77 · 103
5.05 · 103ν̄µ3

5.26 · 10
4.87 · 103
5.16 · 103
4.89 · 103
4.76 · 103
5.09 · 103
5.08 · 103
5.49 · 103
5.81 · 103
5.29 · 103
5.32 · 103
5.21 · 103
4.82 · 103
4.89 · 103
4.87 · 103
4.62 · 103
4.68 · 103
5.00 · 103
5.26 · 103
5.68 · 103
5.09 · 103
5.01 · 103
5.09 · 103
4.92 · 103
4.90 · 103
4.93 · 103
4.85 · 103
4.69 · 103
5.56 · 103
4.69 · 103
5.11 · 103
4.98 · 103
5.04 · 103
5.67 · 103
5.06 · 103
5.55 · 103

UPMU
νe

3

9.03 · 10
8.43 · 103
8.68 · 103
8.73 · 103
8.34 · 103
8.76 · 103
8.77 · 103
7.66 · 103
8.71 · 103
8.96 · 103
8.58 · 103
8.73 · 103
8.59 · 103
8.34 · 103
7.43 · 103
9.97 · 103
1.03 · 104
8.76 · 103
8.84 · 103
6.22 · 103
8.71 · 103
8.68 · 103
8.59 · 103
8.56 · 103
8.31 · 103
8.42 · 103
7.68 · 103
8.83 · 103
8.24 · 103
8.55 · 103
8.74 · 103
8.66 · 103
8.78 · 103
7.11 · 103
8.66 · 103
1.10 · 104ν̄e4

1.13 · 10
3.76 · 104
1.27 · 104
1.31 · 104
3.02 · 104
3.68 · 104
3.61 · 104
1.29 · 104
4.31 · 104
1.14 · 104
1.18 · 104
1.30 · 104
1.71 · 104
4.00 · 104
3.37 · 104
3.22 · 104
3.30 · 104
3.61 · 104
3.66 · 104
1.58 · 104
3.67 · 104
3.37 · 104
3.59 · 104
3.01 · 104
4.12 · 104
4.04 · 104
3.49 · 104
3.11 · 104
3.40 · 104
1.17 · 104
4.10 · 104
2.77 · 104
3.68 · 104
2.31 · 104
4.06 · 104
4.97 · 104

45.98 · 10
9.24 · 104
6.57 · 104
5.51 · 104
8.06 · 104
1.15 · 105
1.15 · 105
1.15 · 105
2.23 · 105
6.35 · 104
7.29 · 104
6.71 · 104
5.19 · 104
1.02 · 105
1.17 · 105
7.41 · 104
7.23 · 104
1.07 · 105
1.50 · 105
1.70 · 105
7.51 · 104
6.94 · 104
1.17 · 105
6.40 · 104
9.52 · 104
9.14 · 104
1.18 · 105
8.31 · 104
1.99 · 105
4.68 · 104
8.06 · 104
6.41 · 104
1.11 · 105
2.06 · 105
7.99 · 104
1.51 · 105

νµ

ν̄µ

(unseen)
3.18 · 101(unseen)(unseen)
3.47 · 101(unseen)(unseen)(unseen)
4.37 · 101(unseen)(unseen)(unseen)(unseen)3.75 · 101
3.58 · 101
3.50 · 101
3.53 · 101(unseen)(unseen)(unseen)(unseen)(unseen)(unseen)(unseen)
3.43 · 101
3.64 · 101
3.57 · 101
3.48 · 101(unseen)(unseen)(unseen)(unseen)(unseen)(unseen)(unseen)(unseen)(unseen)
4.42 · 101(unseen)(unseen)
4.56 · 101(unseen)(unseen)(unseen)
9.20 · 101(unseen)(unseen)(unseen)(unseen)
4.80 · 101
4.09 · 101
4.52 · 101
4.71 · 101(unseen)(unseen)(unseen)(unseen)(unseen)(unseen)(unseen)
4.78 · 101
5.10 · 101
4.00 · 101
4.50 · 101(unseen)(unseen)(unseen)(unseen)(unseen)(unseen)(unseen)(unseen)Results can also be summarised in plots, with horizontal=trigger and vertical=limits, like presented in Figure 51.
Also, we compute the fluence limit for Low energy sample.Plot Figure 52 shows all triggers in low energy region.

8.3

Total energy limits

As discussed in section 5.7, we can compute upper limits on the isotropic energy emitted in neutrinos by the source.
They are presented in the Table 18 and in the Figure 53.
Trigger
GW190408 181802
GW190412
GW190413 052954
GW190413 134308
GW190421 213856
GW190424 180648
GW190425
GW190426 152155
GW190503 185404
GW190512 180714
GW190513 205428
GW190514 065416
GW190519 153544
GW190521
GW190521 074359
GW190527 092055
GW190602 175927
GW190620 030421
GW190630 185205
GW190701 203306
GW190706 222641
GW190707 093326
GW190708 232457
GW190719 215514
GW190720 000836
GW190727 060333
GW190728 064510
GW190731 140936
GW190803 022701
GW190814
GW190828 063405
GW190909 114149
GW190910 112807
GW190924 021846
GW190929 012149
GW190930 133541

dGW [Mpc]
1548
734
4190
5182
3166
2568
157
377
1527
1462
2190
4988
3024
4567
1244
3563
3138
3211
956
2152
5184
791
888
4786
906
3609
858
4034
3750
241
2160
4924
1670
572
3902
786

νµ

ν̄µ
58

1.85 · 10
5.25 · 1055
1.41 · 1059
2.11 · 1059
1.25 · 1057
7.45 · 1058
1.92 · 1056
1.21 · 1057
3.45 · 1056
1.69 · 1058
6.21 · 1058
2.15 · 1059
6.59 · 1058
2.77 · 1057
1.94 · 1056
6.21 · 1060
2.36 · 1057
7.50 · 1058
9.44 · 1057
3.66 · 1058
1.11 · 1059
4.60 · 1057
5.44 · 1057
2.02 · 1059
1.34 · 1056
1.63 · 1057
9.10 · 1055
2.27 · 1057
1.21 · 1059
3.82 · 1056
1.82 · 1058
2.40 · 1059
2.38 · 1058
2.80 · 1057
8.18 · 1058
5.22 · 1057

νe58

3.63 · 10
7.33 · 1055
2.79 · 1059
4.36 · 1059
1.64 · 1057
1.31 · 1059
3.71 · 1056
2.11 · 1057
7.31 · 1056
3.32 · 1058
1.12 · 1059
4.25 · 1059
1.40 · 1059
3.67 · 1057
2.16 · 1056
1.04 · 1061
3.15 · 1057
1.49 · 1059
1.58 · 1058
5.84 · 1058
1.83 · 1059
9.26 · 1057
1.08 · 1058
2.94 · 1059
1.72 · 1056
2.28 · 1057
9.83 · 1055
2.94 · 1057
2.13 · 1059
8.12 · 1056
2.34 · 1058
4.87 · 1060
4.38 · 1058
4.61 · 1057
2.26 · 1059
1.17 · 1058ν̄e58

1.44 · 10
3.54 · 1057
1.16 · 1059
1.75 · 1059
8.18 · 1058
9.08 · 1058
1.97 · 1056
9.96 · 1056
1.70 · 1058
1.33 · 1058
4.96 · 1058
1.73 · 1059
6.38 · 1058
1.69 · 1059
1.23 · 1058
1.34 · 1061
1.79 · 1059
8.65 · 1058
8.08 · 1057
3.31 · 1058
2.09 · 1059
5.30 · 1057
6.29 · 1057
1.82 · 1059
8.87 · 1057
1.00 · 1059
8.40 · 1057
1.48 · 1059
1.13 · 1059
3.23 · 1056
3.58 · 1058
7.96 · 1060
3.46 · 1058
2.60 · 1057
2.00 · 1059
4.47 · 1057

All-flavours
58

2.48 · 10
7.80 · 1057
2.07 · 1059
3.36 · 1059
1.77 · 1059
1.61 · 1059
3.83 · 1056
1.92 · 1057
3.37 · 1058
2.28 · 1058
8.41 · 1058
3.04 · 1059
1.11 · 1059
3.47 · 1059
2.67 · 1058
1.37 · 1061
3.80 · 1059
1.70 · 1059
1.56 · 1058
7.42 · 1058
3.06 · 1059
8.16 · 1057
1.17 · 1058
2.89 · 1059
1.92 · 1058
2.09 · 1059
1.76 · 1058
3.16 · 1059
2.59 · 1059
6.43 · 1056
5.58 · 1058
4.99 · 1060
6.50 · 1058
5.72 · 1057
2.08 · 1059
1.07 · 1058

3.13 · 1058
1.81 · 1056
2.51 · 1059
3.81 · 1059
4.25 · 1057
1.30 · 1059
3.40 · 1056
2.10 · 1057
1.38 · 1057
2.87 · 1058
1.05 · 1059
3.70 · 1059
1.21 · 1059
8.86 · 1057
6.09 · 1056
1.43 · 1061
8.22 · 1057
1.39 · 1059
1.60 · 1058
6.79 · 1058
1.97 · 1059
8.20 · 1057
9.95 · 1057
3.15 · 1059
4.34 · 1056
5.62 · 1057
2.84 · 1056
7.51 · 1057
2.30 · 1059
7.05 · 1056
3.13 · 1058
3.78 · 1059
4.39 · 1058
5.27 · 1057
1.43 · 1059
1.00 · 1058

Table 18: 90& upper limits on Eiso [in erg] either flavour-by-flavour or by combining flavours and putting limits on
the all-flavour emission, assuming equipartition.

66

GW190620_030421

GW190602_175927

GW190527_092055

triggers

67
PC
GW190728_064510
GW190731_140936
GW190803_022701
GW190814
GW190828_063405
GW190909_114149
GW190910_112807
GW190924_021846
GW190929_012149
GW190930_133541

GW190728_064510

GW190731_140936

GW190803_022701

GW190814

GW190828_063405

GW190909_114149

GW190910_112807

GW190924_021846

GW190929_012149

GW190930_133541

103
GW190727_060333

104

GW190727_060333

UPMU
GW190720_000836(c) νe

GW190720_000836

GW190719_215514

PC

GW190708_232457(d) ν̄eFigure 51: Summary of all flux limits.
GW190930_133541

GW190929_012149

GW190924_021846

GW190910_112807

GW190909_114149

GW190828_063405

GW190814

GW190803_022701

GW190731_140936

GW190728_064510

GW190727_060333

GW190720_000836

GW190719_215514

GW190708_232457

GW190707_093326

PC

GW190719_215514

FC
GW190706_222641

GW190930_133541

GW190929_012149

GW190924_021846

GW190910_112807

GW190909_114149

GW190828_063405

GW190814

GW190803_022701

GW190731_140936

GW190728_064510

GW190727_060333

GW190720_000836

GW190719_215514

GW190708_232457

GW190707_093326

GW190706_222641

GW190701_203306

GW190630_185205

PC

GW190708_232457

Combined
GW190707_093326

FC

GW190706_222641

FC

GW190707_093326

triggers
GW190701_203306

2E 2dn
dE [GeV cm ]FC

GW190706_222641

Combined

GW190701_203306

GW190620_030421

GW190602_175927

GW190527_092055

GW190521_074359

GW190521

GW190519_153544

GW190514_065416

GW190513_205428

GW190512_180714

GW190503_185404

GW190426_152155

GW190425

GW190424_180648

GW190421_213856

GW190413_134308

GW190413_052954

GW190412

GW190408_181802

triggers

GW190701_203306

triggers
GW190630_185205

GW190620_030421

GW190602_175927

GW190527_092055

Combined

GW190630_185205

GW190620_030421

GW190602_175927

GW190527_092055

GW190521_074359

GW190521

GW190519_153544

GW190514_065416

GW190513_205428

GW190512_180714

GW190503_185404

GW190426_152155

GW190425

GW190424_180648

GW190421_213856

GW190413_134308

GW190413_052954

GW190412

GW190408_181802

2E 2dn
dE [GeV cm ]Combined

GW190630_185205

105
GW190521_074359

GW190521

GW190519_153544

GW190514_065416

GW190513_205428

GW190512_180714

GW190503_185404

GW190426_152155

GW190425

GW190424_180648

GW190421_213856

GW190413_134308

GW190413_052954

GW190412

GW190408_181802

2E 2dn
dE [GeV cm ]105

GW190521_074359

GW190521

GW190519_153544

GW190514_065416

GW190513_205428

GW190512_180714

GW190503_185404

GW190426_152155

GW190425

GW190424_180648

GW190421_213856

GW190413_134308

GW190413_052954

GW190412

GW190408_181802

2E 2dn
dE [GeV cm ]

104
UPMU

103

102

101

(a) νµ

104
UPMU

103

102

101

(b) ν̄µ

UPMU

104

103

Figure 52: LOW energy triggers fluence limit

9

Summary and Outlook

9.1

Summary

In this thesis, neutrinos related to gravitational waves are searched for in SK detectors. Gravitational waves canbe produced by the merger of binary neutron stars and black holes. Among them, the merger of neutron stars is
believed to produce light, gamma rays, neutrinos, etc. in addition to gravitational waves. According to different
models at present, the energy of the neutrinos produced is different, and detecting neutrinos will help us to study
the models of binary neutron star merger.
Until now, no clear neutrino signal is observed with GW. In the previous study of SK, neutrinos for GW170817
in SK were searched for and no significant signal found. In this study, 56 GWs for O3(from Apr. 2019 to Mar.
2020) and 39 GWs for O3a (from Apr. 2019 to Sep. 2019) are the target. After reduction process in SK, 46 events
remain in the search window(for O3) and 36 events for O3a, ±500 sec around each GWs trigger time, while thebackground rate is 0.111 events/1000 sec for FC, 0.007 events/1000 sec for PC, 0.016 events/1000 sec for UPMU
and 0.729 events/1000 sec for LOW energy sample(energy < 100M eV ). No obvious excess was seen in the timedevelopment of the number of events and energy distribution of the combined data.
Although we did not find a significant signal in this study, we can compute Upper-limit for each flavors in bothHigh energy and Low energy. For High energy, we assume E −2 neutrinos spectrum to compute the flux limit. For
Low energy, we assume Fermi-Dirac Distribution(average 20 MeV), Flat and mono-energetic spectrum to calculate
Fluence-Limit. In both cases, the detailed results have been presented. We also compute the isotropic total neutrinoenergy limits for each neutrino flavors. Assuming isotropic emissions among the different flavors, upper limits on
the total energy as Eiso of each types of neutrinos were derived. The limits on Eiso for individual GW events, andcombined all GW events at each neutrino flavor are shown in Table 18 and Figure 53.

9.2

Outlook

The focus of this research is on the O3 and O3a phases, but as this paper presented in section 3.3, LIGO will start
the O4 phase (See section 3.2) at next year. In O4 phase, the expected events for GWs are 1 per day. For SK
detector, SK-Gd has started. In SK-Gd, 0.2% gadolinium sulfate Gds(SO4 )3 is added to the SK water in order
to identify ν e events. The identify of ν e is useful for the detection of supernova relic neutrino and direction ofsupernova bursts(Which is the background of GW research). Gadolinium captures a neutrino from inverse beta
decay and emits gamma rays of 8M eV total energy (Figure 54)
In O4 phase, the background of Low energy case will be changed because of SK-Gd. Also more data can beused and significant neutrino signal is expected be found.
In our current automatic system, realtime analysis is not possible, however, development a real-time system iscritical to get neutrino signals more quickly.

68

1061
1060

1060

90% [erg]
Eiso

1059
90% [erg]
Eiso

BNS
NSBH
BBH

1061

BNS
NSBH
BBH

1058

1059
1058

1057

1057

1056

1056
102

103
distance [Mpc]

102

103
distance [Mpc](a) νµ(b) ν̄µ

BNS
NSBH
BBH

1061
1060

1060

1059

90% [erg]
Eiso

90% [erg]
Eiso

BNS
NSBH
BBH

1061

1058

1059
1058

1057

1057
102

103
distance [Mpc]102(c) νe103
distance [Mpc](d) ν̄e
BNS
NSBH
BBH

1061

90% [erg]
Eiso

1060
1059
1058
1057

102

103
distance [Mpc](e) All-flavours

Figure 53: 90& upper limits on Eiso [in erg] (either flavour-by-flavour or by combining flavours and putting limits
on the all-flavour emission, assuming equipartition), as a function of source distance and separating by source type.

69

Figure 54: Nertron tagging with gadolinium.[28]Acknowledgement
Thanks As an international student, the study life of a half-year graduate student and a two-year master in Japan
is very fulfilling and rich. When I first came to Japan, communication and communication seemed a little difficult
due to language problems. With the careful guidance of teacher Koshio-san, Harada-san of D1, and the help of
Okada-san and Sakai-san who are both M2, I gradually got used to life in Japan. Thank you teachers and everyone
for your concern for me. After I officially became a master student, the tutor Koshio-san was very dedicated to
guiding my research. At the same time, I also communicated with the members of SK, which made me progress. I
would like to thank my tutors for their guidance in learning, and also thank all SK members for their help in myresearch. At the same time, in the research room, we will hold a meeting once a week. During the meeting, everyone
will introduce the content of their research and also study the paper. This also allows me to learn more than thecontent of my own research. know how. Thanks to the members of the research room, Special Researcher Ito-san,Harada-san of D1, Sakai-san and Okada-san of M2, Kitagawa-san of M1, Nakanishi-san of B4, and Shiraishi-san,Tano-san of B3. Finally, I would like to thank Yamamoto-san and Koeda-san, the clerks of the research office, foryour support for making my research life smooth.
Thanks to the teachers and everyone who have always supported and encouraged me. It was your support andhelp that allowed me to complete my two-year research life as a master student.

70References[1] M. G. Aartsen et al. A Search for MeV to TeV Neutrinos from Fast Radio Bursts with IceCube. 2019.
[2] M. G. Aartsen et al. IceCube Search for Neutrinos Coincident with Compact Binary Mergers from LIGO-Virgo’s
First Gravitational-Wave Transient Catalog IceCube Collaboration. Technical report, 2020.
[3] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso,
R. Adhikari, et al. Gw150914: The advanced ligo detectors in the era of first discoveries. Physical review letters,
116(13):131103, 2016.
[4] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso,
R. Adhikari, et al. Observation of gravitational waves from a binary black hole merger. Physical review letters,
116(6):061102, 2016.
[5] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso,
R. Adhikari, et al. Observing gravitational-wave transient gw150914 with minimal assumptions. PhysicalReview D, 93(12):122004, 2016.
[6] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso,
R. Adhikari, et al. Properties of the binary black hole merger gw150914. Physical review letters, 116(24):241102,
2016.
[7] R. Abbott et al. GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Halfof the Third Observing Run. 10 2020.
[8] K. Abe et al. Search for Neutrinos in Super-Kamiokande associated with Gravitational Wave Events GW150914
and GW151226. 08 2016.
[9] K. Abe et al. Search for Neutrinos in Super-Kamiokande associated with the GW170817 neutron-star merger.
The Astrophysical Journal, 857(1):L4, 02 2018.
[10] B. Baret et al. Multimessenger Science Reach and Analysis Method for Common Sources of Gravitational
Waves and High-energy Neutrinos. Technical report, 12 2011.
[11] S. Barthelmy. Gcn: The gamma-ray coordinates network. https://gcn.gsfc.nasa.gov/.
[12] L. Blanchet, T. Damour, B. R. Iyer, C. M. Will, and A. G. Wiseman. Gravitational-radiation damping ofcompact binary systems to second post-newtonian order. Physical Review Letters, 74(18):3515, 1995.
[13] J. Braun, J. Dumm, F. De Palma, C. Finley, A. Karle, and T. Montaruli. Methods for point source analysis
in high energy neutrino telescopes. Astroparticle Physics, 29(4):299–305, 2008.
[14] O. Caballero. Neutrino emission, equation of state and the role of strong gravity. In AIP Conference Proceedings,
volume 1753, page 040002. AIP Publishing LLC, 2016.
[15] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlochower. Accurate evolutions of orbiting black-holebinaries without excision. Physical Review Letters, 96(11):111101, 2006.
[16] S. Chatterji, L. Blackburn, G. Martin, and E. Katsavounidis. Multiresolution techniques for the detection ofgravitational-wave bursts. Classical and Quantum Gravity, 21(20):S1809–S1818, sep 2004.
[17] N. J. Cornish and T. B. Littenberg. Bayeswave: Bayesian inference for gravitational wave bursts and instrument
glitches. Classical and Quantum Gravity, 32(13):135012, 2015.
[18] S. Fukuda et al. Search for Neutrinos from Gamma-Ray Bursts Using Super-Kamiokande. The AstrophysicalJournal, 578(1):317–324, 10 2002.
[19] S. Fukuda, Y. Fukuda, T. Hayakawa, E. Ichihara, M. Ishitsuka, Y. Itow, T. Kajita, J. Kameda, K. Kaneyuki,
S. Kasuga, et al. The super-kamiokande detector. Nuclear Instruments and Methods in Physics Research
Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 501(2-3):418–462, 2003.
[20] É. Gourgoulhon. Special relativity in general frames. Springer, 2016.

71

[21] K. Hagiwara et al. Search for Astronomical Neutrinos from Blazar TXS 0506+056 in Super-Kamiokande. The
Astrophysical Journal Letters, 887:L6, 2019.
[22] R. Hussain, J. Vandenbroucke, and J. Wood. A Search for IceCube Neutrinos from the First 33 Detected
Gravitational Wave Events. 2019.
[23] R. P. Kraft, D. N. Burrows, and J. A. Nousek. Determination of Confidence Limits for Experiments with Low
Numbers of Counts. Astrophysical Journal, 374:344, June 1991.
[24] A. Le Tiec and J. Novak. Theory of gravitational waves. In An Overview of Gravitational Waves: Theory,
Sources and Detection, pages 1–41. World Scientific, 2017.
[25] LSC. Gracedb gravitational-wave candidate event database. https://gracedb.ligo.org/.
[26] A. H. Mroue, M. A. Scheel, B. Szilagyi, H. P. Pfeiffer, M. Boyle, D. A. Hemberger, L. E. Kidder, G. Lovelace,
S. Ossokine, N. W. Taylor, et al. Catalog of 174 binary black hole simulations for gravitational wave astronomy.
Physical Review Letters, 111(24):241104, 2013.
[27] H. Nishino, K. Awai, Y. Hayato, S. Nakayama, K. Okumura, M. Shiozawa, A. Takeda, K. Ishikawa,
A. Minegishi, and Y. Arai. High-speed charge-to-time converter asic for the super-kamiokande detector. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors andAssociated Equipment, 610(3):710–717, 2009.
[28] A. Orii. Search for Neutrinos associated with Gamma-ray Bursts in Super-Kamiokande. PhD thesis, School of
Science, The University of Tokyo, 2018.
[29] A. Suzuki, M. Mori, K. Kaneyuki, T. Tanimori, J. Takeuchi, H. Kyushima, and Y. Ohashi. Improvementof 20 in. diameter photomultiplier tubes. Nuclear Instruments and Methods in Physics Research Section A:
Accelerators, Spectrometers, Detectors and Associated Equipment, 329(1-2):299–313, 1993.
[30] E. Thrane et al. Search for neutrinos from GRB 080319B at super-kamiokande. Astrophysical Journal,697(1):730–734, 03 2009.
[31] D. Veske, Z. Márka, I. Bartos, and S. Márka. Neutrino emission upper limits with maximum likelihood
estimators for joint astrophysical neutrino searches with large sky localizations. JCAP, 05:016, 2020.
[32] S. Yamada, K. Awai, Y. Hayato, K. Kaneyuki, Y. Kouzuma, S. Nakayama, H. Nishino, K. Okumura,Y. Obayashi, Y. Shimizu, et al. Commissioning of the new electronics and online system for the superkamiokande experiment. IEEE Transactions on Nuclear Science, 57(2):428–432, 2010.

72