Astroparticle Physics 81 (2016) 39–48

Contents lists available at ScienceDirect

Astroparticle Physics
journal homepage: www.elsevier.com/locate/astropartphys

Real-time supernova neutrino burst monitor at Super-Kamiokande
K. Abe a,C, Y. Haga a, Y. Hayato a,C, M. Ikeda a,C, K. Iyogi a, J. Kameda a,C, Y. Kishimoto a,C,
M. Miura a,C, S. Moriyama a,C, M. Nakahata a,C, Y. Nakano a, S. Nakayama a,C, H. Sekiya a,C,
M. Shiozawa a,C, Y. Suzuki a,C, A. Takeda a,C, H. Tanaka a,C, T. Tomura a,C, K. Ueno a,
R.A. Wendell a,C, T. Yokozawa a, T. Irvine b, T. Kajita b,C, I. Kametani b, K. Kaneyuki b,C,1,
K.P. Lee b, T. McLachlan b, Y. Nishimura b, E. Richard b, K. Okumura b,C, L. Labarga c,
P. Fernandez c, S. Berkman e, H.A. Tanaka e, S. Tobayama e, J. Gustafson d, E. Kearns d,C,
J.L. Raaf d, J.L. Stone d,C, L.R. Sulak d, M. Goldhaber f,1, G. Carminati g, W.R. Kropp g, S. Mine g,
P. Weatherly g, A. Renshaw g, M.B. Smy g,C, H.W. Sobel g,C, V. Takhistov g, K.S. Ganezer h,
B.L. Hartﬁel h, J. Hill h, W.E. Keig h, N. Hong i, J.Y. Kim i, I.T. Lim i, T. Akiri j, A. Himmel j,
K. Scholberg j,C, C.W. Walter j,C, T. Wongjirad j, T. Ishizuka k, S. Tasaka l, J.S. Jang m,
J.G. Learned n, S. Matsuno n, S.N. Smith n, T. Hasegawa o, T. Ishida o, T. Ishii o, T. Kobayashi o,T. Nakadaira o, K. Nakamura o,C, Y. Oyama o, K. Sakashita o, T. Sekiguchi o, T. Tsukamoto o,
A.T. Suzuki p, Y. Takeuchi p,C, C. Bronner q, S. Hirota q, K. Huang q, K. Ieki q, T. Kikawa q,
A. Minamino q, A. Murakami q, T. Nakaya q,C, K. Suzuki q, S. Takahashi q, K. Tateishi q,
Y. Fukuda r, K. Choi s, Y. Itow s, G. Mitsuka s, P. Mijakowski H, J. Hignight t, J. Imber t,
C.K. Jung t, C. Yanagisawa t, M.J. Wilking t, H. Ishino u,∗, A. Kibayashi u, Y. Koshio u,C, T. Mori u,
M. Sakuda u, R. Yamaguchi u, T. Yano u, Y. Kuno v, R. Tacik w,E, S.B. Kim x, H. Okazawa y,
Y. Choi z, K. Nishijima A, M. Koshiba B, Y. Suda B, Y. Totsuka B,1, M. Yokoyama B,C, K. Martens C,
Ll. Marti C, M.R. Vagins C,g, J.F. Martin D, P. de Perio D, A. Konaka E, S. Chen F, Y. Zhang F,
K. Connolly G, R.J. Wilkes GaKamioka Observatory, Institute for Cosmic Ray Research, University of Tokyo, Kamioka, Gifu 506-1205, Japan
Research Center for Cosmic Neutrinos, Institute for Cosmic Ray Research, University of Tokyo, Kashiwa, Chiba 277-8582, Japanc
Department of Theoretical Physics, University Autonoma Madrid, Madrid 28049, Spaind
Department of Physics, Boston University, Boston, MA 02215, USAe
Department of Physics and Astronomy, University of British Columbia, Vancouver, BC V6T1Z4, CanadafPhysics Department, Brookhaven National Laboratory, Upton, NY 11973, USAg
Department of Physics and Astronomy, University of California, Irvine, Irvine, CA 92697-4575, USAhDepartment of Physics, California State University, Dominguez Hills, Carson, CA 90747, USAi
Department of Physics, Chonnam National University, Kwangju 500-757, Republic of Koreaj
Department of Physics, Duke University, Durham, NC, 27708 USA
kJunior College, Fukuoka Institute of Technology, Fukuoka, Fukuoka 811-0295, Japanl
Department of Physics, Gifu University, Gifu 501-1193, Japanm
GIST College, Gwangju Institute of Science and Technology, Gwangju 500-712, Korean
Department of Physics and Astronomy, University of Hawaii, Honolulu, HI 96822, USAoHigh Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japanp
Department of Physics, Kobe University, Kobe, Hyogo 657-8501, Japanq
Department of Physics, Kyoto University, Kyoto 606-8502, Japanr
Department of Physics, Miyagi University of Education, Sendai, Miyagi 980-0845, Japans
Solar Terrestrial Environment Laboratory, Nagoya University, Nagoya, Aichi 464-8602, Japant
Department of Physics and Astronomy, State University of New York at Stony Brook, NY 11794-3800, USAu
Department of Physics, Okayama University, Okayama 700-8530, Japanv
Department of Physics, Osaka University, Toyonaka, Osaka 560-0043, Japanw
Department of Physics, University of Regina, 3737 Wascana Parkway, Regina, SK S4SOA2, Canadax
Department of Physics, Seoul National University, Seoul 151-742, Republic of Koreay
Department of Informatics in Social Welfare, Shizuoka University of Welfare, Yaizu, Shizuoka 425-8611, Japanb∗

Corresponding author. Tel.: +81 862517818.
E-mail address: scishino@s.okayama-u.ac.jp (H. Ishino).
1
Deceased.
http://dx.doi.org/10.1016/j.astropartphys.2016.04.003
0927-6505/© 2016 Published by Elsevier B.V.

40

K. Abe et al. / Astroparticle Physics 81 (2016) 39–48

zDepartment of Physics, Sungkyunkwan University, Suwon 440-746, Republic of Korea
Department of Physics, Tokai University, Hiratsuka, Kanagawa 259-1292, JapanB
The University of Tokyo, Bunkyo, Tokyo 113-0033, JapanC
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, University of Tokyo, Kashiwa,
Chiba 277-8583, JapanD
Department of Physics, University of Torront, 60 St., Torront, Ontario M5S1A7, CanadaE
TRIUMF, 4004 Wesbrook Mall, Vancouver, BC V6T 2A3, CanadaF
Department of Engineering Physics, Tsinghua University, Beijing 100084, China
G
Department of Physics, University of Washington, Seattle, WA 98195-1560, USA
H
National Centre For Nuclear Research, Warsaw 00-681, PolandA

