The University of Tokyo
Graduate School of Science
Department of Physics

Study of Neutron-Tagging Algorithm for SK-Gd(SK-GdkºQã-'P⌧˙¢Î¥Í∫‡nv)

Seungho Han

Supervisor: Prof. Kimihiro Okumura
Submitted in part fulfillment of the requirements for the degree ofMaster of Science in Physics of the University of Tokyo, January 2021

Abstract
Since July 2020, the SK-Gd phase of the Super-Kamiokande (SK) experiment has started bydissolving 0.0110 wt.% of Gadolinium (Gd) in the pure-water detector volume. The large
cross section (48,700 barns) and de-excitation energy (⇠8 MeV) of the Gd(n, ) reaction are
believed to greatly enhance the neutron detection efficiency. If the Gd concentration reaches
the target value of 0.1 wt.%, the efficiency is expected to improve up to ⇠90 %, while thesame efficiency for the pure-water detector volume has been only ⇠25 %. With the enhanced
detection efficiency of neutrino-induced neutrons, we hope to boost sensitivity in the search forSupernova Relic Neutrinos (SRN), neutrino oscillation parameters, and proton decays.
To evaluate the neutron detection efficiency with the current Gd concentration of 0.0110 wt.%,
we have inherited and tested the so-called “neutron-tagging” algorithm that has been usedwithin SK to detect p(n, ) signals in pure-water data. The algorithm searches for a clusterof hits on photomultiplier tubes (PMTs) that are localized in time, and applies a pre-trained
Multilayered Perceptron (MLP) to classify the found clusters into signal and noise. The author has parameterized the factors that a↵ect the algorithm’s detection efficiency, and built asoftware that is able to transparently tune each parameter and create/train a custom MLP. Using the software and the SK-Gd simulator under development, the parameters were optimized
to improve detection efficiency on both simulated p(n, ) and Gd(n, ) events. In a realistic
search setting, the neutron capture detection efficiency of the optimized algorithm to simulated
neutron events was ⇠92 % for Gd(n, ), ⇠26 % for p(n, ), and ⇠61 % overall.
The tuned algorithm was applied to AmBe neutron source data taken after 0.0110 wt.% of
Gd was fully loaded. For comparison, the same algorithm was applied to a simplified neutron
source simulation without source geometry and detector triggers. (The full simulation is underdevelopment.) The simplified simulation has proved to be e↵ective in reproducing the feature
distributions of the neutron source data, although the simulation physics model appears to havea larger Gd(n, )-to-p(n, ) ratio. As a consequence, the algorithm shows a higher detection
efficiency on the simulation (⇠60 %) than on the neutron source data (⇠50 %). While wesuspect the neutron transport/capture physics model used in the simulation to be the major
origin of the apparent data-simulation discrepancy, we suggest that the lack of source geometry
and trigger in simulation can bias the simulated neutron detection efficiency.

iAcknowledgements
I am truly grateful to my thesis advisor Prof. Kimihiro Okumura, for all the research opportunities he has provided me with and his insightful guidance that has kept me focused andmotivated whenever needed, especially during the rather unfortunate year of 2020 possessed bya hopeless series of self-quarantine.
Dr. Takatomi Yano has led the AmBe neutron calibration campaign. The Automated Source
Deployment System has made the calibration process so much easier, and a large amount ofdata could be taken regularly without much of a hassle, all thanks to the system that he hasdeveloped. I also thank him for convening the AmBe calibration meeting, and giving me tonsof advice and insights that have helped me progress.
Prof. Yoshinari Hayato has guided me both as a mentor and an expert. I have learned muchabout the detector, simulation, and analysis software from the Okayama group, namely, Prof.
Yusuke Koshio, Dr. Kaito Hagiwara, and Masayuki Harada. Masayuki Harada has been themost inspiring colleauge with whom I discussed most of the topics covered in this thesis.
The SK simulation tool skdetsim has been developed by the SK collaborators for more than
20 years. Prof. Kimihiro Okumura, Prof. Yusuke Koshio, and Dr. Patrick de Perio has helped
me learn how to use it. skdetsim-skgd, the Geant4 neutron capture model attachment toskdetsim, has been developed by Dr. Tsui Ka Ming, Dr. Yosuke Ashida, Dr. Ryosuke Akutsu,
and Masayuki Harada. SKG4, the Geant4-based SK simulator that I used to estimate theamount of neutron-induced BGO scintillation within the triggered data is a work of Dr. Kaito
Hagiwara, Masayuki Harada, Dr. Guillaume Pronost, Seiya Sakai, Shotaro Yamamoto, Prof.
Yusuke Koshio, and myself. The neutron tagging software that I have developed for this study
is based on the toy tagging program made by Dr. Tsui Ka Ming.
Most of the analysis presented in this thesis is either preceded (for the pure-water phase) or
inspired by the discussions given in the works of Dr. Tristan James Irvine [1], Dr. Tsui Ka Ming[2], Toshiki Mochizuki [3], and Dr. Ryosuke Akutsu [4].
I would like to thank Prof. Hiroyuki Sekiya, Dr. Hiroshi Ito, Dr. Yasuhiro Nakajima, Dr.
Kiseki Nakamura, Dr. Linyan Wan, Akira Takenaka, Shintaro Miki, and Shuehei Watanabe fordiscussion and comments relevant to this work.
ii

iii
I would like to express my gratitude to the Gd-loading team/shifts and all SK collaborators,
for working hard day and night to successfully load Gd and keep the detector running. I deeply
appreciate the Japanese Ministry of Education and Culture, Sports, Science and Technology
for financial support, without which I could not have been a part of this brilliant collaboration.
Special thanks goes to the fellow members of the ICRR neutrino group and the Kamioka
Observatory, without whom my student life would have been miserable: Dr. Guillermo Daniel
Megı́as Vázquez, Dr. Ryosuke Akutsu, Junjie Xia, Dr. Hiroshi Ito, Dr. Kiseki Nakamura, Yutaro
Sonoda, Akira Takenaka, Kohei Okamoto, Yu-chin Chen, Shuichi Imaizumi, Takeshi Okada,
Yoshiki Nagao, Minato Shibata, Yuki Kanemura, Shintaro Miki, and Shuhei Watanabe.
Final thanks goes to my family, friends, and everyone who has supported me along the way.

Contents

Abstracti

Acknowledgementsii

1 Introduction2

1.1

Why neutrons? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.2

Neutron production in SK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.3

Neutron reactions in Gd-loaded water . . . . . . . . . . . . . . . . . . . . . . . .

8

1.4

Outline of this study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2 The Super-Kamiokande detector
2.1

2.2

14

Specification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.1.1

ID PMTs and electronics . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.1.2

Water purification

2.1.3

Software trigger . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

SK-Gd . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
iv

CONTENTS

v

3 The neutron detection algorithm

24

3.1

Signal properties of physics simulation models . . . . . . . . . . . . . . . . . . . 26

3.2

Stage 1: Hit-cluster search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.3Stage 2: Binary classification with MLP . . . . . . . . . . . . . . . . . . . . . . 38

3.4

3.3.1

Features (input variables to MLP) . . . . . . . . . . . . . . . . . . . . . . 39

3.3.2

Training the MLP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

4 Testing with a neutron source
4.1

49

Neutron source data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.1.1

The source: AmBe + BGO . . . . . . . . . . . . . . . . . . . . . . . . . 50

4.1.2

Data specification and reduction . . . . . . . . . . . . . . . . . . . . . . . 56

4.2

Neutron source simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

4.3

Data and MC comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

4.4

4.3.1

Features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

4.3.2

Signal likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.3.3

Reconstructed capture time . . . . . . . . . . . . . . . . . . . . . . . . . 72

4.3.4

Estimated neutron detection efficiencies . . . . . . . . . . . . . . . . . . . 74

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.4.1

Sources of systematic uncertainties . . . . . . . . . . . . . . . . . . . . . 75

4.4.2

The data-MC discrepancy . . . . . . . . . . . . . . . . . . . . . . . . . . 77

4.4.3

Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

vi

CONTENTS

5 Summary and prospects

81

Appendices

83

A Neutrino oscillation

83

B Neutron production mechanism of AmBe neutron source

88

C Feature distributions of neutron source data and MC for nonzero heights

90

D Hypothesis test on NHits

95

Bibliography99

List of Tables

2.1

Summary of the software trigger types. . . . . . . . . . . . . . . . . . . . . . . . 21

3.1

The number density n, the thermal neutron absorption cross section

abs , and

their product n abs for each of the top three neutron-absorbing elements: Gd, H,
and O. The number densities are for 0.0110 wt.% Gd-loaded water and the crosssections are taken from the ENDF/B-VII.1 data [20]. The n abs is equivalent tothe capture rate per neutron path length. . . . . . . . . . . . . . . . . . . . . . . 24
3.2

The summary of the performance test of di↵erent NHits thresholds. All test
sets 1-4 share the same MC configuration as explained in Section 3.1 and 3.2.
Test sets 1-3 consist of 3,000 neutron MC events each. The Stage 2 selection
efficiencies are obtained with Keras [29] using the same MLP architecture as
explained in Section 3.3 for reference. The signal likelihood thresholds (Thr.)and the combined selection efficiencies are the ones that give 2% noise ratio
within the final selected candidates. Test set 4 with 300,000 neutron MC eventsis included for comparison with test set 3. All parameters other than the NHits
threshold are as in Table 3.3.

3.3

. . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

The list of parameter values used for Stage 1 (NHits -based candidate search) of
the neutron detection algorithm. For reference, the primary particle generation
time is 1 µs in the global simulation time frame, and the time range of each data
event is [-5, -535] µs. The residual hit time range of [1, 535] µs corresponds tothe full time range following neutron emission. . . . . . . . . . . . . . . . . . . . 37
vii

viii

LIST OF TABLES
3.4

TMVA’s MLP options used for training. . . . . . . . . . . . . . . . . . . . . . . 46

3.5

Summary of the capture selection efficiencies ✏(n, ) of the neutron detection algorithm used in this study. rNoise indicates the noise ratio within the tagged
candidates.

4.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

The estimated gamma-correlated and uncorrelated neutron intensity of the source.
The assumed uncertainty is ⇠22% for the estimation in this study, and 20% for
Watanabe et al. [38]. The 4.4 MeV gamma-ray intensity should be the same as
the correlated neutron intensity. . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

4.2

Relevant properties of BGO at room temperature. [43] . . . . . . . . . . . . . . 53

4.3

The input t and the fit parameters used in Equation 4.1. . . . . . . . . . . . . . 55

4.4

Neutron source data runs used in the analysis. . . . . . . . . . . . . . . . . . . . 57

4.5

Summary of the event cuts applied to data events. The removed ratios (among
the events passing the previous cut) are obtained for Run 85386. . . . . . . . . . 62

4.6

Summary of the candidate cut applied to data events. . . . . . . . . . . . . . . . 62

4.7

Summary of MC neutron capture detection efficiencies, for each source z position.
Most of the di↵erence in the noise rate rNoise between source and sourcelessdummy MCs should be random uncorrelated neutrons from the source. . . . . . 71

4.8

The estimated neutron detection efficiencies with statistical uncertainties. (Dataaveraged from August 12th to December 1st , 2020.) . . . . . . . . . . . . . . . . 75

4.9

Summary of the possible sources of systematic errors and their rough fractional
uncertainty estimates. All errors are assumed to be normal and independent. . . 77

4.10 The measured neutron detection efficiencies and the rough systematic estimates.
The statistical uncertainties are neglected. The tuned MC efficiencies are takenfrom the MC with its Gd cocentration lowered to ⇠0.0103 wt.% to match thefitted neutron capture time constant to the one measured in data (⇠116 µs). . . 79

LIST OF TABLES

ix

A.1 Recent (circa 2017) best-fit values, 3 ranges, and errors for the oscillation parameters, assuming normal mass hierarchy (m1 < m2 < m3 ). The unit ofis 10 3 eV2 , and the unit ofm232

m221 is 10 5 eV2 [50] . . . . . . . . . . . . . . . . . 86

List of Figures

1.1(a) number flux of SRNs compared with other background neutrinos, and (b)event rate of SRNs and invisible muon decay products. Reprinted from Ref. [9] .

1.2

6

Simulated ⌫/¯⌫ CCQE neutron multiplicity distributions from primary ⌫N interaction (left), from ⌫N interaction including FSIs (center), and for ⌫N interaction
including both FSIs and SIs (right). Reprinted from Ref. [4]. . . . . . . . . . . .

1.3

8Neutron cross sections on 1 H, 16 O, and 155 Gd from the JENDL-4.0 library [12].
Neutrons are more likely to be captured with lower kinetic energy, qualitatively as
they tend to spend more time in the vicinity of nuclei with lower speed. The large
fluctuation in the intermediate energy range can be understood as manifestationof resonances in the distinguishable discrete energy levels. A formal approach
can be found in Ref. [10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.4

Neutron reaction cross sections on 16 O and 155 Gd, from the JENDL-4.0 library[12]. Note the small fraction of channels either removing or multiplying thenumber of neutrons. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.5

Capture ratio on Gd, by Gd mass fraction in water. Reprinted from Ref. [13]. . 12

1.6

The de-excitation schematic of the excited state 157 Gd + n ! 158 Gd to the⇤158

Gd ground state. Multiple gamma-rays from the 157 Gd(n, )158 Gd reaction

can share 7.937 MeV. Reprinted from Ref. [14]. . . . . . . . . . . . . . . . . . . 12x

LIST OF FIGURES
2.1

xi

Schematics of Cherenkov radiation (left) and how it forms a ring-shaped image
on a photodetector plane (right). Reprinted from Ref. [16]. . . . . . . . . . . . . 14

2.2

⌫e /⌫µ event identification amongst beam (T2K) neutrino data taken in SK, by
the clarity of Cherenkov rings. Reprinted from Ref. [18]. . . . . . . . . . . . . . 16

2.3

The SK detector. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.4

The schematic of the SK detector. . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.5

Sample event display of 0.6 GeV muon-neutrino. . . . . . . . . . . . . . . . . . . 19

2.6

ID PMT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.7

The schematic of gates issued after a PMT hit. The charging and discharging
gates opens for 400 ns and 350 ns, respectively. All signal currents are ignored
while one of the discharging gate and VETO is open. The total amount of timerequired to process a single input current is 900 ns. The timing of each hit isdetermined by the comparator threshold. . . . . . . . . . . . . . . . . . . . . . . 20

2.8

The Gd2 (SO4 )3 · 8 H2 O powder (left) and the amount of the Gd-sulfate loaded
into the detector by date (right). Reprinted from Ref. [19]. . . . . . . . . . . . . 22

2.9

A block diagram of the Gd-dissolving system. Reprinted from Ref. [19]. . . . . . 22

2.10 Measured Gd concentration by height and date. Reprinted from Ref. [19]. . . . . 23

3.1

The number of hits (left) and the raw PMT hit time distributions (right) obtained
from simulations for p(n, ) and Gd(n, ) reactions. (random neutron vertex,
area-normalized) The simulation configuration is as explained in Section 3.1. . . 25

3.2

Sample MC event displays for a p(n, ) (left) and a Gd(n, ) event (right). The
top figures are the 3D drawings of the inner tank dimensions with all neutron
signal hits (red dots) and gamma-rays (red arrows) in each event. The bottom
histograms are the raw PMT hit time distributions of each event with bin width
being 1 µs and 10 ns. The simulation configuration is as explained in Section 3.1. 27

xii

LIST OF FIGURES
3.3

The gamma-ray multiplicity (left) and the energy sum of all emitted gammarays (right) obtained from simulation with the GCALOR and ANNRI-Gd model.
(area-normalized) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3.4

The energy distribution of the individual gamma-rays obtained from simulation
with the GCALOR and ANNRI-Gd model (left) and the same for Gd(n, ) reaction with four di↵erent physics models available in the simulator, the default
Geant4 [26], GLG4Sim [27], GGarnet [24], and ANNRI-Gd [24]. (area-normalized) 28

3.5

The distribution of distance from neutron’s initial vertex to the capture vertex(left, area-normalized) and the mean of the same distance by neutron momentum(right) obtained from simulation with the GCALOR and ANNRI-Gd model. . . 29

3.6

The distribution of the distances from neutron’s initial vertex to scattered electron vertices obtained from simulation with the GCALOR and ANNRI-Gd model.
(area-normalized) The red dotted line is positioned at the 1 value, which is ⇠48
cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

LIST OF FIGURES
3.7

xiii

Sample distributions of raw PMT hit times of signal (top) and the residual,i.e., ToF-corrected, hit times of signal (center) and noise (bottom). The ToFcorrected hit times show a much sharper peak for signal hits that are generated
close to a specific vertex within the tank volume, whereas they don’t for noise
hits which mostly occur randomly near the tank wall, and thus lack position
dependence. For this reason, a sliding time window searches for any cluster ofhits more than a pre-defined threshold (e.g., the red dotted line) in the sorted
array of ToF-corrected hit times. Any found set of hits are saved as a neutron
capture “candidate”, and most of the candidates are accidental noise selected
as the example shown in the bottom figure. With the search conditions statedin Table 3.3, the signal-to-noise ratio is ⇠1:30 as shown in Figure 3.12. To
prevent double-counting a single signal, a parameter TMINPEAKSEP is set as
the minimum time di↵erence between any two consecutive candidates. If twoneighboring candidates are separated within TMINPEAKSEP, the candidate
with more hits replace the other. TMINPEAKSEP should be significantly larger
than the resolution (⇠20 ns) shown in Figure 3.9. . . . . . . . . . . . . . . . . . 31

3.8

The coordinate resolution (left) and the distance resolution (right) of the smeared
neutron vertices. (area-normalized) The maximum distance from the smeared
vertex to the true vertex was limited to 150 cm. . . . . . . . . . . . . . . . . . . 32

3.9

The linear-scale (left) and log-scale (right) residual (ToF-corrected) hit time
distribution for signal hits. (area-normallized) The solid lines are the resultsof ideal vertex resolution (fit neutron vertex = true neutron vertex) and thedotted lines are those of realistic resolution (1 resolution being ⇠18 cm). Bothideal and realistic resolutions show almost identical peak widths of ⇠14 ns. The
robustness to the vertex resolution within 50 cm is as expected by Figure 3.6.
The tail following the peak should represent photons traveling longer than the
shortest possible path due to reflection on the tank wall and the glass bulbs
covering PMTs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

xiv

LIST OF FIGURESp
3.10 The estimated significance S/ S + B versus the search window width, for p(n, )signals (top) and Gd(n, ) signals (bottom). The solid line is for ideal and thedotted line is for realistic resolution similar to fiTQun (⇠18 cm). The noise hits
were assumed to be perfectly uniform with the rate of 0.08 hits/ns, which is agood approximation of the noise hit rate within SK6. The maximum significancewas achieved with the search window width of ⇠13 ns for p(n, ) reaction and⇠30 ns for Gd(n, ) reaction. The overall significance was maximized with thewidth of ⇠19 ns. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.11 The raw (left) and area-normalized (right) NHits distributions for noise, p(n, ),and Gd(n, ) candidates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.12 Mean number of candidate class per event. . . . . . . . . . . . . . . . . . . . . . 34

3.13 The MLP input feature distributions of noise candidates (gray), true Gd(n, )candidates (pink), and true p(n, ) candidates (cyan) (area-normalized). 5,000
candidates each are randomly sampled candidates a cold (p =1 keV/c) neutron
MC with the MC configuration explained in Section 3.1. . . . . . . . . . . . . . 43

3.14 Correlation matrices for noise (top), p(n, ) (center), and Gd(n, ) candidates(bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3.15 BCE versus epochs for test and train data sets. . . . . . . . . . . . . . . . . . . 45

3.16 Feature importance.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

3.17 The signal likelihood distributions (top) and the MLP performance (bottom). . . 47

4.1

The pictures of the deployed source, reprinted from [1, 4]. The AmBe compoundis encapsulated in a pellet-sized container, with 8 BGO crystals surrounding it. . 51

LIST OF FIGURES
4.2

xv

Reference AmBe neutron kinetic energy spectra for all neutrons (cited in [40],blue), and for gamma-correlated neutrons ([40], orange). The approximated flat
2-8 MeV input spectrum for correlated neutron simulation is shown in green.
The orange and green histograms are area-normalized, and the blue histogram
is normalized in a way that the peak value around 3 MeV match the one in theorange histogram. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

4.3

The QISMSK distribution of the neutron source data events for Run 85386.
Events with QISMSK > 2,000 are mostly triggered by cosmics. Events with
QISMSK ⇠ 1,000 are likely to have been triggered by BGO scintillation. . . . . 53

4.4

The averaged ToF-corrected hit time distributions near the primary trigger
time range, for events with 800 < QISMSK < 1250 (top) and events with
QISMSK < 400 (bottom). The top histogram is scaled to make the fitted
scintillation histogram (orange) in the bottom figure. The scale factor of 0.11
minimizes the variance of the subtracted (green) bin values with reconstructed
hit time larger than 1,200 ns (i.e., gives the “flattest” subtracted histogram). . . 54

4.5

The FLUKA-simulated AmBe neutron’s energy deposit in BGO (left) and themeasured energy deposit of both source neutrons and gamma-rays (right). Reprinted
from [47]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

4.6

The NHITAC distributions for data and MC before Cut 2 (left) and the same
distributions after Cut 2 (right). The red rectangle represents the selection rangeof Cut 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

4.7

The QISMSK distributions before and after Cut 2. Cut 2 apparently removes
high QISMSK events triggered by cosmic rays. . . . . . . . . . . . . . . . . . . . 60

4.8

The QISMSK distributions of data events before applying Cut 3 (blue), after
applying Cut 3 (orange), and removed by Cut 3 (green). . . . . . . . . . . . . . 60

xvi

LIST OF FIGURES
4.9

The QISMSK distributions for three source z positions. Data for z = ±12 m
show slightly larger QISMSK peak values, as the source is closer to the tankwall. Cut 3 for z = ±12 m leaves some of the events with moderate QISMSK
unlike for z = 0 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

4.10 The MaxNonTrgN1300 distributions for data and MC (explained in Section 4.2
including source/sourceless dummy) (left) and the same distributions within therange of MaxNonTrgN1300 2 [0, 200] (right). The red rectangles represent the
selection range of Cut 5. The long tail in the data distribution should be comingfrom the source itself, and the two distinct scintillation peaks are observable. . . 61

4.11 The TDi↵ distributions for all SHE events (i.e., no cut applied) (blue) and for
AFT events only (i.e., Cut 1 applied) (orange). The exponential fit to the blue
histogram gives the decay time of ⇠12 ms. This manifests the Poissonian natureof the SHE triggers, and thus the reciprocal of the decay time directly relates tothe trigger rate of ⇠83 Hz. For some unknown reason, the AFT flag seems to beomitted more frequently for events close to the previous. . . . . . . . . . . . . . 61

4.12 ReconCT distributions for all candidates after event cuts (left) and after the final
candidate cut (right). The red rectangle in the left figure represents the selection
range of ReconCT > 4,000 ns. The exponential fit to the right figure gives thetime constant of 84 µs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

4.13 SHE-triggered event rates (top) and cut-passing event rates (bottom) by date. . 65

4.14 Mean number of candidates per event, after all event and candidate cuts. The
reason why sourceless dummy MCs have less number of candidates is explainedin Figure 4.15. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

LIST OF FIGURES

xvii

4.15 Mean number of candidates per event by date. The increase of the number
by date has been reported after Gd-loading, and the suspected is the possible
bacterial growth within the tank. Asymmetry in z has also been reported for thenoise hits [48], which may be a reason why data with z = 12 have more candidates
than the other two z positions. The di↵erence in the number of candidates per
event shown in Figure 4.14 should come from the di↵erence in the dates thedummy trigger data were taken. Sourceless dummy data were taken from lateAugust to early September, 2020, whereas the source dummy data was taken on
October 27th , 2020. Data taken on October 27th have slightly more candidates
than the source dummy, which should be due to the BGO’s phosphorescence. . . 65

4.16 Area-normalized feature distributions of sourceless dummy MC and data taken
at z = 0 (Run 85386), for all selected candidates. (blue: data, pink: MC (Gd),cyan: MC (H), gray: MC (noise)) . . . . . . . . . . . . . . . . . . . . . . . . . . 66

4.17 Area-normalized feature distributions of source dummy MC and data taken atz = 0 (Run 85386), for tagged candidates. (blue: data, pink: MC (Gd), cyan:
MC (H), gray: MC (noise)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

4.18 Area-normalized feature distributions of source dummy MC and data taken atz = 0 (Run 85386), for all selected candidates. (blue: data, pink: MC (Gd),cyan: MC (H), gray: MC (noise)) . . . . . . . . . . . . . . . . . . . . . . . . . . 68

4.19 Area-normalized feature distributions of source dummy MC and data taken atz = 0 (Run 85386), for tagged candidates. (blue: data, pink: MC (Gd), cyan:
MC (H), gray: MC (noise)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

4.20 Signal likelihood distribution of source dummy MC and data taken at z = 0 m(area-normalized). The red dotted line indicates the taggable likelihood threshold of 0.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

xviii

LIST OF FIGURES

4.21 Signal likelihood distributions of sourceless dummy MC and data taken at z = 12m (top), z = 0 (center), z =12 m (bottom) (area-normalized). The red dottedline indicates the taggable likelihood threshold of 0.7. . . . . . . . . . . . . . . . 71
4.22 Reconstructed capture time distributions of data (left) and MC (right). The MC
histogram for each z was scaled by (# of data events of the same z)/(# of MC
events). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.23 Fitted capture time constant ⌧ by date. . . . . . . . . . . . . . . . . . . . . . . . 73
4.24 Fitted background B by date. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.25 Estimated neutron detection efficiency ✏est.
by date. . . . . . . . . . . . . . . . . 75
n4.26 Scaled QISMSK distributions by interaction channels, for gamma-correlated neutrons (left) and for prompt gamma-rays (right). (area-normalized) The probability that neutrons will mimic the 4.4-MeV scintillation (i.e., QISMSK⇠1,000p.e.) is O(10 2 ) lower than the prompt gamma-rays. . . . . . . . . . . . . . . . 77

A.1 Feynman diagrams for ⌫l n (left) and ⌫¯l p (right) CCQE interactions. . . . . . . . 86
A.2 The theoretical and experimental values of neutrino-nucleon (left) and antineutrinonucleon CC interaction cross sections (right). Reprinted from Ref. [51]. . . . . . 87

B.1 The neutron yield per alpha from measurement. 10% errors assumed for all data
points. [52] Exponential ae bx + c fit gives a good approximation. . . . . . . . . 89

C.1 Feature distributions of source dummy MC and data taken at z = 12 m (Run85390), for all candidates. (blue: data, pink: MC (Gd), cyan: MC (H), gray:
MC (noise)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
C.2 Feature distributions of source dummy MC and data taken at z = 12 m (Run85390), for tagged candidates. (blue: data, pink: MC (Gd), cyan: MC (H), gray:
MC (noise)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

LIST OF FIGURES
C.3 Feature distributions of source dummy MC and data taken at z =1
12 m (Run85384), for all candidates. (blue: data, pink: MC (Gd), cyan: MC (H), gray:
MC (noise)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
C.4 Feature distributions of source dummy MC and data taken at z =12 m (Run85384), for tagged candidates. (blue: data, pink: MC (Gd), cyan: MC (H), gray:
MC (noise)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

D.1 The disintegrated histograms used for parameter estimation (left) and the fitted
NHits distributions with the base histograms (right). . . . . . . . . . . . . . . . 98
D.2 Parameter fit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

Chapter 1
Introduction
The Super-Kamiokande (SK) experiment was designed to explore realms of physics that havenever been known to mankind, by searching for hints of phenomena that should be either
non-existent or extremely faint, with a gigantic water Cherenkov detector. By detecting supernova burst neutrinos [5] and finding evidence of atmospheric neutrino oscillation [6], the water
Cherenkov detector technique employed by SK has been proven more than just cost-e↵ective,
still running for over 20 years while constantly updating the borderline of physics. Having
cleared two Nobel-prize-winning problems already, SK now aims to search for the remaining
hidden secrets of the universe: the Supernova Relic Neutrinos (SRNs), neutrino CP-violation,and proton decays. All of these unknowns hold keys to the origin of the matter and the universe.
The following is a brief summary of the above three physics programs:• Supernova Relic Neutrinos (SRNs)
Detectable supernova bursts within our galaxy are fairly rare (O(1)/century), but supernovae themselves are not (O(1)/s). Considering the burst neutrino flux and the extremely
low interaction cross section of neutrinos, our universe should be su↵used with the remnant neutrinos from all supernova explosions that has happened so far. By measuring the
SRN flux, one can investigate the mechanisms of nucleosynthesis and stellar evolution,and test the well-established supernova neutrino emission model as well.
2

3• Neutrino CP-violation (CPV)CPV represents the matter-antimmatter asymmetry, which should exist to explain thematter-dominant universe. It was discovered that quarks indeed violate CP-symmetry
via weak interaction, but the amount of quark CPV is insufficient to explain the baryogenesis. The remaining source of CPV is the neutrino, whose three generations and thenon-zero/non-degenerate masses confirmed by neutrino oscillation should enable CPV.
The neutrino CPV in experiments manifests itself as a di↵erence between neutrino andantineutrino oscillation probability. The neutrino oscillation is detailed in Appendix A.

• Proton decay
The Standard Model conserves the baryon number, so that a free proton, the lightest
baryon to exist, can never decay on its own to form a non-baryonic matter. However,the baryon number violation is one of the so-called “Sakharov conditions” [7] that arenecessary to explain the baryogenesis. Most grand unified theories (GUTs) explicitly
break baryon number conservation to predict a proton decay process, e.g., p ! e+ + ⇡ 0with lifetime of O(1030 ) to O(1036 ) years.

The SK detector is explained in Chapter 2.

Unfortunately, with the current detector sensitivity, it may take another 20 years to gather
statistics to an amount required to prove the existence of the above-mentioned phenomena.
To boost sensitivity to its target physics, SK has chosen to load Gadolinium (Gd), the mostneutron-absorbing element known, to its pure-water volume. Why? To detect neutrons, ofcourse. Detecting neutrons, one of the most studied particles in the universe, to identify therarest types of events in physics history, may sound too simple; but even the first discoverers[8] of neutrinos have employed neutrons to identify neutrinos! Besides, in liquid scintillatorsfor neutrino physics, loading Gd has become quite an established techinque. So, why neutrons
all of a sudden? What do they have to do with the rare physics that we pursue?

4

Chapter 1. Introduction

1.1

Why neutrons?

Some of the neutron’s importance in the search for rare physics comes from the characteristics of the water Cherenkov detector. With water Cherenkov detectors, it is difficult to detectheavy charged particles due to high Cherenkov threshold (1.4 GeV for protons) that increases
with particle mass. Neutrons are the only exception, as they are never observed via Cherenkov
radiation, but are detected via capture (i.e., absorption, or (n, )) reactions that produce observable gamma-rays. Hence in water Cherenkov detectors, neutrons are practically the only
long-lived hadrons that can be detected regardless of energy.
The physical importance of neutrons stems from the fact that they are copiously produced via
neutrino interactions in matter. Apart from protons that are rarely visible in SK, neutrons are
the most abundantly produced nucleons in neutrino interaction. Yet, the detection efficiency
is low in pure water, because of the low binding energy (2.2 MeV) of deuterons produced in
the middle of the p(n, ) reaction. By counting the number of detected neutrons produced in
one interaction, we could better classify and reconstruct the event topology, thus improving
sensitivity to the rare physics.
Here are some benefits we expect from detecting neutrons in SK:• SRN search
The dominant interaction of electron neutrinos in the SRN energy range (10 - 60 MeV) isthe so-called “inverse beta decay” (IBD, ⌫¯e +p ! e+ +n) on a free proton. Unfortunately,
there are several sources of backgrounds, especially the solar and atmospheric neutrinos
that can mimic the IBD in a conventional pure-water Cherenkov detector with minimal
neutron detection capability. For example, an “invisible muon” event, where a lowenergy (below Cherenkov threshold) neutrino-induced muon decays into an electron withitself undetected, has the same observable particles and thus becomes indistinguishable
from signal (IBD) events. Another possibility of fake SRN events is a scenario where ⌫eelastically scatters a bound electron. The SRNs are expected to be far outnumbered bysuch backgrounds. The estimated flux distributions of signal SRN and its backgrounds

1.1. Why neutrons?

5

are shown in Figure 1.1. If we could efficiently detect neutrons from data, then the
aforementioned most dominant backgrounds that do not necessarily produce neutrons
should be largely removable from analysis by requiring a neutron. We could also improve
the supernova direction reconstruction, as we could purify the SRN scattering events
where forward-scattered electrons keep the direction of the SRNs, from the SRN IBD
events where neutrons obscure the original SRN direction.
• Neutrino CPV search and oscillation analysis
Neutrino CPV search is a part of the neutrino oscillation analysis, which is a processthat estimates the neutrino oscillation parameters by fitting a parameterized “oscillated”(i.e., predicted by the oscillation probability with a given set of oscillation parameters,
the current knowledge on neutrino cross sections, and the estimated initial neutrino flux)
neutrino energy spectrum to the one measured in data. For a brief introduction to theoscillation probability and parameters, refer to Appendix A.
The neutrino sources used in the SK oscillation analysis are atmospheric, accelerator,and solar neutrinos, each with distinct ranges of energy and path length. The main
target interaction in the oscillation analysis is the charged-current quasi-elastic (CCQE,⌫l +n ! l +p for neutrinos and ⌫¯l +p ! l+ +n for antineutrinos, with flavor l 2 {e, µ, ⌧ })interaction, which gives out the interacting neutrino’s flavor by producing a same-flavored
lepton. The SK detector is capable of separating electrons from muons as illustrated in
Figure 2.2, so that we can estimate the counts of e-like and µ-like CCQE events within
the detector. Neutral-current (NC) neutrino scattering on a bound electron can be asource of background, since the scattered electron does not carry any flavor information.
Requiring a neutron from an event sample can help remove this background, which isespecially dominant for low-energy solar antineutrinos.
Another necessity of neutron detection in SK is supported by the fact that in SK, detecting final-state neutrons is practically the only way of distinguishing antineutrinos fromneutrinos. (Recall that only antineutrinos emit neutrons via CCQE interaction.) The⌫/¯⌫ separability is especially important for atmospheric neutrino analysis, in which we
expect both neutrinos and antineutrinos in comparable amounts. With greater neutron

6

Chapter 1. Introduction
detection efficiency, we expect improved sensitivity to neutrino mass hierarchy (MH) inatmospheric analysis [1]. MH is important, since it is correlated to the overall oscillation
analysis including the CPV search. Refer to Appendix A for MH. For accelerator neutrino analysis, the advantage of the separability is relatively small compared to that in
the atmospheric analysis because the initial ⌫/¯⌫ flux can be controlled. Still, accidental
wrong-sign backgrounds can be efficiently removed via neutron detection to improve theaccelerator analysis’ sensitivity to CPV.
Another expected advantage lies in the neutrino energy reconstruction. The initial neutrino energy spectrum in typical oscillation analysis is quite wide (O(0.1) GeV for accelerator experiments, to say the least), while the oscillation probability depends on theincident neutrino energy. For this reason, precision in the energy reconstruction shouldbe critical to the oscillation parameter estimation. Multiple neutrons produced by the
intranuclear cascade of hadronic interaction should carry a part of the neutrino energy,
which should be possible to retrieve if we could efficiently detect neutrons.
• Proton decay search
If proton decays were to occur via p ! e+ + ⇡ 0 in the SK detector, there would be ideallyan event with three electron-like rings, one created by the positron and the other two
imprinted by the two gamma-rays from the instantaneously decaying ⇡ 0 . Atmospheric
neutrinos can fake the event topology by producing a charged lepton and a ⇡ 0 (via either
primary ⌫N interactions, FSI, or SI explained in Section 1.3) within the detector, andrequiring zero detected neutron in the sample can efficiently reject this background.

Figure 1.1: (a) number flux of SRNs compared with other background neutrinos, and (b) event
rate of SRNs and invisible muon decay products. Reprinted from Ref. [9]1.2. Neutron production in SK

1.2

7

Neutron production in SK

The background neutrons to the SK detector are mainly produced via spallation, where highenergy cosmic rays inelastically reacts with matter nuclei to scatter o↵ a large number ofnucleons. Hadronic particles, such as neutrons themselves, protons, alpha particles, and pions
can also interact with nuclei to free bound neutrons. Another background radioactivity is
spontaneous fission, in which a heavy nuclide spontaneously undergoes nuclear fission to emit
neutrons. A small number of 238 U contained in Gd salt is subject to this process, although
its contribution to the background neutrons should be negligible due to its very long half-life(O(109 ) years).
Neutrons can also be produced via neutrino-nucleus interaction. The following is the summaryof the three stages in which neutrons can be produced:• Primary ⌫N interaction:
Antineutrinos can primarily undergo a CC or NC exchange with nuclei to knock out aneutron, where neutrinos can knock out neutrons via NC exchange only. Both types witha higher energy can also excite the nucleon to a baryonic resonance state (RES), which
immediately de-excite to a nucleon and produce pions, kaons, and gamma-rays. At evenhigher energies, neutrinos can scatter o↵ a single quark to produce more hadrons (DIS).
• FSI (Final-state interactions)
The final-states of the primary ⌫N interaction (mostly nucleons and hadrons) are very
likely to undergo strong interaction within the nucleus, to produce even more hadrons.
• SI (Secondary interactions)
Knockout nucleons and other hadrons (mostly pions) produced within the nucleus canfurther interact with the other surrounding nuclei even after they exit the primary nucleus.

Figure 1.2 shows how the FSI/SI can obscure the final neutron multiplicity in CCQE samples.
Nonetheless, a non-zero multiplicity proves to be an e↵ective indicator of ⌫¯-CCQE.

8Chapter 1. Introduction

Figure 1.2: Simulated ⌫/¯⌫ CCQE neutron multiplicity distributions from primary ⌫N interaction (left), from ⌫N interaction including FSIs (center), and for ⌫N interaction including both
FSIs and SIs (right). Reprinted from Ref. [4].

1.3

Neutron reactions in Gd-loaded water

Neutrons are charge-neutral, so that they do not undergo Coulomb interaction, but only strong
interaction with nuclei in extreme proximity (⇠ O(1) fm). Due to such property, neutrons are
more penetrable than other charged particles of similar mass. Moreover, due to their gas/wavelike properties at lower energies, neutrons are often described using optical and thermodynamic
models. It is conventional to classify neutrons by their kinetic energy as follows:• Cold neutrons (K.E. < 0.025 eV)Cold neutrons are neutrons slower than most neutrons at room temperature (290 K).
Ultracold (< O(100) neV) neutrons are subject to total reflection upon matter nuclei
with high optical potential, so that they can be treated as ideal gas. Cold neutrons arewidely used in the study of internal crystal structure of solids, as well as in precision
measurements of free neutron properties such as electron dipole moment. Despite their
physical versatility, their preponderance within SK (⇠285 K) should be negligible.
• Thermal neutrons (K.E. ⇠0.025 eV)
Thermal neutrons are neutrons with kinetic energy close to the mode (⇠25 meV) of theMaxwell–Boltzmann distribution at around room temperature (⇠290 K). For neutrons
within SK, this energy range is practically the lower limit they can reach. Thermal
neutrons have wavelengths of O(0.1) nm.

1.3. Neutron reactions in Gd-loaded water9• Fast neutrons (1 < K.E. < 20 MeV)
Fast neutrons with kinetic energy of O(1) MeV are often products of nuclear fission andcosmic spallation. Many of the neutrons produced by high energy hadronic scattering(including ⌫N interaction) are expected to fall into this range.

We expect most of the neutrons produced in SK to have kinetic energy ranging from fast to
thermal. The neutron-nuclear reactions can be categorized into two types, according to theirinteraction times [10]:• Direct reaction
A direct reaction can occur when a projectile neutron has enough kinetic energy to have awavelength of O(1) fm (i.e., ⇠20 MeV nucleon) and a relatively large impact parameter,
so that it can interact with a single valence nucleon of the target nucleus. The unstablebound state of the incident neutron and the valence nucleon typically lives for about
10 22 s, which is in the order of time required to make several round trips (in classicalterms) around the nucleus, before one of the projectile neutron and the nucleon finds a
way to escape. The ejected particle is generally forwardly-directed with strong angularcorrelation, following the energy-momentum conservation.
• Compound Nucleus (CN) reaction
When the projectile neutron collides with a nucleus with a small impact parameter, itinteracts many times inside the nucleus, boosting individual nucleons into excited states,
until it comes to rest inside the nucleus. The thermal equilibrium model is used to modelsuch nuclear reactions. The excess energy introduced to the “compound nucleus” by
the incident neutron is shared by many of the bound nucleons via multiple intranuclear
scatterings, so the de-excitation probability of one specific nucleon to gain a large enoughshare of the energy to escape becomes statistically small. Such statistical constraint
makes the overall lifetime of the CN to be O(10 18 ) to O(10 16 ) s, which is substantiallylonger than that of direct reactions. The angular distribution of the products are nearly
isotropic, owing to the random intranuclear interactions, especially for light projectiles.

10

Chapter 1. Introduction

Neutrons produced within the SK detector loaded with Gd can undergo both direct and CN
reactions with matter nuclei. In pure water, the most dominant direct reaction channel forneutrons with kinetic energy below 20 MeV is the elastic scattering with free protons (i.e.,p(n, n)p). Refer to Figure 1.3 for the cross sections. The elastic scattering of an energetic
neutron within the sea of free protons at rest quickly moderates the neutron down to thethermal regime within a few µs.1 A moderated neutron is mostly captured by a proton to form
a deuteron2 , emitting the excess 2.2 MeV binding energy as a form of a monoenergetic gammaray. This absorption process can be viewed as a CN reaction where the CN is the deuteron andthe lifetime is infinite (stable), so that the radiation is almost isotropic.
In case a small amount of Gd is dissolved in water, Gd can also capture a fraction of moderated
neutrons, as Gd has a few order of magnitude higher absorption cross section than 1 H or 16 O.
The fraction of captures on Gd with varying Gd mass fraction in water is shown in Figure1.5. The radiative neutron capture on Gd (i.e., Gd(n, )) is an example of a CN reaction,
where the CN in thermal equilibrium de-excites via emitting multiple gamma-rays. In caseof the most naturally abundant and most strongly absorbing 157 Gd, the capture formula is157⇤

Gd + n ! 158 Gd ! 158 Gd +(7.94 MeV). The total radiated energy of 7.94 MeV is themass di↵erence between the initial and final states, and is shared by multiple gamma-rays
that are successively emitted in the process of multi-level de-excitation. The schematic of
the process is shown in Figure 1.6. The energy distribution of the gamma-rays is random innature, although there is a high chance of producing multiple low-energy gamma-rays due tothe existence of quasi-continuum levels as shown in Figure 1.6.
Although most neutrons are expected to be thermally captured via either p(n, ) or Gd(n, ),a small fraction (⇠2 % or below according to the simulation explained in Section 3.1) may
undergo inelastic interactions before they thermalize, sometimes multiplying or removing theincident neutron. The cross sections are shown in Figure 1.4. While quantitative evaluation of
those channels is necessary, we focus on p(n, ) and Gd(n, ) reactions of moderated neutrons.
11

H is actually the most e↵ective nuclide in moderating energetic neutrons. The mean logarithmic reduction2
1
2
of neutron energy per collision is described as < ⇠ >=< ln EE0 >= 1 + (A2A1) ln A
A+1 ⇡ A [11], where A is themass number of the target nuclide.
2
The thermal neutron capture cross section on 16 O is negligible compared to that on 1 H.

1.3. Neutron reactions in Gd-loaded water

11

Figure 1.3: Neutron cross sections on 1 H, 16 O, and 155 Gd from the JENDL-4.0 library [12].
Neutrons are more likely to be captured with lower kinetic energy, qualitatively as they tend to
spend more time in the vicinity of nuclei with lower speed. The large fluctuation in the intermediate energy range can be understood as manifestation of resonances in the distinguishable
discrete energy levels. A formal approach can be found in Ref. [10].

12Chapter 1. Introduction

Figure 1.4: Neutron reaction cross sections on 16 O and 155 Gd, from the JENDL-4.0 library [12].
Note the small fraction of channels either removing or multiplying the number of neutrons.

Figure 1.5: Capture ratio on Gd, by Gd mass fraction in water. Reprinted from Ref. [13].

⇤Figure 1.6: The de-excitation schematic of the excited state 157 Gd + n ! 158 Gd to the 158 Gd
ground state. Multiple gamma-rays from the 157 Gd(n, )158 Gd reaction can share 7.937 MeV.
Reprinted from Ref. [14].

1.4. Outline of this study

1.4

13

Outline of this study

With high neutron detection efficiency brought forth by dissolving Gd into the SK’s detector
volume, we expect to improve background rejection as explained in Section 1.1. The SK-Gd
phase has officially started since July 2020 with Gd mass fraction of 0.0110 wt.%, and we needto evaluate the change in neutron detection efficiency.
SK already has the so-called “neutron-tagging” algorithm that has been used to detect neutronsin pure-water data. The workflow of this algorithm involves a multivariate analysis technique
where we evaluate the signal likelihood of unseen data with an artificial neural network trainedwith labeled data extracted from the detector simulation. The process is detailed in Section
3, and the method has been proven valid for p(n, ) detection in pure water [1, 3, 4]. In thisstudy, we take the same approach to the Gd-loaded data as well, but this time, we optimize
the relevant parameters to the SK-Gd simulation that has become ready for use and validation.
The tuned algorithm was applied to a neutron source data, to evaluate the neutron detection
efficiency of the Gd-loaded detector and the reproducibility of the data with the simulation.
The analysis process and the results are discussed in Section 4.

Chapter 2
The Super-Kamiokande detector
The Super-Kamiokande (SK) detector is an underground water Cherenkov detector, which
as the name suggests, detects Cherenkov light emitted in water. When a charged particle
passes through a dielectric (polarizable) medium, it can electromagnetically interact with thesurrounding molecules to excite them and make them radiate. In case the particle travels fasterthan the phase velocity of electromagnetic wave in the medium, simultaneous radiation fromdi↵erent positions can share the wavefront to form a sort of shock-wave of angle ✓ = cos 1 ( /n).
The phenomenon itself and its anisotropy was first discovered by P. A. Cherenkov [15] for whom
the phenomenon was named. The schematic is shown in Figure 2.1.

Figure 2.1: Schematics of Cherenkov radiation (left) and how it forms a ring-shaped image ona photodetector plane (right). Reprinted from Ref. [16].

1415The energy radiated per track length is modeled by the Frank-Tamm formula (in cgs) [17]
dE
(ze)2
= 2
dx
c

where ze and

Z

✓

! 1
✏(!)>1/ 2

12 ✏(!)

◆

d!

(2.1)

are each the charge and the speed (v/c) of the moving particle, ✏ is thepermittivity of the medium, ! is the frequency of the radiation (assuming magnetic permeabilitypµ(!) = 1 so that n(!) = ✏(!)). For an ultrarelativistic particle of ⇡ 1, the energy loss perpath length dE/dx in Equation 2.1 and the Cherenkov angle ✓ = cos 1 ( /n) shown in Figure
2.1 are almost fixed, given the particle charge and the dielectric properties of the medium. For
instance, the energy loss due to an ultrarelativistic particle with unit charge passing throughwater where the Cherenkov angle ✓ is ⇠42 is O(1) keV/cm.
Since the Cherenkov radiation stops as the particle is slowed down below the Cherenkov threshold speed ofmin = 1/n, the radiation from most charged particles with moderate energies are

ring-shaped, as illustrated in Figure 2.1. Hence, if we create a photodetector plane as in Figure
2.1 to measure the radius of the ring and estimate the amount of radiated photons, we cantechnically reconstruct the position and the energy of the incoming charged particle. Although
the energy deposited to the detector system via Cherenkov radiation is extremely small (lessthan 1 %) compared to the energy lost by ionization1 , reconstructing an event from detected
Cherenkov light is a practical option for a gigantic detector with zero net electromagnetic field.
A disadvantage to this technique is that the Cherenkov light from heavy particles is often
barely visible at best due to the energy threshold of Cherenkov radiation increasing with particle mass.2 This is an inevitable sacrifice in this technique whose priority is set to increasing
neutrino event statistics. Fortunately, electrons and muons are detectable, and even separable with the water Cherenkov technique. Electrons are more likely to induce electromagnetic
interactions including pair production and bremmstrahlung3 , so that the Cherenkov rings generated by electrons are often “fuzzier” than their muon counterparts. The “fuzziness” can bequantified, and be used as a particle identification discriminator as illustrated in Figure 2.2.
1

For example, the total stopping power of ⇠2 MeV/cm for minimum ionizing particle in water.
0.8 MeV for e, 160 MeV for µ, 210 MeV for ⇡, 750 MeV for K, and 1.4 GeV for p.
3
The critical energy of electrons in water is ⇠70 MeV, whereas the same for muons are ⇠1 TeV.
2

16

Chapter 2. The Super-Kamiokande detector

Figure 2.2: ⌫e /⌫µ event identification amongst beam (T2K) neutrino data taken in SK, by theclarity of Cherenkov rings. Reprinted from Ref. [18].

2.1

Specification

The detector is as illustrated in Figure 2.3. The detector is located at 36 25’32.6” N, 137 18’37.1”
E, about 1,000 m underground in the Mozumi Mine in Hida’s Kamioka area with 2,700 mwater-equivalent (m.w.e.) mean overburden, which decreases the comsic ray flux by the factorof O(10 5 ) compared to the sea level intensity.4 The detectable background gamma-ray intensity (E > 0.5 MeV) is O(0.1) cm2 ·s 1 ·sr 1 , and the background neutron rate is O(10 5 )
cm2 ·s 1 .5 The average atmospheric neutrino event rate is ⇠8 per day.
The detector is cylinder-shaped, 39 m in diameter and 42 m in height. The inner photomultipliertube (PMT) support structure divides the tank into two distinct, optically isolated volumes.
The inner detector (ID) volume has dimensions of 33.8 m (diameter) by 36.2 m (height),
containing 32 ktons of water and 11,146 inward-facing 50 cm PMTs (Hamamatsu R3600) on
its inner wall. The inner wall except the PMTs is covered with black sheet to reduce reflection.
The fiducial volume (FV) of the ID is defined as the cylindrical volume 2 m inward from the
ID wall, and its total mass is 22.5 ktons. The mean distance from the inner tank wall to theouter wall is about 2.5 m.
4
5The cosmic muon flux at SK is ⇠ O(10 7 ) cm2 ·s 1 ·sr 1 . (⇠2 Hz event rate at the detector)
The measurement was done before 1994.

2.1. Specification

17

The outer detector (OD) volume is instrumented with 1,885 outward-facing 20 cm PMTs
(Hamamatsu R1408), serving as an active veto counter against incoming particles as well as apassive shield for neutrons and gamma-rays from the surrounding rocks. The OD region has
a width of about ⇠2 m with a total water mass of about 18 kilotons. The OD wall is covered
with white Tyvek sheet to maximize reflection and photon detection efficiency.

The two detector volumes are isolated from each other with light-proof sheets covering bothinner and outer surfaces of the PMT support structure. The ID PMTs are mounted on a 70
cm grid, with 7,650 on the barrel (side walls), 1,748 on the top, and the remaining 1,748 on thebottom. The e↵ective photocoverage of the ID surface is ⇠40 %. The schematic of the detector
volumes is shown in Figure 2.4. The Cherenkov radiation induced by a charged particle moving
inside the detector can trigger the detector (if it’s strong enough) to imprint a ring-shaped
image on the detector surface. In case the particle is created within the ID, there would be
no simultaneously detected signal in the OD. Thus, neutrino-induced particles are without any
OD signature, and a sample muon-neutrino CC event display is illustrated in Figure 2.5.

The detector should be calibrated every once in a while to understand the detector capability
and keep it from changing over time. The optical parameters of water within the detector
volume are calibrated using laser. The relative PMT gain and timing distributions is calibratedusing a di↵user ball, which is a point light source with uniform intensity over all solid angles.
The absolute energy scale of the PMTs are calibrated using various physical sources, including low-energy electrons from a linear accelerator and neutrons and the accompanied capture
gamma-rays from a point source. Cosmic ray muons have a well-known energy distribution sothat they are used to cross-check the calibration. The calibrated parameters are introduced tothe detector simulation to maximize the reproducibility.

To minimize the signal distortion caused by the geomagnetic field (⇠450 mG), 26 Helmholtz
coils are lined up along the tank wall, reducing the net magnetic field to 50 mG.

18

Chapter 2. The Super-Kamiokande detector

Figure 2.3: The SK detector.

Figure 2.4: The schematic of the SK detector.

2.1. Specification

19

Figure 2.5: Sample event display of 0.6 GeV muon-neutrino.

2.1.1

ID PMTs and electronics

The ID PMTs are 50 cm Hamamatsu R3600 PMTs. The design is as shown in Figure 2.6. The
photocathode of the ID PMT is coated with a bialkali (Sb–K–Cs) that matches the wavelengthof the Cherenkov light. Its sensitivity range is 280 to 660 nm. The maximum quantum efficiency(QE) is about 21% at its peak ( ⇠380 nm). The average collection efficiency at the first dynodeis 70 %, which is uniform within ±7 % across the whole PMT surface. The dynode structure
follows the Venetian Blind style, and the gain is 1 ⇥ 107 for 2,000-V high voltage. The relativetransit time spread for the single photoelectron signal is about 2.2 ns (1 ). A PMT hit isdefined as an individual signal pulse with over 0.25-photoelectron-equivalent charge deposit.
The average dark-hit rate at the 0.25-photoelectron threshold is known to be 7 kHz.
The current signal from each ID PMT is sent to a QTC6 -based Electronics with Ethernet(QBEE) system. The maximum detectable integrated charge without saturation for each PMT
6Charge-to-time converter.

20

Chapter 2. The Super-Kamiokande detector

is 2500 pC (equivalent to ⇠1,300 photoelectrons). A time-to-digital converter (TDC) includedin the system digitizes the time and deposited charge of each PMT hit, and send the information
to online frontend PCs. Technically, all recorded hits are saved o✏ine. A software trigger isapplied to the recorded/sorted hits to issue an event.

Figure 2.6: ID PMT.

Figure 2.7: The schematic of gates issued after a PMT hit. The charging and discharging
gates opens for 400 ns and 350 ns, respectively. All signal currents are ignored while one ofthe discharging gate and VETO is open. The total amount of time required to process a single
input current is 900 ns. The timing of each hit is determined by the comparator threshold.

2.1.2

Water purification

The water that flows into the SK detector is taken from the mine, and it is critical to keep
the water as pure and transparent as possible. Radioisotopes, notably Radon and Radium,
are naturally dissolved in air and water. Such unwanted radioactivity increases noise hits and2.1. Specification

21

mimics low-energy PMT hits such as neutron (n, ) signals. For these reasons, the water within
the detector is constantly circulated through a purification system at a flow rate of ⇠50 tons/h.
The purification system has many stages, but they all boil down to three main objectives:filtering out particulates of finite size, controlling growth of bacteria, and reducing Radoncontents through degasification. The purity of water is represented by the measured resisitivity.
While the value before applying the purification system is typically about 11 M⌦·cm, the value
after applying the system is 18.2 M⌦·cm, reaching the chemical limit of the ultra-pure water.
The photon attenuation length at this state is ⇠90 m.

2.1.3

Software trigger

Since all hits are saved o✏ine, events are issued by “software” triggers. If the number of events
within a certain 200-ns window is over a pre-defined threshold, records of PMT hits within theevent window is set as one event. Table 2.1 summarizes the thresholds and the window widths
for each event trigger type.
Trigger TypeSpecial High Energy (SHE)High Energy (HE)Low Energy (LE)Super Low Energy (SLE)
Outer Detector (OD)

Threshold [hits/200 ns]
58
50
47
31
22 (in OD)

Event window [µs][-5, +35][-5, +35][-5, +35][-0.5, +1.0][-5, +35]

Table 2.1: Summary of the software trigger types.
For SHE triggers only, a special attachment flag named AFT (after-trigger) is raised in order tocontain (n, ) reactions that happen with a time constant of O(100) µs, and additional 500-µs
record of PMT hits is attached to the SHE-triggered event. The AFT attachment is not done
for all SHE events, however, due to additional triggers near the primary SHE events.
Another special type of trigger is a random trigger, also known as the “dummy” trigger within
SK. A “dummy” event is triggered by nothing but a clock with fixed frequency7 , so that it istypically used to model dark noise hits in the detector simulation.
7

Usually this “clock” is set to the T2K beam frequency, but it can also be used as a “dummy” trigger incase the beam is o↵! Dummy-triggered events have an event window of [-500, +500] µs.

22

2.2

Chapter 2. The Super-Kamiokande detector

SK-Gd

13 tons of Gd2 (SO4 )3 · 8 H2 O was put into the tank via a dedicated dissolving system, startingfrom July 14th , 2020 until August 17th , 2020. The loaded amount corresponds to a 0.026 %concentration of Gd2 (SO4 )3 · 8 H2 O by weight dissolved in 50 ktons of total pure water in thedetector system, and the concentration is equivalent to the Gd mass fraction of 0.0110 wt.%.
The loaded Gd-sulfate and the amount loaded by date are shown in Figure 2.8.

Figure 2.8: The Gd2 (SO4 )3 · 8 H2 O powder (left) and the amount of the Gd-sulfate loaded into
the detector by date (right). Reprinted from Ref. [19].
A new Gd-loaded water purification and circulation system was developed to load Gd while
maintaining the pure water transparency as is. To achieve this, a special ion exchange resin
was introduced to remove all ionic impurities in water except Gd3+ and SO4 2 .

Figure 2.9: A block diagram of the Gd-dissolving system. Reprinted from Ref. [19].
Pure water was sent from the SK tank to the Gd-loading system at a flow rate of 60 tons/h.

2.2. SK-Gd

23

The flow was separated into two streams, one at a rate of 48 tons/h, and the other at a rate of
12 tons/h into which the Gd-sulfate powder was dissolved. The Gd2 (SO4 )3 · 8 H2 O powder was
put into the the powder feeder, and then sent to a pump for dissolving. The compound wascirculated at a high speed within the dissolution tank, resulting in a 0.13 % Gd2 (SO4 )3 · 8 H2 O
solution. The solution was purified in the “pre-treatment” section of the system, and mergedwith the 48 tons/h pure water stream to make a 0.026 % Gd2 (SO4 )3 · 8 H2 O solution. The final
solution was sent to the SK tank using the “recirculation” section of the system, which continues
to recirculate the SK water, even after the Gd loading was completed. The temperature of the
supplied Gd water was set to about 0.3 C lower than that of the tank water just before the Gdloading, to make it gradually fill the tank from the bottom. The overall system is as illustrated
in Figure 2.9.
Since the Gd-loaded water was filled from the bottom, the concentration by height was measured
using the Atomic Absorption Spectroscopy (AAS) technique8 and the conductivity measurement. A sample of Gd-loaded water was taken from the tank with varying heights, and the
Gd concentration was investigated daily during the Gd-loading phase. The results are shownin Figure 2.10. The Gd concentration at z = 12 m reaches ⇠250 ppm by August 14th , 2020.

Figure 2.10: Measured Gd concentration by height and date. Reprinted from Ref. [19].
The neutron data used in this study was mainly taken during the 3-month commissioning phase(September - December, 2020) after the Gd-loading was completed.

8This technique makes use of the wavelengths of light specifically absorbed by a vaporized element, in ourcase, Gd.

Chapter 3
The neutron detection algorithm
Most neutrons produced within the SK detector eventually lose energy to be absorbed by anatomic nucleus within the tank volume. The compound nucleus that has absorbed a neutron
immediately de-excites and emits one or several gamma-rays with a fixed sum of energy. If
these gamma-rays scatter o↵ a bound electron hard enough (ve > 0.75c, KEe > 0.8 MeV) togenerate Cherenkov photons, the PMTs on the inner surface of the tank will be able to detectthem. We call such neutron-induced PMT hits a neutron capture signal, or in short, a neutron
signal. A neutron capture signal might have only a couple of hits, or sometimes, more thana few dozen hits depending on the capturing nucleus and its position within the tank. The
detector size of O(10) m limits the hit time spread to be O(100) ns. Figure 3.1 shows thenumber of PMT hits and the time spread distributions obtained from simulation for p(n, )and Gd(n, ) reactions, together which are expected to take up more than 99.98% of the total
captures within ⇠0.11 wt.% Gd-loaded water according to Table 3.1.
GdHOn [/cm3 ]4.22 ⇥ 1014
6.72 ⇥ 1019
3.36 ⇥ 1019

abs [barns]
44.87 ⇥ 10
3.32 ⇥ 10 1
1.90 ⇥ 10 4

n abs [10 5 /cm]
2.06
2.23
6.38 ⇥ 10 4

Table 3.1: The number density n, the thermal neutron absorption cross section abs , andtheir product n abs for each of the top three neutron-absorbing elements: Gd, H, and O. The
number densities are for 0.0110 wt.% Gd-loaded water and the cross sections are taken fromthe ENDF/B-VII.1 data [20]. The n abs is equivalent to the capture rate per neutron path
length.
24

25

Figure 3.1: The number of hits (left) and the raw PMT hit time distributions (right) obtained
from simulations for p(n, ) and Gd(n, ) reactions. (random neutron vertex, area-normalized)
The simulation configuration is as explained in Section 3.1.
Figuring out whether a set of PMT hits is coming from a neutron, given their number andtime spread, is not an easy task, as the total of ⇠12,000 PMTs would record on average 8 hits
per 100 ns even without any external source. Such hits that originate from the detector itself
or the surrounding environment are called the “dark” hits, and their primary source is thethermionic emissions (i.e., dark current) within each PMT. Dark hits would sometimes total
up to 30 hits in a 100 ns interval, making it challenging for us to pick out neutron signals by
their hit counts. If we were to do so, we would have to compromise between signal purity anddetection efficiency. To maximize both, we need to find other distinctive features of signal hits.
Multivariate statistics has studied viable solutions to classifying unseen data by statisticallyanalyzing the known feature distributions of di↵erent categories. The same problem is now
widely addressed in machine learning and experimental particle physics as well, where classifying signals of physical interest is increasingly becoming difficult. For the past few years, SK has
also tested a couple of multivariate analysis techniques such as a Multilayer Perceptron (MLP)[1] and Boosted Decision Trees (BDT) [2, 21] in identifying neutron captures within the purewater data and the simulated SK-Gd events. The methods have proven valid for neutron signal
classification, and thus the same methodology was inherited and improved.
In this study, we target both p(n, ) and Gd(n, ) signals, as the Gd-loading has just finishedits commissioning phase with Gd concentration being ⇠10% of the target value of 0.1 wt.%.

26

Chapter 3. The neutron detection algorithm

The author would like to clarify that designing the algorithm from scratch was not part of theauthor’s original work. The author’s contribution to this algorithm lies in developing an easily
tunable program and optimizing it with simulation of the ⇠0.0110 wt.% Gd-loaded detector.
The improvements made by this study will be specified. The neutron detection program1 has
been fully re-written by the author.
The neutron detection algorithm used in SK neutron analysis consists of two stages. In the firststage, we look for a cluster of hits whose number and hit time distribution resemble those of
signal hits. In the second stage, we statistically analyze each selected hit-cluster’s features to
examine its signal likelihood, using a multivariate analysis technique. The following subsections
will discuss the specific design and the tuning process for each stage of the algorithm.

3.1

Signal properties of physics simulation models

Before testing and optimizing the neutron detection algorithm, the properties of signal hits
were investigated with Monte Carlo (MC) simulation. The simulator skdetsim-skgd [22],whose framework other than the neutron capture models fully inherits skdetsim, was used
to simulate neutron capture events. For low-energy (KEn < 20 MeV) neutron transport,
GCALOR (Geant-CALOR interface) is called [23]. For neutron captures, the same GCALOR
is called for p(n, ) reaction, and the ANNRI-Gd model [24] based on the results of the ANNRI
experiment [14, 25] at J-PARC is used for Gd(n, ) reaction. There are a couple of caveats
worth mentioning; the water parameters such as transparency and absorption lengths are tuned
for pure water, and captures on nuclides other than Gd or H are not simulated. Noise hits are
not simulated as well, but extracted from randomly sampled dummy trigger events taken at
the end of the SK6 commissioning phase (Run 85200 - 85245), for more realistic reproductionof data events. A brief explanation on the dummy trigger can be found in Subsection 2.1.3.

1

The code is maintained on a GitHub repository at https://github.com/seunghophysics/NTag.
SK libraries and ROOT are required. Basic instructions can be found at https://ntag.readthedocs.io.
For the API documentation, please refer to https://www-sk.icrr.u-tokyo.ac.jp/~han/NTag/index.html.
(access to SK’s internal network required)3.1. Signal properties of physics simulation models

27

Figure 3.2: Sample MC event displays for a p(n, ) (left) and a Gd(n, ) event (right). The top
figures are the 3D drawings of the inner tank dimensions with all neutron signal hits (red dots)
and gamma-rays (red arrows) in each event. The bottom histograms are the raw PMT hit time
distributions of each event with bin width being 1 µs and 10 ns. The simulation configuration
is as explained in Section 3.1.

28

Chapter 3. The neutron detection algorithm

Figure 3.2 shows sample displays of the simulated p(n, ) and Gd(n, ) events. As can beinferred from the figure, the p(n, ) and Gd(n, ) signal hits should have distinctive properties.
Figure 3.3 shows the gamma-ray multiplicity and total energy distributions for captures on both
nuclides. The p(n, ) reaction emits a single monoenergetic 2.22 MeV gamma-ray, whereas theGd(n, ) reaction emits 4 to 8 gamma-rays with their energies totaling 7.98 MeV for captures
on 157 Gd and 8.54 MeV for captures on 155 Gd. Figure 3.4 shows the energy distribution ofindividual simulated gamma-rays, which are mostly sub-MeV for the Gd(n, ) reaction. These
individual gamma-rays are generated with uniform angular distribution for 4⇡ solid angle.

Figure 3.3: The gamma-ray multiplicity (left) and the energy sum of all emitted gamma-rays(right) obtained from simulation with the GCALOR and ANNRI-Gd model. (area-normalized)

Figure 3.4: The energy distribution of the individual gamma-rays obtained from simulation
with the GCALOR and ANNRI-Gd model (left) and the same for Gd(n, ) reaction with fourdi↵erent physics models available in the simulator, the default Geant4 [26], GLG4Sim [27],GGarnet [24], and ANNRI-Gd [24]. (area-normalized)

3.1. Signal properties of physics simulation models

29

Neutrons travel through a random path before they are thermalized and finally captured bya surrounding nucleus. The distances from a neutron’s initial vertex to its capture position
with its momentum between 60 to 120 MeV/c were obtained from simulation, and Figure 3.5
shows their distribution. The left figure suggests that low-energy neutrons with kinetic energy
below 10 MeV are likely to be captured ⇠20 cm away from the initial vertex. Clear momentum
dependence of the distance can be observed from the right figure. The figures also imply thatthe capture vertex resolution of any neutrino vertex fitter, which only works on the recoil
leptons, should be at least ⇠20 cm.

Figure 3.5: The distribution of distance from neutron’s initial vertex to the capture vertex (left,area-normalized) and the mean of the same distance by neutron momentum (right) obtained
from simulation with the GCALOR and ANNRI-Gd model.
However, what physically hit the PMTs are the Cherenkov photons generated by Comptonscattered electrons, rather than the capture gamma-rays themselves. These gamma-rays should
travel another few tens of centimeters before they scatter o↵ an electron. Since the scatteredelectrons are usually barely over the Cherenkov threshold, they can only radiate instantaneously.
Hence, it is safe to assume that all Cherenkov radiation has originated from the electron vertex.
Figure 3.6 displays the distribution of the distances from the neutron’s initial vertex to verticesof scattered electrons. According to this figure, signal hits should be coming from coordinates
20 to 60 cm (1 : 48 cm) away from the neutron’s production point. This implies that theresolution of the neutrino vertex fitter, if we were to use it to estimate the neutron vertex, would
not significantly a↵ect the overall performance of neutron detection as long as the resolution
stays within 50 cm.

30

Chapter 3. The neutron detection algorithm

Figure 3.6: The distribution of the distances from neutron’s initial vertex to scattered electron
vertices obtained from simulation with the GCALOR and ANNRI-Gd model. (area-normalized)
The red dotted line is positioned at the 1 value, which is ⇠48 cm.

3.2

Stage 1: Hit-cluster search

The first step of the neutron detection algorithm is to pick out neutron signal “candidates” basedon their hit counts. However, both signal and dark hits have comparable hit counts in O(100)ns, so we might end up picking up an enormous amount of dark hits as candidates, if we wereto use the detector’s global time. Therefore for each event, we estimate the photon production
time instead, by subtracting from the PMT hit time the time-of-flight(ToF) a photon takes totravel from its creation vertex to the tank wall. The estimated photon production time is also
called the ToF-corrected hit time or the residual hit time. With vertex resolution good enough,the signal hits have residual hit time spread of O(10) ns, much smaller than the global hit time
spread of O(100) ns. Hence working with the ToF-corrected hit time reduces the number ofnoise hits, most of which have no positional dependence. If a particular set of hits recorded
within O(10) ns of residual hit time has more hits than a pre-defined threshold, then the hits
are taken as a neutron signal candidate. This process is illustrated in Figure 3.7.

3.2. Stage 1: Hit-cluster search

31

Figure 3.7: Sample distributions of raw PMT hit times of signal (top) and the residual, i.e.,
ToF-corrected, hit times of signal (center) and noise (bottom). The ToF-corrected hit times
show a much sharper peak for signal hits that are generated close to a specific vertex within
the tank volume, whereas they don’t for noise hits which mostly occur randomly near the tankwall, and thus lack position dependence. For this reason, a sliding time window searches for any
cluster of hits more than a pre-defined threshold (e.g., the red dotted line) in the sorted arrayof ToF-corrected hit times. Any found set of hits are saved as a neutron capture “candidate”,
and most of the candidates are accidental noise selected as the example shown in the bottom
figure. With the search conditions stated in Table 3.3, the signal-to-noise ratio is ⇠1:30 asshown in Figure 3.12.
To prevent double-counting a single signal, a parameter TMINPEAKSEP is set as the minimum time di↵erence between any two consecutive candidates. If two neighboring candidates
are separated within TMINPEAKSEP, the candidate with more hits replace the other. TMINPEAKSEP should be significantly larger than the resolution (⇠20 ns) shown in Figure 3.9.

32

Chapter 3. The neutron detection algorithm

The performance of the hit-cluster search was optimized and evaluated with the same simulation
configuration as in Section 3.1. Since applying a neutrino vertex fitter requires a charged lepton
for a primary trigger, for simplicity we smeared the true neutron vertex to simulate the vertex
fit, rather than actually applying a dedicated fitter. Each of the true neutron vertex coordinates
was smeared by a Breit-Wigner distribution of= 7. Smearing with the given

value resultsin vertex resolution of ⇠18 cm, which is comparable to that of fiTQun [28], whose 1 vertex
resolution for sub-GeV 1-ring muon neutrino events is ⇠18 cm. This is well within the ⇠50 cm
vertex resolution of Cherenkov-radiating electrons shown in Figure 3.6. The smeared neutron
vertices are as illustrated in Figure 3.8.

Figure 3.8: The coordinate resolution (left) and the distance resolution (right) of the smeared
neutron vertices. (area-normalized) The maximum distance from the smeared vertex to thetrue vertex was limited to 150 cm.
There are several tunable parameters within this stage. One of the most critical ones that
directly a↵ects the cluster search is the time width of the search window that slides throughthe sorted array of ToF-corrected hit times. Let us call this parameter the “search window
width.” With a wider width, more signal hits as well as noise hits are likely to be included
within the search window. From Figure 3.9, we expect the optimal value to be within O(10) ns,
however, it had not been optimized quantitatively ever since the algorithm was designed. So
far, the width of 10 ns had been used since the value is optimal for containing electron events,
and thus commonly used within the SK collaboration. In this study, the search window width
was tuned to 14 ns, which is close to the peak widths shown in Figure 3.9 and shows higher
signal-to-noise significance for p(n, ) signals according to Figure 3.10.

3.2. Stage 1: Hit-cluster search

33

Figure 3.9: The linear-scale (left) and log-scale (right) residual (ToF-corrected) hit time distribution for signal hits. (area-normallized) The solid lines are the results of ideal vertex resolution(fit neutron vertex = true neutron vertex) and the dotted lines are those of realistic resolution(1 resolution being ⇠18 cm). Both ideal and realistic resolutions show almost identical peakwidths of ⇠14 ns. The robustness to the vertex resolution within 50 cm is as expected by Figure
3.6. The tail following the peak should represent photons traveling longer than the shortest
possible path due to reflection on the tank wall and the glass bulbs covering PMTs.

p
Figure 3.10: The estimated significance S/ S + B versus the search window width, for p(n, )signals (top) and Gd(n, ) signals (bottom). The solid line is for ideal and the dotted line is forrealistic resolution similar to fiTQun (⇠18 cm). The noise hits were assumed to be perfectly
uniform with the rate of 0.08 hits/ns, which is a good approximation of the noise hit rate
within SK6. The maximum significance was achieved with the search window width of ⇠13 ns
for p(n, ) reaction and ⇠30 ns for Gd(n, ) reaction. The overall significance was maximized
with the width of ⇠19 ns.

34

Chapter 3. The neutron detection algorithm

Figure 3.11: The raw (left) and area-normalized (right) NHits distributions for noise, p(n, ),and Gd(n, ) candidates.

Figure 3.12: Mean number of candidate class per event.

3.2. Stage 1: Hit-cluster search

35

The second parameter of importance is the number-of-hits (NHits ) threshold which should becorrelated with the search window width. The NHits distributions for the three candidate
types are shown in Figure 3.11. Determining the threshold is again a compromise between
signal purity and detection efficiency. The threshold that has been used so far within SK is 7.
With the search window width fixed at 14 ns, three threshold values, 5, 6, and 7 were tested
in hopes of getting higher efficiency, and the summary of the test results is shown in Table
3.2. Lower thresholds show dramatically higher selection efficiency for Stage 1, especially forp(n, ) captures. However, the number of noise candidates, which directly translates to thecomputational load, increases significantly as well. Besides, signal candidates with a too small
number of hits are often indistinguishable from noise candidates, and therefore they are likely
to be rejected by the MLP that is used to separate signals and noise in the following Stage 2. If
most of the signal candidates with a small number of hits were to be rejected in the end, there
would be no reason to lower the threshold after all, which would only sacrifice signal purity andincrease the computational load exponentially.
To verify the validity of lowering the NHits threshold, a sample MLP with the same architecture
as the one used in Stage 2 (explained in Section 3.3) was made with Keras [29], and was trainedwith candidates selected with di↵erent threshold values. The top three rows in Table 3.2 show
the selection results of the threshold values 5, 6, and 7 from a 3,000-event single neutron MC.
The selection efficiencies of both Stage 1 and 2 combined listed in Table 3.2 are the ones thatgive the noise ratio of 2% within the selected candidates. The overall selection efficiency ishigher for lower threshold, so we can conclude that lowering the threshold definitely enhances
detection efficiency, especially for p(n, ) signals. However, whether it is worth computing ⇠15times longer amount of time to get a few percent increased amount of signals boils down to
one’s personal preference. Also, the enhanced efficiency with lower thresholds may merely bea result of larger input statistics, i.e., the much larger amount of noise candidates taken in one
event.
To check the e↵ect that the amount of statistics has on the final selection efficiency, an additional
test set of candidates were taken from a 300,000-event single neutron MC with threshold valueof 7. The results are listed in the final row of Table 3.2, and they are almost comparable to the36

Chapter 3. The neutron detection algorithm

results with threshold 6. Hence, it is evident that the lower selection efficiency that comes witha higher threshold value can be compensated by increasing the quantity of input statistics. In
fact, event-wise parallel computing enables us to process 300,000-event MC as fast as 3,000event MC, as long as we can split the job into 100 pieces and run them in parallel. For thisreason, the threshold was kept as 7 in this study. However if the detection efficiency for thep(n, ) channel is crucial as with pure water, one would have to consider lowering the thresholdby the expense of the computation time.
No.
1
2
3
4

NHits cut
NHits 5
NHits 6
NHits 7
NHits 7

Stage 1 e↵. (%)
Gd
HAll
99.4 82.9 91.8
99.3 71.4 86.4
99.3 58.8 80.6
99.0 57.5 79.3

NCandidates
Noise
Signal
1,477,453
2,753
382,762
2,592
96,1602,418
9,660,836 237,913

Thr.
0.767
0.740
0.425
0.740

Stage 1&2 e↵. (%)
Gd
HAll
97.4 26.4 66.0
96.9 20.6 62.2
95.8 17.0 58.7
96.7 24.4 62.6

Table 3.2: The summary of the performance test of di↵erent NHits thresholds. All test sets 1-4
share the same MC configuration as explained in Section 3.1 and 3.2. Test sets 1-3 consist of
3,000 neutron MC events each. The Stage 2 selection efficiencies are obtained with Keras [29]using the same MLP architecture as explained in Section 3.3 for reference. The signal likelihood
thresholds (Thr.) and the combined selection efficiencies are the ones that give 2% noise ratio
within the final selected candidates. Test set 4 with 300,000 neutron MC events is included forcomparison with test set 3. All parameters other than the NHits threshold are as in Table 3.3.
Another important parameter that needs to be tuned is the PMT dead-time. It is to reduce
the PMTs’ random burst noise which is known to occur by a higher chance to a PMT with aphoton hit rather than a PMT without any. This rather peculiar phenomenon was reported
quite recently within the SK collaboration [3]. This noise can be reduced by imposing hypothetical dead-time and removing all hits within the dead-time from the selected candidates.
The previous studies [3, 4] suggest 3-6 µs as appropriate values. However in this study, thePMT dead-time was neglected and set to 0, as the search window width had been optimized
beforehand without considering this parameter. The full optimization is left for future work.
Table 3.3 lists the parameter values used for the evaluation of the algorithm. For evaluation,
the algorithm was applied to 300,000 cold (p = 1 keV/c) neutron MC events with random
neutron vertex within the detector’s fiducial volume. With the parameter values listed in Table
3.3, the signal selection efficiency of Stage 1, the hit-cluster search, was ⇠99% for Gd(n, )signals, ⇠58% for p(n, ) signals, and ⇠80% overall, as shown in Table 3.2.

3.2. Stage 1: Hit-cluster search

No.
1
2
3
4
5
6
7
8
9

Parameter
TWIDTHNHITSTHNHITSMX
N200MX
T0TH
T0MX
TMINPEAKSEPPVXRESTRBNWIDTH

Description
Search window width for hit-clusters
Lower threshold for NHits within 14 ns
Upper threshold for NHits within 14 ns
Upper threshold for NHits within 200 ns
Lower threshold for residual hit time
Upper threshold for residual hit time
Minimum time di↵erence between candidates
in Breit-Wigner vertex smearing
PMT dead-time

37

Unit
ns
ns
µs
µs
ns
cm
µs

Value
14
7
70
200
1
535
60
7
0

Table 3.3: The list of parameter values used for Stage 1 (NHits -based candidate search) of theneutron detection algorithm. For reference, the primary particle generation time is 1 µs in the
global simulation time frame, and the time range of each data event is [-5, -535] µs. The residual
hit time range of [1, 535] µs corresponds to the full time range following neutron emission.

38

Chapter 3. The neutron detection algorithm

3.3

Stage 2: Binary classification with MLP

The second step is to extract features from the selected candidates and statistically find eachcandidate’s signal likelihood based on its features. If the candidate’s signal likelihood is greaterthan the pre-defined threshold, we say that the candidate is “tagged.” As for the multivariate analysis framework, this study used the Multilayer Perceptron (MLP) method, commonly
known as the “artificial neural network,” which is also currently utilized within the SK experiment to detect neutrons from atmospheric and beam neutrino events. If the user inputs
distinguishing features between signal and noise to an MLP, then following the pre-defined algorithm (in this study, backpropagation), the network will optimize itself to output the predictionthat gives the highest efficiency or the least loss. This process is called supervised learning.
The working principles of the MLP method can be found in Ref. [30].
The Neut-Fit algorithm [1] was used to estimate the origin vertex of signal photons. For eachcandidate, the algorithm searches for a vertex that minimizes the variance of the ToF-corrected
hit times by looping over a fixed number of grid points that evolve increasingly finer after
each loop. Specifically for a p(n, ) signal, this algorithm should give nearly ideal electron
vertex resolution in theory, if the candidate does not include any noise hit.2 For this reason,
past studies [1, 3, 4] on pure water neutron tagging have utilized Neut-Fit vertices as well
as neutron fit vertices in ToF correction. However for Gd(n, ) signals, the algorithm is notexpected to work as well, since the Gd capture signals mostly originate from multiple electron
vertices. (refer to Figure 3.2 or 3.3) Therefore, the number of grid points per each coordinate
axis was reduced by 80 % (the initial grid width was changed from 100 cm to 500 cm), and thetarget grid width was changed from 0.5 cm to 20 cm. This change has reduced the computation
time by the factor of O(10 2 ) without altering the overall detection performance. The 1 vertex
resolution of Neut-fit for 2.2 MeV gamma-ray was ⇠270 cm for both before and after the change.
2

In reality, most candidates have a couple of noise hits as in the center plot of Figure 3.7. The Neut-Fit
algorithm is highly sensitive to noise hits, and the vertex resolution is in general much worse (O(100) cm)
than that of neutrino vertex fitters (O(10) cm). Nonetheless, Neut-fit is currently the only option available for
capture vertex estimation of p(n, ) candidates, as the target NHits of low energy neutrino vertex fitter is muchlarger (⇠60 Hits for ⇠10 MeV ⌫e events) than typical NHits of p(n, ) candidates. For Gd(n, ) candidates, theviability of applying the likelihood-based low energy neutrino vertex fitter BONSAI [31] has been demonstrated
in a past fake-data study with 0.1 wt.% Gd concentration [2].

3.3. Stage 2: Binary classification with MLP

3.3.1

39

Features (input variables to MLP)

16 features are selected based on the past studies [1, 3, 4], and are extracted from each candidate
to be used as inputs to the MLP. Distributions are shown in Figure 3.13, and their correlation
matrices are as illustrated in Figure 3.14. The candidates are classified into noise, p(n, ),and Gd(n, ) candidates by comparing the reconstructed capture time (ReconCT) and thetrue neutron capture time within a simulated event. The reconstructed capture time of acandidate is defined as the average of the maximum and minimum ToF-corrected hit times of
the candidate. Any candidate whose ReconCT is within 50 ns from the true neutron capture
time, or the capture gamma-ray production time in the simulation time frame, is taken as signal
candidates.3 The following is the description and categorization of the features according totheir definitions:• NHits -related:– NHits: The number of hits within the residual search window width (14 ns).
NHits should be one of the most strongly separating features. Most of the Gd(n, )signals and some of the p(n, ) signals could be separated by a simple cut on NHits,for example, NHits > 15. NHits for noise candidates decreases exponentially.
– N200: The number of hits within (-100, +100) ns window centered at ReconCT.
N200 has separation power on candidates that are selected due to accidental increase of noise hits (radioactivity, spallation, etc.) within proximity. N200 for noisecandidates show much weaker correlation with NHits than for signal candidates.
• Timing-related:– TRMS: Standard deviation of residual hit times (within 14 ns).
Noise hits are mostly random in time, and thus they generally have a larger TRMS
than signal hits. Gd(n, ) candidates show a sharper distribution, since they havemuch larger signal-to-noise ratio than p(n, ) candidates.
3

This strategy that has been used so far within SK is convenient, however, is a source of error as noisecandidates coincident with a neutron capture can be classified as signals. It should better be replaced with a
hit-based alternative in the future.

40

Chapter 3. The neutron detection algorithm• Vertex-related:– prompt nfit: Distance from the (smeared) neutron fit vertex to the Neut-Fit vertex.
Neut-fit vertices of signal candidates should be mostly within 500 cm distance fromthe neutron vertex, considering the resolution of Neut-fit (1 : 270 cm for 2.2 MeV
gamma-ray). Candidates with smaller prompt nfit are more likely to be signals, andthose with very large prompt nfit (> 1,000 cm) are most likely to be noise.
– DWall: (Shortest path) distance from the neutron fit vertex to the tank wall.
A neutron fit vertex is technically a vertex that is used for ToF subtraction. DWall
itself has no separation power whatsoever, as a neutrino vertex fitter, unlike NeutFit, gives out only one fixed vertex per event. However, DWall combined withother vertex-related features should provide a better picture of the neutron vertex
to the MLP. For example, if a candidate has a small DWall and a relatively large
DWallMeanDir (explained below), the candidate is more likely to be a noise since thephotons that have hit the detector should have traveled further against attenuation
than DWall. Including DWall as a feature has indeed increased the overall detection
efficiency by ⇠1% (abs).
– DWall n: Distance from the Neut-fit vertex to the tank wall.
DWall n has high separation power for noise candidates, many of which have NeutFit vertices on the tank wall (DWall n = 0). DWall n for signal candidates shows
stronger correlation with DWall than for noise candidates. Unfortunately, DWall ncannot replace DWall, since the resolution of Neut-fit is much worse (O(100) cm)
than that of dedicated lepton vertex fitters (O(10) cm).
– DWallMeanDir: DWall in the averaged hit direction.
DWallMeanDir should be smaller for signal candidates whose PMT-hitting photons
have to travel against attenuation. Many noise candidates, however, are a set of
spontaneous dark hits coming from the PMTs themselves, and thus they tend tohave a larger DWallMeanDir. DWallMeanDir is negatively correlated with PMT
acceptance (refer to [1, 4] for details), or the probability of PMTs to be hit.

3.3. Stage 2: Binary classification with MLP

41– ThetaMeanDir: Mean of the angles made by a vector from the neutron fit vertex
to each hit PMT and a unit vector in the averaged hit direction.
ThetaMeanDir of p(n, ) candidates should have its peak close to 42 , which is the
Cherenkov angle of water. ThetaMeanDir for Gd(n, ) candidates does not sharethe same physical meaning due to the gamma-ray multiplicity of Gd(n, ) reaction.
Noise hits are often clustered locally, and thus noise candidates have a tail in thelower end of the distribution. However, the 42 peak of noise candidates had notbeen identified in the past study [1], and further investigation is needed.
• Angularity-related:– AngleMean, AngleStdev, AngleSkew: Mean, sample standard deviation, andskewness of the opening angles of all possible 3-hit combinations.
Any 3-hit combination has an opening angle, and there exist

NHits
3

possible 3-hit

combinations for a candidate with NHits hits. The features AngleMean, AngleStdev,
and AngleSkew are all statistics of the opening angle distributions of a candidate.
The introduction of these features is inspired by one of the past studies [4] that
has introduced Cherenkov-likelihood of PMT hits. Its separability is based on thediscrepancy between the opening angle distributions of signal and noise candidates.
Noise hits are more likely to be either clustered locally or distributed randomly, whichmakes tails at ends. AngleMean shows very strong correlation with ThetaMeanDir.
– Beta(l): Angular isotropy parameters of hits.
The Beta(l) parameter is defined asl =2NHits (NHits

X1) i,jPl (cos ✓ij )

(3.1)

i6=j

where Pl is the Legendre polynomial of degree l (1  l  5), and ✓ij is the angle
between two vectors from the neutron fit vertex to the i-th and j-th hit PMTs. In
general, lower values of the Beta(l) parameters indicate greater angular isotropy ofhits. These parameters up to scaling are the estimated sum of the spherical harmon-

42

Chapter 3. The neutron detection algorithm
ics expansion coefficients, which must be rotationally invariant. They characterize
the angular distribution of the PMT hits, and have shown separation power between
signal and non-signal hits in the SNO experiment [32]. For that reason, they havebeen adopted to neutron detection in the SK experiment by the past studies [2, 4].
These parameters show distinct distributions between signal and noise candidates,
however, it seems like the number of signal-like noise candidates has increased compared to the previous work [4].
The correlations that Beta(l) parameters have with other features are interesting.
Beta1 has very strong correlation with ThetaMeanDir, since Beta1 is equivalent
to a cosine of ThetaMeanDir. All Beta(l) parameters are mutually correlated fornoise candidates, but not as much for signal candidates. Beta(l) parameters are alsostrongly correlated with the opening angle statistics, especially with AngleMean andAngleStdev, and moderately correlated with prompt nfit. The correlation is smaller
for signal candidates. As in the SNO experiment [32], a linear combination of a
selection of Beta(l), such as

14 :=1 + 4 4 , might be enough for an input.

3.3. Stage 2: Binary classification with MLP

43

Figure 3.13: The MLP input feature distributions of noise candidates (gray), true Gd(n, )candidates (pink), and true p(n, ) candidates (cyan) (area-normalized). 5,000 candidates eachare randomly sampled candidates a cold (p =1 keV/c) neutron MC with the MC configuration
explained in Section 3.1.

44

Chapter 3. The neutron detection algorithm

NHits

TRMS

N200

prompt_nfit

DWall

DWall_n

DWallMeanDir

ThetaMeanDir

AngleMean

AngleStdev

AngleSkew

Beta1

Beta2

Beta3

Beta4

NHits

1.00

-0.08

0.34

0.08

-0.03

-0.08

-0.08

-0.21

-0.20

0.26

0.03

0.21

0.28

0.29

0.29

Beta5

0.28

TRMS

-0.08

1.00

0.01

-0.03

0.01

0.02

0.02

0.06

0.05

-0.07

0.00

-0.06

-0.08

-0.09

-0.09

-0.10
0.21

N200

0.34

0.01

1.00

0.12

0.00

-0.05

-0.11

-0.19

-0.15

0.19

-0.01

0.18

0.24

0.23

0.22

prompt_nfit

0.08

-0.03

0.12

1.00

-0.02

-0.27

-0.03

-0.51

-0.42

0.36

0.13

0.50

0.57

0.57

0.56

0.56

DWall

-0.03

0.01

0.00

-0.02

1.00

0.71

-0.43

0.27

0.29

-0.04

-0.23

-0.30

-0.07

-0.02

-0.02

-0.02

DWall_n

-0.08

0.02

-0.05

-0.27

0.71

1.00

-0.27

0.36

0.32

-0.18

-0.17

-0.37

-0.24

-0.21

-0.20

-0.20

DWallMeanDir

-0.08

0.02

-0.11

-0.03

-0.43

-0.27

1.00

-0.16

-0.22

-0.10

0.30

0.19

-0.12

-0.17

-0.16

-0.13

ThetaMeanDir

-0.21

0.06

-0.19

-0.51

0.27

0.36

-0.16

1.00

0.81

-0.55

-0.40

-0.98

-0.76

-0.67

-0.66

-0.66

AngleMean

-0.20

0.05

-0.15

-0.42

0.29

0.32

-0.22

0.81

1.00

-0.38

-0.69

-0.83

-0.60

-0.54

-0.54

-0.54

AngleStdev

0.26

-0.07

0.19

0.36

-0.04

-0.18

-0.10

-0.55

-0.38

1.00

-0.08

0.49

0.66

0.59

0.59

0.58

AngleSkew

0.03

0.00

-0.01

0.13

-0.23

-0.17

0.30

-0.40

-0.69

-0.08

1.00

0.46

0.13

0.06

0.08

0.09

Beta1

0.21

-0.06

0.18

0.50

-0.30

-0.37

0.19

-0.98

-0.83

0.49

0.46

1.00

0.73

0.63

0.62

0.62

Beta2

0.28

-0.08

0.24

0.57

-0.07

-0.24

-0.12

-0.76

-0.60

0.66

0.13

0.73

1.00

0.87

0.82

0.81

Beta3

0.29

-0.09

0.23

0.57

-0.02

-0.21

-0.17

-0.67

-0.54

0.59

0.06

0.63

0.87

1.00

0.90

0.86

Beta4

0.29

-0.09

0.22

0.56

-0.02

-0.20

-0.16

-0.66

-0.54

0.59

0.08

0.62

0.82

0.90

1.00

0.92

Beta5

0.28

-0.10

0.21

0.56

-0.02

-0.20

-0.13

-0.66

-0.54

0.58

0.09

0.62

0.81

0.86

0.92

1.00

Beta5

NHits

TRMS

N200

prompt_nfit

DWall

DWall_n

DWallMeanDir

ThetaMeanDir

AngleMean

AngleStdev

AngleSkew

Beta1

Beta2

Beta3

Beta4

NHits

1.00

-0.16

0.55

-0.27

-0.10

-0.06

-0.10

0.00

-0.12

0.16

0.07

0.08

0.02

-0.01

0.01

0.01

TRMS

-0.16

1.00

-0.03

0.04

0.04

0.04

0.15

0.15

0.11

-0.09

-0.08

-0.18

-0.07

0.01

0.00

-0.01

N200

0.55

-0.03

1.00

-0.09

-0.05

-0.03

-0.01

0.00

-0.06

0.08

0.04

0.04

0.01

0.01

0.02

0.03

prompt_nfit

-0.27

0.04

-0.09

1.00

0.00

-0.16

0.04

-0.26

-0.18

0.04

0.10

0.26

0.19

0.12

0.09

0.06

DWall

-0.10

0.04

-0.05

0.00

1.00

0.88

0.10

0.02

0.05

-0.03

-0.03

-0.05

-0.03

0.00

0.02

0.01

DWall_n

-0.06

0.04

-0.03

-0.16

0.88

1.00

0.12

0.05

0.06

-0.04

-0.03

-0.07

-0.06

-0.02

0.00

0.00

DWallMeanDir

-0.10

0.15

-0.01

0.04

0.10

0.12

1.00

-0.02

-0.15

-0.13

0.16

0.03

-0.11

0.00

0.02

0.01

ThetaMeanDir

0.00

0.15

0.00

-0.26

0.02

0.05

-0.02

1.00

0.76

-0.38

-0.55

-0.97

-0.41

-0.11

-0.05

-0.07

AngleMean

-0.12

0.11

-0.06

-0.18

0.05

0.06

-0.15

0.76

1.00

-0.20

-0.82

-0.79

-0.27

-0.05

-0.03

-0.07

AngleStdev

0.16

-0.09

0.08

0.04

-0.03

-0.04

-0.13

-0.38

-0.20

1.00

-0.07

0.35

0.52

0.12

0.05

0.02

AngleSkew

0.07

-0.08

0.04

0.10

-0.03

-0.03

0.16

-0.55

-0.82

-0.07

1.00

0.60

-0.02

-0.19

-0.04

0.01

Beta1

0.08

-0.18

0.04

0.26

-0.05

-0.07

0.03

-0.97

-0.79

0.35

0.60

1.00

0.41

0.08

0.05

0.08

Beta2

0.02

-0.07

0.01

0.19

-0.03

-0.06

-0.11

-0.41

-0.27

0.52

-0.02

0.41

1.00

0.37

0.11

0.07

Beta3

-0.01

0.01

0.01

0.12

0.00

-0.02

0.00

-0.11

-0.05

0.12

-0.19

0.08

0.37

1.00

0.40

0.15

Beta4

0.01

0.00

0.02

0.09

0.02

0.00

0.02

-0.05

-0.03

0.05

-0.04

0.05

0.11

0.40

1.00

0.45

Beta5

0.01

-0.01

0.03

0.06

0.01

0.00

0.01

-0.07

-0.07

0.02

0.01

0.08

0.07

0.15

0.45

1.00

Beta5

NHits

TRMS

N200

prompt_nfit

DWall

DWall_n

DWallMeanDir

ThetaMeanDir

AngleMean

AngleStdev

AngleSkew

Beta1

Beta2

Beta3

Beta4

NHits

1.00

-0.20

0.89

-0.35

-0.15

-0.14

-0.10

-0.10

-0.22

0.21

0.17

0.20

0.07

0.00

0.06

0.07

TRMS

-0.20

1.00

-0.14

0.14

0.06

0.06

0.12

0.05

0.07

-0.07

-0.05

-0.08

-0.03

0.01

-0.01

-0.02

N200

0.89

-0.14

1.00

-0.26

-0.12

-0.12

-0.06

-0.10

-0.21

0.19

0.17

0.19

0.07

0.00

0.07

0.07

prompt_nfit

-0.35

0.14

-0.26

1.00

0.04

-0.01

-0.01

-0.22

-0.16

0.05

0.13

0.19

0.13

0.05

0.01

0.00

DWall

-0.15

0.06

-0.12

0.04

1.00

0.97

0.12

-0.02

0.02

-0.05

0.00

-0.02

-0.03

-0.01

0.02

0.02

DWall_n

-0.14

0.06

-0.12

-0.01

0.97

1.00

0.15

-0.01

0.02

-0.05

0.01

-0.02

-0.04

-0.01

0.02

0.02

DWallMeanDir

-0.10

0.12

-0.06

-0.01

0.12

0.15

1.00

0.09

0.01

-0.11

0.02

-0.08

-0.07

0.00

0.03

0.02

ThetaMeanDir

-0.10

0.05

-0.10

-0.22

-0.02

-0.01

0.09

1.00

0.90

-0.61

-0.78

-0.96

-0.25

-0.04

-0.06

-0.08

AngleMean

-0.22

0.07

-0.21

-0.16

0.02

0.02

0.01

0.90

1.00

-0.54

-0.91

-0.94

-0.22

-0.03

-0.06

-0.10

AngleStdev

0.21

-0.07

0.19

0.05

-0.05

-0.05

-0.11

-0.61

-0.54

1.00

0.29

0.58

0.60

0.02

0.07

0.03

AngleSkew

0.17

-0.05

0.17

0.13

0.00

0.01

0.02

-0.78

-0.91

0.29

1.00

0.83

-0.03

-0.14

0.06

0.07

Beta1

0.20

-0.08

0.19

0.19

-0.02

-0.02

-0.08

-0.96

-0.94

0.58

0.83

1.00

0.29

0.04

0.07

0.10

Beta2

0.07

-0.03

0.07

0.13

-0.03

-0.04

-0.07

-0.25

-0.22

0.60

-0.03

0.29

1.00

0.20

0.05

0.06

Beta3

0.00

0.01

0.00

0.05

-0.01

-0.01

0.00

-0.04

-0.03

0.02

-0.14

0.04

0.20

1.00

0.21

0.06

Beta4

0.06

-0.01

0.07

0.01

0.02

0.02

0.03

-0.06

-0.06

0.07

0.06

0.07

0.05

0.21

1.00

0.28

Beta5

0.07

-0.02

0.07

0.00

0.02

0.02

0.02

-0.08

-0.10

0.03

0.07

0.10

0.06

0.06

0.28

1.00

Figure 3.14: Correlation matrices for noise (top), p(n, ) (center), and Gd(n, ) candidates(bottom).

3.3. Stage 2: Binary classification with MLP

3.3.2

45

Training the MLP

Modeling and executing the binary classification with MLP was done with TMVA [33] (Toolkit
for Multivariate Data Analysis with ROOT) v4.2.0 [34]. This is a package developed by agroup of physicists during the early 2000s. It is very convenient if the data is in ROOT format.
The version 4.2.0 has been the latest standalone version that is available since the program
was integrated into ROOT in 2013. Although the convergence is quite slower than the latest
optimizers such as Adam [35], it shows similar performance as Keras [29] or the latest TMVA
that is maintained by ROOT. The MVA::MLP class of TMVA was used to model the MLP. The
options used are specified in Table 3.4. The hyperparameters have been tuned heuristicallyrather than systematically.
The binary cross entropy (BCE) is used as the loss function. The BCE is defined as
N

BCE =

1 X trueylog p(yi ) + 1
N i=1 iyitrue log (1p(yi ))(3.2)

where N is the number of candidates, yitrue is the true class (label), and p(yi ) is the output
node value, i.e., the signal likelihood of the i-th candidate. The output signal likelihood p(yi )must be restricted within the range [0, 1], as we use the sigmoid f (x) = 1/(1 + e x ) 2 [0, 1] as
the activation function of all nodes. The change of loss versus epochs is shown in Figure 3.15.
The convergence within 1% loss variance was achieved after ⇠200 epochs with the learning anddecay rates specified in Table 3.4.

Figure 3.15: BCE versus epochs for test and train data sets.

46

Chapter 3. The neutron detection algorithm
OptionWeighting
Variable transformation
Activation function
EpochsArchitecture
Loss function
Optimizer
Learning rate
Decay rate
Train/test split
UseRegulator

Used values
None
Normalization
Sigmoid500
16:17:15:1
Binary cross entropy
Backpropagation
0.02
0.01
50:50 (randomly sampled)
True

Table 3.4: TMVA’s MLP options used for training.
The importance of the features is illustrated in Figure 3.16. The importance of the i-th input
variable Ii in TMVA [34] is defined as the sum of the squared weights between the variable andthe nodes in the first hidden layer:(1)

Ii = x̄2i

nnodes ⇣

X
j

(1)

!ij

⌘2

(3.3)

(1)

where x̄i is the sample mean of the i-th input variable, nnodes is the number of nodes in the first(1)

hidden layer, and !ij is the weight of the connection between the i-th input variable and thej-th node in the first hidden layer. The opening angle statistics AngleMean, AngleStdev, andAngleSkew which were introduced to replace the Cherenkov-likelihood used in the past study[4] showed the weakest importance.

Figure 3.16: Feature importance.

3.4. Evaluation

3.4

47

Evaluation

Evaluation of the combined (Stage 1&2) neutron detection algorithm was estimated using single
cold (pn = 1 keV/c) neutron MC (without any other particle generation). The test data set
was made by randomly sampling 50% of signal and noise candidates each from 300,000 MC
events. The other half is used for training MLP, which is explained in Subsection 3.3.2. The
numbers of total noise candidates and total signal candidates in the test data set were 4,748,992
and 116,677, respectively. The neutron fit vertex, i.e., the estimated neutron emission vertex
was smeared from a true neutron vertex using the smearing condition explained in Section 3.2.
Noise hits were taken from randomly sampled events taken at the end of the SK6 commissioning
phase (Run 85200 - 85245). The performance of the algorithm is summarized in Figure 3.17.

Figure 3.17: The signal likelihood distributions (top) and the MLP performance (bottom).
The results are as expected; most noise candidates have signal likelihood close to 0, while
signal candidates have larger signal likelihood in general. The p(n, ) candidates show smaller
signal likelihood than the Gd(n, ) candidates. The bottom figure represents the performanceof the trained MLP according to the threshold values set on signal likelihoods. For example, if

48

Chapter 3. The neutron detection algorithm

we tag candidates with signal likelihood larger than 0.5, we reject about a half of the p(n, )candidates selected in Stage 1, while almost all Gd(n, ) candidates from Stage 1 are tagged.
The noise ratio within the tagged candidates, i.e., candidates with signal likelihood larger than
the threshold, should be close to 4 %.
To evaluate the performance of Stage 1 and 2 of the detection algorithm, we define the neutron
capture detection efficiency ✏(n, ) as:✏(n, ) =(# of true tagged neutron captures)(# of true neutron captures)(3.4)

The capture detection efficiency ✏(n, ) in principle is di↵erent from the neutron detection efficiency ✏n , whose denominator should be all neutrons including those that are not captured to
produce gamma-rays. Since the fraction of such uncaptured neutrons is very low (less than 2% according to the models used in this study), ✏(n, ) is essentially ✏n . All efficiencies on MC
data set shown in this chapter are ✏(n, ) values.
The summary of the trained MLP, for the signal likelihood threshold of 0.7, is given in Table
3.5. The two rows represent di↵erent ReconCT time range used for the signal search. Of course
it is better to have a wider search time range, but in actual neutron data the first few µs areoften contaminated with triggering light coming from recoil leptons or scintillation. Hence, theneutron search time range of [1, 535] µs is practically applicable to MC events only.
The overall neutron detection efficiency of the combined detection algorithm was estimated tobe ⇠60% (with noise ratio rNoise ⇠2 %) for 0.0110 wt.% Gd concentration.
Time range[1, 535] µs[4, 535] µs

Stage 1 e↵. (%)
Gd
HAll
99.0 57.6 79.4
95.3 56.8 77.0

Stage 2 e↵. (%)
Gd
HAll
98.0 45.0 78.8
96.6 44.9 78.6

Combined (%)
GdHAll rNoise95.6 25.9 62.6 2.01
92.1 25.5 60.5 2.06

Table 3.5: Summary of the capture selection efficiencies ✏(n, ) of the neutron detection algorithm
used in this study. rNoise indicates the noise ratio within the tagged candidates.

Chapter 4

Testing with a neutron source

The neutron detection efficiency of a detector should be a convolution between the detection
algorithm’s performance and the detector’s physical capability.1 In the performance evaluation
explained in the previous chapter, we have replaced the former with the semi-realistic simulationequipped with recently developed physics models and simulation parameters tuned to the purewater phase. We have to demonstrate that it actually reproduces the detector’s response.

A traditional approach to this has been putting a point neutron source inside the tank andcompare the detector’s response between the actual data and its simulation. We have followed
this approach and recorded data since the beginning of the Gd phase. Then, the neutron
detection algorithm was applied to both data and MC to compare the detection efficiencies.
The neutron calibration campaign is still ongoing.

The following subsections will explain the details of the neutron source data and its simulation,compare various distributions of the two, and discuss the sources of systematic uncertainties.
1

Apart from the simulation parameters that are irrelevant with Gd-loading, the neutron capture ratio on Gd
and the optical parameters of the Gd-loaded water should be ideally the ones that characterize the detector.
The former is correlated with the capture cross sections, and both are correlated with the Gd concentration.
Tuning these parameters is important, however, is beyond the scope of this study.

49

50

Chapter 4. Testing with a neutron source

4.1

Neutron source data

The source we have employed is the so-called AmBe (Americium-Beryllium) neutron source. It
is one of the alpha neutron sources where decay alpha particles from one element (e.g., 241 Am)
are bombarded onto the other element (e.g., 9 Be) to produce neutrons. The design of this sourceis very simple; it is typically made by encapsulating a mixture of AmO2 and Be powders in athin stainless steel container. Since the alpha particles can only travel a few dozen micrometers
within the compound, the source can be very compact in size (⇠ O(1) cm), and thus be a goodneutron point source for detector calibration. The source was deployed inside the tank using
the Automated Source Deployment System [36].
When the source emits a neutron, it is known to simultaneously emit a 4.4 MeV gammaray, also called the “prompt” gamma-ray, with ⇠60% chance [37]. These gamma-rays let us
know the exact neutron production time, in the same way a recoil lepton does in the inverse
beta decay. The prompt gamma-rays are what make the efficiency measurement possible, andwithout those, we would never be certain if a neutron was actually emitted or not. In order totake advantage of the prompt gamma-rays, the source was put inside a transparent plastic casewith 8 symmetrically positioned Bismuth Germanate (Bi4 Ge3 O12 , BGO) crystals encompassing
the source. The arrangement is shown in Figure 4.1. BGO is an inorganic scintillator commonlyused to detect gamma-rays. If a gamma-ray hits one of the BGO crystals and scatters o↵ boundelectrons, the crystal will scintillate. If the scintillation is strong enough to trigger the detector,
then the detector response within a certain time range will be recorded as one event. The
recorded data should have several sources of contamination, and they must be removed with aproper set of cuts before we compare the data with MC simulation.

4.1.1
241

The source: AmBe + BGO

Am is an alpha-decaying isotope with half-life of ⇠432 years. The intensity of 241 Am of

our source is 97 µCi, or 3.6 MBq [4, 38], and the neutron emission rate was estimated bymultiplying the alpha intensity by the neutron yield per alpha particle. The latter, however,4.1. Neutron source data

51

Figure 4.1: The pictures of the deployed source, reprinted from [1, 4]. The AmBe compoundis encapsulated in a pellet-sized container, with 8 BGO crystals surrounding it.
is highly dependent on the source design, which unfortunately is unknown either. Hence, theneutron yield per alpha particle was assumed to be 6.0 ⇥ 10 5 , based on the other sources [39–41]. The uncertainty was assumed to be 20% according to the reported maximum experimental
yield [42]. The neutron emission rate of the source, estimated using the alpha intensity and theassumed neutron yield per alpha particle, is 216 ± 48 neutrons per second. It is larger thanthe estimate of ⇠163 Hz from a previous study [38] by ⇠1 . The neutron intensity could bebetter estimated with a dedicated measurement.
Table 4.1 shows the estimated gamma-correlated and uncorrelated source neutron intensity,
based on the recommended prompt gamma-ray emission probability of 0.575 ± 4.8% [37], incomparison with the same estimates on the same source by Watanabe et al. [38].
Neutron typeGamma-correlated
Uncorrelated

Estimated intensity
124 ± 28 Bq
92 ± 21 Bq

Watanabe et al. [38]87 ± 17 Bq
76 ± 15 Bq

Table 4.1: The estimated gamma-correlated and uncorrelated neutron intensity of the source.
The assumed uncertainty is ⇠22% for the estimation in this study, and 20% for Watanabe et al.
[38]. The 4.4 MeV gamma-ray intensity should be the same as the correlated neutron intensity.
Figure 4.2 illustrates the measured kinetic energy spectra of the source neutrons. The maximum kinetic energy is ⇠12 MeV for uncorrelated neutrons, and ⇠7 MeV for gamma-correlated
neutrons whose 4.4 MeV energy should be taken away by the prompt gamma-ray.
The details of the neutron emission mechanism is explained in Appendix B.

52

Chapter 4. Testing with a neutron source

Figure 4.2: Reference AmBe neutron kinetic energy spectra for all neutrons (cited in [40],blue), and for gamma-correlated neutrons ([40], orange). The approximated flat 2-8 MeV input
spectrum for correlated neutron simulation is shown in green. The orange and green histograms
are area-normalized, and the blue histogram is normalized in a way that the peak value around
3 MeV match the one in the orange histogram.

4.1. Neutron source data

53

Some of the relevant properties of BGO are summarized in Table 4.2. We can imagine thata prompt gamma-ray hitting a ⇠3-cm-wide BGO crystal would lose most of its energy within
the crystal to scintillate enough amount of light to trigger the detector. The amount of photoelectrons deposited in the ID PMTs within the primary trigger time range (within 1.3 µs after
trigger) is called QISMSK within SK, and its distribution is shown in Figure 4.3.
Light yield
Emission peak
Rise time
Decay time
Radiation length

8,400 photons / MeV
480 nm< 2 ns
300 ns1.1 cmTable 4.2: Relevant properties of BGO at room temperature. [43]Figure 4.3: The QISMSK distribution of the neutron source data events for Run 85386. Eventswith QISMSK > 2,000 are mostly triggered by cosmics. Events with QISMSK ⇠ 1,000 are
likely to have been triggered by BGO scintillation.

At the moment, we lack a proper BGO scintillation simulator which enables us to identify
and analyze the composition of the QISMSK peaks. We can only infer from the data itself
to identify the peaks qualitatively. Figure 4.4 shows the average pulse shape of data eventswith selected peaks. The pulse shapes of the residual hit time distributions were fitted with a
two-component pulse modeling function

f (t) = ⇥(t

T0 )

X

i2{fast,slow}



ai e

(t T0 )/⌧d,i✓

1
1 + e (t T0 5⌧r )/⌧r

1

◆

+b

(4.1)

54

Chapter 4. Testing with a neutron source

Figure 4.4: The averaged ToF-corrected hit time distributions near the primary trigger time
range, for events with 800 < QISMSK < 1250 (top) and events with QISMSK < 400 (bottom).
The top histogram is scaled to make the fitted scintillation histogram (orange) in the bottom
figure. The scale factor of 0.11 minimizes the variance of the subtracted (green) bin values withreconstructed hit time larger than 1,200 ns (i.e., gives the “flattest” subtracted histogram).

4.1. Neutron source data55

where the parameters are as described in Table 4.3. Each term in the square bracket represents
the pulse decay and rise for each of the fast and slow pulse components. Since the typical
exponential pulse rise function is replaced with a smoother sigmoid function of the form 1/(1 +e t ), fitting a pulse shape with the function in Equation 4.1 should give a larger fit value forthe rise time ⌧r than the typical exponential fit.
t
T0
⇥
ai
⌧r
⌧d,i
b

ToF-subtracted hit time [ns]Time of first hit on BGO [ns]Heaviside step function
Coefficient for the i-th component
Rise time [ns]Decay time of the i-th component [ns]Baseline number of hits in one bin

Table 4.3: The input t and the fit parameters used in Equation 4.1.

For 800 < QISMSK < 1250, the fitted decay time of ⇠370 ns is in reasonable agreement with
300 to 400 ns decay time of BGO with temperature in between 10 and 20 C [44, 45]. Thus,
we can conclude that scintillation is the dominant trigger for the events with high QISMSK.
However for QISMSK < 400, the same fit with the timing parameters fails to fit the pulse shapeproperly. Rather, the tail fits quite well with the high QISMSK pulse shape scaled by 0.11,
which implies BGO scintillating with either much lower intensity or lower frequency.
The subtracted histogram in the bottom of Figure 4.4 also demonstrates that there may be PMT
hits earlier than scintillation, for events with low QISMSK. We can imagine an event scenario
where a leaking prompt gamma-ray or a few of the Gd-capture gamma-rays are detected by
the PMTs without triggering the detector. If a correlated neutron or the remaining capture
gamma-rays hit the BGOs to boost up the earlier hits and eventually trigger the detector,
such a scenario can mimic the pulse shape of the low QISMSK events, imprinting a sharp
peak with a following signature of scintillation with long decay time. As with the possibilityof scintillation due to neutrons, Vincke et al. [46] reports the AmBe neutron’s energy deposit
within a BGO crystal ⇠4 cm in size, as shown in Figure 4.5. A small amount of neutrons
mimicking the 4.4-MeV scintillation can be a source of systematic uncertainty, which will be
discussed in Subsection 4.4.1.

56

Chapter 4. Testing with a neutron source

Figure 4.5: The FLUKA-simulated AmBe neutron’s energy deposit in BGO (left) and themeasured energy deposit of both source neutrons and gamma-rays (right). Reprinted from [47].

4.1.2

Data specification and reduction

The neutron source data has been taken since the beginning of the Gd phase (July 2020), either
weekly or bi-weekly with di↵erent source positions aligned along the z-axis that passes near
the tank center. A dataset for each run has a ⇠30-minute record of the detector’s response,the most important of which is the hit time and charge information of individual PMTs. Hits
more than 60 within a 200 ns time window would trigger the SHE flag of the detector, and thedetector response within [-5, 35] µs time window around the trigger would be recorded as one
event. The AFT trigger is forced for all SHE triggers. However, for some technical reasons thathave not been clarified yet, only ⇠50% of the SHE-triggered events accompany a subsequent
AFT-recorded event. An AFT-recorded event has an additional 500-µs record of the detector’s
response following the end of the SHE-triggered event. In other words, about a half of thedata events will have ⇠535 µs length record of the detector’s response following the trigger,
while the other half will have only ⇠35 µs length of data. In this study, runs taken after the
Gd-loading with z 2 { 12, 0, +12} [m] were selected and analyzed. The specifications of thetarget runs are summarized in Table 4.4.

4.1. Neutron source data

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

Run
85122
85128
85134
85149
85155
85161
85178
85184
85189
85210
85212
85214
85228
85230
85232
85247
85249
85251
85272
85274
85276
85299
85301
85303
85336
85338
85340
85384
85386
85390
85435
85437
85439
85480
85482
85484

57

Date
x [m]
2020-08-05-3.89
2020-08-05-3.89
2020-08-05-3.89
2020-08-12-3.89
2020-08-12-3.89
2020-08-12-3.89
2020-08-19-3.89
2020-08-19-3.89
2020-08-19-3.89
2020-08-26-3.89
2020-08-26-3.89
2020-08-26-3.89
2020-09-01-3.89
2020-09-01-3.89
2020-09-01-3.89
2020-09-08-3.89
2020-09-08-3.89
2020-09-08-3.89
2020-09-15-3.89
2020-09-15-3.89
2020-09-15-3.89
2020-09-29-3.89
2020-09-29-3.89
2020-09-29-3.89
2020-10-15-3.89
2020-10-15-3.89
2020-10-15-3.89
2020-10-27-3.89
2020-10-27-3.89
2020-10-27-3.89
2020-11-10-3.89
2020-11-10-3.89
2020-11-10-3.89
2020-12-01-3.89
2020-12-01-3.89
2020-12-01
-3.89

y [m]
-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707-0.707

z [m]
-12
0
12
-12
0
12
-12
0
12
-12
0
12
-12
0
12
-12
0
12
-12
0
12
-12
0
12
-12
0
12
-12
0
12
-12
0
12
-12
0
12

Duration [min]
33
32
33
37
31
30
30
30
30
31
31
37
30
31
30
30
30
30
31
39
36
32
30
38
30
33
34
32
32
32
31
30
32
32
31
30

Table 4.4: Neutron source data runs used in the analysis.

58

Chapter 4. Testing with a neutron source

Here, we present the details of the reduction process, specifically for Run 85386 (z = 0). It
consists of 6 cuts that are supposed to reduce the following possibly contaminated events:(1) Non-AFT events
Events without AFT flag (i.e., TrgType 6= 2) are discarded, which removes about a halfof all SHE-triggered events. This cut ensures all selected events have a [-5, 535] µs range.
(2) Cosmic ray events
Cosmic rays can also trigger the SHE flag of the detector. Since cosmic rays are supposed
to hit OD PMTs as well as ID PMTs, a cut on the number of OD hits is necessary. The
distribution of NHITAC, the number of OD hits within the primary trigger time range,
for all SHE-triggered data events and the AmBe neutron MC (explained in Section 4.2)
are shown in Figure 4.6. Events with NHITAC > 10 are removed.
(3) Events possibly triggered by neutron capture gamma-rays
Multiple capture gamma-rays from a Gd(n, ) reaction may trigger the detector andsome of them can even cause scintillation. Since we are not supposed to find a neutron
signal from events triggered by capture gamma-rays, they must be removed. As can beobserved in Figure 4.4, most of the triggering scintillation starts to occur at 860 ns in thereconstructed capture time scale. Thus, events with TrgTagged — the number of tagged
candidates with reconstructed capture time < 900 ns — larger than 0 are regarded as
possibly triggered by neutron capture gamma-rays and thus removed. Unfortunately, hits
induced by a prompt gamma-ray or dim scintillation light are likely to be selected as atagged candidate, so this cut removes a huge portion of events with moderate QISMSK,
as illustrated in Figure 4.8.
(4) Events not likely to have been triggered by 4.4-MeV scintillation
Only the events in the higher QISMSK peak are selected. Events whose QISMSK is notin the range [850, 1200] for z = 0 m and [900, 1250] for z = ±12 m are removed. The
QISMSK distributions for three z positions (Run 85384, 85386, 85390) after applying theabove cuts are shown in Figure 4.9, along with the QISMSK cut range limits.

4.1. Neutron source data

59(5) Events with a large number of noise hits within the neutron search time range
An event variable MaxNonTrgN1300 was defined to represent the event’s maximum number of hits in 1.3 µs within the neutron search time range, which is reconstructed capture time 2 [4, 535] µs. Figure 4.10 shows the distributions of MaxNonTrgN1300 for
data and MC. Any event with MaxNonTrgN1300 > 160 was regarded as a scintillationcontaminated event, and thus discarded.
(6) Events that may have a neutron mixed in from the preceding trigger
If some two consecutive events are too close in time (i.e.,t . O(100) µs), a neutron fromthe earlier event may mix into the subsequent event. To prevent such a set of events,
an event variable TDi↵ [ms] was defined to represent the time di↵erence between two
consecutive SHE events. The TDi↵ distribution of data events are shown in Figure 4.11.
Any event with TDi↵ < 1 ms is removed.

Figure 4.6: The NHITAC distributions for data and MC before Cut 2 (left) and the same
distributions after Cut 2 (right). The red rectangle represents the selection range of Cut 2.

60

Chapter 4. Testing with a neutron source

Figure 4.7: The QISMSK distributions before and after Cut 2. Cut 2 apparently removes high
QISMSK events triggered by cosmic rays.

Figure 4.8: The QISMSK distributions of data events before applying Cut 3 (blue), after
applying Cut 3 (orange), and removed by Cut 3 (green).

Figure 4.9: The QISMSK distributions for three source z positions. Data for z = ±12 m show
slightly larger QISMSK peak values, as the source is closer to the tank wall. Cut 3 for z = ±12m leaves some of the events with moderate QISMSK unlike for z = 0 m.

4.1. Neutron source data

61

Figure 4.10: The MaxNonTrgN1300 distributions for data and MC (explained in Section
4.2 including source/sourceless dummy) (left) and the same distributions within the rangeof MaxNonTrgN1300 2 [0, 200] (right). The red rectangles represent the selection range of Cut
5. The long tail in the data distribution should be coming from the source itself, and the two
distinct scintillation peaks are observable.

Figure 4.11: The TDi↵ distributions for all SHE events (i.e., no cut applied) (blue) and for AFT
events only (i.e., Cut 1 applied) (orange). The exponential fit to the blue histogram gives the
decay time of ⇠12 ms. This manifests the Poissonian nature of the SHE triggers, and thus thereciprocal of the decay time directly relates to the trigger rate of ⇠83 Hz. For some unknown
reason, the AFT flag seems to be omitted more frequently for events close to the previous.

62

Chapter 4. Testing with a neutron source

Table 4.5 summarizes the applied event cuts. For Run 85386, the number of total SHE-triggered
events was 162,900, and the number of events passing Cuts 1-6 was 22,077 (13.55% pass ratio).
Neutron candidates were sought for within the reconstructed capture time (i.e., ReconCT [ns])range of [4, 535] µs, which corresponds to the data recorded ⇠3 µs after the presumed neutron
emission. The lower limit of 4 µs is to avoid the scintillation light, as shown in Figure 4.12.
Cut
1
2
3
4
5
6

Removed events
Non-AFT
Cosmic ray
Possibly triggered by neutronsNon-4.4-MeV scintillation
Noisy within neutron search time range
Too close to the previous event

Removed if
TrgType 6= 2
NHITAC > 10
TrgTagged > 0⇠(850 < QISMSK < 1200)
MaxNonTrgN1300 > 160
TDi↵ > 1

Removed ratio
0.54
0.08
0.34
0.51
< 0.01
0.01

Table 4.5: Summary of the event cuts applied to data events. The removed ratios (among theevents passing the previous cut) are obtained for Run 85386.
Cut
1

Removed candidates
Too close to scintillation

Removed if
ReconCT < 4000

Table 4.6: Summary of the candidate cut applied to data events.

Figure 4.12: ReconCT distributions for all candidates after event cuts (left) and after the final
candidate cut (right). The red rectangle in the left figure represents the selection range ofReconCT > 4,000 ns. The exponential fit to the right figure gives the time constant of 84 µs.
Event cuts 2, 5 and Candidate cut 1 were applied to MC dataset as well, as the noise hits for
MC are taken from data. The signal likelihood threshold for candidate-tagging was set to 0.7.

4.2. Neutron source simulation

4.2

63

Neutron source simulation

Since it was rather challenging to implement an inorganic scintillator with a complex geometry
in the Geant3-based simulator, the simulation of BGO and the trigger system was completely
neglected in this study. A single neutron with momentum range [60, 120] MeV/c (close to
kinetic energy range [2, 8] MeV) without any prompt gamma-ray or uncorrelated neutrons was
simulated with skdetsim-skgd. No other source elements were implemented in the simulator.
As will be explained later, the lack of BGO and trigger simulation should be the core reason
behind the discrepancy between data and MC.
The same physics models — G4NDL neutron cross section library for low energy neutrontransport and ANNRI-Gd model for Gd(n, ) reaction modeling — were used as in generationof MLP training dataset (refer to Section 3.1). The Gd concentration was set to ⇠0.0110 wt.%,
which is within 1% agreement with the current estimate. The neutron production vertex was
fixed at each calibration position, and no vertex smearing was applied. The neutron production
time and the trigger time was fixed at t =1,000 ns in the simulated detector’s global raw hit
time. Since the neutrons in data events are expected to be produced at t ⇠ 850 ns in raw PMT
hit time, searching for neutrons in the same time frame for both data and MC would causesome systematic error to the search result. However, this factor is completely neglected in thefollowing analyses, as the neutron’s capture time constant (⇠ O(100) µs) is much longer thanthe delay. Less than 0.2 % error is expected with the current Gd concentration.
The noise hits were extracted from the same dummy trigger data as in Section 3.1. One thing tonote is that the source itself also randomly generates unwanted noise hits, for example, uncorrelated neutrons and scintillation mixing into the neutron search time range. To compensate forsuch source-induced noise hits, another set of random trigger data was taken with the source at
z = 0 m for ⇠30 minutes (Run 85388). Let us call the dummy trigger data without the source“sourceless dummy”, and the ones with the source inside the tank “source dummy.” These twodi↵erent dummies serve as background measurement. 300,000 single neutron MC events weregenerated for each z position with sourceless dummy hits attached, and another 30,000 eventswere generated for z = 0 with source dummy hits attached.

64

4.3

Chapter 4. Testing with a neutron source

Data and MC comparison

Figures 4.13, 4.14, and 4.15 should give a general idea on the number of events and candidates
in data and MC. One thing to note is that there has been a reported increase of noise hits
throughout the entire run range used in the analysis. The reason for this is currently underinvestigation, while one of the suspected is fluorescent bacterial growth within the tank [48].
The current analysis should be a↵ected by this unresolved phenomenon, and it is discussed
qualitatively in Section 4.4. Further analysis should be done after the noise level is stabilized.
From this point and on, we assume that each data event is on an equal footing with a neutron MC event. The following subsections will compare distributions of the features, signal
likelihood, and reconstructed capture time.

4.3.1

Features

The feature distributions of all candidates from data and MC with sourceless dummy (z = 0
only, distributions for other z positions are shown in Appendix C) are shown in Figures 4.16.
Although there are more candidates per data events than per MC events as shown in Figure 4.14,
the normalized feature distributions are generally in good agreement, proving that modelingdata events with a single neutron MC plus noise hits from sourceless dummy data is an e↵ective
approach. Still, data events seem to have slightly more signal-like candidates than MC events,
as in Beta(l) distributions of Figure 4.16. Such tendency is washed out in Figure 4.18, wheredata events are compared with source dummy MC, suggesting that the afterglow of BGO shouldhave similar angular characteristics as signal hits.
The tagged-candidates distributions in Figures 4.17 and 4.19 also show reasonable agreement
between data and MC, except that data events show more p(n, )-like candidates and less
Gd(n, )-like candidates than MC events. Such tendency is observed in both sourceless andsource dummy MCs, and it should stem from the simulation settings such as the neutrontransport/capture models, cross sections, photon attenuation within Gd-doped water, quantum
efficiency of PMTs, etc. Those settings should be further tuned in the future.

4.3. Data and MC comparison

65

Figure 4.13: SHE-triggered event rates (top) and cut-passing event rates (bottom) by date.

Figure 4.14: Mean number of candidates per event, after all event and candidate cuts. The
reason why sourceless dummy MCs have less number of candidates is explained in Figure 4.15.

Figure 4.15: Mean number of candidates per event by date. The increase of the number bydate has been reported after Gd-loading, and the suspected is the possible bacterial growthwithin the tank. Asymmetry in z has also been reported for the noise hits [48], which may be
a reason why data with z = 12 have more candidates than the other two z positions.
The di↵erence in the number of candidates per event shown in Figure 4.14 should come fromthe di↵erence in the dates the dummy trigger data were taken. Sourceless dummy data weretaken from late August to early September, 2020, whereas the source dummy data was taken
on October 27th , 2020. Data taken on October 27th have slightly more candidates than thesource dummy, which should be due to the BGO’s phosphorescence.

66

Chapter 4. Testing with a neutron source

Figure 4.16: Area-normalized feature distributions of sourceless dummy MC and data taken atz = 0 (Run 85386), for all selected candidates. (blue: data, pink: MC (Gd), cyan: MC (H),
gray: MC (noise))4.3. Data and MC comparison

67

Figure 4.17: Area-normalized feature distributions of source dummy MC and data taken atz = 0 (Run 85386), for tagged candidates. (blue: data, pink: MC (Gd), cyan: MC (H), gray:
MC (noise))

68

Chapter 4. Testing with a neutron source

Figure 4.18: Area-normalized feature distributions of source dummy MC and data taken atz = 0 (Run 85386), for all selected candidates. (blue: data, pink: MC (Gd), cyan: MC (H),
gray: MC (noise))4.3. Data and MC comparison

69

Figure 4.19: Area-normalized feature distributions of source dummy MC and data taken atz = 0 (Run 85386), for tagged candidates. (blue: data, pink: MC (Gd), cyan: MC (H), gray:
MC (noise))

70

4.3.2

Chapter 4. Testing with a neutron source

Signal likelihood

Figure 4.21 shows the signal likelihood distributions of all candidates selected in data andsourceless dummy MC, and Figure 4.20 compares the same distributions of data and source
dummy MC, for three source z positions. Overall, the distributions for all three positions seemto have a similar shape, so both the trained MLP and the simplified source simulation withsourceless dummy prove to be e↵ective. One caveat to the source dummy is that it not only
includes phosphorescent noise hits, but also random neutron capture gamma-rays which shouldbe very signal-like but definitely should not be classified as signals in MC. A huge increase ofcandidates with signal likelihood close to 1 between the sourceless dummy MC (center, Figure4.21) and the source dummy MC (Figure 4.20) can be observed. Another thing to note is thatin Figure 4.21, the noise distribution for MC at z = 12 m, especially within the range of output
likelihood 2 [0.8, 1.0], is di↵erent from the ones of other two z positions. A reason suspected
is the reported increase of low energy noise hits centered around the top of the tank after the
Gd-loading was completed [48]. Further investigation is needed.
The neutron capture detection efficiencies ✏(n, ) for signal likelihood threshold of 0.7 are listed
in Table 4.7. The expected overall neutron capture detection efficiency with sourceless dummy
MC was ⇠60% for z = 12 m and ⇠56% for z = 0. The z-asymmetry in the noise ratio rNoise isnoticeable, which may be due to the reported increase of noise hits near the top of the tank.

Figure 4.20: Signal likelihood distribution of source dummy MC and data taken at z = 0 m(area-normalized). The red dotted line indicates the taggable likelihood threshold of 0.7.

4.3. Data and MC comparison

71

Figure 4.21: Signal likelihood distributions of sourceless dummy MC and data taken at z = 12m (top), z = 0 (center), z = 12 m (bottom) (area-normalized). The red dotted line indicates
the taggable likelihood threshold of 0.7.

z [m]
12
0
0
-12

Dummy
Sourceless
Sourceless
SourceSourceless

Efficiency (%)
Threshold GdH0.7
95.8 24.0
0.7
94.7 18.0
0.7
94.9 17.9
0.7
94.1 24.5

All
60.3
56.8
56.9
59.8

rNoise
3.30
1.68
11.07
1.54

Table 4.7: Summary of MC neutron capture detection efficiencies, for each source z position.
Most of the di↵erence in the noise rate rNoise between source and sourceless dummy MCs shouldbe random uncorrelated neutrons from the source.

72

4.3.3

Chapter 4. Testing with a neutron source

Reconstructed capture time

The reconstructed capture time (ReconCT) of a candidate is defined as the mean residual
hit time of the earliest and latest hits (in residual time) within a candidate. Considering the
maximum residual time width (14 ns) of a candidate and the time resolution of the residual
hit time illustrated in Figure 3.9, the reconstructed capture time should give a small enoughresolution to fit the exponential distribution of ⌧ ⇠ O(100) µs. The reconstructed capture time
distributions of data and MC are shown in Figure 4.22. The histograms were fitted with afunction of the form Ae t/⌧ + B within ReconCT 2 [15, 535] µs. The lower limit of 15 µs is
the expected time that neutrons require to thermalize. The mean of the fitted capture time
constant ⌧ was 116.7 ± 0.2 µs for data and 112.5 ± 0.4 µs for MC. The two fitted ⌧ s disagree,which implies that the neutron capture ratio on Gd is higher for MC than in data, as expectedby Figures 4.17 and 4.19. Since the error of the Gd concentration is ⇠1%, the Gd capture crosssections or the neutron thermalization process within Geant3 should be checked in the future.

Figure 4.22: Reconstructed capture time distributions of data (left) and MC (right). The MC
histogram for each z was scaled by (# of data events of the same z)/(# of MC events).
Figure 4.23 shows the time evolution of the fitted time constant ⌧ for each of the three source z
positions. The Gd-loaded water was injected from the bottom of the tank, and we can observethat the fitted ⌧ stabilizes for z = 12 m at around September 1st , 2020. The fitted ⌧ of data isslightly larger than that of MC for the same Gd concentration of 0.019 wt.%, as in Figure 4.22.
B, the number of fitted background candidates, is much larger for data events than MC events,
as data events include random uncorrelated neutrons from the source. For z = 0, the values offitted B are 142 ± 15, 21 ± 4, and 109 ± 26 tagged candidates per second for data, sourceless4.3. Data and MC comparison

73

dummy MC, and source dummy MC, respectively. Although the noise levels between the data
(Run 85386, Oct. 27, 2020) and the sourceless dummy MC (noise hits taken from late August to
early September, 2020) are clearly di↵erent as suggested in Figure 4.15, the number of tagged(i.e., signal-like) candidates should not be significantly a↵ected. Thus, the di↵erence between
the fitted B for sourceless and source dummy MCs should provide a rough yield of sourceinduced tagged background candidates. The estimated yield was 88 ± 26 tagged candidates
per second. Assuming roughly the neutron detection efficiency of 55 ± 10 % for z = 0 (basedon Table 4.7), we get the estimated source intensity of 160 ± 57 neutrons per second. This
number is in agreement with the source neutron intensity estimates 216 ± 48 (this study) and
163 ± 33 (Watanabe et al. [38]), although plagued by huge systematic uncertainties. Also, thenumber of leaking prompt gamma-rays, that misses the BGO but directly hits the detector to
be tagged as signal candidates, as well as the amount of tagged background candidates removed
by Event cut 5 (see Table 4.5) is neglected in this estimation. Further study of scintillation andlonger measurement of source dummy is needed to fully understand the connection between
the source particle intensity and the fitted background B.

Figure 4.24 shows the time evolution of the fitted background B for each of the three sourcez positions. The fitted Bs seem constant throughout the campaign, except for z = 12 m data
which show a slight increase of B over time. This may be due to the reported increase of noise
hits near the top of the tank [48].

Figure 4.23: Fitted capture time constant ⌧ by date.

74

Chapter 4. Testing with a neutron source

Figure 4.24: Fitted background B by date.

4.3.4

Estimated neutron detection efficiencies

The estimated neutron detection efficiency ✏est.
of both data and MC was defined as
n
✏est.
n =

(# of tagged candidates) (# of estimated background candidates)(# of cut passing events)

where the number of estimated background candidates is defined as B T /w = 53.1B.

(4.2)

T

is the neutron search time range of 531 µs and w is the bin width 10 µs of the reconstructedcapture time histogram. The definition of ✏est.
in Equation 4.2 assumes that each cut passing
nevent has a single emitted source neutron as in each MC event and that all time-correlated
neutrons are signals. The true neutron detection efficiency ✏n for MC was defined as

✏true
=
n

(# of tagged signal candidates)(# of cut passing events)(4.3)

For MC, the number of cut passing events should be the same as the number of emitted
neutrons. The estimated neutron detection efficiencies by date for both data and MC are as
illustrated in Figure 4.25. We can observe that the neutron detection efficiency for data looks
systematically lower than that of MC. The efficiencies of source and sourceless dummy MCs did
not show significant discrepancy. The detection efficiencies of data and MC are summarized in
est.
Figure 4.25 and Table 4.8. The ✏true
n ’s for MC agree with ✏n ’s within statistical uncertainties.

4.4. Discussion

75

by date.
Figure 4.25: Estimated neutron detection efficiency ✏est.
n
z [m]
12

Data est. (%)
49.96 ± 0.44

046.03 ± 0.33

-1250.79 ± 0.46

MC est. (%) Data est./MC est.
59.33 ± 0.33
0.84 ± 0.01
55.76 ± 0.24
0.83 ± 0.02
56.51 ± 1.45
0.81 ± 0.03
59.02 ± 0.38
0.86 ± 0.01

MC true (%)
59.40
55.91
56.07
58.91

MC dummy
sourcelesssourcelesssourcesourcelessTable 4.8: The estimated neutron detection efficiencies with statistical uncertainties.
(Data averaged from August 12th to December 1st , 2020.)

4.4

Discussion

The following subsections will discuss the possible explanations for the apparent data-MC
discrepancy shown in Figure 4.25 and Table 4.8, as well as the currently known sources ofsystematic bias in the neutron detection efficiency measurement.

4.4.1

Sources of systematic uncertainties

Apart from the data-MC discrepancy, there also exist several sources of systematic uncertainties
to the data ✏n measurement. The following are the possible sources and the rough estimates of
the related uncertainties:• Noise rate by date
Apparently the noise level had increased over time within the selected runs, and it cana↵ect the overall tagging efficiency. To verify the e↵ect of the noise level on the detection
efficiency ✏est.
n , another set of noise hits were extracted from the dummy trigger data taken

76

Chapter 4. Testing with a neutron source
on October 27th , 2020, and attached to the z = 0 MC. The detection efficiency ✏MCfornthe sourceless dummy MC was 53.9 ± 1.1 %. Comparing this value with the same ✏MCofnthe z = 0 sourceless dummy MC in Table 4.8, maximum ± 5 % uncertainty is assumed.
• Neutron-induced scintillation
Neutrons can hit the BGO crystals to make them scintillate, and a few of them can be as
bright as scintillation induced by the 4.4-MeV prompt gamma-rays. Moreover, they canbe captured close to BGO, or sometimes within BGO to cause scintillation.
To estimate the amount of such fake triggers within the selected data events, a Geant4based gamma-correlated neutron MC simulation was carried out along with the BGO
geometry and the default Geant4 scintillation model implemented in the simulator. The
amount of neutron-induced scintillation was estimated by tracking all interaction products within BGO. The simulated QISMSK distribution was linearly scaled to fit the two
QISMSK peaks in data events shown in Figure 4.9. The scaled QISMSK distributions foreach interaction channel of neutrons and prompt gamma-rays are given in Figure 4.26.
Among 56,748 cut passing MC events, 3,864 events had neutron-induced scintillation. Using the neutron intensity estimates shown in Table 4.1, the neutron-induced scintillation
rate within the cut-passing events was estimated to be 12 ± 3 % [49].
Since neutrons can both scintillate the crystal and leak into the fiducial volume to be
detected as signals, the actual percentage of “fake” events without any signal neutron is
expected to be lower than 12 %. We assume +10 % uncertainty.
• PMT gain
About ± 3% uncertainty due to the fluctuation of PMTs’ gain is reported by a past study[4], over ⇠8 years of time. We disregard the uncertainty in this analysis which was done
over a four-months range of runs.
• MLP classification reliability
About ± 5% uncertainty on the MLP classification was imposed by a past study [4]. We
impose the same for this study.

Table 4.9 summarizes the above listed sources of systematics. According to Tables 3.5 and 4.7,

4.4. Discussion

77

Figure 4.26: Scaled QISMSK distributions by interaction channels, for gamma-correlated neutrons (left) and for prompt gamma-rays (right). (area-normalized) The probability that neutrons will mimic the 4.4-MeV scintillation (i.e., QISMSK⇠1,000 p.e.) is O(10 2 ) lower thanthe prompt gamma-rays.
the true neutron capture detection efficiency ✏(n, ) of sourceless dummy MCs for z = ±12 m(60.3 % for z = 12 m, 59.8 % for z =12 m) are closer to that of random vertex MC (60.5 %).

We estimate the overall neutron detection efficiency of the SK detector by taking the mean of
the two z = ±12 m results, as 50 +86 (sys.) ±1 (stat.) %.
SourceNoise rate
Fake triggers
MLP classification
Total

Fractional uncertainty (%)( 5, +5)( 0, +10)( 5, +5)( 7, +12)

Table 4.9: Summary of the possible sources of systematic errors and their rough fractional
uncertainty estimates. All errors are assumed to be normal and independent.

4.4.2

The data-MC discrepancy

There are several factors within our MC, that are likely to cause the data-MC discrepancy:• Gd concentration in MC
Although rather unlikely, the input Gd concentration to the MC, which is based on the
total input Gd-sulfate mass to the detector system, might simply have been overestimatedcompared to the true concentration within the tank after purification, resulting in a higher
detection efficiency and a lower capture time constant (⌧MC ⇡ 113 µs, ⌧Data ⇡ 116 µs).

78

Chapter 4. Testing with a neutron source• Lack of BGO in MC
A few source neutrons can mimic the 4.4-MeV scintillation of prompt gamma-rays. Events
triggered by such neutron-induced scintillation are expected to lack detectable neutrons,
and thus become a source of systematic bias, lowering the estimated detection efficiency.
A simulation with the BGO geometry and a scintillation model tuned to the current
measurement is required to estimate the bias.
• Neutron physics (neutron transport/capture models and cross sections)
The fraction of neutron captures on Gd is apparently higher in MC than in data, according
to the feature distributions and the fitted capture time constants. Di↵erent neutron
physics models and cross sections should be tested to estimate the model-dependent bias.
• Untuned detector characteristics (photon attenuation in Gd-loaded water, PMT QE, etc.)
Almost two years have passed since the detector characteristics applied in the SK simulator were last tuned during the pure water phase (SK4). Those parameters can affect the horizontal placement of the NHits distributions of signal candidates, increasing/decreasing overall NHits of candidates.

Table 4.10 compares the measured neutron detection efficiencies plus their rough systematic
estimates with the same efficiencies of the MC event sets, including one with a lowered Gd
mass fraction (⇠0.0103 wt.%) to match the (n, ) capture time constant of ⇠116 µs of data.
While the tuned MC efficiencies for z = ±12 m are within the bounds of the systematics,
the two efficiencies still disagree for z = 0. It suggests that the Gd concentration in MC
alone cannot be a sole cause of the data-MC discrepancy. The possible relative significanceof the aforementioned factors are discussed in Appendix D, by fitting a parameterized NHits
distribution based on a hypothesis that the gamma-ray production model and the detector
characteristics are fully compatible with data.

4.4. Discussion

z [m]
12
0
-12

79

Data e↵. (%)
50
46
51

Sys. (%)
+6
4
+6
3
+6
4

Tuned MC e↵. (%)
54
55
57

Untuned MC e↵. (%)
59
56
59

Table 4.10: The measured neutron detection efficiencies and the rough systematic estimates.
The statistical uncertainties are neglected. The tuned MC efficiencies are taken from the MCwith its Gd cocentration lowered to ⇠0.0103 wt.% to match the fitted neutron capture time
constant to the one measured in data (⇠116 µs).

80

Chapter 4. Testing with a neutron source

4.4.3

Future work

Below is the list of future work:• BGO and trigger simulation
Geant4-based simulation is under development. The BGO-geometry is fully implemented,
and scintillation physics is yet to be tuned. If the BGO simulation combined with SHE
trigger simulation can reproduce the data QISMSK peak, we would be able to decompose
the peak and estimate the amount of fake triggers within the cut-passing events.
• Neutron physics models
Di↵erent models should be tested to estimate the model-dependent systematic.
• Tuning detector simulation parameters, stabilizing the noise level
The systematic coming from the untuned detector characteristics in simulation should beruled out to better analyze the significance of the fake triggers and the physics model.
• Longer source/dummy measurement with di↵erent positions

Chapter 5
Summary and prospects
To improve visibility of neutrino-induced neutrons, the SK detector has fully loaded Gd — themost powerful element known in capturing thermal neutrons — with 0.0110 wt.% concentration
since the end of July, 2020. Although some of the parameters still remain tuned to the purewater phase, a new simulation framework with the recent p(n, )/Gd(n, ) models and crosssections implemented has become ready for use since September, 2020.
To evaluate the detector’s achievable neutron detection efficiency, we have inherited the sameneutron detection algorithm [1, 3, 4] that has been used so far within SK for p(n, ) detection.
The algorithm consists of two stages; in Stage 1, we pick out any cluster of PMT hits thatare tightly localized in the ToF-corrected time frame, and in Stage 2, we categorize those
selected clusters into signal and noise according to the signal likelihood output of the pretrained MLP. The parameters that characterize Stage 1 were tuned in this study to maximize
neutron detection efficiency in the Gd-loaded (0.0110 wt.%) SK’s single neutron MC (detailedin Section 3.1) with the Geant4 neutron cross section library and the ANNRI-Gd model as
neutron capture physics models. A new MLP was designed and trained on 300,000 single
neutron (with random vertex) MC events generated by the same simulation framework. With
this detection algorithm and the realistic neutron search time range of [T0 +3, T0 +534] µs where
the neutron production time T0 = 1 µs, the neutron capture detection efficiency obtained ona test MC dataset was ⇠92 % for Gd(n, ) captures, ⇠26 % for p(n, ) captures, and ⇠61 %81

82

Chapter 5. Summary and prospects

overall. The parameter optimization and the training of MLP to Gd-loaded detector MC havenot a↵ected the p(n, ) neutron detection efficiency in pure water ( only) which was reported
to be ⇠25 % by the past studies [1, 3, 4].
The same algorithm was applied to both AmBe neutron source data and single neutron MC.
The signal likelihood and feature distributions showed reasonable agreement between data andMC, proving the e↵ectiveness of our primitive approach that neglects complex source geometry
and simplifies the detector trigger system. However, the MC showed higher neutron capture
ratio on Gd than in data, which should be one of the causes for the discrepancy in the estimatedneutron detection efficiencies between data (⇠50 %) and MC (⇠60 %). Another source of theapparent discrepancy should be fake neutron event triggers due to scintillation induced bysource neutrons, which are not at all accounted for in MC. The neutron physics models, crosssections, untuned detector parameters, and fluctuation in the detector dark noise level shouldcause bias as well. Roughly estimating the systematics, we set 50+64 %, 46+63 %, 51+64 % foreach z being 12, 0, -12 m, as the neutron detection efficiencies of the algorithm and the current
SK detector of 0.0110 wt.% Gd.
We hope to lower the uncertainties in the neutron detection efficiency measurement, whichis the major systematic in the neutron multiplicity of neutrino-nucleus interaction. Enhanced
neutron detection efficiency with smaller uncertainty would definitely accelerate mass hierarchy
and CP violation search with atmospheric/beam neutrinos and supernova relic neutrino search.

Appendix A
Neutrino oscillation
The three neutrino weak flavor eigenstates that have been discovered so far form a complete,orthonormal basis for the three lepton generations of the Standard Model. If the neutrinos
were massive, we may define another complete, orthonormal basis of three mass eigenstates,
and this basis should not necessarily coincide with the flavor basis. Let the flavor eigenstate bea coherent superposition of the mass eigenstates, related by a 3 ⇥ 3 unitary matrix U :|⌫l i =

Xk

Ulk⇤ |⌫k i(A.1)

where |⌫l i is the flavor eigenstate of flavor l 2 {e, µ, ⌧ }, and |⌫k i is the k-th mass eigenstate wherek 2 {1, 2, 3}. If the neutrino mass eigenvalues are non-degenerate, an interesting phenomenon
can happen. Let |⌫i denote a neutrino state that has been produced via weak interaction as|⌫l i at time x0 = 0. Then, each mass eigenstate should gain a distinct phase as it propagatesin (vacuum) spacetime, so that|⌫(x)i =
=

X
k

X
l0 ,k

where

Ulk⇤ |⌫k i e i k (x)(A.2)
Ulk⇤ Ul0 k |⌫l0 i e i k (x)k (x) := pk · x is the phase gain of |⌫k i after propagation x.

83

84

Chapter A. Neutrino oscillation

Naturally, the probability that the state |⌫i that has started out with flavor l will be detected
to have another flavor l0 is

P (⌫l ! ⌫l0 )(x) ==X

2Ulk⇤ Ul0 k e i k (x)kll0 + 2ReX⇥

Ulk⇤ Ul0 k Ulk0 Ul⇤0 k0 e i

kk0 (x)k0 6=kwhere

kk0 (x) ⌘

k (x)

(A.3)

⇤

k0 (x). If all mass eigenstates propagate in-phase, i.e., the neutrino

masses are degenerate so thatk (x) is the same for all k, P (⌫l ! ⌫l0 ) =ll0 by the unitarityof U , and thus there is zero chance that a certain flavor eigenstate will change its flavor as
it propagates. But if the neutrino masses are non-degenerate so that the mass eigenstates
become out-of-phase, the probability in Equation A.3 should be non-zero for l 6= l0 . Such flavor
transition allowed by the non-degeneracy of neutrino masses is called “oscillation,” due to thesinusoidal form of the transition probability in Equation A.3. The oscillation of neutrinos havebeen confirmed by experiments for various sources: solar, atmospheric, and artificial.

In the ultrarelativistic limit,

k (x) reduces to2
k (x) ⇡ mk L/2E, where mk is the non-degeneratemass of the k-th mass eigenstate |⌫k i, L = |~x|2 is the distance that the neutrino has traveled,and E = p0k = E⌫ is assumed to be the energy of the neutrino, same for all k. Hence, theoscillation probability P (⌫l ! ⌫l0 ) depends on the components of the so-called neutrino mixing(or the PMNS) matrix U , the squared di↵erences

m2kk0 ⌘ m2k

m2k0 of two distinct neutrino

masses, and the neutrino’s kinematics: the traveled distance L and the neutrino energy E.

In a scenario where only the Dirac masses are concerned, the unitary 3 ⇥ 3 neutrino mixing
matrix U is parameterizable with four real free parameters, considering the neutrino and lepton
spinors of three generations absorbing the relative phases and removing

3
2

1 degrees of

freedom from U , so that only the remaining four degrees of freedom can be actually physical.
The standard way of parameterizing the mixing matrix U is:

85

2 3 2
⌫e
1
6 7 6
4⌫µ 5 = 40
⌫⌧
0
2

3
0
0
7
cos ✓23 sin ✓23 5
sin ✓23 cos ✓23

3
cos ✓130 sin ✓13 e i CP
6
7
0
1
0
4
5
i CP
sin ✓13 e0cos ✓132
32 3
cos ✓12 sin ✓12 0
⌫1
6
76 7
4 sin ✓12 cos ✓12 05 4⌫2 5
0
0
1
⌫3

(A.4)

where ✓ij for i, j 2 {1, 2, 3} are the mixing angles that represent the internal mixing of mass
eigenstates, andCPrepresents the overall complex phase of U . Note that the complex phase

is not allowed for the eigenspace dimension smaller than 3. The complex phase of the mixing
matrix U allows for the CP violation (CPV), i.e., a non-zero di↵erence between the oscillation
probabilities of neutrinos and antineutrinos, according to Equation A.3. The three mixing
angles ✓ij , the (Dirac) CP phase CP , and the two independent neutrino mass-squared di↵erences
m221 andm232 are the fundamental parameters that fully characterize the oscillation with agiven neutrino energy and propagation. These are called the neutrino oscillation parameters.
Most of the oscillation parameters except

CP are relatively precisely measured as in Table A.1.

CP remains largely unknown, while CP conservation ( CP = 0 or ⇡) is disfavored according to

the recent data analysis [18]. A non-zero and non-⇡

CP , if precisely measured, may accountfor the matter-dominant universe, as the same possible CP violation currently observed in thequark sector is known to be either non-existent or extremely small. Another pending problem
is the sign ofm232 . The sign is not known yet, as the leading order (vacuum) oscillation

probability is not sensitive to its sign. This problem is called the neutrino mass hierarchy(MH) problem. If(NH), and if

m232 > 0, i.e., m3 > m2 > m1 , we call this scenario the normal hierarchym232 < 0, i.e., m2 > m1 > m3 , we call this scenario the inverted hierarchy (IH).

The mass hierarchy problem is correlated with the CPV search to some extent, as in matter
unlike in vacuum, the sign ofantineutrinos.

m232 can fake CPV to a di↵erent extent between neutrinos and86

Chapter A. Neutrino oscillation

Parameter

Best-fit

3 range

3 error| m232 |2.525

2.411 – 2.646

9.31%

m221

7.37

6.93 – 7.96

13.8%

sin2 ✓13

0.02150.0190 – 0.0240

23.3%

sin2 ✓120.297

0.250 – 0.354

35.0%

sin2 ✓23

0.425

0.381 – 0.615

55.1%

CP /⇡

1.38

0 – 0.17

0.76 – 2

102%

Table A.1: Recent (circa 2017) best-fit values, 3 ranges, and errors for the oscillation parameters, assuming normal mass hierarchy (m1 < m2 < m3 ). The unit of m232 is 10 3 eV2 , andthe unit of m221 is 10 5 eV2 [50]
⌫l

l

⌫¯l

W+
n

l+

W
p

p

n

Figure A.1: Feynman diagrams for ⌫l n (left) and ⌫¯l p (right) CCQE interactions.

In practice, we indirectly assume the presence of neutrinos by detecting the products of neutrinos’ weak interactions with matter, i.e., matter nuclei and bound electrons. The weak interaction between neutrinos and the other fermions can occur either through charged current(CC) exchange mediated by one of the two W bosons (W + , W ), or through neutral current(NC) exchange mediated by the neutral Z boson. All neutrinos can elastically scatter o↵ abound electron via NC exchange, and ⌫e can do the same via CC exchange as well. Although
the neutrino-electron scattering cross section is theoretically well-estimated, the interaction is
far rarer than the neutrino-nucleon scattering by O(10 3 ) for moderate neutrino energies (⇠1
GeV). While neutrinos can go through NC exchange to interact with and scatter nucleons, amore dominant and visible interaction channel with nucleons of sub-GeV neutrinos is the CC
quasi-elastic (CCQE) interaction, where neutrinos exchange a W boson with a nucleon to produce a clearly detectable charged lepton (whose nonzero mass explains the name quasi -elastic)
and a recoil nucleon which is invisible more than often in experiments. The simplified feynman
diagrams of ⌫n and ⌫¯p CCQE interactions are shown in Figure A.1.

87

Figure A.2: The theoretical and experimental values of neutrino-nucleon (left) and antineutrinonucleon CC interaction cross sections (right). Reprinted from Ref. [51].
Due to the simplicity of the interaction, relatively high cross section, and the detectability of arecoil charged lepton, most of the detected neutrino event statistics that are used to constrain
the oscillation parameters are dominated by CCQE. In reality, the neutrino targets are not freenucleons, and the neutrinos technically interact with each nucleon bound to a nucleus. Although
the possible intranuclear interactions of the neutrino-nucleon final states — also known as the
final state interactions (FSI) — complicates the event topology and makes it challenging to⇤Zmodel the dynamics, in the ideal 2-body CCQE scenarios ⌫l + A
X ! l + Z A1 X + p and⇤

Z 1Z⌫¯l + A
X ! l+ + A
1 Y + n with one outgoing lepton l and one nucleon n/p each, the energy of

the incoming neutrino E⌫ can be kinematically deterimined as

E⌫ =

2
(min
n/p )

e↵ 2
(mout,)n/pe↵2(mout,n/pe↵
m2l + 2(mout,)Eln/p

El + pl cos ✓l )out,e↵
where min⇡ moutn/p is the mass of the incoming nucleon, mn/p
n/p

(A.5)

Eb is the e↵ective mass ofthe outgoing nucleon considering the nuclear binding energy Eb , and ml , El , pl , and ✓l are eachthe mass, energy, momentum, and scattering angle of the outgoing lepton.

Appendix B
Neutron production mechanism of
AmBe neutron source
The 241 Am is an alpha-decaying isotope with a half-life of ⇠432 years:
241Am

432 years 237!

Np + 4 ↵(5.5 MeV) + (59.5409 keV)

(B.1)

The long half-life of 241 Am ensures stable intensity throughout a long period of time, unlike
spontaneous fission neutron source like 252 Cf. Since the second most frequent decay channel
has a branching ratio of O(10 12 , we can regard all 241 Am isotopes undergo the alpha decay
described in Equation B.1.
The AmBe source is a mixture of 241 Am, the alpha emitter, and 9 Be, the alpha absorber. The
source itself is designed to be a stainless steel casing encapsulating the Americium-oxide andBe metallic powder. When bombarded with alpha particles, 9 Be absorbs one of those to form⇤⇤

a compound nucleus of 13 C , which instantaneously decomposes into 12 C and a neutron. The
Q-value is 5.7 MeV and thus exothermic, and kinematically alpha with zero kinetic energy⇤

can make this reaction happen. The excited 12 C can be in one of the few possible energy
levels including the ground state, and if it is in the first excited state (4.44 MeV, spin-parity
2+ ), it de-excites with emitting the same energy gamma-ray. In the second excited state (7.66
88

89

Figure B.1: The neutron yield per alpha from measurement. 10% errors assumed for all data
points. [52] Exponential ae bx + c fit gives a good approximation.
MeV) that happens with alpha threshold energy of ⇠2 MeV, it de-excites with emitting analpha particle and rarely a gamma-ray. For the third and higher excited states, it undergoes
a breakup reaction of emitting 3 alpha particles. This should be negligible for AmBe. The
neutron emission process correlated with 4.4 MeV gamma-ray is described in the equations
below.

9⇤⇤

Be + ↵ ! 13 C ! 12 C + n + Q(5.704 MeV)
12C

⇤ ⇠60% 12!

C+

1 (4.438 MeV)

(B.2)
(B.3)

Considering the measured cross sections of each channel of 9 Be(↵, n)12 C reaction with the
kinetic energy of ↵ being ⇠ 5.5 MeV, about 78% of the fast neutrons emitted via Equation B.2
should be correlated with 4.4 MeV gamma-rays. However in reality, the alpha particles may
slow down before they are absorbed. Z. Liu et al [37] recommends the R-value, or the 4.438
MeV gamma-ray to the total neutron ratio, of 0.575 ± 4.8% which should be dependent on the
AmBe capsule, oxygen component, etc.

Appendix C
Feature distributions of neutron source
data and MC for nonzero heights
The following are the feature distributions of neutron source data and MC for z = ±12 m:

90

91

Figure C.1: Feature distributions of source dummy MC and data taken at z = 12 m (Run85390), for all candidates. (blue: data, pink: MC (Gd), cyan: MC (H), gray: MC (noise))

92

Chapter C. Feature distributions of neutron source data and MC for nonzero heights

Figure C.2: Feature distributions of source dummy MC and data taken at z = 12 m (Run85390), for tagged candidates. (blue: data, pink: MC (Gd), cyan: MC (H), gray: MC (noise))

93

Figure C.3: Feature distributions of source dummy MC and data taken at z = 12 m (Run85384), for all candidates. (blue: data, pink: MC (Gd), cyan: MC (H), gray: MC (noise))

94

Chapter C. Feature distributions of neutron source data and MC for nonzero heights

Figure C.4: Feature distributions of source dummy MC and data taken at z = 12 m (Run85384), for tagged candidates. (blue: data, pink: MC (Gd), cyan: MC (H), gray: MC (noise))Appendix D
Hypothesis test on NHits
It is difficult to evaluate the factors mentioned in Subsection 4.4.1 at the current stage, however,
we can test a hypothesis to check which factor is the most significant. The hypothesis is thatthe detector characteristics of the simulator and the neutron capture models (i.e., the gammaray production models) are fully compatible with data. The cross sections are allowed to di↵er
between data and MC, so that the capture ratio on Gd can be a free parameter. Under thisassumption, we can model the tagged data distribution of feature x with individual disintegrated
histograms we obtain from either MC or data.
First, let h1n
MC (x) denote the bin value of the x histogram for single neutron source dummy MC.
Then, h1n
MC (x) can be expanded with a set of individual histograms as in Equation D.1:

h1n
MC (x) ' cMCnrMC hGd
MC (x) + (1rMC )hHMC (x)

o

B,all
+ hB
dark (x) + hsource (x)(D.1)

Here, hiMC (i 2 {Gd, H}) represents the tagged x histogram of a single neutron MC where thesimulated neutron is completely captured by the ith nuclide. hiMC should be a convolution of thedetector’s response (e.g., simulation physics models, water parameters, PMT gain, etc.) andthe algorithm’s detection efficiency (e.g., neural network, likelihood threshold, cut conditions,etc.). hB
dark denotes the tagged x histogram of completely dark detector (i.e., sourceless dummy
configuration), whereas hB,all
source denotes the same histogram with the source dummy configuration95

96

Chapter D. Hypothesis test on NHits

subtracted by hB
dark . rMC is the MC’s captured neutron ratio to all generated neutrons, andrMC is the MC’s neutron capture ratio on Gd. The values of cMC and rMC are 0.984 and 0.507
for ⇠0.0110 wt.% Gd concentration.
In a similar way, we can also model the x distribution hcorr
data of a gamma-correlated data event
where we expect a single neutron emitted from the source to be detectable, as in Equation
D.2. Values and distributions that might di↵er from MC are subscripted “data.” hB
dark shouldbe the same between data and MC, since we take those from dummy trigger data rather fromsimulation. Another thing to note is that the gamma-correlated data events should have somepart of its radiation as noise, while the rest (i.e., gamma-correlated neutrons) should be included
B,allwithin hidata (i 2 {Gd, H}). Hence, we denote the former as hB,corrsource := bhsource , where the factor

b should be close to the AmBe R-value (ratio of correlated neutrons to all neutrons) of ⇠0.58
or higher in case a fraction of prompt gamma-rays are tagged by the algorithm.

n
Gd
hcorr
'
c
data rdata hdata (x) + (1dataoB,corrrdata )hH
(x)
+ hB
data
dark (x) + hsource (x)(D.2)

Now, let us introduce the parameter s that denotes the ratio of correlated neutron events to all
event-cut passing events. Then, we can model the x distribution of all cut passing data eventswith Equation D.3.

B,allhdata (x) ' shcorrs) hBdata + (1dark (x) + hsource (x)noHB,corrB= s cdata rdata hGd
(x)
+
c
(1
r
)h
(x)
+
h
(x)
+ (1 s)hB,alldata
data data
data
sourcesource (x) + hdark (x)
n
o
H
B,all
B
= s cdata rdata hGd
(x)
+
c
(1
r
)h
(x)
+
bh
(x)
+ (1 s)hB,alldata
data data
data
sourcesource (x) + hdark (x)= scdata rdata hGddata (x) + scdata (1rdata )hH
data (x) + {1(1Bb)s}hB,all
source (x) + hdark (x)(D.3)

Assuming hiMC ' hidata for i 2 {Gd, H}, i.e., assuming that the simulated detector response andthe algorithm efficiency are the same in both data and MC, we can model hdata with MC signal
distributions hiMC , the dark distribution hdark from sourceless dummy, and the source-induced

97
noise distribution hB,all
source from source dummy, along with four fittable parameters s, c(= cdata ),r(= rdata ), and b, each within the range [0,1]. This is expressed in Equation D.4.

Gdhfit
MC (x) := scrhMC (x) + sc(1r)hH
MC (x) + {1(1Bb)s}hB,all
source (x) + hdark (x)(D.4)

Below is the summary of the four free parameters:• s: fraction of events with one signal neutron (sMC = 1)
• c: fraction of captured neutrons to all emitted neutrons (cMC = 0.984)• r: fraction of neutron captures on Gd (rMC = 0.507)• b: scaling parameter between gamma-uncorrelated source-induced background and all
source-induced background

With the above hypothesis assuming Gaussian, we have estimated the best-fit parameter valuesof s, c, r, and b using a simple

2-test on the NHits distribution for Run 85386 (z = 0). hB
dark

was taken from Run 85392, the closest sourceless dummy run to Run 85386. hB,all
source was taken
as the NHits distribution the only source dummy run available (Run 85388, z = 0) subtractediby hB
dark . hMC are taken from the source dummy MC used in Section 4.3. For all datasets, thesignal likelihood threshold was set to 0.546, which gives 56.61 % capture detection efficiency
with 1.99 % noise ratio for within sourceless dummy MC whose noise hits are taken from Run
85392. The corresponding capture detection efficiency is 94.62 % for Gd(n, ) and 17.52 % forp(n, ) reaction. Each h(NHits) is shown in the left figure of Figure D.1.
2=

P

x (hdata (x)

2
fit
hfit
MC (x)) /hMC (x) was calculated for bins with more than 5 candidates inthe data NHits distribution. The

2-minimizing values of parameters s, c, r, and b were 1.000,

1.000, 0.423, and 0.585, respectively, and the fitted hfit
MC along with hdata is shown in the rightof Figure D.1. Since the assumption hiMC ' hidata for i 2 {Gd, H} is only approximate as thewater parameters have not been tuned for Gd-loaded water, there is no reason to expect a
perfect fit. Nonetheless, the fitted MC histogram shows reasonable agreement with data. The

98

Chapter D. Hypothesis test on NHits

Figure D.1: The disintegrated histograms used for parameter estimation (left) and the fitted
NHits distributions with the base histograms (right).
rMC value of 0.507 and the fitted parameter values of s = 1.00 and r(= rdata ) = 0.43 indicatethat neutronless triggers may be a negligible source of uncertainty, compared to the neutron
capture ratio on Gd. The fitted b is 0.585, which is close to the AmBe R-value of ⇠0.575.
With the Gaussian assumption, the parameter likelihood is assumed to be proportional toe

2 /2, and the normalized single parameter likelihood distributions are as illustrated in Figure

D.2. The likelihood peaks of s and b are very localized so that their peak widths are within 2 %of the peak values, backing up that the neutron capture ratio on Gd should be a more serious
source of uncertainty than neutronless triggers under the current fit model.

Figure D.2: Parameter fit.
fit
The ratio rMC /rdatais ⇠0.84. The Gd-neutron capture ratio r itself should ber⇡

Z EmaxEmin

RE
nGd 0 n Gd (E)/(dE/dx) dE⇢(En ) P
dEn
R En
X (E)/(dE/dx) dE
X nX 0(D.5)

99
where Emax and Emin are the maximum and minimum kinetic energy of AmBe neutrons, ⇢(En )is the fraction of AmBe neutrons with kinetic energy En , nX is the number density of element
X,

X (E) is the capture cross section of element X for neutron kinetic energy E, and dE/dx

is the stopping power of neutrons in Gd-loaded water. For simplicity in MC, X is either Gd
or H. Since the number densities are estimated within 1 % precision, ⇢(E),

X (E), and dE/dx

might be causing errors. To verify the source of uncertainty, a cross-check simulation with a
precise neutron energy spectrum and di↵erent physics models should be carried out.
For reference, tuning the Gd concentration in MC to fit the capture time constant ⌧ to ⌧data ⇡116 µs has yielded ✏est.
⇠55 % for sourceless dummy MC and ✏est.
⇠52 % for source dummynnMC. The Gd concentration was set to ⇠0.0103 wt.%, and the corresponding rMC was 0.4953.
The discrepancy remains to persist even with lowered Gd concentration.

Bibliography[1] T. J. Irvine. “Development of Neutron-Tagging Techniques and Application to Atmospheric Neutrino Oscillation Analysis in Super-Kamiokande”. PhD thesis. The Universityof Tokyo, 2014.
[2] K. M. Tsui. Super-K Gd project: Neutron-tagging algorithm. 2017 Autumn Meeting, The
Physical Society of Japan. Sept. 12, 2017. url: http://www- sk.icrr.u- tokyo.ac.
jp / ~nakahata _ s / tokusui / talk / ntag _ jps _ autumn _ 2017 _ kaming . pdf (visited on
12/14/2020).
[3] T. Mochizuki. “Development of 50-cm Diameter Photomultiplier Tubes and a Neutron
Tagging Algorithm for Hyper-Kamiokande”. MA thesis. The University of Tokyo, 2019.
(Japanese).
[4] R. Akutsu. “A Study of Neutrons Associated with Neutrino and Antineutrino Interactions
on the Water Target at the T2K Far Detector”. PhD thesis. The University of Tokyo,
2019.
[5] K. Hirata et al. “Observation of a neutrino burst from the supernova SN1987A”. In: Phys.
Rev. Lett. 58 (14 Apr. 1987), pp. 1490–1493. doi: 10.1103/PhysRevLett.58.1490. url:
https://link.aps.org/doi/10.1103/PhysRevLett.58.1490.
[6] Y. Fukuda et al. “Evidence for Oscillation of Atmospheric Neutrinos”. In: Phys. Rev.
Lett. 81 (8 Aug. 1998), pp. 1562–1567. doi: 10 . 1103 / PhysRevLett . 81 . 1562. url:
https://link.aps.org/doi/10.1103/PhysRevLett.81.1562.

100

BIBLIOGRAPHY

101[7] A. D. Sakharov. “SPECIAL ISSUE: Violation of CP in variance, C asymmetry, andbaryon asymmetry of the universe”. In: Soviet Physics Uspekhi 34.5 (May 1991), pp. 392–
393. doi: 10.1070/PU1991v034n05ABEH002497.
[8] C. L. Cowan et al. “Detection of the Free Neutrino: a Confirmation”. In: Science 124.3212(1956), pp. 103–104. issn: 0036-8075. doi: 10 . 1126 / science . 124 . 3212 . 103. url:
https://science.sciencemag.org/content/124/3212/103.
[9] S. Ando and K. Sato. “Relic neutrino background from cosmological supernovae”. In: New
J. Phys. 6 (2004), p. 170. doi: 10.1088/1367-2630/6/1/170. arXiv: astro-ph/0410061.
[10] K. S. Krane. Introductory nuclear physics. New York, NY: Wiley, 1988.
[11] “Neutron Nuclear Reactions”. In: Nuclear Reactor Physics. John Wiley & Sons, Ltd,
2007. Chap. 1, pp. 1–32. isbn: 9783527611041. doi: https : / / doi . org / 10 . 1002 /9783527611041.ch1. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/
9783527611041.ch1.
[12] K. Shibata et al. “JENDL-4.0: A New Library for Nuclear Science and Engineering”.
In: Journal of Nuclear Science and Technology 48.1 (2011), pp. 1–30. doi: 10 . 1080 /
18811248.2011.9711675. url: https://doi.org/10.1080/18811248.2011.9711675.
[13] L. Marti Magro et al. “Evaluation of gadolinium’s action on water Cherenkov detector
systems with EGADS”. In: Nuclear Instruments and Methods in Physics Research Section
A: Accelerators, Spectrometers, Detectors and Associated Equipment 959 (Feb. 2020),
p. 163549. doi: 10.1016/j.nima.2020.163549.
[14] K. Hagiwara et al. “Gamma-ray spectrum from thermal neutron capture on gadolinium157”. In: Progress of Theoretical and Experimental Physics 2019.2 (2019). issn: 2050-3911.
doi: 10.1093/ptep/ptz002. url: https://doi.org/10.1093/ptep/ptz002.
[15] P. A. Cherenkov. “Visible emission of clean liquids by action ofradiation”. In: Dokl.

Akad. Nauk. SSSR 2, 451 (1934).
[16] T. Lou. “Development of 50 cm Diameter Hybrid Photo-Detector for Hyper-Kamiokande”.
MA thesis. The University of Tokyo, 2017.

102

BIBLIOGRAPHY[17] J. D. Jackson. Classical electrodynamics. Third edition. New York : Wiley, [1999] ©1999.
url: https://search.library.wisc.edu/catalog/999849741702121.
[18] K. Abe et al. “Constraint on the matter–antimatter symmetry-violating phase in neutrino
oscillations”. In: Nature 580.7803 (2020), pp. 339–344. issn: 1476-4687. doi: 10.1038/
s41586-020-2177-0. url: https://doi.org/10.1038/s41586-020-2177-0.
[19] Super-Kamiokande Collaboration. Introduction of Gadolinium into Super-Kamiokande
and the Start of New Observations. Aug. 21, 2020. url: http : / / www - sk . icrr . u tokyo.ac.jp/sk/news/2020/08/sk-gd-detail-e.html.
[20] M. B. Chadwick et al. “ENDF/B-VII.1 Nuclear Data for Science and Technology: Cross
Sections, Covariances, Fission Product Yields and Decay Data”. In: Nuclear Data Sheets
112.12 (2011), pp. 2887–2996. issn: 0090-3752. doi: https://doi.org/10.1016/j.
nds.2011.11.002. url: http://www.sciencedirect.com/science/article/pii/
S009037521100113X.
[21] C. Simpson et al. “Sensitivity of Super-Kamiokande with Gadolinium to Low Energy
Antineutrinos from Pre-supernova Emission”. In: The Astrophysical Journal 885.2 (Nov.
2019), p. 133. doi: 10.3847/1538-4357/ab4883. url: https://doi.org/10.3847%7B%
5C%%7D2F1538-4357%7B%5C%%7D2Fab4883.
[22] K. M. Tsui et al. skdetsim-skgd. https://github.com/kmtsui/skdetsim-skgd. 2020.
commit a8caebc.
[23] C. Zeitnitz and T. A. Gabriel. “The GEANT-CALOR interface and benchmark calculations of ZEUS test calorimeters”. In: Nuclear Instruments and Methods in Physics
Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment
349.1 (1994), pp. 106–111. issn: 0168-9002. doi: https://doi.org/10.1016/01689002(94)90613- 0. url: http://www.sciencedirect.com/science/article/pii/
0168900294906130.
[24] M. Sakuda et al. Website of the ANNRI-Gd model. http://www.physics.okayamau.ac.jp/~sakuda/ANNRI-Gd_ver1.html. (Visited on 12/14/2020).

BIBLIOGRAPHY

103[25] T. Tanaka et al. “Gamma-ray spectra from thermal neutron capture on gadolinium-155
and natural gadolinium”. In: Progress of Theoretical and Experimental Physics 2020.4(Nov. 2020). issn: 2050-3911. doi: 10.1093/ptep/ptaa015. url: https://doi.org/10.
1093/ptep/ptaa015.
[26] E. Mendoza and D. Cano-Ott. “Update of the Evaluated Neutron Cross Section Libraries
for the Geant4 Code”. In: ().
[27] G. Horton-Smith. GLG4sim, Generic liquid-scintillator anti-neutrino detector Geant4
simulation. url: http://neutrino.phys.ksu.edu/~GLG4sim/ (visited on 01/06/2021).
[28] S. Tobayama. “An analysis of the oscillation of atmospheric neutrinos”. PhD thesis. University of British Columbia, 2016.
[29] F. Chollet et al. Keras. https://keras.io. 2015.
[30] I. Goodfellow, Y. Bengio, and A. Courville. Deep Learning. MIT Press, 2016.
[31] M. Smy. “Low Energy Event Reconstruction and Selection in Super-Kamiokande-III”.
In: 30th International Cosmic Ray Conference. Vol. 5. July 2007, pp. 1279–1282.
[32] A. Bellerive et al. “The Sudbury Neutrino Observatory”. In: Nuclear Physics B 908 (July
2016), pp. 30–51. issn: 0550-3213. doi: 10 . 1016 / j . nuclphysb . 2016 . 04 . 035. url:
http://dx.doi.org/10.1016/j.nuclphysb.2016.04.035.
[33] A. Hoecker et al. “TMVA: Toolkit for Multivariate Data Analysis”. In: PoS ACAT (2007),
p. 040. arXiv: physics/0703039.
[34] A. Hoecker et al. User’s Guide. Oct. 4, 2013. url: http://tmva.sourceforge.net/
old_site/docu/TMVAUsersGuide.pdf (visited on 12/22/2020).
[35] D. P. Kingma and J. Ba. Adam: A Method for Stochastic Optimization. 2017. arXiv:
1412.6980 [cs.LG].
[36] T. Yano, Y. Isobe, and A. Suzuki. Installation of Automated Source Deployment System.
Nov. 8, 2018.

104

BIBLIOGRAPHY[37] Z. Liu et al. “The 4.438MeV gamma to neutron ratio for the Am–Be neutron source”. In:
Applied Radiation and Isotopes 65.12 (2007), pp. 1318–1321. issn: 0969-8043. doi: https://doi.org/10.1016/j.apradiso.2007.04.007. url: http://www.sciencedirect.
com/science/article/pii/S0969804307001200.
[38] H. Watanabe et al. “First study of neutron tagging with a water Cherenkov detector”.
In: Astroparticle Physics 31.4 (2009), pp. 320–328. issn: 0927-6505. doi: https://doi.
org/10.1016/j.astropartphys.2009.03.002. url: http://www.sciencedirect.
com/science/article/pii/S0927650509000401.
[39] Eckert & Ziegler Nuclitec GmbH. Sealed Radiation Sources Product Information. url:
http://www.hanshin-eng.co.kr/index/img/5Nuclitec.pdf (visited on 12/26/2020).
[40] J. Scherzinger et al. “A comparison of untagged gamma-ray and tagged-neutron yields
from 241AmBe and 238PuBe sources”. In: Applied Radiation and Isotopes 127 (2017),
pp. 98–102. issn: 0969-8043. doi: https://doi.org/10.1016/j.apradiso.2017.05.
014. url: http://www.sciencedirect.com/science/article/pii/S0969804316309861.
[41] Neutron sources:Americium-241/Beryllium and Californium-252. 1976. url: https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-UR-11-01808.
The source is re-referenced from the URL.
[42] K. W. Geiger and L. Van Der Zwan. “Radioactive neutron source spectra from 9Be(↵,n) cross section data”. In: Nuclear Instruments and Methods 131.2 (1975), pp. 315–321.
issn: 0029-554X. doi: https : / / doi . org / 10 . 1016 / 0029 - 554X(75 ) 90336 - 5. url:
http://www.sciencedirect.com/science/article/pii/0029554X75903365.
[43] BGO Bismuth Germanate. url: https://www.crystals.saint-gobain.com/products/bgo (visited on 12/27/2020).
[44] W. Wolszczak et al. “Temperature Properties of Scintillators for PET Detectors: a Comparative Study”. In: Nov. 2014. doi: 10.1109/NSSMIC.2014.7430768.
[45] C. L. Melcher et al. “Temperature Dependence of Fluorescence Decay Time and Emission
Spectrum of Bismuth Germanate”. In: IEEE Transactions on Nuclear Science 32 (1985),
pp. 529–532.

BIBLIOGRAPHY

105[46] H. Vincke et al. “Response of a BGO detector to photon and neutron sources: simulations
and measurements”. In: Nuclear Instruments and Methods in Physics Research Section
A Accelerators Spectrometers Detectors and Associated Equipment 484 (May 2002). doi:
10.1016/S0168-9002(01)01966-0.
[47] E. Gschwendtner. “Benchmarking the Particle Background in the LHC Experiments”.
Presented on 17 Sep 2000. 2000. url: https://cds.cern.ch/record/532546.
[48] Y. Nakajima. Status of water system and water in the tank. Dec. 2, 2020. (internal).
[49] S. Han. Estimating the ratio of neutron-less scintillation with SKG4. url: http : / /indico-sk.icrr.u-tokyo.ac.jp/event/5538/contributions/15861/attachments/
15854/18637/han_neutronless_scint.pdf (visited on 01/06/2020). (SK internal).
[50] F. Capozzi et al. “Global constraints on absolute neutrino masses and their ordering”.
In: Phys. Rev. D 95 (9 May 2017), p. 096014. doi: 10.1103/PhysRevD.95.096014. url:
https://link.aps.org/doi/10.1103/PhysRevD.95.096014.
[51] J. L. Hewett et al. Fundamental Physics at the Intensity Frontier. 2012. arXiv: 1205.2671[hep-ex].
[52] G. N. Vlaskin, Y. S. Khomyakov, and V. I. Bulanenko. “Neutron Yield of the Reaction(↵, n) on Thick Targets Comprised of Light Elements”. In: Atomic Energy 117.5 (2015),
pp. 357–365. issn: 1573-8205. doi: 10.1007/s10512-015-9933-5. url: https://doi.
org/10.1007/s10512-015-9933-5.