a r t i c l e

i n f oArticle history:
Received 18 January 2016
Revised 15 March 2016
Accepted 1 April 2016
Available online 19 April 2016
Keywords:
Supernova
NeutrinosSuper-Kamiokandea b s t r a c t
We present a real-time supernova neutrino burst monitor at Super-Kamiokande (SK). Detecting supernova
explosions by neutrinos in real time is crucial for giving a clear picture of the explosion mechanism.
Since the neutrinos are expected to come earlier than light, a fast broadcasting of the detection may
give astronomers a chance to make electromagnetic radiation observations of the explosions right at theonset. The role of the monitor includes a fast announcement of the neutrino burst detection to the world
and a determination of the supernova direction. We present the online neutrino burst detection systemand studies of the direction determination accuracy based on simulations at SK.

1. Introduction
The detection of neutrinos from SN1987A opened a new era ofneutrino astronomy [1]. Although the number of the detected neutrino events [2] was small, signiﬁcant information about the supernova (SN) explosion and neutrino properties was obtained [3]. The
current generation of detectors are waiting for the next SN neutrino burst to accumulate event statistics much larger than thoseof SN1987A.
The SN burst neutrinos arrive at the Earth earlier than the electromagnetic radiation, since the neutrinos generated at the core of
the explosion and emitted from the surface of the neutrinospheretravel at nearly the speed of light, while the shock waves propagating to the outside with velocity much slower than the neutrino
velocity result in the emission of electromagnetic radiation [4]. The
delay may depend on the structure of the envelope of the core aswell as the surrounding stellar environment, and is expected torange between tens of minutes and tens of hours [5]. Therefore,
the detection of the neutrino burst can generate a warning able toallow the astronomers to observe the radiation from the onset ofthe explosion. Such warning systems have been developed by several neutrino observatories [6,7] as well as by the supernova early
warning system (SNEWS) [8].
It is also important to determine the SN direction using theneutrino signal: the direction information can guide optical instruments toward the SN explosion and enable observation of
the onset of radiation. Among the neutrino detectors operating atpresent, Super-Kamiokande (SK) is the only detector able to determine the SN direction using neutrino events. We have developed a
SN direction determination method by applying a maximum likelihood ﬁt.
SK is the world’s largest water Cherenkov detector located
10 0 0 m underground, inside a mountain in Kamioka, Gifu, Japan.
The detector consists of 50,0 0 0 tons of ultra-pure water and about
13,0 0 0 photomultiplier tubes (PMTs). Based on the information of
the yields and the arrival timing of Cherenkov photons for individual PMTs as well as the PMT locations, SK is able to measure the position, direction and energy of a neutrino event inreal time. Details of the SK detector and its performance are described in [9]. In 2008, SK upgraded its readout system. The system has improved the data processing speed signiﬁcantly, lowering the trigger energy threshold and dead time for the SN burst
events [10].

© 2016 Published by Elsevier B.V.

We will describe the SN neutrino burst monitor at SK. In
Section 2, we will describe details of the monitor system and its
performance. We ﬁrst describe SN models used in this report for
performance evaluations. A Monte-Carlo (MC) simulation with the
SN models is utilized for the estimation of the detection eﬃciency.
We will explain the selection criteria to discriminate SN bursts andthe main background in SK: radioactive decays caused by cosmic
ray muon spallation [11]. In Section 3, we will describe a method
to reconstruct the SN direction and studies of its performance.
2. Real time supernova neutrino burst monitor
2.1. Supernova models
We ﬁrst describe the SN models we use in this report. We employ two models for the SN neutrino burst: the Wilson model[12] and the Nakazato model [13]. For the Nakazato model, wechoose two parameter sets: M = 20, trevive = 200 ms and Z = 0.02
(NK1), and M = 13, trevive = 100 ms and Z = 0.004 (NK2), where
M is the progenitor mass in units of the solar mass, trevive is theshock revival time and Z is the metallicity, respectively. We choose
the ﬁrst one for the SN1987A progenitor mass, which is about 20
times the solar mass. The model with the latter parameters gives
the smallest neutrino ﬂuxes in the Nakazato model. Both models provide time dependences of neutrino luminosities and energy
spectra for ν e , ν̄e and ν x , for 18 s (20 s) for the Wilson (NK) model,
where ν x refers to the muon and tau types of neutrinos and antineutrinos. The anti-electron neutrino ﬂuences of the models of
Wilson, NK1 and NK2 are 16.0, 9.8 and 9.4, respectively, in unitsof 1010 /cm2 in the energy range of 7 to 50 MeV at the distance of
10 kpc without neutrino oscillation. We also take into account neutrino oscillations based on [14]. We assume PH = 0 in the parameterization of [14], which implies adiabatic transitions between electron and tau (anti-) neutrinos due to sin2 2θ13 = 0.095 ± 0.010 [15].
In this report, we do not take into account collective effects. We
use the cross sections [16] for the inverse beta decays, [17,18] forthe charged current interactions to oxygen, and [19] for the electron elastic scatterings.
Monte-Carlo simulation samples are generated for the three SN
models and three neutrino oscillation hypotheses by making use of
the full SK detector MC simulator based on Geant3 in the SK innerdetector volume (32.5-kton) in the energy range of 3 to 60 MeV
taking into account the trigger threshold curve. The calibration of

K. Abe et al. / Astroparticle Physics 81 (2016) 39–48

41

Table 1
Numbers of expected events at SK in the 22.5-kton ﬁducial volume with the 7 MeV total energy threshold fora SN burst with a distance of 10 kpc. We estimated these numbers using SK MC: we generate 30 0 0 ensemblesof the MC samples, reconstructed the events with the SK standard reconstruction tool, applied the selection
criteria, and then calculated the average numbers.
Wilsonν̄e + p → e+ + n
νe + e− → νe + e−ν̄e + e− → ν̄e + e−
νx + e− → νx + e−ν̄x + e− → ν̄x + e−
νe +16 O → e− + X
ν̄e +16 O → e+ + X
Total

NK1

NK2

No osc.

NH

IH

No osc.

NH

IH

No osc.

NH

IH

4923
74
25
41
34
8
64
5169

5667
130
29
33
33
662
196
6750

7587
114
37
34
29
479
531
8811

2076
43
10
17
14
22
27
2209

2399
56
12
19
14
78
48
2626

2745
56
14
18
14
74
70
2991

1878
39
9
17
13
16
20
1992

2252
54
11
17
13
72
41
2460

2652
54
13
17
14
68
64
2882

the detector and its simulation are described in [20]. We apply the
SK standard reconstruction program to the generated MC events
to obtain the vertex position, the direction and the total energyof each event. In the SN monitor, we use events with total energy
greater than 7 MeV in the 22.5-kton ﬁducial volume, where theﬁducial volume is deﬁned as the volume whose surface is located
2 m inside from the surface of the SK inner detector volume.
We generate MC samples for the three SN models for the threeneutrino oscillation hypotheses: no oscillation, normal hierarchy(NH) and inverted hierarchy (IH). Table 1 shows the expected numbers of events of the three SN models at SK in the 22.5-kton ﬁducial volume with the total energy threshold of 7 MeV, obtained byaveraging the 30 0 0 MC ensembles at the distance of 10 kpc. Fig. 1
shows the reconstructed energy distributions at SK for the Wilson
and NK1 models with and without neutrino oscillations to display
the effect of the neutrino oscillations. Fig. 1 also shows the energy
spectrum of the spallation events found in the silent warnings described in Section 2.2.

2.2. The monitor system
In this section, we will describe details of the SN neutrino burst
monitor system. Fig. 2 shows a ﬂow diagram of the system. The SKdata collected by the data acquisition system are sent to the event
builder. At the event builder, the event data are packed and stored
in a data ﬁle which we call a sub-run ﬁle. Each sub-run ﬁle contains about one minute of event data. The sub-run ﬁles are sent
to both the oﬄine process and the SN burst monitor. In the ofﬂine process, the data ﬁles are converted to an oﬄine data format
that is used for various physics analyses and detector calibrations.
The SN monitor system is running on a single computer on whicha control process operating continuously handles all the processes
and the data ﬁles.
For each sub-run ﬁle sent to the SN monitor, two processes areautomatically executed by the control process: a reformat process(the ﬁrst process) and an event reconstruction process (the second
process). The reformat process converts the online data format tothe oﬄine data format. Using the oﬄine format data, the event
reconstruction process reconstructs the vertex position, directionand energy for each event. It takes about two minutes to ﬁnish
the reformat and event reconstruction for one sub-run ﬁle. Eventswith total energy greater than 7 MeV and vertex position within
the 22.5-kton ﬁducial volume in SK are selected. We remove cosmic ray muons and their subsequent decay electron events. Afterthe reconstruction of each selected event, a 20-s time window is
opened backwards in time from the event, and the number of selected events in the window (Ncluster ) is counted. If there is a subrun ﬁle boundary, the time window extends to the previous subrun ﬁle.

We also compute a variable D that characterizes the vertex distribution. The variable D identiﬁes the dimension of the vertex
distribution and is an integer number from 0 to 3, corresponding to point-, line-, plane- and volume-like distributions, respectively. The variable D is determined by comparing χ 2 values obtained from the lengths of the major and minor axes that correspond to the eigenvalues of a correlation matrix of the vertex distribution. The correlation matrix is a 3 × 3 matrix whose elements
are deﬁned as (xi − xi  )(x j − x j  ), where i, j = 1, 2, 3 identify
the vertex position axes and x is the mean value of a variableNcluster 
x. We construct a χ 2 = n=1|dn − d(s )|2 , where dn is the nthevent vertex position and d (s ) is a position closest to dn on eithera point, a line or a plane with parameters s that determine thegeometry of the three cases. The three eigenvalues λi (i = 1, 2, 3,
and λ1 ≤ λ2 ≤ λ3 ) are used to construct the minimum χ 2 values that are (λ1 + λ2 + λ3 )/3, (λ1 + λ2 )/2 and λ1 , computed byassuming the vertex distribution is point-, line- and plane-like, respectively, with the condition of ∂ χ 2 /∂ s = 0. The comparison ofthe χ 2 values to determine a D value is tuned using MC simulations so that the calculated D value reproduces the input one.
In case of a real SN burst, the vertex distribution should be uniform in SK, and we would have D = 3, depending on the numberof burst events, which is conﬁrmed by a simulation. In contrast,
for the case of a background burst mainly originating from spallation events, the vertex positions distribute along the parent muontracks, and we would have D = 2, 1 or 0, where the spallation
events are the radioactivities created by both high-energy cosmic
ray muons and by constituents of the resulting hadronic showers.
When the process ﬁnds Ncluster ≥ 60 events and D = 3, it generates a prompt SN burst warning which initiates phone-callings andemails sent to experts in the SK collaboration within a few minutes
after the SN burst occurs. We call such a warning a “golden” warning. Subsequent to a golden warning, the experts start a meeting inorder to make a world-wide announcement within one hour. The
threshold of Ncluster is determined so that we would have 100% SN
detection eﬃciency at the Large Magellanic Cloud (LMC) assuming
the three SN models described in Section 2.1.
The third process in Fig. 2 combines all the sub-run data anddetermines the SN direction by a ﬁt. All the SN burst event information is summarized and sent to the experts by e-mail, whichis also used as the input to the discussions. Following these discussions the announcement containing the information about the
number of observed neutrinos, the burst time duration, the universal time the burst happens and the estimated direction of the SN
in the equatorial coordinate system is broadcast to the ATEL [21],GCN [22], IAU-CBAT [23] and SNEWS [8]. The universal time is determined using 1 pps (pulse per second) signals from the global
positioning system and a local time clock system consisting of acommercial rubidium clock [24]. No golden prompt warning has
been sent so far.

K. Abe et al. / Astroparticle Physics 81 (2016) 39–48

Events/MeV

103

Spallation events

a

102

103

Events/MeV

42

102

10

10

1

1

10

Events/MeV

103

15

20

25
30
35
40
45
50
Reconstructed Energy (MeV)

Spallation events

c

102

10-1

103

Events/MeV

10-1

10

1

1

15

Events/MeV

20

25
30
35
40
45
50
Reconstructed Energy (MeV)

Spallation events

e

102

10-1

15

20

25
30
35
40
45
50
Reconstructed Energy (MeV)

Spallation events

d

10

15

103

Events/MeV

10

103

20

25
30
35
40
45
50
Reconstructed Energy (MeV)

Spallation events

f

102

10

10

1

1

10-1

10

102

10

10-1

Spallation events

b

10

15

20

25
30
35
40
45
50
Reconstructed Energy (MeV)

10-1

10

15

20

25
30
35
40
45
50
Reconstructed Energy (MeV)

Fig. 1. Reconstructed energy distributions at SK for (a) the Wilson model and (b) the NK1 model at 10 kpc for the four neutrino interaction channels without neutrino
oscillation, (c) the Wilson model and (d) the NK1 model with the NH hypothesis for neutrino oscillations, and (e) the Wilson model and (f) the NK1 model with the IH
hypothesis for neutrino oscillations. The purple dotted line is the energy spectrum of the spallation candidates found in the silent warnings. The spallation event histogram
is normalized to the averaged number of accidental background events per 20 s multiplied by a factor of 104 . (For interpretation of the references to color in this ﬁgurelegend, the reader is referred to the web version of this article).

When the SN burst has less than 60 events, the golden warning will not be generated. Instead we set another threshold. The
threshold for generating the warning is determined so that thebackgrounds are suppressed: we set the threshold of Ncluster ≥ 25
and require D = 3. The warning generated with this condition iscalled a normal warning. The normal warning is sent to the experts only by e-mail, without phone-calls, and to SNEWS. Conveners among the experts check the event cluster found by this warning and make a decision about whether to have a meeting for theannouncement. The normal warning threshold is set so that wewould have 100% SN detection eﬃciency at the Small Magellanic
Cloud (SMC) assuming the three SN models. The details of the detection eﬃciency will be described in Section 2.3. The reason toprovide the normal warning is to avoid any fake warnings caused

by unexpected software and hardware troubles. We have had nonormal warning so far.
In summary, the SN monitor reformats the data and determines
the vertex position, direction, and energy of events within a few
minutes of the data being collected. It then searches for burstsof events - a cluster occurring within 20 s – with energies above
7 MeV and whose vertices fall within SK’s 22.5-kton ﬁducial volume. Within an hour of a warning being issued, the experts gather
and hold a meeting to determine the appropriate public announcement to make, if any, based on the quality and nature of the detected burst.
We also provide a lower threshold such that we require morethan 13 events in 10 s. We call a warning generated with this condition a “silent” warning. The conditions were tuned so that weK. Abe et al. / Astroparticle Physics 81 (2016) 39–48

tions, i.e., the vertex distribution is concentrated around the parent
muon track, the energies of the events have the typical spallationenergy spectrum up to 20 MeV, and the time distribution is anexponential decay consistent with the lifetimes of the spallation
products. In the ﬁgure, we generate MC simulation samples in therange of 60 ≤ Ncluster ≤ 100 uniformly, and plot the distributions
for the samples. The probability to have D ≤ 2 for SN MC with
the Ncluster range is 8 × 10−4 . No SN MC sample having D ≤ 2 isfound for 930,0 0 0 samples with 10 0 < Ncluster < 1, 0 0 0. For a normal warning condition, i.e., 25 ≤ Ncluster < 60, the probability tohave D ≤ 2 is 1.3%. This demonstrates that the variable D can discriminate between the SN-like clusters and spallation backgroundclusters.
Fig. 4 shows the SN detection eﬃciency as a function of a distance to a SN for the normal and golden warnings for the three SN
models without neutrino oscillation and with neutrino oscillations
for normal and inverted mass hierarchy hypotheses. It is found thatthe system has 100% detection eﬃciency up to the LMC locatedat 50 kpc away for all three models for the golden warning withthe three hypotheses. For the SNe at the SMC, about 64 kpc away,
the eﬃciency depends on the hypotheses for the Nakazato model,
and is 100% for the Wilson model. The normal warning has almost
100% detection eﬃciency for the three models. The eﬃciency isbasically determined by the number of inverse beta decay events.
The difference between detection eﬃciencies among the three hypotheses of the neutrino oscillations is caused by the difference
between the ν̄e energy spectra. The average energy of ν̄e is smaller
than that of ν̄x when those neutrinos are emitted from the neutrinosphere. With neutrino oscillations, ν̄x are converted to ν̄e andtherefore the average energy of ν̄e at SK increases, resulting in ahigher event rate of the inverse beta decays. Due to this effect, theFig. 2. Block diagram of the real time SN neutrino burst monitor. The arrows show
paths of the data ﬂow, and boxes indicate processes. Details are explained in the
text.

would have a few warnings per day from spallation events. The
silent warnings are sent to only a few experts of the monitor system operation and detector condition, and are not used as the fastalert for a SN burst.
2.3. Performance study of the SN burst monitor with simulations

Entries / bin

Data (Spallation events)
SN MC (Wilson)

102

103

10

1

1

3

10

3

10
100 200 300 400 500 600 700 800 900 1000
2 for point distributionc

Data (Spallation events)
SN MC (Wilson)

102bData (Spallation events)
SN MC (Wilson)

102

10

10-1
0

Entries / bin

a

10-1
0

Entries / bin

Entries / bin

Fig. 3 shows the minimum χ 2 distributions for the three geometrical assumptions and D distributions for SN MC and spallation data triggered by the silent warnings. The spallation event
cluster is identiﬁed from the vertex, energy and time distribu-

103

43

105

3

10
100 200 300 400 500 600 700 800 900 1000
2 for line distributiondData (Spallation events)
SN MC (Wilson)

4

10

103
102

10

10
1

1
10-1
10-1
0

3

10
100 200 300 400 500 600 700 800 900 1000
2 for plane distribution-2
10-0.5

0

0.5

1
1.5
2
2.5
3
3.5
Dimension of vertex distribution

Fig. 3. Minimum χ 2 value distributions assuming (a) a point distribution, (b) a line distribution and (c) a plane distribution, and (d) distribution of D (dimension of vertex
distribution) for SN MC (red broken line) and spallation data (black solid line) found as silent warnings. For SN MC clusters, we plot the distributions for the MC samples
with 60 ≤ Ncluster ≤ 100. The histogram of SN MC is normalized to the number of data entries. (For interpretation of the references to color in this ﬁgure legend, the reader
is referred to the web version of this article).

K. Abe et al. / Astroparticle Physics 81 (2016) 39–48

1

No osc.

0.9

Golden: Wilson
Normal: Wilson
Golden: NKI
Normal: NK1
Golden: NK2
Normal: NK2

0.6
0.5
0.4

LMC

0.7

SMC

0.8
Milky Way (edge)

Detection efficiency

44

0.3
0.2
0.1
50

100

150

200
250
Distance to SN (kpc)

1

NH

0.9

Golden: Wilson
Normal: Wilson
Golden: NK1
Normal: NK1
Golden: NK2
Normal: NK2

0.6
0.5
0.4

LMC

0.7

SMC

0.8
Milky Way (edge)

Detection efficiency

0
0

0.3
0.2
0.1
50

100

150

200
250
Distance to SN (kpc)

1

IH

0.9

Golden: Wilson
Normal: Wilson
Golden: NK1
Normal: NK1
Golden: NK2
Normal: NK2

0.6
0.5
0.4

SMC

0.7

LMC

0.8
Milky Way (edge)

Detection efficiency

0
0

0.3
0.2
0.1
0
0

50

100

150

200
250
Distance to SN (kpc)

Fig. 4. Detection eﬃciency of a SN burst as a function of a distance for the normal
and golden warnings with the Wilson and Nakazato models. From top to bottom,
without neutrino oscillation, with neutrino oscillations for the normal hierarchy andwith neutrino oscillations for the inverted hierarchy.

detection eﬃciency of SNe at the SMC increases for the case ofneutrino oscillations.
2.4. Operation of the SN burst monitor
We have operated the SN burst monitor system for about 20
years, since the beginning of SK data-taking in 1996. The SN burst
selection criteria and operation scheme have been changed, updated, and improved throughout this period. The SN monitor system scheme described in this report came into service in April of2013. Before that time, earlier versions of the monitor system had
been running as one of the oﬄine processes.
In Fig. 5(a), we show the silent warning rates per 24 h as afunction of the elapsed days from Jan. 1st, 2010. The rate hastrends as a consequence of – and which track – variations of thewater transparency in SK. The energy scale used in the energy reconstruction program in the SN monitor process has been continually adjusted to compensate for these transparency ﬂuctuations.
Despite these ﬂuctuations, the warning rate has been relatively stable over the last six years, with an average rate of 2.4 warnings per
24 h.
The expected number of accidental background events satisfying the event selection is 0.121 events per 10 s with a root meansquare of 0.007 events. We estimate this by counting the number of events in the SK ﬁducial volume with a total energy greaterthan 7 MeV for one day, and scale this number to a rate per 10 s.
Fig. 5(b) shows the estimated averaged background event rate fora recent period of 434 days. The background events are considered
to be spallation products, since there should be negligible contamination from known radioactivities other than spallation products
given the applied energy threshold.
Fig. 5(c) shows the data processing time distribution for thesilent warnings found. The average time to ﬁnish the processing isabout 170 s; ﬂuctuations are caused by the reconstruction process
and the condition of the network through which the data sample
ﬁles are copied from the SK data acquisition system. The oﬄine SN
monitor that had run before April 2013 took about ﬁve minutes toﬁnish the reformat and reconstruction processes, as the oﬄine reconstruction program was tuned for physics analysis and calibration. We have optimized the reconstruction program for the online
SN monitor to increase the processing speed without degrading its
performance.
Fig. 5(d) shows the averaged monthly duty cycle of the SN monitor system over a recent 34-month period; the SN monitor operates with a duty cycle of about 97%. The monitor searches for
SN event bursts during normal SK running, but it does not operate during SK detector calibration runs, particularly during thosecalibration runs employing artiﬁcial sources that intentionally generate event bursts. Most of the 3% loss of the SN monitor duty cycle comes from planned calibration. Note that even when the SN
monitor is off, SK still has a non-realtime capability to detect a SN
burst during these calibration runs. This is achieved via dedicatedoﬄine analyses which remove likely source events based on theirvertex positions and event timings.
We use the spallation events found as silent warnings to estimate the false alarm rate by assuming constant rate Poissonprocesses rather than generating simulation samples of the backgrounds. We combine multiple silent warnings randomly and form
a combined cluster to estimate a probability of having a golden(normal) warning. Using 2551 silent warnings, we combined twoof them for all combinations of two spallation bursts, to form
2551C2 = 3, 252, 525 patterns, and estimate the probability for the
combined burst to pass the criteria for a golden (normal) warning to be 0 (4.3 × 10−6 ). For three combinations with 2551 C3 patterns, the probability is estimated to be 4.1 × 10−5 (1.9 × 10−4 )for a golden (normal) warning. One silent warning happens every 10 hours. The probability to have two (three) spallation clusters coincident within 20 s is 5.6 × 10−4 (3.1 × 10−7 ). Therefore
the probability to have a golden (normal) warning due to accidentally coincident spallation bursts is 3.1 × 10−7 × 4.1 × 10−5 =1.3 × 10−11 (5.6 × 10−4 × 4.3 × 10−6 = 2.4 × 10−9 ). The false alarm
rates are calculated to be once per 9.0 × 107 years for a golden
warning and once per 4.7 × 105 years for a normal warning.
The processes in the SN burst monitor are kept under observation by a web-based monitor running on a dedicated PC. Any
problems are displayed on the web monitor immediately after theyK. Abe et al. / Astroparticle Physics 81 (2016) 39–48

a

b

c

d

45

Fig. 5. Stability of the SK operation: (a) the silent warning rate per 24 h as a function of elapsed days from Jan. 1st of 2010, (b) averaged number of accidental background
events per 10 s with one histogram entry for one day, only 0.02% of the 10 s windows have more than 13 as shown in the inset histogram, (c) process time to ﬁnish thereformat and reconstruction to issue the silent warnings found, and (d) SN monitor duty cycle averaged over one month for each entry.

are found, and the SK shift takers are notiﬁed by visible and audio
alerts.
3. Determination of the SN direction
The determination of the SN direction is crucial since the direction information provided by neutrinos is useful for astronomers to
observe the SN explosion process from the onset via electromagnetic waves. At present, SK is the only operating experiment with
suﬃcient detector mass to determine the neutrino direction fromelastic scattering events which preserve the SN direction. Though
the inverse beta decay events also have a correlation with the
SN direction, the elastic scattering events mainly dominate the direction determination power. A study of the determination of the
SN direction using neutrinos was performed in [18,25]. Here, we
present a method we have developed. We will explain the algorithm to determine the SN direction, and then will show its performance obtained using SK MC.
3.1. Algorithm
We determine the SN direction based on a maximum likelihood
method. A likelihood function Li for i-th event is deﬁned as:

Li =


r

Nrk pr (Ei , dˆi ; dˆSN ),(1)

where the index r indicates one of the four neutrino interaction
channels: inverse beta decay (ν̄e p), electron elastic scattering ofanti-electron neutrino (ν̄e e), other elastic scatterings (ν e) and thecharged-current interactions on oxygen (ν 16 O). The index k indicates the energy bin, running from 1 to 5 for the energy rangesof 7 < E < 10, 10 < E < 15, 15 < E < 22, 22 < E < 35 and35 < E < 50, respectively, where E is the measured total electron energy in MeV. Nrk is the number of events of the interaction r in the k-th energy bin. Ei is the i-th event total electron energy, which uniquely determines the index k, dˆi is thei-th event direction and dˆSN is the SN direction we want to determine. The pr (Ei , dˆi ; dˆSN ) function is a probability density function (PDF) for interaction r as a function of the energy Ei and an
inner-product of dˆi · dˆSN = cos θSN . The PDF is determined using SK
MC. The number of ν̄e elastic scattering events can be inferredfrom the number of inverse beta decay events with the relation
Nν̄e e,k = m Akm Nν̄e p,m , where the matrix Akm is calculated from aratio of the total cross sections of the two interactions. We determine PDFs for elastic scatterings with the following procedure.
We divide the SN MC sample elastic scattering events generated
with the Wilson model into one-MeV bins from 7 to 35 MeV. Forenergies greater than 35 MeV, we combined all events into one
bin. Then we ﬁt the cos θ SN distribution with the known SN direction in MC using a model function that is the superpositionof four exponential functions and containing eight parameters foreach energy bin. For a given energy value, we compute the eightK. Abe et al. / Astroparticle Physics 81 (2016) 39–48

100
90
80
70
60
50
40
30
20
10
0
-1
20
18
16
14
12
10
8
6
4
2
0
-1

1

100
90
80
70
60
50
40
30
20
10
0
-1

Entries/bin
1

100
90
80
70
60
50
40
30
20
10
0
-1

a

-0.8 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

cos SN

Entries/bin

50
45
40
35
30
25
20
15
10
5
0
-1

c

-0.8 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

cos SN

350

Entries/bin

Entries/bin

Entries/bin

Entries/bin

46

e

300
250

b

-0.8 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

-0.4

-0.2

0

0.2

0.4

0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

0.8

1

0.8

1

cos SNd-0.8 -0.6

cos SN

f

200
150
100
50

-0.8 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

cos SN0-1-0.8 -0.6

cos SN

Fig. 6. Distributions of cos θ SN with the Wilson model at 10 kpc for the ﬁve energy bins: (a) 7 < E < 10 MeV, (b) 10 < E < 15 MeV, (c) 15 < E < 22 MeV, (d) 22 <E < 35 MeV and (e) 35 < E < 50 MeV, where E is the measured total electron energy, and (f) all energies combined. The superimposed solid lines are the ﬁtted likelihood
functions.

parameter values by interpolating the parameter values of neighboring two energy bins and applying those to the model function
to obtain the PDF value. A similar procedure is applied to the PDFs
for inverse beta decays and interactions on oxygen to determine
the PDF values. We construct a likelihood L = exp( k,r Nrk ) i Li ,and maximize it so that:

∂L
∂L
=
= 0,
∂ Nrk ∂ dˆSN(2)where for Nrk we vary r = ν̄e p, ν e and ν 16 O. For r = ν 16 O, we assume the cos θ SN is same for neutrino and anti-neutrino interactions. We set Nrk = 0 for r = ν 16 O with k = 1, 2, 3, as the expected
number of charged current interactions on oxygen is negligible inthose energy ranges, as shown in Fig. 1. The SN direction dˆSN contains two parameters: zenith and azimuth angles, that are translated to the direction in the equatorial coordinate system with thetime the burst is found. Hence, we vary 14 parameters of Nrk anddˆSN in total.
When we perform a ﬁt with the likelihood method, we ﬁrst determine the initial value of the direction based on a grid search:
we scan all dˆSN to the 4π directions with a coarse grid step andcount the number of events that satisfy cos θ SN > 0.8 at each step,
and we set the initial value that gives the largest number of events.
3.2. PerformanceFig. 6 demonstrates cos θ SN distributions of a ﬁt to a MC sample of the Wilson model at 10 kpc for the ﬁve energy bins andcombined one with the superpositions of the ﬁtted likelihood functions. Fig. 7 shows the corresponding direction distribution on asky map in the equatorial system. The red (blue) points are the reconstructed directions of each elastic scattering event (inverse beta
decay or charged current reaction on oxygen), and the star markshows the reconstructed SN direction. The elastic scattering events
concentrate around the reconstructed SN direction, while the distribution of inverse beta decays and charged current events is almost uniform across the entire sky.

Fig. 7. Reconstructed direction on a sky map in the equatorial system obtained bya MC simulation with the Wilson model at 10 kpc. Red points are the directions ofelastic scattering events, blue points are event directions of inverse beta decay andcharged currents on oxygen, and the star point is the reconstructed SN direction.
(For interpretation of the references to color in this ﬁgure legend, the reader isreferred to the web version of this article).

Fig. 8 shows θ distributions of the three models (Wilson, NK1
and NK2) for 30 0 0 MC samples at 10 kpc without neutrino oscillation, where θ is the angle between the input SN direction andthe ﬁtted direction. The solid lines are ﬁt results using the vonMises–Fisher (MF) distribution [26]:f (θ ; κ ) =κ

2 sinh κeκ cos θ sin θ ,(3)

where κ determines the sharpness of the distribution concentration on a sphere.
We estimate the angular resolutions of the SN direction determination using an ensemble estimation. In order to cope with any
possible combinations of the elastic scatterings and inverse beta
decays, we employ the following method. We generate a numberof MC samples for various combinations of ﬁtted yields of the elastic scatterings and inverse beta decays in the ranges up to 1500 for
the former and 60,0 0 0 for the latter. We divide each range into 15

400
350

Wilson
NK1
NK2

300
250
200
150
100
50
0
0

2

4

6

8

10

12

14

16

18
20
( deg.)

Fig. 8. Distributions of θ for three models at a distance of 10 kpc without neutrino oscillation. The superimposed solid lines are the ﬁtted MF functions.

30

a
25

Wilson model w/o

osc.

Ensemble estimation
Likelihood curvature
Grid search

20
15
10

at 68% coverage (deg.)

σ=

5

c
20

10

15

where the second derivative is the curvature calculated at the
point on the sphere at which the likelihood function L becomesmaximal. We calculate σ values along four planes that include themaximal point and have different azimuth angles with respect tothe point, and employ the maximum among the four σ values.
Then we obtain the resolution θ σ (q) that covers an area with afraction of q = 1 − p of the MF distribution with a p Value:θσ (q ) = arccos

1
κ



ln(1 − q + qe−2κ ) + 1 ,(5)and κ = 1/σ 2 .

30

b
25

Wilson model
Without osc.
With osc. (Inverted Hierarchy)
With osc. (Normal Hierarchy)

20
15
10

NK1 model w/o

0
0

20
25
30
Distance to SN (kpc)

osc.

Ensemble estimation
Likelihood curvature
Grid search

15
10

5

10

15

20
25
30
Distance to SN (kpc)

30

d
25
20

NK1 model
Without osc.
With osc. (Inverted Hierarchy)
With osc. (Normal Hierarchy)

15
10
5

5
0
0

(4)

5

30
25

1
− ∂ 2 ln L ,
2
∂θSN

5
0
0

en or



en at 68% coverage (deg.)en orat 68% coverage (deg.)

to obtain a 15 × 15 matrix. Each matrix element contains about
30 0 0 MC samples. For each element, we determine the angle θ enthat covers 68.2%, 90% and 95% of the MC samples. We also provide probabilities to have the true SN position in 2°, 5° and 10°
with respect to the ﬁtted direction. We generate the matrices for
the Wilson and NK1 models with (NH) and without neutrino oscillations. The dependence of the θ en on the models has about a 10%
variation. We employ the largest values of θ en and the smallest
value of the probabilities for each matrix element conservatively.
When we ﬁnd a SN neutrino burst, we apply the ﬁt to the burst
events to obtain the SN direction and the yields of the elastic scatterings and inverse beta decays. With the ﬁtted yields, we iden-

47

tify the matrix element and obtain the angular resolutions and theprobabilities that are to be announced to the public. For example,
we ﬁnd θen = 3.1 − 3.8◦ (4.3 − 5.9◦ ) at 68.2% coverage for the Wilson (NK1) model at 10 kpc, where the range covers various neutrino oscillation scenarios in Table 1.
The Nakazato model also provides a SN model with a black holeformation for M = 30 solar mass [13] . The model predicts thatneutrino emission suddenly stops 842 ms after the core bounce.
If the SN monitor observes an abrupt cutoff of the supernova neutrino ﬂux, this could be the signature of the birth of a black hole.
Based on this model, we generate MC samples and apply the ﬁt.
We obtain an angular resolution of 2.3° at 10 kpc with 68.2% conﬁdence level. Therefore this may help the identiﬁcation of a positionof a disappeared massive star as proposed by [27].
In order to understand the behavior of the estimated angularresolution, we make use of the curvature of the likelihood at its
maximal position: we deﬁne a value σ as:

en at 68% coverage (deg.)

Entries/bin

K. Abe et al. / Astroparticle Physics 81 (2016) 39–48

5

10

15

20
25
30
Distance to SN (kpc)

0
0

5

10

15

20
25
30
Distance to SN (kpc)

Fig. 9. θ that covers 68% of SNe (the area) as a function of distance. In (a) and (c), the solid, dashed and dotted lines correspond to θ en obtained by the ensemble
estimation, θ σ (0.682) obtained by the likelihood curvature with Eq. (5), and that obtained by the grid search, respectively, for the Wilson (a) and NK1 (c) models. Figure (b)and (d) show θ en as a function of distance for the three neutrino oscillation hypotheses for Wilson (b) and NK1 (d) models.

48

K. Abe et al. / Astroparticle Physics 81 (2016) 39–48

In Fig. 9(a) and (c), we show the obtained angular resolution at
the 68.2% coverage for the ensemble estimation of the ﬁt, the likelihood curvature method and the ensemble estimation using a gridsearch as a function of the distance for the Wilson and NK1 models. We ﬁnd that the θ distributions are well modeled by the MF
function and the likelihood curvature estimation θσ (q )(q = 0.682 )is consistent with that of the ensemble estimation for the statistics
larger than those of SNe at 10 kpc (7 kpc) for the Wilson (NK1)model. However, θ σ (q) is found to be an underestimate for thesmaller statistics. The degradation of the angular resolution for the
distant SNe is likely due to failure in giving a proper initial value ofthe SN direction by the grid search, indicated by the fact that theθ en value approaches that of the grid search for larger distances(smaller statistics). The small statistics produce a large ﬂuctuationthat sometimes makes a fake peak on the grid direction search.
That makes the initial value a wrong direction and the likelihood
ﬁt ﬁnds a local minimum around the direction. Fig. 9(b) and (d)show the ensemble estimation of the angular resolution as a function of the SN distance for three neutrino oscillation hypotheses.
The angular resolutions with the two neutrino oscillation hypotheses are smaller than those without neutrino oscillation. This is due
to an increase in elastic scattering events as shown in Table 1.
We estimate the precision of the angular uncertainty using
Eq. (5) under the assumptions of the SN models used; the precision
 of θ σ is determined using θσ = dθσ /dq × q = dθσ /dq ×
q/N, where N = 30 0 0 is the number of samples in the ensemble. We ﬁnd θσ = 0.09° with q = 0.682 for the NK1 model with
the NH neutrino oscillation at a distance of 10 kpc, which is muchsmaller than the angular resolution variation of various neutrino
oscillation scenarios.
4. Summary
We describe a real-time monitor of a SN neutrino burst at SK.
The monitor is able to provide a fast warning to the world within
one hour. The system is operating on a dedicated computer independent of the oﬄine processes. The SN neutrino burst selection
criteria are determined so that fake event bursts mainly causedby spallation events are rejected. Using MC simulations, we ﬁndthat the system has 100% detection eﬃciency up to the LMA for
the three SN models with the golden warning criteria. The expected total number of neutrino events with the selection criteria
is about 5200 (2200) for the Wilson (NK1) model at 10 kpc without neutrino oscillation. Neutrino oscillations increase the number of inverse beta decay events and enhance the detection eﬃciency for SNe at the SMC for both normal and inverted hierarchy
hypotheses.
The SN direction pointing is of importance as it provides achance for astronomers to observe the SN explosion from its onsetwith electromagnetic waves. SK is the only detector that enables us
to reconstruct the SN direction using only neutrinos among the existing neutrino observatories. We have developed an algorithm toidentify the SN direction and its error using a maximum likelihood
method. The pointing accuracy estimated by the ensemble study
is found to be 3.1–3.8◦ (4.3–5.9◦ ) at 68.2% coverage for the Wilson(NK1) model at 10 kpc, where the range covers various neutrino
oscillation scenarios.

Acknowledgments
We gratefully acknowledge the cooperation of the Kamioka
Mining and Smelting Company. The Super-Kamiokande experimenthas been built and operated from funding by the Japanese Ministry of Education, Culture, Sports, Science and Technology, theU.S. Department of Energy, and the U.S. National Science Foundation. Some of us have been supported by funds from the Research
Foundation of Korea (BK21 and KNRC), the Korean Ministry of Science and Technology, the National Research Foundation of Korea(NRF- 20110 0240 09), the European Union (H2020 RISE-GA641540SKPLUS), the Japan Society for the Promotion of Science, theNational Natural Science Foundation of China under Grant no.
11235006, the Natural Sciences and Engineering Research Council(NSERC) of Canada, and the Scinet and Westgrid consortia of Compute Canada. This work was partly supported by the Grant-in-Aid
for Scientiﬁc Research on Innovative Areas [JSPS No. 26104006].
References[1] K. Hirata, et al., Phys. Rev. Lett. 58 (1987) 1490; R.M. Bionta, et al., Phys.
Rev. Lett. 58 (1987) 1494; E.N. Alekseev, L.N. Alekseeva, V.I. Volchenko,
I.V. Krivosheina, JETP Lett. 45 (1987) 589; M. Aglietta, et al., Europhys. Lett.
3 (1987) 1315.
[2] In this paper, according to particle-physics nomenclature, “event” refers to aneutrino interaction recorded by Super-Kamiokande.
[3] D.N. Schramm, J.W. Truran, Phys. Rep. 189 (1990) 89; M. Koshiba, Phys. Rep.
220 (1992) 229.
[4] S.M. Adams, et al., Astrophys. J. 778 (2013) 164.
[5] T. Shigeyama, K. Nomoto, M. Hashimoto, D. Sugimoto, Nature 328 (1987) 320;
S. Blinnikov, et al., Astrophys. J. 532 (20 0 0) 1132.
[6] R.C. Svoboda, NASA-CR-194367 (1992); Y. Oyama, et al., Nucl. Instr. Meth. A
340 (1994) 612; M. Ambrosio, et al., Astropart. Phys. 8 (1998) 123.
[7] E.N. Alexeyev, L.N. Alexeyeva, J. Exp. and Theo. Phys. 95 (2002) 5; M.E. Monzani, et al., Nuovo Cim. C 29 (2006) 269; N. Yu Agafonova, et al., Astropart.
Phys. 28 (2008) 516; A.A. Aguilar-Arevalo, et al., Phys. Rev. D 81 (2010) 032001;
R. Abbasi, et al., A&A A535 (2001) 109; H. Wei for the Daya Bay Collaboration,
Phys. Proc. 61 (2015) 802; K. Asakura, et al., Astrophys. J. 818 (2016) 91; K. Zuber, Nucl. Part. Phys. Proc. 265-266 (2015) 233.
[8] P. Antonioli, et al., New J. Phys. 6 (2004) 114; K. Scholberg, Astron. Nachr. 329(2008) 337.
[9] S. Fukuda, et al., Nucl. Instrum. Methods A 501 (2003) 418.
[10] H. Nishino, et al., Nucl. Instrum. Methods A 610 (2009) 710.
[11] S.W. Li, J.F. Beacom, Phys. Rev. C 89 (2014) 045801.
[12] T. Totani, et al., Astrophys. J. 496 (1998) 216.
[13] K. Nakazato, K. Sumiyoshi, H. Suzuki, T. Totani, H. Umeda, S. Yamada, Astrophys. J. Suppl. 205 (2013) 2.
[14] A.S. Dighe, A.Y. Smirnov, Phys. Rev. D 62 (20 0 0) 0330 07.
[15] K.A. Olive, et al. (Particle Data Group), Chin. Phys. C 38 (2014) 090 0 01; http://pdg.lbl.gov.
[16] P. Vogel, J.F. Beacom, Phys. Rev. D 60 (1999) 053003.
[17] E. Kolbe, K. Langanke, P. Vogel, Phys. Rev. D 66 (2002) 013007.
[18] R. Tomas, D. Semikoz, G.G. Raffelt, M. Kachelrieß, A.S. Dighe, Phys. Rev. D 68
(2003) 093013.
[19] J.N. Bahcall, M. Kamionkowski, A. Sirlin, Phys. Rev. D 51 (1995) 6146.
[20] K. Abe, et al., Nucl. Instrum. Methods A737C (2014) 253.
[21] R.E. Rutledge, Astron. Soc. Pac. 110 (1998) 754; The Astronomer’s Telegram,
http://www.astronomerstelegram.org/.
[22] S.D. Barthelmy, T.L. Cline, P. Butterworth, AIP Conf. Proc. 587 (2001) 213; The
Gamma-ray Coordinates Network, http://gcn.gsfc.nasa.gov.
[23] IAU: Central Bureau for Astronomical Telegrams, http://www.cbat.eps.harvard.
edu/.
[24] K. Abe, et al., Nucl. Instr. Meth. A 659 (2011) 106.
[25] J.F. Beacom, P. Vogel, Phys. Rev. D 60 (1999) 033007; S. Ando, K. Sato, Prog.
Theo. Phys. 107 (2002) 957.
[26] R.A. Fisher, 1953 Proc. R. Soc. Lond. Ser. A, 217 295–305; K.V. Mardia, P.E. Jupp,
Directional Statistics, Wiley, 1972.
[27] C.S. Kochanek, et al., Astrophys. J. 684 (2008) 1336.