Measurement of Tau Neutrino Appearance andCharged-current Tau Neutrino Cross Section with
Atmospheric Neutrinos in Super-Kamiokande
byZepeng Li
Department of Physics
Duke University
Date:
Approved:

Christopher Walter, Supervisor

Kate Scholberg

Robert Behringer

Shailesh Chandrasekharan

Ying Wu

Dissertation submitted in partial fulfillment of the requirements for the degree ofDoctor of Philosophy in the Department of Physics
in the Graduate School of Duke University
2017Abstract
Measurement of Tau Neutrino Appearance andCharged-current Tau Neutrino Cross Section with
Atmospheric Neutrinos in Super-Kamiokande
byZepeng Li
Department of Physics
Duke University
Date:
Approved:

Christopher Walter, Supervisor

Kate Scholberg

Robert Behringer

Shailesh Chandrasekharan

Ying Wu
An abstract of a dissertation submitted in partial fulfillment of the requirements forthe degree of Doctor of Philosophy in the Department of Physics
in the Graduate School of Duke University
2017Copyright c 2017 by Zepeng Li
All rights reserved except the rights granted by theCreative Commons Attribution-Noncommercial LicenceAbstract
Super-Kamiokande is a 50 kiloton water-Cherenkov detector in Japan, which has
been collecting atmospheric neutrino data for more than 20 years. Tau neutrino
appearance is expected in atmospheric neutrinos due to neutrino oscillations. The
wide span of neutrino energies and neutrino path lengths in atmospheric neutrinos
and the large target mass of Super-Kamiokande allow a detection of charged-current
tau neutrino interactions in Super-Kamiokande.
This thesis describes a search for atmospheric tau neutrino appearance and a measurement of the charged-current tau neutrino cross section in Super-Kamiokande. A
neural network is applied to identify charged-current tau neutrino interactions contained in the Super-Kamiokande detector. Using 5,326 days of atmospheric neutrino
data, Super-K measures the tau normalization to be 1.47˘0.32 under the assumptionof the normal hierarchy, relative to the expectation of unity for nominal oscillation
parameters and assumed charged-current tau neutrino cross section. The result excludes the hypothesis of no-tau-appearance with a significance level of 4.6σ, thus
giving a direct evidence of neutrino flavor change due to neutrino oscillations. By
scaling the cross sections in the simulations to match Super-K data, a flux-averaged
charged-current tau neutrino cross section is measured to be p0.94 ˘ 0.20q ˆ 10´38
cm2 , compared with the theoretical flux-averaged charged-current tau neutrino crosssection of 0.64 ˆ 10´38 cm2 . This is the second reported measurement of this process.
The measured cross section is consistent with the Standard Model prediction.

iv

To my family and my girlfriend.

v

Contents
Abstractiv

List of Tables

xi

List of Figures

xiii

List of Abbreviations and Symbols

xxiv

Acknowledgements

xxvi

1 Introduction1

1.1

Neutrino History . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.2

The Standard Model Electroweak Theory . . . . . . . . . . . . . . . .

4

1.3

Atmospheric Neutrinos . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.4

Neutrino Oscillations . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

1.4.1

Discovery of Neutrino Oscillations . . . . . . . . . . . . . . . .

9

1.4.2

Neutrino Mixing . . . . . . . . . . . . . . . . . . . . . . . . .

10

1.4.3

Neutrino Oscillations in Vacuum . . . . . . . . . . . . . . . .

11

1.4.4

Neutrino Oscillations in Matter . . . . . . . . . . . . . . . . .

13

1.4.5

Experimental Status . . . . . . . . . . . . . . . . . . . . . . .

16

1.4.6

Remaining Questions . . . . . . . . . . . . . . . . . . . . . . .

20

Outline of Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

1.5

2 The Tau Neutrino
2.1

23

Direct Observation of the Tau Neutrino . . . . . . . . . . . . . . . . .
vi

23

2.2

Previous Search for Tau Neutrinos from Neutrino Oscillations . . . .

26

2.2.1

Previous Search for Tau Neutrino Appearance in Super-K . .

27

2.2.2

OPERA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

3 Super-Kamiokande

31

3.1

Cherenkov Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

3.2

Super-Kamiokande Detector . . . . . . . . . . . . . . . . . . . . . . .

32

3.2.1

Inner Detector (ID) . . . . . . . . . . . . . . . . . . . . . . . .

33

3.2.2

Outer Detector (OD) . . . . . . . . . . . . . . . . . . . . . . .

34

3.3

Photomultiplier Tubes (PMT) . . . . . . . . . . . . . . . . . . . . . .

35

3.4

Electronics and Data Acquisition (DAQ) . . . . . . . . . . . . . . . .

38

3.4.1

Inner Detector Electronics and DAQ (SK I-III ) . . . . . . . .

38

3.4.2

Outer Detector Electronics and DAQ (SK I-III) . . . . . . . .

39

3.4.3

Triggering (SK I-III) . . . . . . . . . . . . . . . . . . . . . . .

40

3.4.4

SK IV Electronics and DAQ Upgrade . . . . . . . . . . . . . .

40

Background Reduction . . . . . . . . . . . . . . . . . . . . . . . . . .

41

3.5.1

Water Purification . . . . . . . . . . . . . . . . . . . . . . . .

41

3.5.2

Air Purification System. . . . . . . . . . . . . . . . . . . . .

41

Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

3.6.1

PMT Calibration . . . . . . . . . . . . . . . . . . . . . . . . .

42

3.6.2

Water Calibration . . . . . . . . . . . . . . . . . . . . . . . . .

46

3.6.3

Energy Scale . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

3.5

3.6

4 Atmospheric Neutrino Simulations

53

4.1

Atmospheric Neutrino Flux . . . . . . . . . . . . . . . . . . . . . . .

53

4.2

Simulation of Neutrino Interactions . . . . . . . . . . . . . . . . . . .

56

4.2.1

57

Quasi-elastic Scattering . . . . . . . . . . . . . . . . . . . . . .

vii

4.2.2

Meson Exchange Current . . . . . . . . . . . . . . . . . . . . .

59

4.2.3

Single Meson Production . . . . . . . . . . . . . . . . . . . . .

59

4.2.4

Deep Inelastic Scattering . . . . . . . . . . . . . . . . . . . . .

61

Simulation of Charged-current Tau Neutrino Interactions . . . . . . .

62

4.3.1

Simulation of Tau Neutrino Interactions in Water with NEUT

62

4.3.2

Modeling the Polarization of Tau Lepton . . . . . . . . . . . .

64

4.3.3

Simulation of Tau Lepton Decay with Tauola . . . . . . . . .

64

4.4

Final State Interaction . . . . . . . . . . . . . . . . . . . . . . . . . .

66

4.5

Detector Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

4.3

5 Atmospheric Neutrino Event Reduction and Reconstruction
5.1

68

Full Contained Sample . . . . . . . . . . . . . . . . . . . . . . . . . .

69

5.1.1

FC1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

5.1.2

FC2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

5.1.3

FC3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

5.1.4

FC4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

5.1.5

FC5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

5.1.6

Final Reduction . . . . . . . . . . . . . . . . . . . . . . . . . .

80

Event reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . .

81

5.2.1

Vertex Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

5.2.2

Ring Counting . . . . . . . . . . . . . . . . . . . . . . . . . .

83

5.2.3

Particle Identification . . . . . . . . . . . . . . . . . . . . . . .

85

5.2.4

Momentum Reconstruction . . . . . . . . . . . . . . . . . . . .

88

5.3

Decay Electron Search . . . . . . . . . . . . . . . . . . . . . . . . . .

89

5.4

Summary of Reduction and Reconstruction . . . . . . . . . . . . . . .

90

5.2

viii

6 Atmospheric Tau Neutrino Appearance in Super-K

91

6.1

Atmospheric Tau Neutrino Appearance from Neutrino Oscillations . .

91

6.2

Charged-current Tau Neutrino Events in Super-Kamiokande . . . . .

94

6.3

A Neural Network Algorithm for Tau Neutrino Identification . . . . .

95

6.3.1

Input Variables to the Neural Network . . . . . . . . . . . . .

96

6.3.2

Training and Testing the Neural Network . . . . . . . . . . . . 100

7 Search for Tau Neutrino Appearance

104

7.1

Fit the Tau and Background Normalization. . . . . . . . . . . . . . 104

7.2

Systematic Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
7.2.1

Individually Estimate Systematic Uncertainties . . . . . . . . 119

7.2.2

Simultaneous Fitting of the Tau Normalization and Systematic
Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

8 Charged-current Tau Neutrino Cross Section Measurement
8.1

8.2

Charged-current Tau Neutrino Cross Section Measurement . . . . . . 134
8.1.1

Charged-current Tau Neutrino Cross Section Measurement Basedon Two-dimensional Fit of Tau Normalization . . . . . . . . . 135

8.1.2

Comparing the Number of Measured Tau Events with Expectation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137

Comparisons with Previously Reported Results . . . . . . . . . . . . 144
8.2.1

Flux-averaged Tau Neutrino Cross Section . . . . . . . . . . . 144

8.2.2

Comparison of σ{E with DONUT Result . . . . . . . . . . . . 146

9 Three-flavor Oscillation Analysis with Tau Identification
9.1

9.2

134

148

Three-flavor Atmospheric Neutrino Oscillation . . . . . . . . . . . . . 149
9.1.1

Electron Neutrino Resonance Induced by Matter Effect . . . . 149

9.1.2

Oscillation Analysis Framework . . . . . . . . . . . . . . . . . 149

Constraining the Tau Background with Tau Identification . . . . . . . 153ix

10 Conclusion

157

A Ratio of Charge in the Most Energetic Ring

159

B Neural Networks and Their Application in Particle Physics

161

C Partially Contained Reduction

163

C.1 Partially Contained Reduction . . . . . . . . . . . . . . . . . . . . . . 163
C.1.1 PC1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
C.1.2 PC2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164
C.1.3 PC3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
C.1.4 PC4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
C.1.5 PC5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
C.1.6 PC Reduction Summary . . . . . . . . . . . . . . . . . . . . . 167
Bibliography

169

Biography

179xList of Tables
1.1

Oscillation parameter values as listed in (Patrignani et al., 2016). For∆m223 , the sign is unknown and the uncertainty is on the absolute value. 20

3.1

Specification of the 20-inch PMT used in Super-Kamiokande. . . . . .

37

3.2

Specifications of the ATM. . . . . . . . . . . . . . . . . . . . . . . . .

39

4.1

Decay modes of tau leptons with branching ratio larger than 1%.
Adapted from (Patrignani et al., 2016). . . . . . . . . . . . . . . . . .

66

Summary of selection efficiencies in each stages and event rates in the
final FC samples for each Super-K period. . . . . . . . . . . . . . . .

81

Estimated background contaminations in each Super-K period. SubGeV means events with Evis ă 1.3 GeV and Multi-GeV means eventswith Evis ą 1.3 GeV. . . . . . . . . . . . . . . . . . . . . . . . . . . .

81

5.1
5.2

6.1

MC sample sizes for each stage of the MLP. The same statistics areused for SKI-IV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

6.2

The breakdown of interaction modes of both background and expectedsignal for the SK-I period. By cutting the NN output at 0.5, each modeis separated into tau-like and non-tau-like. . . . . . . . . . . . . . . . 103

7.1

Data fit results of tau normalization, background normalization and
DIS background normalization in each Super-K period. . . . . . . . . 107

7.2

Systematic errors used in tau neutrino appearance search that arecommon to all Super-K run periods. 1σ size of the systematic error isshown as a percentage. . . . . . . . . . . . . . . . . . . . . . . . . . . 119

7.3

Systematic errors used in tau neutrino appearance search that are
independent between Super-K run periods. 1σ size of the systematic
error is shown as a percentage. . . . . . . . . . . . . . . . . . . . . . . 121

xi

7.4

Change in fitted tau normalization after shifting systematic error by`1σ{ ´ 1σ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

7.5

Summary of errors of oscillation parameters. . . . . . . . . . . . . . . 123

7.6

Summary of systematic errors for both the expected and measured
number of ντ events. The errors for each category have been added inquadrature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

7.7

Fit results for common systematic errors in all Super-K run periods
in simultaneous fit of tau normalization and systematic errors. . . . . 131

7.8

Fit result of systematic errors that are independent between Super-K
run periods in the simultaneous fit tau normalization and systematic
errors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

8.1

Number of background, tau signal and data passing the event selectionin each Super-K period. . . . . . . . . . . . . . . . . . . . . . . . . . 139

8.2

exp
after shifting one systematic error by 1σ for
Change in Nτexp and NBG
common systematic errors in all Super-K run periods after selection
cuts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

8.3

exp
Change in Nτexp and NBG
after shifting one systematic error by 1σ forsystematic errors that are independent between Super-K run periods. 140

9.1

Summary of best-fit information for fits assuming the normal or inverted hierarchy with the Super-K neutrino data only (517 d.o.f.). . . 152

C.1 Estimated contamination by cosmic ray muons for PC. . . . . . . . . 168

xii

List of Figures
1.1

1.2

1.3
1.4

1.5

1.6
1.7

Feynman diagrams of neutrino interactions in the Standard Model.
(a) Feynman diagram vertex for charged-current interactions where aneutrino changes to a charged lepton (or vice versa) via the exchangeof the W ˘ vector boson with another particle. (b) Feynman diagram
vertex for neutral-current interactions where a neutrino interacts via
the exchanges of the Z 0 vector boson. . . . . . . . . . . . . . . . . . .

6

Measurements of the hadron production cross-section around the Z
resonance. The curves indicate the predicted cross-section for two,
three and four neutrino species with SM couplings and negligible mass.
Taken from (Schael et al., 2006). . . . . . . . . . . . . . . . . . . . . .

7

A schematic representation of atmospheric neutrino production in cosmic ray interactions, taken from (Kajita, 2006). . . . . . . . . . . . .

8The zenith angle distributions of µ´like and e´ like events in SuperKamiokande. Upward-going particles have cos Θ ą0 and downwardgoing ones have cos Θ ă0. The hatched regions show the Monte Carlo
expectations for no oscillations normalized to the data live-time withstatistical errors, and the solid line is the expectations with the best-fit
oscillation. Taken from (Fukuda et al., 1998). . . . . . . . . . . . . .

10Fluxes of 8 B solar neutrinos in flavor of µ or τ vs electrons deductedfrom the three neutrino reactions in SNO. The diagonal bands show
the total flux from the Standard Solar Model prediction (dashed) andthe measurement with NC reaction in SNO. Taken from (Ahmad et al.,2002). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

Neutrino CC forward elastic scattering in matter. This reaction adds
a potential to the electron neutrinos. . . . . . . . . . . . . . . . . . .

14

Ratio of the observed ν̄e to the expectation for no oscillation as afunction of L/E. Taken from (Araki et al., 2005). . . . . . . . . . . .

17

xiii

Allowed contours of sin2 θ12 and ∆m221 for Super-K solar neutrino data
(green solid), Super-K+SNO (green dashed), KamLAND (blue), andcombined results (red). Taken from (Sekiya, 2013). . . . . . . . . . .

18

Allowed contours of sin2 θ23 and |∆m232 | assuming the normal hierarchfrom T2K, NOvA, MINOS+, Super-K, and IceCube. Taken from (Abeet al., 2017). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

1.10 Left: Layout of the Daya Bay experiment. Right: Ratio of measured
versus expected signal in each detector, assuming no oscillation. The
oscillation survival probability at the best-fit value is given by thesmooth curve. Taken from (An et al., 2012). . . . . . . . . . . . . . .

19

1.11 The ordering of the neutrino mass eigenstates. For the normal hierarchy, m3 is the heaviest state and ∆m232 ą0, while for the inverted
hierarchy m3 is the lightest and ∆m232 ă0. Taken from (King andLuhn, 2013). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

1.8

1.9

2.1
2.2
2.3

2.4

2.5

2.6

Schematic view the the experimental setup of DONUT experiment.
Taken from (Kodama et al., 2001). . . . . . . . . . . . . . . . . . . .

24

Four examples of ντ CC interaction events in DONUT. The neutrinos
are incident from the left. Taken from (Kodama et al., 2001). . . . . .

25

The suppression of CC tau neutrino cross section compared to CCνµ /ν̄µ cross sections due to the tau lepton mass. Taken from (Kodama
et al., 2008). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

The likelihood (top) and neural network output (bottom) distributions of downward (left) and upward (right) events for data (dots),
atmospheric neutrino background MC (dashed histogram), and tau
and background MC (solid histogram). The events for likelihood ą
0 and neural network ą0.5 are defined as tau-like. Taken from (Abeet al., 2006). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

Zenith angle distribution for the likelihood (top) and neural network(bottom) methods. A fitted excess of tau-like events in the upwardgoing direction is shown in the shaded area. Taken from (Abe et al.,2006). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

Fit results showing projections in the NN output and zenith angle
distribution. The simulations and data sets have been combined from
SK-I through SK-III. The fitted tau signal is shown in gray. Taken
from (Abe et al., 2013). . . . . . . . . . . . . . . . . . . . . . . . . . .

29

xiv

2.7

Event display of a candidate ντ events in OPERA, in which tau decay
quickly into secondary particles after production. Taken from (Agafonova et al., 2015). . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

3.1

Schematic of the Cherenkov wavefront. . . . . . . . . . . . . . . . . .

33

3.2

A sketch of Super-K detector, and the location of the detector under
Mt. Ikenoyama. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

Schematic view of the Tyvek sheet installed in the OD that divides
the OD into top, barrel and bottom regions. . . . . . . . . . . . . . .

35

A schematic of the 20-inch PMT used in the ID. Taken from (Fukuda
et al., 2003). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

3.3
3.4
3.5

Quantum efficiency of the ID PMT, with relative spectrum of Cherenkov
light in pure water overlaid. Taken from (Nishino, 2009). . . . . . . . 36

3.6

Pulse height distribution for a single p.e. in a 20-in PMT. The peaknear zero is from dark current (Fukuda et al., 2003). . . . . . . . . . .

37

3.7

Single photoelectron distribution of a typical ID PMT. The bump at
about 2.0 pC corresponds to on p.e.. Taken from (Fukuda et al., 2003) 43

3.8

Schematic view of the relative gain measurement system. Taken from
(Lee, 2012) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

Schematic view of the relative timing measurement using a laser.
Taken from (Ishihara, 2010). . . . . . . . . . . . . . . . . . . . . . . .

45

3.10 A characteristic TQ map for an ID PMT (Abe et al., 2014). . . . . .

45

3.11 Charge linearity of all ID PMTs. The points shows the average foreach expected number of photoelectrons region. Taken from (Abeet al., 2014). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

3.12 Typical fitted water coefficient functions used in the SK-MC. The
points are the data obtained in April 2009. Taken from (Abe et al.,2014). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

3.9

3.13 Effective observed charge for through-going muons as a function ofphoton pathlengths (l), where Q is the detected charge and f pθq is the
PMT acceptance (top).Time variation of water attenuation lengths for
SK-I, SK-II, SK-III and SK-IV (bottom). Taken from (Abe et al., 2014). 48

xv

3.14 The averaged Pp.e. /Pθ where Pθ is the momentum reconstructed with
θC for SK-IV data (solid) and MC (dashed) (left), and the double ratio
between MC and data (right). Taken from (Lee, 2012). . . . . . . . .

50

3.15 Invariant mass distribution of neutrino-induced π 0 events for data(dot) and MC (boxes). MC is normalized to the live time of data. . .

50

3.16 Momentum of decay electrons for data (crosses) and MC (solid histogram). Taken from (Abe et al., 2014). . . . . . . . . . . . . . . . .

51

3.17 Energy scale uncertainty as a function of momentum for calibration
measurements in SK-IV, decay electrons (solid circle), single π 0 (solid
triangles), low-energy stopping muon (open squares) and high-energy
stopping muons (solid squares). Taken from (Ishihara, 2010). . . . . .

52

4.1

4.2

4.3

4.4

The primary cosmic ray protons at solar minimum as calculated in(Honda et al., 2004). Overlaid data points (see (Honda et al., 2004)
for references) are taken from Webber (crosses), LEAP (upward triangles), MASS1 (open circles), CAPRICE (vertical diamonds), IMAX
(downward triangles), BESS98 (circles), AMS (squares), Ryan (horizontal diamonds), JACEE (downward open triangles), Ivanenko (upward open triangles), Kawamura (open squares) and Runjob (open
diamonds). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

Zenith angle dependence of atmospheric neutrino flux at Super-K averaged over the azimuthal at 0.32 GeV (left), 1.0 GeV (middle) and
3.2 GeV (right). cos θ=1 represents vertically downward going andcos θ=-1 represents vertically upward going neutrinos. Taken from
(Honda et al., 2015). . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

The direction-averaged atmospheric neutrino flux at Super-K as afunction of neutrino energy. The left plot is absolute flux and theright plot is the flux ratios. Honda flux model is shown in red. Other
model predictions are compared with Honda model to estimate thesystematic uncertainties on the flux. Taken from (Honda et al., 2015).

56

CCQE cross sections of νµ (top) and ν̄µ from NEUT calculation, andoverlaid measurement data from ANL (Barish et al., 1977), Gargamelle(Bonetti et al., 1977; Pohl et al., 1978), BNL (Ahrens et al., 1986), Serpukhov (Belikov et al., 1985) and SKAT (Brunner et al., 1990). Solidlines represent scattering off a free proton, and dashed lines represent
scattering off bound nucleons in 16 O. Taken from (Hayato, 2009). . .

58

xvi

4.5

Cross section for ν̄µ charged current single pion productions. Solidlines are NEUT calculations, points are experimental data. Taken
from (Hayato, 2009). . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

4.6

Charged-current tau neutrino cross sections in NEUT.

. . . . . . . .

63

4.7

The contour plot of the DIS τ ´ polarization in the plane of pτ cos θand pτ sin θ for ντ charged-current DIS process at Eν =10 GeV. The
thick gray line represents the kinematic boundary. The QE scattering
contributed along the boundary, and the RES process contributes justinside of the boundary. The length of the arrows shows the degree ofpolarization, and the direction gives the τ ´ spin in the τ ´ rest frame.
Taken from (Hagiwara et al., 2003). . . . . . . . . . . . . . . . . . .

65

4.8

Feynman diagram of a tau lepton decay. . . . . . . . . . . . . . . . .

65

4.9

π ` ´12 C scattering cross sections, where the cross sections for eachinteraction process are shown separately by the colors described in thelegend. The data pions are the measurement from π ` ´ C scatteringexperiments and the solid lines are the cross sections from the tuned
NEUT pion cascade model. The dashed lines are for the previous
NEUT pion FSI model. Taken from (de Perio, 2011). . . . . . . . . .

67

5.1

Schematic view of three types of neutrino events in Super-K. . . . . .

69

5.2

Distributions of the atmospheric MC and data before the FC2 cuts 2a(left) and 2b (right). MC events are high-energy FC neutrino events
in FV, and data represent the events that pass the FC1 reduction.
Taken from (Lee, 2012). . . . . . . . . . . . . . . . . . . . . . . . . .

71

Two-dimensional scatter plot of NHITAin (x-axis) and NHITAout forthrough-going muons selection in FC3 reduction for SK-III. (a) shows
the data after the FC2 reduction. (b) shows MC events after the FC2
reduction. (c) shows the final sample for data (bottom left) and Monte
Carlo (bottom right) events. Cut values are indicated by solid lines.
Taken from (Lee, 2012). . . . . . . . . . . . . . . . . . . . . . . . . .

72

Number of hit OD PMTs around the entrance point of muons forevents in FC3 stopping muon reduction in SK-IV. (a) shows the data
after FC2 reduction. (b) shows FC MC events after FC2 reduction.
(c) shows the final samples for data and MC after FC reduction. Taken
from (Lee, 2012). . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

74

5.3

5.4

xvii

5.5

Schematic view of the positions of the twelve cable holes on the stainless steel support. Four holes in red circles are installed with vetocounter. Taken from (Fukuda et al., 2003). . . . . . . . . . . . . . .

75

Schematic view of a cable hole muon and a veto counter. Taken from
(Fukuda et al., 2003). . . . . . . . . . . . . . . . . . . . . . . . . . .

75

The timing distribution of a flasher (top) and an atmospheric neutrino
MC event (bottom). The arrows indicate the time window for counting
NMIN100 . Taken from (Lee, 2012). . . . . . . . . . . . . . . . . . . .

76

NMIN100 distributions for events with more than 800 hit ID PMTs
(left) and less than 800 hit ID PMTs (right). (a) shows the data after
the FC 2nd reduction. (b) shows FC Monte Carlo events after the
FC 2nd reduction. (c) shows the final samples for data and the Monte
Carlo events, in which the number of Monte Carlo events is normalized
to that of the data. Selection cuts are indicated by arrows. Taken from
(Lee, 2012) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

The Hough transform method. Virtual rings with a 42-degree opening
angle are drawn around each hit PMT, and the ring center is foundat the point at which most virtual rings overlap. . . . . . . . . . . .

84

5.10 The ring likelihood distribution for sub-GeV events (left) and multiGeV events (right) for SK-I. The data is shown in the black dots,
and the MC with neutrino oscillations applied at ∆m2 =2.1ˆ10´3 eV2and sin2 p2θq=1.0 is shown as the histogram. A cut in the likelihood
separates single-ring from multi-ring events. . . . . . . . . . . . . . .

84

5.11 Simulated e´ and µ´ , both with momentum 1000 MeV/c. The ringpattern is sharper for µ´ and more diffuse for e´ . The color scale
indicated the amount of charge deposited on the PMTs. The main
display shows the ID, while the small display shows on the top rightshows the OD, in which there is negligible activity. . . . . . . . . . .

86

5.12 The distribution of particle identification likelihood for sub-GeV (top)and multi-GeV (bottom) single-ring FC events. The hatched histogram shows the MC contribution due to CC νµ interactions, theblack dots show the data and the empty histogram shows the total
MC with neutrino oscillations applied as in Fig. 5.10 . . . . . . . . .

88

5.6
5.7

5.8

5.9

6.1

Probability of tau neutrino appearance from muon neutrinos (left) andelectron neutrinos (right) as a function of neutrino energy and zenith
angle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xviii

93

6.2

Flux of atmospheric ντ (red) and ν̄τ (blue) from neutrino oscillations
as a function of neutrino energy at Super-Kamiokande for the threeflavor oscillation parameters given in the text. . . . . . . . . . . . . .

93

6.3

Feynman diagrams of weak interactions that include a tau neutrino. .

94

6.4

Total charged-current cross section as a function of neutrino energy
for ντ (red) and ν̄τ (blue) for neutrino energy between 3.5 GeV and
70 GeV in NEUT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

95

6.5

Feynman diagram of a tau lepton decay. . . . . . . . . . . . . . . . .

96

6.6

Simulation of a background (left) event with 2.2 GeV visible energy
and a CC ντ event (right) with 2.3 GeV visible energy in Super-K,
which produces multiple rings in the ID. . . . . . . . . . . . . . . . .

97

Simulation of a background event that has a deep-inelastic interaction,
which produces multiple rings in the detector. . . . . . . . . . . . . .

97

6.7
6.8

The input variables to the neural network for the downward MC (black
histogram) and data (black dots), along with tau MC (gray shaded
histogram). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

6.9

Energy weights for the training sample. . . . . . . . . . . . . . . . . . 102

6.10 Neural network output of training (filled histogram) and test (dots)
sample for background (red) and signal (blue) simulations. . . . . . . 102
6.11 Efficiency of signal and background rejection by cutting the NN output
in SK-IV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
7.1

2-dimensional probability distribution likelihood as a function of zenith
angle and neural network output for tau (left) and background (right)
built with MC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

7.2

The breakdown of background simulation into DIS (orange) and NonDIS (red) as a function of neural network output. . . . . . . . . . . . 105

7.3

Tau normalization fit results with 2D unbinned likelihood fit of toy
MC generated with α “ 1. The distribution has a σ of 0.24 with aGaussian fit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

7.4

Fit of tau normalization in toy MC study with tau normalizations set
to 1.37 in the toy MC generation. The distribution is centered at 1.37,
and the fit of the distribution with a Gaussian function found a σ of
0.26. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

xix

7.5

Fit results showing the projections in NN output and zenith angle
distribution for tau-like (NN ą 0.5), upward-going (cospθq ą 0.1), nontau-like (NN ă 0.5) and downward-going (cospθq ă 0.1) events for boththe two-dimensional PDFs and data. The PDFs have been rescaledby the fit results. The PDFs and data sets have been combined from
SK-I through SK-IV. The fitted tau signal is shown in gray. . . . . . . 109

7.6

Uncertainty on the absolute flux normalization as a function of neutrino energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

7.7

Flavor ratio (νµ ` ν̄µ )/(νe ` ν̄e ) in FLUKA model (solid) and Bartol
model (dashed) normalized by Honda flux model. . . . . . . . . . . . 111

7.8

Anti-neutrino/neutrino ratio normalized by ratio in Honda flux modelfor ν e {νe (left) and ν µ {νµ (right) for FLUKA (solid) and Bartol (dashed)models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

7.9

Neutrino path length as a function of the zenith angle calculated withthe atmosphere density in Honda model (solid) and the atmospheredensity compressed by 10% (dashed). Taken from (Lee, 2012) . . . . 113

7.10 The Hernandez model normalized by the base model in NEUT forνµ p Ñ µ´ pπ ` (solid line), νµ n Ñ µ´ pπ 0 (dashed line) and νµ n Ñµ´ pπ ` (dotted line) are shown in for the theoretical cross sections(left) and the double ratio for ν̄{ν ratio (right). Figure is taken from
(Lee, 2012). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
7.11 Re-weighted simulation of π ` scattering on carbon target with NEUT
and experimental data in: quasi-elastic channel (top left), absorption
channel (top right), total reaction channel (bottom left) and single
charge exchange channel (bottom right). Taken from (Lee, 2012). . . 117
7.12 NHITAC distributions for FC and PC final samples for SK-I (top left),
SK-II (top right), SK-III (bottom left) and SK-IV (bottom right).
Data are shown by dots with error bar, MC events assuming no oscillation are shown with dashed lines and MC with νµ Ñ ντ 2-flavor
oscillation with (sin2 2θ, ∆m2 ) = (1.00, 2.5 ˆ 10´3 eV2 ) are shownwith solid lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
7.13 An asymmetric Gaussian function is constructed with the fitted tau
normalization and it uncertainty. The significance of tau appearance
is calculated by integrating the negative part of the function. . . . . . 124

xx

7.14 An example histogram of the PDF of K/π ratio uncertainty for SK-I
background. The vertical axis is the output of the NN; the horizontal
axis is the cosine of the event zenith distribution. . . . . . . . . . . . 126
7.15 Tau normalization in toy MC generated and fitted with all atmospheric
neutrino analysis systematic errors, which is centered at 1.0 and has
a p value of 3.6ˆ10´4 . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
7.16 Histograms of the PDFs for uncertainty on absolute normalization ofneutrino flux below 1 GeV that causes small change in event distributions in SK-I. The color demonstrates the change of events number inthe bin after shifting the systematic error by 1σ. . . . . . . . . . . . 128
7.17 Tau normalization in toy MC generated and fitted with selected atmospheric neutrino analysis systematic errors, which is centered at 1.0
and has a p value of 4.4ˆ10´4 . . . . . . . . . . . . . . . . . . . . . . . 129
7.18 Fitted values of systematic error in simultaneous fit of tau normalization and systematic errors. . . . . . . . . . . . . . . . . . . . . . . . . 130
7.19 Fit results showing projections in the NN output and zenith angle
distribution for tau-like (NNą0.5), upward-going [cospθq ą0.1], nontau-like (NN ă0.5) and downward-going [cospθq ă0.1] events for boththe two-dimensional PDFs and data. The PDFs and data sets havebeen combined from SK-I through SK-IV. The fitted tau signal isshown in gray. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
7.20 Distribution of the fitted tau normalization in toy MC generated with
1.47 times of the nominal expected tau events. . . . . . . . . . . . . . 133
8.1

Energy spectrum of tau neutrino events after selection in SK-I through
SK-IV simulation after weighting by oscillation probabilities and live
time (left), and integral of number of events in the left plot with energy
above the energy value in the bin (right). . . . . . . . . . . . . . . . . 136

8.2

Theoretical charged-current tau neutrino cross sections (solid lines) as
a function of neutrino energy for ντ (blue) and ν̄τ (red), atmospheric
tau neutrino flux (dashed lines) from neutrino oscillations by oscillating the electron and muons neutrinos in HONDA model prediction to
tau neutrinos for ντ (blue) and ν̄τ (red) and flux-averaged tau neutrino
cross section (gray dashed line). . . . . . . . . . . . . . . . . . . . . 137

xxi

8.3

Measured flux-averaged charged-current tau neutrino cross section(black), together with theoretical differential cross sections (ντ in blueand ν̄τ in red), flux-averaged theoretical cross section (dashed gray)and tau events after selection in MC simulations (gray histogram). . . 138

8.4

Visible energy of events passing the cuts in Table 8.1. Expectedbackground is subtracted from data, shown in dots, with statisticaluncertainties from data (orange) and systematic uncertainties frombackground simulations. Expected tau events are shown in the blue
histogram, with statistical (blue) and systematic (shaded green) uncertainties. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141

8.5

Sτ distribution in toy MC experiments generated with expected signal,
background and data after selections. The solid black line shows thecentral value of Sτ as 1.5. The dashed line denotes the boundaries of
68% C.L. region, which is 0.91-2.15. . . . . . . . . . . . . . . . . . . . 142

8.6

Flux-averaged tau neutrino cross section between 3 GeV and 70 GeV. 143

8.7

Tau normalization in toy MC generated and fitted with one-dimensional
PDFs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143

8.8

Suppression of tau neutrino cross section relative to electron and muon
neutrinos due to the tau lepton mass. Taken from (Kodama et al., 2001).145

8.9

Comparison of the Super-K measured flux-averaged CC tau neutrino
cross section with DONUT result reweighted by the Super-K tau
neutrino flux. DONUT cross section digitized from (Kodama et al.,2008), and weighted by the tau neutrino from at Super-Kamiokande,ă σ ą“ p0.37 ˘ 0.18q ˆ 10´38 cm2 . . . . . . . . . . . . . . . . . . . . . 146

9.1

Probabilities of νµ Ñ νe (left) and ν̄µ Ñ νe under the assumptionof the normal hierarchy (top) and inverted hierarchy (bottom). It
is noted that resonance in the upward multi-GeV region shows up
in different neutrino types under the assumptions of the normal or
inverted hierarchy. Oscillation probabilities are calculated assuming
the oscillation parameters: ∆m212 “ 7.6 ˆ 10´5 eV 2 , ∆m232 “ 2.4 ˆ
10´3 eV 2 , sin2 θ12 “ 0.3, sin2 θ23 “ 0.5, sin2 θ13 “ 0.0238 and δCP “ 0.

9.2

150

Definition of the binning scheme in the three-flavor oscillation analysis.
410 bins in FC, 60 bins in PC, 50 bins in UPMU. White boxes indicatethat momentum bin is uniformly divided into 10 zenith angle bins
between -1 and 1 for FC and PC, between -1 and 0 for UPMU. Shadedboxes indicate that the momentum bin is not further divided into
zenith angle bins. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
xxii

9.3

Constraints on neutrino oscillation parameters from the SKI-IV atmospheric neutrino data. The orange line denotes the result in the IH
assumption, which has been offset from the NH result, shown in blue,by the difference of their ∆χ2 . . . . . . . . . . . . . . . . . . . . . . . 153

9.4

Zenith angle distribution of five samples that include electron appearance signal in the Super-K three-flavor neutrino oscillation analysis.
The MC expected under the assumption of the normal (inverted) hierarchy is shown in blue (orange) histogram, without taking tau events
into account. The tau events are shown in the green filled histogram. 154

9.5

Distribution of neural network output for events with 3.0 ă log Evis ă
3.4 and ă ´1 cos θ ă-0.8 in the multi-GeV e-like νe sample, corresponding to one bin the standard three-flavor oscillation analysis. The
simulations are scaled to the live time of 5326 days in SK-I through
SK-IV. The νe and νµ events are shown in the black histogram, while
the tau events are shown in red filled histogram on top of the otherevents. The single bin is divided into two bins at the neural network
equal to 0.25. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

9.6∆χ2 distribution as a function of sin2 θ23 for standard binning (red)and with additional binning in neural network output (yellow). Under
the assumption of true (inverted) normal hierarchy, the ∆χ2 between
normal and inverted hierarchies is shown in solid (dashed) lines. . . . 156

A.1 The fraction of energy carried by the most energetic ring for data
(dots) and MC (histogram) in SK-I (left) and SK-II (right). SK-III
and SK-IV have similar data and MC agreement as in SK-I. . . . . . 160
A.2 Ratio of charge carried by most energetic ring calculated with PMT
charge in each SK period. . . . . . . . . . . . . . . . . . . . . . . . . 160
B.1 Structure of a multilayer perceptron with seven input nodes in the
input layer, one output node in the output layer, and ten nodes in thehidden layer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
C.1 Scatter plot of Distsegclust and NHITAsegclust in MC (left) and data(right) before the last corner clipper cut in PC5. The data eventswith NHITAsegclust ą30 and Evis ă 1500 MeV are mostly background
events tagged in the eye-scan. . . . . . . . . . . . . . . . . . . . . . . 167

xxiii

List of Abbreviations and Symbols
Abbreviations
All abbreviations are defined in the text, but this section is designed to be a helpfulreference for some of the more common ones.
CC
CCQE

Charged-current
Charged-current Quasi-Elastic

DIS

Deep Inelastic Scattering

ES

Elastic Scattering

Evis

Visible Energy

FC

Fully-contained

FSI

Final State Interaction

FV

Fiducial Volume

ID

Inner Detector

MC

Monter Carlo

MEC

Meson Exchange Current

MLP

Multilayer Perceptron

NC

Neutral-current

NN

Neural Network

OD

Outer Detector

PC

Partially-contained
xxiv

PID

Particle Identification

PMT

Photomultiplier Tubep.e.

photoelectron

SK

Super-Kamiokande

UPMU

Upward-going Muon

VME

Versa Module Europaxxv

Acknowledgements
I want to first express my gratitude to my advisors, Professor Chris Walter and Kate
Scholberg. They have guided me with their knowledge, patience and encouragementin the last six years. They are always ready to answer any of my questions on research.
This thesis would not have been completed without their precious guidance.
I also want to thank the great people in the Duke neutrino group: Dr. Josh
Albert, Dr. Tarek Akiri, Dr. Alex Himmel, Dr. Erin O’sullivan, Dr. Taritree
Wongjirad, Gleb Sinev and Justin Raybern. I learned a lot during the many talks
with you. I would like to especially thank Dr. Alex Himmel for guiding me in researchand taking me to Japanese restaurants. Alex can always provide valuable help incoding, physics, and academic writing. I began to enjoy sushi after he took me tosushi restaurant when we were in Japan. It has been a great experience working
with these great people.
It is really enjoyable working on the Super-Kamiokande experiment, especiallyin the ATMPD group. Prof. Ed Kearns and Prof. Masato Shiozawa provided
valuable advise in this thesis project. Prof. Masato Shiozawa also helped to prepare
document for my Japanese visa application for my every travel to Japan. Prof. Roger
Wendell can always answer my questions on the ATMPD analysis software. Prof.
Yoshinari Hayato taught me a lot on neutrino cross sections. There are many othercollaborators from Japan and US that I want to thank. The analysis in this thesis
relies on the great work from these collaborators.

xxvi

I would like to thank my friends and colleagues at Duke. You helped me to getused to Durham and enjoy the life in Durham when I first came to US. You took
me to many good Chines restaurants. Among these people, I owe Dong Wang a big
thank because he drove me to the airport at 5 am and picked me up from the airport
at 12 am in most of my trips to Japan or China. I forgot how many times you havedone it. But I still remember you are ready to give a ride every time I turn to you
for help. I feel lucky to meet so many good friends at Duke.
Finally, I want to thank my family. They have supported me all the way to my
PhD degree. My parents do not know particle physics, but they always encourage me
to chase my dream. My brother has given me a lot of help in life and study. He used
his own experience to show me the wonderful life in college when I was still in middle
school. It was especially important for a child from a rural area. He encouraged me
to study abroad when I was in college. I could not have reached today without him.
I also want to thank my girlfriend, Guirong Liu, for her love.

xxvii

1
Introduction

The neutrino is one of the fundamental particles that make up the universe. It is thesecond most abundant particle in the Standard Model, after only photons, with tensof billions of neutrinos passing through every cubic centimeter of the Earth every
second. Ths abundant particle is mysterious because it only interacts with othermatter via weak interactions. The discovery of non-zero neutrino mass in the late
1990s opened the new field of neutrino oscillation. The study of neutrino oscillations
has broadened our view of nature.
In this dissertation, I will talk about the study of one type of neutrino oscillations, atmospheric tau neutrino appearance in the Super-Kamiokande experiment.
Super-Kamiokande is a large water Cherenkov detector in Gifu, Japan, designed tosearch for proton decay and to study neutrinos from various sources. It is credited forthe discovery of the neutrino oscillations by measuring a deficit in the atmosphericmuon neutrinos in 1998. This thesis describes a search for atmospheric tau neutrino
appearance in Super-Kamiokande. Tau neutrino appearance is expected in the atmospheric neutrinos from three-flavor neutrino oscillations. Benefiting from the wide
energy range of atmospheric neutrinos, Super-K is expected to detect charged-current
1tau neutrino interactions in its large target mass. By measuring charged-current
tau neutrino interactions, a measurement of tau neutrino appearance in Super-K
is performed. This measurement is a direct confirmation of three-flavor neutrino
oscillations. Additionally, the large expected sample of tau neutrinos in Super-K
can shed light on the understanding of tau neutrino interactions. There is only one
previously reported measurement of the charged-current tau neutrino cross section.
This thesis reports a measurement of the cross section with a smaller uncertainty.
The tools developed in the analysis are capable of identifying tau events statistically.
Including this sample will increase the sensitivity to the neutrino mass hierarchy in
the Super-K three-flavor atmospheric neutrino oscillation analysis.

1.1 Neutrino History
The study of β decay provided the first hint of the existence of the neutrino. The
process was thought to be a two-body decay which should produce a discrete energy
spectrum as in alpha radiation. However, the spectrum of the beta particle wasmeasured to be continuous, which was in apparent contradiction to the law of conservation of energy. The neutrino was postulated by Wolfgang Pauli (Pauli, 1930) in
1930 to explain the continuous energy spectrum of β decay. In Pauli’s new theory, anelectrically neutral particle that was not detected in the β decay was responsible forcarrying away the energy. By adding this new particle, energy and momentum areconserved in β decay. Based on the evidence at that time, Pauli also predicted thatthe new neutral particle had a small mass and 1/2 spin, and named it the “neutron”.
An electrically neutral particle was discovered soon by Chadwick and he also
named it a neutron. However, that particle is much more massive than that postulated by Pauli. We now know that Chadwick discovered what is now known as
the neutron. Later, in 1934, Enrico Fermi developed a theory of beta decay (Fermi,
1934). In this theory, the neutron discovered by Chadwick decays into a proton,2an electron and the “neutron” postulated by Pauli. Fermi renamed the particle theneutrino, meaning “little neutral one” in Italian. Fermi’s theory of beta decay posits
four fermions interacting directly with one anothern Ñ p ` e´ ` ν̄e .

(1.1)

Fermi’s theory of beta decay has a four-point vector current, which only describes thebeta decays between nuclear states with identical angular momenta. Later in 1936,
Fermi’s theory was extended by Gamow and Teller to explain some types of beta
decay that involved a change in 1 unit of nuclear spin (Gamow and Teller, 1937).
Two decades after Pauli’s proposal, the neutrino was first discovered by Reines
and Cowan in 1956 (Cowan et al., 1956) via an inverse beta decay process: ν̄e `p Ñ n ` e` . Using CdCl2 -doped water, Reines and Cowan observed this reaction
through the coincidence measurements of gamma rays from positron annihilationand the delayed gamma rays from the neutron’s capture on cadmium. Later in 1962,
Lederman, Schwartz and Steinberger confirmed the existence of the muon neutrino(Danby et al., 1962) using a neutrino beam produced by the decay of charged pions
into a beam of muon neutrinos. The experiment found that the interactions ofneutrinos in the beam produced muons but no electrons, confirming that the muon
neutrino is a distinct particle. The tau neutrino was discovered in 2001 by the
DONuT experiment (Kodama et al., 2001) through the detection of tau leptons inneutrino interactions. Electron neutrino, muon neutrino and tau neutrino make up
three generations of neutrinos, corresponding to three generations of charged leptons(electron, muon and tau).
In 1955, Davis exposed 3,800 liters of CCl4 to the Brookhaven Graphite Research
Reactor to search for the reaction of ν̄e ` 37 Cl Ñ e´ ` 37 Ar (Davis, 1955). By notobserving the process, it was confirmed that νe and its anti-particle ν̄e are different.
The same applies to muon neutrinos, which was later found at CERN because in3νµ interactions only µ´ in the final state were detected, but not µ` (Bienlein et al.,1964).
Parity is defined as a symmetry transformation by an inversion of one spatial
coordinate at the origin. In 1956, Lee and Yang found that parity conservation had
never been tested for weak interactions and suggested a number of tests (Lee andYang, 1956). Wu then measured the beta decays of polarized 60 Co nuclei (Wu et al.,1957). The decay of 60 Co is given by60

Co Ñ 60 Ni˚ ` e´ ` ν̄e .

(1.2)

It was observed that electrons are preferentially emitted in the opposite spin direction to that of the 60 Co nucleus, showing that parity is maximally violated in weakinteractions. The experiment implies that the neutrino exist in only one helicity state
where helicity is defined as the projection of particle’s spin along the direction of its
momentum. Later in the classic experiment by Goldhaber that used the electron
capture of 152 Eu (Goldhaber et al., 1958), the helicity of the neutrino was measured
to be -1.
The discovery of parity violation suggested that a new approach was needed todescribe the true form of weak interactions. In 1957, Robert Marshak and George
Sudarshan (Sudarshan and Marshak, 1994) and Richard Feynman and Murray Gellmann (Feynman and Gell-Mann, 1958) independently proposed a V-A (vector minus
axial vector) Lagrangian for the weak interaction. In this theory, the weak interaction
only acts on left-handed particles and on right-handed antiparticles, and maximizes
parity violation.

1.2 The Standard Model Electroweak Theory
The Standard Model of particle physics is a theory that describes the electroweak(a unification of the weak and electromagnetic interactions) and strong interactions.
4

In the Standard Model, neutrinos only interact via the weak force; thus electroweak
theory is the only part of the Standard Model that is discussed below.
Schwinger was the first to suggest the existence of intermediate vector bosons in
weak interactions. In 1957, he proposed a gauge theory of weak interactions thatare mediated by two charged bosons (which he named Z ` and Z ´ ) (Schwinger,
1957). In the same paper, he also first suggested the unification of electromagnetic
interactions and weak interactions. The unification was first formulated by Glashow
as a gauge field with SU(2)ˆU(1) symmetry (Glashow, 1961), and a neutral gaugeboson Z 0 was added in the new theory. Salam and Ward proposed a similar model
based on SU(2)ˆU(1) symmetry independently (Salam and Ward, 1964). The SU(2)symmetry defines the weak isospin group. The leptons and quarks that participate
in the weak interactions are grouped into isospin doublets, and the leptons that
do not interact via the weak interaction are defined as weak isospin singlets. The
spontaneous breaking of symmetry through the Higgs mechanism leads to a set of
massive vector bosons, W ` , W ´ and Z 0 that mediate the weak interaction (Weinberg, 1967; Salam, 1968). This new theory of electroweak interactions is also known
as Glashow-Weinberg-Salam theory.
In the Standard Model, the neutrinos only interact with other matter throughthe exchange of these bosons. The interactions are classified into charged-current
interactions and neutral-current interactions based on which boson takes part in
the interaction. In the charged-current interaction, a neutrino changes to a charged
lepton of the same flavor (or vice versa) via exchanging of a W ˘ vector boson withanother particle. In the neutral-current interaction, a neutrino exchanges a Z 0 vectorboson with another particle, where the flavor of the neutrino is irrelevant. Fig. 1.1
shows the Feynman diagrams of neutrino interaction vertices.
The Standard Model electroweak theory has withstood all experimental challenges since its inception. In 1973, neutral current interactions were observed in the
5

νl

νl

l

νl

Z0

W˘

(b) Neutrino-current interactions.

(a) Charged-current interactions

Figure 1.1: Feynman diagrams of neutrino interactions in the Standard Model. (a)Feynman diagram vertex for charged-current interactions where a neutrino changes
to a charged lepton (or vice versa) via the exchange of the W ˘ vector boson withanother particle. (b) Feynman diagram vertex for neutral-current interactions wherea neutrino interacts via the exchanges of the Z 0 vector boson.

Gargamelle experiment at CERN (Hasert et al., 1973) and at Fermilab (Benvenuti
et al., 1974), which was proposed in the GWS theory but not yet observed at thetime the theory was proposed. The W ˘ vector bosons were discovered in the UA1and UA2 experiments at CERN (Arnison et al., 1983; Banner et al., 1983), and theZ 0 was discovered later in the same experiments. Both experiments measured themasses of three vector bosons that were in good agreement with the Standard Model
prediction. The existence of the Higgs field which is responsible for the mass of W ˘and Z 0 bosons was also proved by the discovery of the Higgs boson in the ATLAS
and CMS experiments at CERN (Aad et al., 2012; Chatrchyan et al., 2012).
The number of light neutrino types can also be determined by measuring thewidth of the Z 0 boson in collider experiments. Z 0 bosons decay into quarks andcharged leptons as well as into neutrinos which are invisible in collider experiments.
The invisible partial width, Γinv , due to Z 0 decays into neutrinos is determined bysubtracting the measured partial widths, corresponding to Z 0 decays to quarks andcharged leptons, from the total Z 0 width. Fig. 1.2 shows the measurement of thehadron production cross section around the Z 0 resonance from four LEP experiments
compared with predicted cross section for two, three and four types of neutrino. The

6

combined result from four LEP experiments measured the number of neutrino types
to be 2.984˘0.008.

Figure 1.2: Measurements of the hadron production cross-section around the Z
resonance. The curves indicate the predicted cross-section for two, three and four
neutrino species with SM couplings and negligible mass. Taken from (Schael et al.,2006).

1.3 Atmospheric Neutrinos
The analysis in this thesis is performed with the study of atmospheric neutrinos.
Atmospheric neutrinos are produced in the interactions of cosmic rays with nuclei inthe upper atmosphere, as shown in Fig. 1.3. The interaction produces a shower ofhadrons, mostly pions and kaons. The dominant decay (branching ratio ą 99%) ofcharged pions produces a muon and a neutrino, and the muon further decays into aelectron, one electron neutrino and one muon neutrino.

π ` Ñ µ ` ` νµë e` ` νe ` ν̄µ7(1.3a)(1.3b)

Figure 1.3: A schematic representation of atmospheric neutrino production incosmic ray interactions, taken from (Kajita, 2006).

π ´ Ñ µ´ ` ν̄µë e´ ` ν̄e ` νµ(1.3c)(1.3d)

Based on this kinematic chain, the ratio of muon neutrinos to electron neutrinos
is expected to be around 2 below 1 GeV. The ratio is modified at higher energies
where neutrinos are produced by decays of more energetic parents. Some muons produced in these decays reach the Earth before decaying due to Lorentz time dilation.
Consequently, no electron neutrinos are produced and the ratio of muon to electron
neutrino increases with energy.
The flux of cosmic rays above about 10 GeV is isotropic about the Earth; thus
the flux of neutrinos is expected to be equal for equal angles with the zenith andnadir. The path length of atmospheric neutrinos varies from Op10q km to Op104 q
km. Due to the wide energy span of cosmic rays, atmospheric neutrinos have energies
spanning from Op10q MeV to above PeV. The wide span in path length and energy
makes atmospheric neutrinos a robust souce for studying neutrino oscillations.

8

1.4 Neutrino Oscillations
1.4.1

Discovery of Neutrino Oscillations

The Sun’s energy comes from nuclear fusion reactions in which neutrinos are also
produced. The flux of neutrinos at the Earth is several ten billion per square centimeter per second. In the late 1960s, the Homestake experiment first measured theflux of neutrinos from the Sun, and found a deficit compared to the prediction of
the Standard Solar Model (Davis, 1994). The deficit was confirmed by subsequent
experiments, such as Kamiokande (Fukuda et al., 1996), SAGE (Abdurashitov et al.,1999) ,GALLEX (Hampel et al., 1999) and Super-Kamiokande (Fukuda et al., 2001).
The discrepancy is known as the solar neutrino anomaly (Bahcall et al., 1996).
Another anomaly was later found in the atmospheric neutrinos. As discussed in
Section 1.3, the ratio of νµ {νe is expected to be around 2 at energies below 1 GeV.
However, atmospheric neutrino experiments measured a value significantly smaller
than 1 (Fukuda et al., 1994; Becker-Szendy et al., 1992; Allison et al., 1999), indicating that either there was less νµ in the data than in the prediction, or there wasmore νe , or both.
The efforts to resolve these two anomalies led to the discovery of neutrino oscillations. The atmospheric neutrino anomaly was first solved by the successorof Kamiokande experiment, Super-Kamiokande, which measured the flux of atmospheric neutrinos with unprecedented precision. The Super-Kamiokande collaboration announced that they had found the evidence for neutrino oscillations in atmospheric neutrinos in 1998 (Fukuda et al., 1998). As shown in Fig. 1.4, the measured
neutrino flux showed dependence on zenith angle for muon neutrinos, but not forelectron neutrinos. The deficit in the observed νµ coming from below, but not fromabove, was consistent with νµ Ñ ντ oscillation.
The solar neutrino anomaly was later solved by the Sudbury Neutrino Observa9

Figure 1.4: The zenith angle distributions of µ´like and e´ like events in SuperKamiokande. Upward-going particles have cos Θ ą0 and downward-going ones havecos Θ ă0. The hatched regions show the Monte Carlo expectations for no oscillations
normalized to the data live-time with statistical errors, and the solid line is the
expectations with the best-fit oscillation. Taken from (Fukuda et al., 1998).

tory (SNO) collaboration in 2002, which found direct evidence for neutrino flavor
transformation in solar neutrinos (Ahmad et al., 2002). SNO used a tank of D2 O
as its target to detect charged-current (CC) interactions, electron elastic scattering(ES), and neutral-current (NC) interactions. The CC interaction is only sensitive toνe therefore only measures νe flux. For ES, electron neutrinos can interact via boththe charged and neutral currents, while νµ and ντ only interact via neutral current at
solar neutrino energies. This channel measures a combination of the electron, muonand tau neutrino flux. The NC interaction is sensitive to all neutrino flavors, thus
measures the total flux. As shown in Fig. 1.5, SNO measured a total flux consistent
with the Standard Solar Model prediction, but a deficit of electron neutrinos. The
result is consistent with νe Ñ νµ{τ oscillation. The solar neutrino anomaly was finally
understood, and neutrino flavor transformation was confirmed.
1.4.2

Neutrino Mixing

Neutrino oscillations occur because neutrino flavor eigenstates are not identical to themass eigenstates. A neutrino flavor eigenstate is composed of a linear combination
10

Figure 1.5: Fluxes of 8 B solar neutrinos in flavor of µ or τ vs electrons deductedfrom the three neutrino reactions in SNO. The diagonal bands show the total flux
from the Standard Solar Model prediction (dashed) and the measurement with NC
reaction in SNO. Taken from (Ahmad et al., 2002).

of mass eigenstates, ν1 , ν2 , ν3 , through a mixing matrix, which is analogous to the
CKM matrix used for quarks. The flavor eigenstates of neutrinos can be expressedby mixing of mass eigenstates as:
¨ ˛ ¨
˛¨ ˛
νe
Ue1 Ue2 Ue3
ν1˝νµ ‚ “ ˝Uµ1 Uµ2 Uµ3 ‚˝ν2 ‚ντUτ 1 Uτ 2 Uτ 3
ν3

(1.4)

the matrix is known as Pontecorvo-Maka-Nakagawa-Sakata (PMNS) matrix, which
is a complex unitary matrix. Therefore, a neutrino of flavor α can be expressed as:
να “ÿ

Uαi νi .

(1.5)

i

1.4.3

Neutrino Oscillations in Vacuum

Suppose a neutrino beam is produced containing flavor eigenstate να at the originof x-axis and time 0 and the neutrinos are aimed along the x-axis in free space, themass eigenstates propagate according to the time-dependent Schrödinger equation

11

with no potentialsiB|νi px, tqy “ H |νi px, tqy
Bt

(1.6)

The time evolution of the wave-function is determined by the time evolution of themass eigenstates and can be written as|νk px, tqy “ e´iφk |νk y ,(1.7)

where φk “ pi ¨ xi is the phase of the plane wave representing the mass eigenstate.
We can transform the Schrödinger equation to be expressed in the flavor basis as:|να px, tqy “ÿUαk e´iφk

k

ÿ

˚
Uβk
|νβ y “βÿÿβ

˚ ´iφkUβke
Uαk |νβ y .

(1.8)

k

The amplitude for detecting a neutrino of flavor β is
Aαβ “

ÿ

˚
e´iφk Uαk
Uβk ,

(1.9)

k

The probability of a neutrino changing to a different flavor is then:P pνα Ñ νβ q “|Aαβ |2“δαβ ´ 4ÿ˚˚RepUαiUβi Uαj Uβjq sin2 p∆m2ij

iąj

ˆ
`2

ÿ

˚
˚
ImpUαiUβi Uαj Uβj
q sin

∆m2ij

iąj

L
qq
4E

(1.10)

˙Lq,
2E

where ∆m2ij “ m2i ´m2j , L is the path length of the neutrino, and the phase difference
between two mass states, φi ´ φj can be expressed as L

m2i ´m2j
2E

since the mass ofthe neutrino (ă 1 eV) is much smaller than the energies (ą 1 MeV) of current
experiments.
A simple case to consider is the oscillation between two neutrino flavors. In this
scenario, there is only one mass-squared difference, and the mixing matrix U reduces
12to a 2ˆ 2 unitary matrix with one degree of freedom.
ˆU“˙cos θ sin θ,´ sin θ cos θ(1.11)

where θ is the mixing angle. Equation 1.10 reduces to a simpler form#
2L
,α “ β1 ´ sin2 2θ sin2 ∆m4E
P pνα Ñ νβ q “2
2 ∆m2 Lsin 2θ sin 4E , α ‰ β,(1.12)

where ∆m2 “ m22 ´ m21 . Note that here and above, we have used natural units. If
converting from natural units to laboratory units, there is the change,∆m2 L
1.27∆m2 L ” eV 2 ¨ km ı
Ñ
.
4E
E
GeV
1.4.4

(1.13)

Neutrino Oscillations in MatterAll the calculations of neutrino oscillations so far have been assuming the neutrinos
propagate in vacuum. The oscillations arise from a difference in phase between themass eigenstates solely due to the mass differences. However, this has to be modified
in matter where the phase difference is determined by the total energy of the mass
eigenstate. Neutral-current interactions are flavor independent, and so produce acommon phase for all components of neutrinos. However, since the Earth contains
a large number of electrons, electron neutrinos may interact with them via the W ˘bosons, while the muon and the tau neutrino do not. Therefore, electron neutrinos
passing through matter are affected by coherent elastic forward scattering, as shownin Fig. 1.6, adding a potential to electron neutrinos. The resulting effective potentialthat electron neutrinos feel is expressed as?
Ve “ ˘ 2GF ne ,(1.14)

where GF is Fermi’s constant and ne is the electron density of the matter, and thecontribution is positive for νe and negative for ν̄e . The influence of matter on neutrino
13

oscillations is often referred to as the MSW effect (Wolfenstein, 1978; Mikheev andSmirnov, 1985).
νe

e´

W`
νe

e´
Figure 1.6: Neutrino CC forward elastic scattering in matter. This reaction adds
a potential to the electron neutrinos.

It is fairly easier to demonstrate the MSW effect in two-neutrino mixing and oscillations where one component is the electron neutrino. The neutrino is then assumed
to pass through a medium with constant density. In this model, the Schrödinger
equation is
d
i
dt

ˆ ˙
ˆ ˙
νe
νe
“H
,
νβ
νβ

(1.15)

where H is the Hamiltonian operator. In a vacuum, this operator isˆ
Hvac “

E1 0
0 E2

˜

˙
« pI `

m21
2E

0

¸
0
m22
2E

.

(1.16)We define the neutrino mixing by
ˆ ˙ ˆ
˙ˆ ˙
νe
cos θ sin θ
ν1
“
.
νβ
´ sin θ cos θν2(1.17)

Then we transform the Schrödinger equation to be expressed in the flavor basis as
d
i
dt

ˆ ˙
ˆ ˙
ˆ
˙ˆ ˙
∆m2
νe
sin2 θ
sin θ cos θνe: νe“ U Hvac U
“
.
2
νβ
νβ
sin
θ
cos
θ
cos
θ
ν
2E
β

14

(1.18)

Note that we have subtracted a component proportional to the identity matrix thatwill not affect the oscillation. Then we have to add an effective potential to the
Schrödinger equation that only appears in the e-e component of the weak basis
ˆ
Hm “

˙
Ve 0
.
0 0

(1.19)The total effective Hamiltonian in the weak basis isˆ∆m2
H “ U Hvac U ` Hm “
2E
:

˙
2E
sin2 θ ` ∆msin θ cos θ
2 Ve.
sin θ cos θ
cos2 θ

(1.20)

Note that the mass matrix is no longer diagonal, meaning that the mass eigenstatesin matter are not identical to the mass eigenstates in the vacuum. In order to get
the correct eigenstates, we need to diagonalize this matrix. The eigenvalues are2

˘∆m2M “ ˘ ∆mf where
4E m
c
fm “

sin2 2θ ` pcos 2θ ´

2E
Ve q2 .
∆m2

(1.21)

And the mixing matrix in matter, UM , isˆUM “

cos θM sin θM´ sin θM cos θM˙,(1.22)

where it turns out thatcos 2θM “1∆m2
pcos 2θ ´
Ve q
fm
2E

sin 2θM “1
sin 2θ.
fm(1.23)The two-neutrino oscillation probability in matter in terms of the variables θM and∆m2M isP pνe Ñ νβ q “ sin2 2θM sin2
15

∆m2M L
.
E(1.24)

Obviously, the vacuum oscillations are restored as the matter density gets smaller.
Additionally, there exists a resonance effect at certain energies and densities whensin2 2θM “ 1, which was first pointed out by Mikheev and Smirnov (Mikheev andSmirnov, 1985).
It should be noted that the matter effect can determine the sign of the mass splitting. In the mass splitting parameter, ∆m2M , there is a term p∆m2 ´ 2EVe { cos θq2 .
Depending on the sign of ∆m2 , the matter potential, Ve , is part of a term that
adds or subtracts from the vacuum mass splitting. In atmospheric neutrinos, thematter effect is a sub-leading factor because neutrinos coming from the bottom of
the detector travel a long distance through matter. Therefore, atmospheric neutrino
oscillations can in principle help to distinguish the sign of the mass splitting.
For three-flavor neutrino oscillations, it is not possible to calculate the matter
effect so elegantly. The change of density in the Earth adds more difficulties to thecalculation in the atmospheric neutrino oscillations. For the analysis described here,
we use a solution in (Barger et al., 1980) to account for the matter effect. For thedensity profile of the Earth, we numerically calculate the probabilities by iterating
the calculation through a series of constant density layers (Wendell et al., 2012).
1.4.5

Experimental Status

Since the discovery of neutrino oscillations, many experiments have been performed
to measure different flavor oscillations using diverse sources of neutrinos, such as
solar neutrinos (Ahmad et al., 2002; Fukuda et al., 2001), atmospheric neutrinos(Abe et al., 2011b), accelerator neutrinos (Ahn et al., 2003; Michael et al., 2006; Abeet al., 2012a) and reactor neutrinos (Abe et al., 2012b; An et al., 2012; Ahn et al.,2012). The experimental status of neutrino oscillations is summarized below.
The current values of sin2 θ12 and ∆m221 come from the measurement of solar
neutrinos and reactor neutrinos. Super-K and SNO measured the oscillations of16

solar neutrinos as described earlier in this chapter (Ahmad et al., 2002; Fukuda
et al., 2001). The oscillations were later confirmed by the KamLAND experiment
which measured ν̄e disappearance using ν̄e coming from reactors (Araki et al., 2005).
The neutrinos measured by KamLAND were from sources at various distances totheir detector; thus they were exposed to sources with a wide range of L{E. This
allowed the experiment to measure a precise sinusoidal L{E dependence of neutrino
oscillations (Araki et al., 2005). Fig. 1.7 shows the number of observed events as afunction of L/E.

Figure 1.7: Ratio of the observed ν̄e to the expectation for no oscillation as afunction of L/E. Taken from (Araki et al., 2005).

The combination of the these three experiments, Super-K, SNO, and KamLAND
put the tightest constraints on the oscillation parameters, sin2 θ12 and ∆m221 . Fig.
1.8 shows the constraints on these parameters from the three experiments. In thisanalysis, the values for the solar parameters come from a global fit in the 2014 PDG(Olive et al., 2014).
As discussed earlier in this chapter, Super-K played an important role in the
discovery of neutrino oscillations. Super-K has been measuring the oscillations ofatmospheric neutrinos since its start, and provides one of the best measurements of17

Figure 1.8: Allowed contours of sin2 θ12 and ∆m221 for Super-K solar neutrino
data (green solid), Super-K+SNO (green dashed), KamLAND (blue), and combined
results (red). Taken from (Sekiya, 2013).

sin2 θ23 and ∆m232 (Abe et al., 2011b). IceCube-DeepCore is another experiment thatmeasured these parameters by measuring the oscillations of atmospheric neutrino
with energies higher than 10 GeV (Aartsen et al., 2016). Recently, the long-baseline
experiments, such as MINOS (Adamson et al., 2011), T2K (Abe et al., 2017) andNOvA (Adamson et al., 2016), have provided the most precise measurements ofthese parameters. These experiments measured oscillations of neutrinos produced
at accelerators. Fig. 1.9 shows the constraints on sin2 θ23 and ∆m232 from both
atmospheric and accelerator neutrino experiments. Note that the figure below is theresult of a very recent analysis. In the oscillation described in this thesis (Abe et al.,2017), the value for the two parameters are from an older Super-K result in (Abeet al., 2011b).
The mixing angle of θ13 was initially expected to be small since the solar andatmospheric neutrino results were consistent with the value of zero. In 2011, T2K
first reported the indication of non-zero sin2 2θ23 via measuring electron neutrino
appearance from muon neutrinos (Abe et al., 2011a). Shortly after the indication
18

Figure 1.9: Allowed contours of sin2 θ23 and |∆m232 | assuming the normal hierarchfrom T2K, NOvA, MINOS+, Super-K, and IceCube. Taken from (Abe et al., 2017).

from T2K, several reactor experiments, such as Daya Bay (An et al., 2012), RENO
(Ahn et al., 2012) and Double-Chooz (Abe et al., 2012b), measured the value ofsin2 2θ13 precisely. These experiments measured the disappearance of ν̄e with one (ora set of detectors) near the reactor about a few hundred meters away and another(or a set of detectors) a couple of kilometers away. The effect of neutrino oscillations
can be measured by comparing the observed neutrinos in the near and far detectors.
Fig. 1.10 shows the layout of the Daya Bay experiment, and the measured sin2 2θ13
via observing ν̄e in their near and far detectors.

Figure 1.10: Left: Layout of the Daya Bay experiment. Right: Ratio of measured
versus expected signal in each detector, assuming no oscillation. The oscillation
survival probability at the best-fit value is given by the smooth curve. Taken from(An et al., 2012).

19

Recently, all these three experiments have updated their results with more data,
making a more precise measurement of sin2 θ13 as recorded in (Patrignani et al.,2016). In the analysis described in this thesis, the value for sin2 θ13 comes from thean older global fit which are recorded in the 2012 PDG (Beringer et al., 2012).
The current best measured values of the oscillation parameters are summarizedin Table 1.1.
Table 1.1: Oscillation parameter values as listed in (Patrignani et al., 2016). For∆m223 , the sign is unknown and the uncertainty is on the absolute value.
parametervalue
2
sin θ120.304˘0.014
2
sin θ130.0219˘0.0012
sin2 θ23
0.51˘0.05
2∆m21(7.53˘0.18)ˆ10´5 eV2∆m232
0.00244˘0.00004 eV2

1.4.6

Remaining Questions

Although most of the oscillation parameters have been measured, there are still
important open questions to be answered in order to fully understand the nature ofneutrinos.
First of all, it is not known whether θ23 mixing is exactly maximal (θ23 =45˝ ).
Current measurements of sin2 θ23 are consistent with 0.5, but have a large uncertainty.
If θ23 mixing is maximal, it could be a hint of some underlying symmetry (Xing andZhao, 2016).
While the squared mass differences between the neutrino mass eigenstates havebeen measured, the sign of ∆m232 (or ∆m231 ) has not yet been determined. In thenormal hierarchy: m3 ą m2 ą m1 , while in the inverted hierarchy: m2 ą m1 ą m3 ,as shown in Fig. 1.11. Due to this ambiguity, current experiments have to analyze
their results under different assumptions of mass hierarchy. The true hierarchy also
20has an important implication in the neutrino-less double beta decay experiments thatsearch for a process involves simultaneous emission of two positrons but no neutrinos.
If the mass hierarchy is inverted, the next generation of neutrino-less double beta
decay experiments have a chance to detect a signal (Beringer et al., 2012).

Figure 1.11: The ordering of the neutrino mass eigenstates. For the normal hierarchy, m3 is the heaviest state and ∆m232 ą0, while for the inverted hierarchy m3 isthe lightest and ∆m232 ă0. Taken from (King and Luhn, 2013).

Another important question is whether CP is violated in neutrino oscillations.
The value of δCP must be non-zero in order for CP violation in neutrino oscillations
to happen. It can be probed by comparing the probabilities of the CP-conjugate
oscillations, Pνµ Ñνe and Pν̄µ Ñν̄e . If CP turns out to be violated in leptons, it is ofimportance for the understanding of the matter and anti-matter asymmetry in theuniverse (Dimopoulos and Susskind, 1978).

1.5 Outline of Thesis
This dissertation is organized as follows: a brief introduction of neutrino history
and neutrino physics is presented in Chapter 1. Further discussion of tau neutrino
21and the study of tau neutrino appearance are presented in Chapter 2. The experimental setup of Super-Kamiokande is described in Chapter 3. The simulation ofneutrino events is described in Chapter 4. The event selection and reconstruction
are presented in Chapter 5. A neural network algorithm for tau event identification
is described in Chapter 6. The atmospheric tau neutrino appearance search is presented in Chapter 7. A measurement of charged-current tau neutrino cross section
is presented in Chapter 8. Future application of the tau identification in three-flavor
neutrino oscillation is discussed in Chapter 9. Finally, in Chapter 10, I summarize
the tau neutrino analysis.
The work described in this dissertation is made possible with efforts from many
Super-Kamiokande collaborators. The works done by Super-K collaborators arepresented here to have a full description of the analysis. My contributions to theanalysis are listed here: I was one of the data reduction experts in charge of thepartially-contained reduction. I used the neural network to select tau events indatasets from SK-I through SK-IV, and developed a new method of search for tau
neutrino appearance. I selected and optimized the input variables to the neutral
network in the analysis. I carried out the analysis of the measurement of the chargedcurrent tau neutrino cross section. I was a member of oscillation analysis group thatdeveloped the three-flavor oscillation analysis with tau neutrino identification.

22

2
The Tau Neutrino

The tau neutrino, together with the tau, forms the third generation of leptons. It
was the last of the leptons in the Standard Model to be discovered. The tau neutrino
is responsible for the deficit of muon neutrinos observed in Super-K because the tau
neutrinos from muon neutrino flavor change are rarely observed in Super-K. In thischapter, I briefly discuss the status of study of the tau neutrino to better set thecontext for the analysis in this thesis.

2.1 Direct Observation of the Tau Neutrino
The tau neutrino was first directly observed in DONUT experiment in 2000, as mentioned in Section 1.1. DONUT was designed to search for charged-current (CC) tau
neutrino interactions by identifying tau as the only lepton produced at the interaction vertex. Fig. 2.1 shows the experimental setup of the DONUT experiment.

The neutrinos were created by directing 800 GeV/c protons from the Fermilab
Tevatron to a 1 meter tungsten target (Kodama et al., 2001). The dominant component in the neutrino beam is electron neutrino and muon neutrino from the decays
23

Figure 2.1: Schematic view the the experimental setup of DONUT experiment.
Taken from (Kodama et al., 2001).

of charmed mesons in the beam dump. The tau neutrinos were primarily producedby the leptonic decays of a Ds meson into τ and ν̄τ , and subsequent decay of τ to ντ .
The muons produced in the proton interactions in the beam dump were absorbed or
swept away from the downstream detector using shielding and magnets.
The DONUT detectors consist of a scintillation counter veto wall, the hybrid
emulsion target, trigger hodoscopes, and the charged particle spectrometer. The
hybrid emulsion target was also used as a tracker which has a spatial resolutionof „1 µm (Kodama et al., 2001). The good spatial resolution is important in theidentification of tau at the interaction vertex because of its extremely short lifetime.
In total, DONUT observed 9 CC ντ events with an estimated background of 1.5
events based on 3.6ˆ1017 protons on target (Kodama et al., 2008). Fig. 2.2 shows
four CC ντ candidates in DONUT.
DONUT also measured the CC tau neutrino cross section based on the direct
observation of CC ντ interactions. The total CC ντ cross section per nucleon isfactored as
σpEq “ σ const EKpEq,(2.1)

where KpEq is the kinematic effect due to the tau lepton mass as shown in Fig. 2.3.
The kinematics of ντ n Ñ τ ´ p (or ν̄τ p Ñ τ ` n) requires that the minimum neutrino
24

Figure 2.2: Four examples of ντ CC interaction events in DONUT. The neutrinos
are incident from the left. Taken from (Kodama et al., 2001).

energy in the rest frame of the target nucleon isEmin “ mτ r1 `mτs “ 3.5GeV.
2mN

(2.2)

The observed number of neutrino interactions can be written as
Nint “ T OT Nνtgt

Nnucl
σ,
area

(2.3)

where T OT is the experimental efficiency for observing the interaction, Nνtgt is thenumber of neutrinos passing through the target area, Nnucl is the number of nucleonsin the target area, and σ is the interaction cross section per nucleon.
Because DONUT can not distinguish ντ and ν̄τ , the measured CC cross section is
actually an average of ντ and ν̄τ cross sections. The flux of neutrinos was measured
to be approximately equal to the flux of antineutrinos (Kodama et al., 2008). σ const
is measured to be (0.39˘0.13˘0.13)ˆ10´38 cm2 GeV´1 . The first error is statistical,
the second one is systematic. The measured cross section has a large statistical error
25

Figure 2.3: The suppression of CC tau neutrino cross section compared to CCνµ /ν̄µ cross sections due to the tau lepton mass. Taken from (Kodama et al., 2008).

because DONUT only observed 9 neutrino interactions. The systematic error is also
large because of the large uncertainty on the production and decay of Ds mesons.

2.2 Previous Search for Tau Neutrinos from Neutrino Oscillations
In three-flavor neutrino oscillations, νµ Ñ ντ oscillation is the reason for the observed
deficit of muons neutrinos in Super-K atmospheric neutrinos and long-baseline accelerator neutrinos. Because charged-current tau neutrino cross sections are suppressedby the kinematic effect of relatively large mass of the tau, tau neutrinos from neutrino
oscillations interact more rarely than electron neutrinos or muon neutrinos in the detectors. The 3.5 GeV energy threshold is higher than the peak energy of neutrinos in
long-baseline experiments that are tuned to maximal disappearance at their baseline.
This results in the deficit of observed neutrinos in Super-K and long-baseline neutrino
experiments. Additionally, the direct observation of tau neutrinos from neutrino oscillations has been challenging due to the difficulties of detecting and identifying
26

tau leptons. Nonetheless, Super-K and OPERA experiments have searched for tau
neutrino appearance in atmospheric neutrinos and high-energy accelerator neutrinos
respectively.
2.2.1

Previous Search for Tau Neutrino Appearance in Super-K

Atmospheric neutrinos have a wide range of energies as described in Section 1.3. Thismeans that the higher energy tau neutrinos from neutrino oscillations have enoughenergy for charged-current interactions to occur. Super-K published a result of asearch for atmospheric tau neutrino appearance based on SK-I data in 2006 (Abeet al., 2006), then updated the result based on SK-I-III data in 2013 (Abe et al.,2013).
Instead of exclusively reconstructing the track of tau at the neutrino interaction
vertex, Super-K performed a search for hadronic decays of tau leptons. The decays
typically produce multiple pions, which are not often seen in electron or muon neutrino events. The separation of tau events from electron or muon neutrino eventswas based on a likelihood method or a neural network algorithm. Fig. 2.4 shows
the likelihood and neural network output distributions for data taken in SK-I. The
events with likelihood ą 0 or neural network output ą0.5 are defined as tau-like.
The selected tau-like events were fitted with a combination of the expected tau
appearance and atmospheric neutrino background in the zenith angle distribution.
Ntotal pcos θq “ αNtau ` βNbkg ,
2

χ “

10
ÿ
pN obs ´ αN tau ´ βN bkg q2
i

i

σi2

i“1

i

,

(2.4)
(2.5)

The best-fit value of tau normalization is 1.76 (1.71) for the likelihood method (neuralnetwork method), corresponding to 138˘48(stat)`15´32 (syst) observed tau events. The
hypothesis of no-tau-appearance was excluded with a significance of 2.4σ. Fig. 2.5
shows the zenith angle distributions for tau-like events after fit.
27

Figure 2.4: The likelihood (top) and neural network output (bottom) distributionsof downward (left) and upward (right) events for data (dots), atmospheric neutrino
background MC (dashed histogram), and tau and background MC (solid histogram).
The events for likelihood ą 0 and neural network ą0.5 are defined as tau-like. Taken
from (Abe et al., 2006).

In 2013, Super-K updated the result using improved neural network for tau separation and more data collected in SK-I through SK-III. Additionally, the data
were fitted using an unbinned likelihood fit with a combination of the expected
tau events and atmospheric neutrino background in two-dimensional distributionsof neural network output and zenith angle. The tau normalization was fitted to`17.8
be 1.42˘0.35(stat)`0.14´0.12 , corresponding to 180.1˘44.3(stat)´15.2 (syst) observed tau

events. The significance of excluding no-tau-appearance was increased to 3.8σ.
2.2.2

OPERA

The OPERA experiment in the underground Gran Sasso Laboratory (LNGS) was
designed to detect tau neutrino appearance using accelerator neutrinos. It used a28

Figure 2.5: Zenith angle distribution for the likelihood (top) and neural network(bottom) methods. A fitted excess of tau-like events in the upward-going direction
is shown in the shaded area. Taken from (Abe et al., 2006).

Figure 2.6: Fit results showing projections in the NN output and zenith angle
distribution. The simulations and data sets have been combined from SK-I through
SK-III. The fitted tau signal is shown in gray. Taken from (Abe et al., 2013).

29

emulsion target similar to that in the DONUT experiment to achieve precise reconstruction of the short track of tau. The detector was placed 730 km downstream
the high energy long-baseline CERN to LNGS neutrino beam. Instead of using aneutrino beam containing tau neutrinos at production in DONUT, OPERA used apure muon neutrino beam that contained almost no tau neutrinos at production. ντappearance is expected after the νµ neutrinos travel 730 km to reach the detector.
Fig. 2.7 shows a candidate ντ event observed in OPERA. At the end of the OPERA

Figure 2.7: Event display of a candidate ντ events in OPERA, in which tau decay
quickly into secondary particles after production. Taken from (Agafonova et al.,2015).

experiment, 5 ντ events were detected with an estimated background of 0.25˘0.05,
rejecting the background-only hypothesis at 5.1σ (Agafonova et al., 2015).
The studies of the tau neutrino and related neutrino oscillations have produced
interesting results in the last two decades. However, the studies are still limitedby the statistics because of the difficulties in both production and detection of tau
neutrinos. The large target mass of Super-Kamiokande, coupling with the wide
energy span of atmospheric neutrinos, offers a good opportunity to study the tau
neutrino and related neutrino oscillations.

30

3Super-Kamiokande

Super-Kamiokande (Super-K) is a large water Cherenkov detector in Japan. It is
located in the Mozumi mine at a depth of 1,000 underground in Gifu, Japan. The
main purpose of the experiment is to search for proton decays, as well as to study
neutrinos from various natural sources, such as neutrinos from the Sun, atmosphere,supernovae, and other astrophysical sources. It was also used as the far detector in
K2K experiment (Ahn et al., 2003) between 1999 and 2004. It has been playing thesame role in T2K long-baseline neutrino experiment (Abe et al., 2011a) since 2009.
Super-K started taking data in 1996, and has been running for more than 20 years.
Datasets have been taken in four run periods defined by the detector configurations,
SK-I (1996–2001), SK-II (2002–2005), SK-III (2006–2008) and SK–IV (2008-now).
The configurations of each phase and detection methods are described here. Moredetailed information can be found in (Fukuda et al., 2003) and (Abe et al., 2014).

31

3.1 Cherenkov Radiation
Cherenkov radiation is electromagnetic radiation produced by a charged particle
passing through a dielectric medium with a speed faster than the speed of light inthat medium. The photons are emitted in a cone along the direction of the particle
with an angle ascos θc “

1
,
nβ

(3.1)

where θc is the angle between the direction of Cherenkov light and the direction of
the charged particle, n is the refractive index of the medium and β is the ratio ofthe velocity of the particle to the speed of light in the vacuum. A common analogy
is the sonic boom of a supersonic object. The sound waves generated by the object
cannot propagate forward from the object, instead forming a shock front. Fig. 3.1
illustrates the production of Cherenkov light and the direction of the light relative tothe direction of charged particle. For the ultra-pure water of Super-K, n=1.33 at 580
nm, so that 1/n«0.75. In water, the energy threshold for Cherenkov radiation is 0.76
MeV for electrons and 158.7 MeV for muons. The number of photons generated perunit track length and per unit wavelength is given by the expression (Mead, 1958),
d2 N
2πα
1
“ 2 p1 ´ 2 2 q
dxdλ
λ
nβ

(3.2)

where α=1/137 is the fine structure constant, x is the track length, and λ is the
wavelength of the photon. It corresponds to about 400 photons per cm at λ=300550 nm.

3.2 Super-Kamiokande Detector
The Super-K detector consists of a large cylindrical water tank filled with 50,000
ton of ultra pure water. The tank has a diameter of 39.3 m and a height of 41.4 m.
32

c
t
n
θ

βct

Figure 3.1: Schematic of the Cherenkov wavefront.

The detector is optically separated into two concentric cylindrical regions, the innerdetector (ID) and the outer detector (OD). The OD is defined as the region extending
2.05 m inward the wall of the barrel, and 2.2 m from its top and bottom. A 55-cm
thick stainless steel support structure, which is opaque to light, is installed between
the ID and the OD. Photo-multipliers are mounted on both sides of the stainless steel
structures to detect Cherenkov radiation produced by charged particles in water. The
layout of the detector is shown in Fig. 3.2.
3.2.1

Inner Detector (ID)The ID has a diameter of 33.8 m and a height of 36.2 m, and contains 32,000 tonsof ultra-pure water. 20 inch inward-facing PMTs are installed uniformly on the
inner surface of the ID. To prevent light leaks into the OD, the space between the
PMTs is lined with black sheet. During the SK-I period, there were 11,146 PMTs,
corresponding to about 40% of the ID surface covered by photocathode. An accident
in 2001 caused by the implosion of PMTs damaged half of the PMTs, the remaining
PMTs were rearranged uniformly in the reconstruction. For SK-II, the detector was
run with 5,182 PMTs in the ID, corresponding to about 19% photocathode coverage.
33

Figure 3.2: A sketch of Super-K detector, and the location of the detector under
Mt. Ikenoyama.

A full reconstruction in 2006 restored the detector with full photocathode coverageof 40% in the ID, with 11129 PMTs in the ID after the restoration. The period
between 2006 and 2008 is known as SK-III. An upgrade of electronics in 2008 marks
the start of the SK-IV period.
3.2.2

Outer Detector (OD)The OD is viewed by 8-inch outward-facing PMTs installed on the OD side of thestainless steel support structure. In all four SK periods, there are 1885 PMTs in
the OD. In order to improve the light detection efficiency, the inner surface of thetank and the space between PMTs on the support structure is covered by reflective
Tyvek R . Tyvek’s reflectivity is „90% at the wavelength of 400 nm.
During the full reconstruction of the detector in 2005, Tyvek sheet was installedin the OD to separate the wall barrel region from the top region and the bottom
region. Fig. 3.3 shows a schematic view of the OD segmentations.

34

OD

ID
OD segmentation

Figure 3.3: Schematic view of the Tyvek sheet installed in the OD that divides the
OD into top, barrel and bottom regions.

3.3 Photomultiplier Tubes (PMT)A photomultiplier tube (PMT) is a kind of sensitive detector of ultraviolet, visible
and near-infrared ranges of light. An incident photon produces a photo-electron on
the photocathode, and the signal current is amplified by multiple dynodes. The
PMT enables the detection of single photons when the incident flux of light is very
low. It has the advantages of high gain, low noise, ultra-fast response and a large
area of collection.
The ID uses a 20-inch diameter Hamamatsu R3600 PMT (Fig. 3.4), developed byHamamatsu Photonics K.K. in cooperation with Kamiokande collaborators (Kumeet al., 1983). The dynode structure and bleeder circuit were modified to improvetiming response and collection efficiency for Super-K (Suzuki et al., 1993). The
photocathode of the PMT is coated with bi-alkali (Sb-K-Cs) that has high sensitivity
to Cherenkov light and low thermal emission. A peak quantum efficiency (22%) isfound in the region of 360–400 nm wavelength, that also corresponds to the peak of
the relative Cherenkov spectrum in pure water as shown in Fig. 3.5. The gain of the35

ID PMT is 107 at a high voltage of 1500 V to 2000 V. Fig. 3.6 shows the pulse height
distribution for a single p.e in the ID PMT. Additional specifications are listed inTable 3.1.

Figure 3.4: A schematic of the 20-inch PMT used in the ID. Taken from (Fukuda
et al., 2003).

Figure 3.5: Quantum efficiency of the ID PMT, with relative spectrum of
Cherenkov light in pure water overlaid. Taken from (Nishino, 2009).

The geomagnetic field at the Super-Kamiokande detector is about 450 mG, which
affects photoelectron trajectories inside the PMTs. In order to compensate the field,
36

distribution (arbitrary unit)

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

0

1

2
3
photo-electron

4

5

Figure 3.6: Pulse height distribution for a single p.e. in a 20-in PMT. The peaknear zero is from dark current (Fukuda et al., 2003).

Table 3.1: Specification of the 20-inch PMT used in Super-Kamiokande.
Shape
Photocathode area
Window material
Quantum efficiency
DynodesGain
Dark current
Dark pulse rate
Cathode non-uniformity
Anode non-uniformity
Transit time
Transit time spread
WeightPressure tolerance

Hemispherical
50 cm diamater
Bialkali (Sb-K-Cs)
22% at λ = 390 nm
11 stage Venetian blind type107 at „ 2000 V
200 nA at 107 gain
3 kHz at 107 gain
ă 10%
ă 40%
90 nsec at 107 gain
2.2 nsec (1 σ) for 1 p.e. equivalent signals
13 kg
6 kg/cm2 water proof

26 sets of horizontal and vertical Helmholtz coils are arranged around the innersurface of tank, reducing the magnetic field to « 50 mG.
Following the accident in November 2001, the ID PMTs were instrumented withacrylic covers. A clear 12-mm thick UV-transparent acrylic dome is installed overthe photocathode of each PMT, and the side of each PMT is protected with FiberReinforced Plastic (FPR), with holes allowing water to flow freely. The transparencyof the cover for photons with normal incident angle in water is ą 96% above 350 nm
37

wavelength.
The OD uses 1885 8-inch PMTs. The photocathode of each PMT is coveredby a 60 cmˆ60 cmˆ1.3 cm wavelength shifter plater (Becker-Szendy et al., 1993),
which increases the light collection efficiency by 60%. The plates degrade the timing
resolution of the OD PMTs from 13 ns to 15 ns. Since the OD is used as a vetofor cosmic ray muons rather than a particle tracker, the light collection efficiency ismore important than the timing resolution.

3.4 Electronics and Data Acquisition (DAQ)The Super-Kamiokande electronics and data acquisition systems (DAQ) were designed to process signals from a large number of channels quickly and efficiently.
The electronics and DAQ were upgraded in 2008, which marked the start of SK-IV.
In this chapter, I briefly describe the features of the Super-K electronics and DAQ
system.
3.4.1

Inner Detector Electronics and DAQ (SK I-III )For the first three running periods, each of the ID PMT was connected to one of 12
input on custom analog timing modules (ATMs) for signal processing. The ATMs
recorded integrated charge and arrival time of all ID PMTs. Each ATM was housed
in a crate containing 20 ATMs, with a total of 48 crates evenly distributed amongfour electronics huts located above the detector tank.
Each input channel of the ATM has two switching channels containing both aCharge to Analog Converter (QAC) and a Time to Analog Converter (TAC), allowing
one channel to be active while the other one is digitizing. This two channel setup
minimizes dead time, in particular for events in which there may be two stages
occurring in rapid succession. Detailed specifications of the ATM are shown in Table
3.2. Each ATM is set at a threshold of 0.23 p.e. equivalent. When a PMT is hit and38

Table 3.2: Specifications of the ATM.
Number of channels
One hit processing time
Charge dynamic range
Timing dynamic range
Charge resolution (LSB)
Charge resolution (RMS)
Timing resolution (LSB)
Timing resolution (LSB)
Temperature dependence (QAC)
Temperature dependence (TAC)
Event number
Data size of one hit
FIFO

12/board„5.5 µs„400-600 pC (12-bit)„1300 ns (12-bit)
0.2 pC
0.2 pC„0.3-0.4 ns
0.4 ns3 Count/deg. Ø0.6 pC/deg.
2 Count/deg. Ø0.8 pC/deg.
8-bit
6 bytes
2 kBytes ( 340 hits)

triggered, a 15 mV 200 ns wide square wave is issued to the ATM front panel and a
900 ns veto is placed on the channel to prevent after-pulsing hits. An analog sum of
all such waves is known as the “HITSUM”. A global trigger is issued at 29 PMT hits,corresponding to roughly 5 MeV of energy. On a global trigger, a stop signal is sent
to TAC, then the information in the QAC and TAC is digitized by the ADC. The
integrated charge and time information from all ATM channels is stored along withthe event number information distributed by GONG (Go/No-Go) modules. Other
header information is sent to Super Memory Partner boards (SMP). Data stored in
the SMPs are later read out by online computers via Versa Module Europa (VME).
More details of the DAQ process can be found in (Fukuda et al., 2003).
3.4.2

Outer Detector Electronics and DAQ (SK I-III)

The paddle cards are used to distribute high-voltage to the OD PMTs. Each paddle
card is connected to 12 PMTs via coaxial cables and is itself powered by a LeCroy
1454 mainframe. The coaxial cable is used to supply high-voltage and transmit
analog signals from the PMTs. The signals are then sent to charge-to-time conversion(QTC) modules that convert the PMT signal to a rectangular pulse with a width
39

proportional to the input charge. This pulse is generated along with a 25 mV high200 ns wide analog signal used to create an OD “hitsum” when a PMT exceeds a
0.25 p.e. threshold in a 200 ns wide time window. When a global trigger is issued,
the logical pulse from the QTC is fed to a time-to-digital converter (TDC) which isread out across all OD channels. OD data are stored 10 µs before and 6 µs after thetrigger and are read out via VME to online computers for analysis.
3.4.3

Triggering (SK I-III)

Three types of Super-K triggers are implemented based on the number of hit PMTs
(HITSUM) in the ID within a 200 ns time window, corresponding roughly to separate
physics analyses of interest. The high energy (HE) trigger is issued when HITSUM
exceeds 31 hits and fires at a rate of „5 Hz. The low energy (LE) trigger corresponds
to 29 hits and has a rate of „11 Hz. Finally, the super low energy (SLE) trigger
was implemented about a year into SK-I period to lower the solar neutrino analysis
threshold. The OD trigger system works similarly, with a trigger threshold of 19
hits within a 200 ns time window. When a trigger is issued, it is sent to a hardware
trigger module TRG that signals the readout of the PMT channels.
3.4.4

SK IV Electronics and DAQ Upgrade

In order to reduce dead time, the Super-K front-end electronics and DAQ systemwere upgraded in 2008 (Yamada et al., 2010; Nishino et al., 2009). Here, I briefly
describe the upgrade.
The new front-end electronics are QTC-Based Electronics with Ethernet (QBEE).
One QBEE board is equipped with eight QTC chips, and each QTC chip has three
input channels. The QTC has three gain ranges (small, medium, large) for PMT
readout that covers a wider dynamic range (0.2 - 2500 pC) than the ATMs. The
timing and charge resolution are similar to those of the ATMs. The QTC outputs

40

a pulse whose leading edge contains the timing information and whose width is
proportional to the total charge. The ethernet cables improves data flow relative tothe old TKO readout.
Instead of hardware triggers in SK-I-III, a software trigger is implemented in
SK-IV that contains all PMT hit information. The software trigger makes selections
based on the same hit threshold described above for SK-I-III trigger system.

3.5 Background Reduction
3.5.1

Water Purification

Since the light collected by the PMTs has to travel through the water, it is importantto keep the water as pure as possible. Water taken from the mine is first passed
through a 1 µm filter to remove large particles. It is then passed through a heatexchanger that cools the water to about 13 ˝ C to limit the growth of bacteria.
Following this, a cartridge polisher removes heavy ions, and a UV sterilizer kills theremaining bacteria. The water then goes through a reverse osmosis and vacuum
degassifier system that removes radon and other dissolved gas. In the final stage
ultra filters remove particles larger than „10 nm in diameter. The water is thenreturned to the tank. Water in the tank circulates through the system at about 60
tons/hour.
3.5.2

Air Purification System

The air in the mine housing Super-K detector has a high concentration (10–5000
Bq/m3 ) of radon. To limit the amount of this background entering Super-K, fresh
air is continuously pumped into the control room and cavern area of the experiment.
Radon concentrations are kept around 20-30 Bq/m3 inside the top of the Super-K
tank and below 100 Bq/m3 at the cavern.

41

3.6 Calibration
3.6.1

PMT Calibration

Absolute gain
For the absolute gain, the single p.e. distribution must be measured in order toconvert the charge measured by the PMTs to the number of p.e.s. The calibration
uses a 252 Cf source surrounded by a spherical mixture of nickel oxide and polyethylene. Spontaneous fission of the 252 Cf produces neutrons that are brought to thermal
equilibrium by the polyethylene, and then captured by the nickel. The neutron
captured by nucleus isotropically emits „9 MeV gamma rays. The nickel source isdeployed at the center of the tank, and generates about 100 photons per capture.
Therefore, most hit PMTs will have only a single p.e. The mean value of the 1 p.e.
charge distribution is taken as a constant to convert charge measured by PMTs to
the number of p.e.s, as shown in Fig. 3.7. The constants in SK-I through SK-IV are
2.044 pC/p.e., 2.297 pC/p.e., 2.243 pC/p.e., and 2.658 pC/p.e. respectively.
Relative PMT gain
The high voltage of each PMT is set so that all PMT gains are approximately
the same. The relative gain describes the response of each PMT relative to the average. Fig. 3.8 shows a schematic view of the relative gain calibration system. Light
from a Xe lamp is passed through an ultraviolet filter and injected into a scintillatorball lowered into the tank. The scintillator ball absorbs the remaining UV light andre-emits photons with wavelength of 440 nm uniformly. The intensity of the Xe lightis monitored by two photodiodes and a 2-inch PMT that also serves as a calibration
trigger. PMT high voltages are adjusted so that the charge measured in each tube
after corrections for light attenuation, acceptance and shape of the scintillation ballagree between all of the PMTs. This measurement is performed at various positionsof the ball within the tank. The relative gains distribution had width of 7.0% at the
42

mean

2.055pC

5

10

3

2

1

0

0

pC

Figure 3.7: Single photoelectron distribution of a typical ID PMT. The bump at
about 2.0 pC corresponds to on p.e.. Taken from (Fukuda et al., 2003)

beginning of SK-I.

Relative timing calibration
The hit timing of PMTs varies across the detector due to several factors. Since
the Super-K tank is very large, the timing difference arises from different PMT cable lengths. Similarly, the timing is a function of the amount of charge collectedby a PMT because larger hits exceed TDC threshold earlier due to slewing in thediscriminators. Therefore, each PMT’s timing response is measured as a function of
charge and stored as a look-up table called the time-charge (TQ) map. TQ maps
are measured with the apparatus in Fig. 3.9. Light from a N2 laser is split into two
beams, one sent to a diffuser ball deployed in the tank and the other to a reference
PMT. The laser produces a 3 ns wide pulse of 384 nm light. The intensity of thelight is varied with an optical fiber in order to measure the timing at different pulse

43Figure 3.8: Schematic view of the relative gain measurement system. Taken from
(Lee, 2012)heights. A TiO2 tip at the center of the diffuser ball works in conjunction with asilica gel to diffuse light without smearing its time profile. A typical TQ map isshown in Fig 3.10.

QBEE charge linearity
For the QBEE modules in SK-IV, a dedicated calibration run is taken to measure the linearity of the readout system. The measurement covers a wide rangeof light intensity, all the way into saturation. Ten ID PMTs, that are read out
with well-understood, high-resolution CAMAC ADC, are selected as reference in themeasurement. The laser system mentioned in the relative timing calibration is used
to inject laser of thirty different intensities into the tank. The relative intensities ofthese 30 sets of different-intensity laser shots were established by the reference PMTs.
The intensity scale for each ID PMT is fixed by the chosen intensity for which its

440.CUGTPORWNUGYKFVJPUGE&[G.CUGTPORWNUGYKFVJPUGE64)/QP2/6

5-+&

8CTKCDNGHKNVGT&KHHWUGTDCNN
Figure 3.9: Schematic view of the relative timing measurement using a laser. Taken
from (Ishihara, 2010).

Figure 3.10: A characteristic TQ map for an ID PMT (Abe et al., 2014).

45

observed charge is closest to 30 p.e.s. Using the relative intensities derived from theten reference PMTs, we can calculate the expected charge for each of the other 29
laser intensities. The charge linearity is measured as the ratio of observed charge
to expected charge. Fig. 3.11 shows the result of the linearity measurement. The
response of the ID PMTs is linear within 1% with less than 200 p.e.s, and increases
to more than 10% at above 1000 p.e.s. This non-linearity is taken into account in
the MC simulations.

Figure 3.11: Charge linearity of all ID PMTs. The points shows the average foreach expected number of photoelectrons region. Taken from (Abe et al., 2014).

3.6.2

Water Calibration

Measurement of light scattering and absorption in water with a laser
The attenuation of light during travel in the water depends on the water transparency. Water transparency is affected by absorption and scattering of light. The
relevant scattering processes in SK are Rayleigh and Mie scattering. The attenuation
length in water can be described by three parameters:
Latten pλq “1αRay pλq ` αM ie pλq ` αabs pλq

.

(3.3)

These three parameters are measured with a laser system that injects light tuned todifferent wavelengths via an optical fiber into the tank from eight different positions.
46

The charge and timing information from the laser light measurement was used totune the water parameters in Monte Carlo simulations (MC) until good agreement
between data and MC was found. The measured parameters are shown in Fig. 3.12

Figure 3.12: Typical fitted water coefficient functions used in the SK-MC. The
points are the data obtained in April 2009. Taken from (Abe et al., 2014).

Water transparency measurement with cosmic ray muonsCosmic ray muons passing through the detector at a rate of a few Hz, making up
the bulk of high energy events. Muons below 100 GeV deposit energy via ionization
at a rate of about 2 MeV/cm. Therefore, they are a natural calibration source
and cover a much wider energy range than the laser sources described above. Only
vertical muons are selected in this measurement. The muon track is reconstructed byconnecting the entrance and exit points in the ID. Assuming the light is not scattered,
the number of photoelectrons detected in a cosmic ray muon event is expressed by:

Qplq “ Q0f pθq ´l{Latten
e
l
47

(3.4)

where Q0 is a constant, f pθq is the PMT angular acceptance with θ as the angle respect to the normal, l is the photon path length, and Latten is the attenuation length.
The water attenuation length is measured continuously with comic ray muons. Fig.
3.13 shows the measurements of Equation 3.4 and time variation of the attenuationlog(Ql/f( θ))length.

500

RUN 3106

P1
P2

5.987
-0.9487E-04

1000

2000

3000

0.4609E-03
0.4444E-06

400
300
200
0

4000

5000

Figure 3.13: Effective observed charge for through-going muons as a function ofphoton pathlengths (l), where Q is the detected charge and f pθq is the PMT acceptance (top).Time variation of water attenuation lengths for SK-I, SK-II, SK-III and
SK-IV (bottom). Taken from (Abe et al., 2014).

3.6.3

Energy Scale

Precise knowledge of the detector energy scale is important for energy reconstruction.
Four independent sources are used to calibrate the energy scale in momentum range48

from a few tens of MeV/c to about 10 GeV/c.
Track length of high energy stopping muonsMomentum of high energy muons is approximately proportional to their track
lengths. Muons entering from the top and stopping inside the detector are selected
for calibration. The track length is reconstructed by connecting the entering point
and the vertex position of the subsequent decay electron. Resolutions of stopping
muon entering position and decay electron vertices are better than 50 cm. The accuracy of the absolute energy scale is checked by comparing calibration data with
Monte Carlo simulations. The agreement between data and MC is within 0.7%,
1.1%, 2.0% and 2.2% for SK-I, SK-II, SK-III, and SK-IV respectively.

Low energy stopping muons
The Cherenkov angles of low energy muons are used for low momentum (ă 400
MeV/c) calibration. As shown in Equation 3.1, the angle of the emitted Cherenkov
photons relative to the charged particle is a function of the particle’s energy. By measuring the Cherenkov angle of stopping muons, their momentum can be determined
as:
m
.
p“ a
pn cos θq2 ´ 1(3.5)

where n is the refraction index of water, m is the mass of muon, and p is the muon
momentum. The Cherenkov angle reaches a limiting value for high momenta, so only
muons with momentum less than 400 MeV/c are used. Fig. 3.14 shows the averaged
Pp.e. /Pθ and the double ratio of MC to data as a function of Pθ . The agreement iswithin 0.7%, 1.3%, 2.1% and 2.1% for SK-I, SK-II, SK-III and SK-IV respectively.
Neutrino induced π 0
Neutral pions are produced by neutral-current interactions of atmospheric neutrinos in the detector, which then decay into two photons. The invariant mass of π 049

SK-IV

SK-IV

1.4
1.04
1.02

MC/DataP(pe)/P(θ)

1.2

1

0.98

0.8

Data
MC

0.6
200

1

250

300

350

400

0.96
200

450

250

300

350

400

450

Reconstructed Momentum in P(θ)

Reconstructed Momentum in P(θ)

Figure 3.14: The averaged Pp.e. /Pθ where Pθ is the momentum reconstructed with
θC for SK-IV data (solid) and MC (dashed) (left), and the double ratio between MC
and data (right). Taken from (Lee, 2012).

can be reconstructed by combining momentum of two photons. The distribution of
the invariant mass is centered around 135 MeV as shown in Fig. 3.15. A Gaussian
fit of the data and MC found that the agreement between data and MC is within
0.7%, 1.3%, 0.3% and 1.7% for SK-I, SK-II, SK-III and SK-IV respectively.
Super Kamiokande IV 2519.9 days : Monitoring

250

200

150

100

50

0
0

50

100

150

200

250

300

Mass(2 e-like rings) MeV/c2

Figure 3.15: Invariant mass distribution of neutrino-induced π 0 events for data(dot) and MC (boxes). MC is normalized to the live time of data.

Decay electrons50Decay electrons from decay of stopping muons are abundant and their energy
spectrum is useful for energy calibration in the tens of MeV region. Muons are
susceptible to muon capture by a nucleon, in which no electrons are produced. These
stopping muons are rejected in the search for decay electrons. Specific selection
criteria for decay electrons will be described in Section 5.2. Fig. 3.16 shows themomentum distribution of decay electrons in data and MC.
SK-IV

Number of events

4000

Mean Values
Data: 37.36 ± 0.05 MeV
MC: 37.11 ± 0.04 MeV

3000

2000

1000

0
0

20

40

60

80

Reconstructed momentum

Figure 3.16: Momentum of decay electrons for data (crosses) and MC (solid histogram). Taken from (Abe et al., 2014).

Summary of energy scaleFig. 3.17 shows the energy scale uncertainty as a function of momentum for thecalibration sources described above. The uncertainty is within 3% for all samples
considered.

51

Figure 3.17: Energy scale uncertainty as a function of momentum for calibration
measurements in SK-IV, decay electrons (solid circle), single π 0 (solid triangles),
low-energy stopping muon (open squares) and high-energy stopping muons (solid
squares). Taken from (Ishihara, 2010).

52

4
Atmospheric Neutrino Simulations

Monte Carlo (MC) simulation is essential for the Super-K atmospheric neutrino
analysis. The neutrino interactions in the detector are simulated, then comparedwith the collected data samples. The MC simulation of atmospheric neutrinos inSuper-K is done in stages - modeling the interactions of neutrinos in the water,subsequent particle tracking inside the detector volume and the detector’s responseto these particles. Each of these stages are described here.

4.1 Atmospheric Neutrino Flux
The analyses of atmospheric neutrinos are essentially done by investigating thechange of the observed neutrino flux at the detector relative to the expectation.
Therefore, it is important to have an accurate model of atmospheric neutrino production. The Honda flux calculation (Honda et al., 2015) is adopted to model theatmospheric neutrino flux at Super-K. It is fully three-dimensional and tracks every
necessary particle from its generation to interaction within the atmosphere. The
model takes into account the primary cosmic ray spectra, the hadronic interactions53

in the atmosphere and the effects of the geomagnetic field in the calculation. The
model is briefly described below.
The primary cosmic rays consist of mostly proton and helium nuclei. The primarycosmic ray flux, which is the input to the neutrino flux calculation, is based on
experimental measurements by AMS (Alcaraz et al., 2000) and BESS (Sanuki et al.,
2000; Haino et al., 2004). The effect of the magnetic field from the solar windis prominent at low cosmic ray momenta („ 1GeV/c), which is included in thecalculation. Fig. 4.1 shows the primary cosmic ray flux in the Honda model.
The US standard atmospheric model (NOAA and Force, 1976) is used as themodel for the density structure of the atmosphere. Geomagnetic field induces a updown and east-west asymmetry of the neutrino flux as well as a low-energy cutoff
in the flux due to the deflection of low energy primary cosmic rays. Its effects arecalculated based on the IGRF2005 model (Thébault et al., 2015).

Figure 4.1: The primary cosmic ray protons at solar minimum as calculated in(Honda et al., 2004). Overlaid data points (see (Honda et al., 2004) for references)
are taken from Webber (crosses), LEAP (upward triangles), MASS1 (open circles),
CAPRICE (vertical diamonds), IMAX (downward triangles), BESS98 (circles), AMS(squares), Ryan (horizontal diamonds), JACEE (downward open triangles), Ivanenko(upward open triangles), Kawamura (open squares) and Runjob (open diamonds).

54

For the production of secondary particles in the atmosphere, DPMJET-III (Roesler
et al., 2000) is used for hadronic interactions above 32 GeV while JAM (Niita et al.,2006) is used for energy below 32 GeV due to its better agreement with low energy
cosmic muon data. The hadronic interaction simulations are fine tuned to match thesimulated cosmic ray muons with experimental measurements from (Motoki et al.,
2003; Achard et al., 2004; Haino et al., 2004; Sanuki et al., 2002).
Fig. 4.2 shows the zenith angle distributions of neutrino flux for three different
neutrino energies simulated with Honda model. The distribution is peaked at thehorizontal because the cosmic ray particles can travel longer distances in the atmosphere horizontally and have more time for decay. The large up-down asymmetry in
the left plot is due to the geomagnetic field.

Figure 4.2: Zenith angle dependence of atmospheric neutrino flux at Super-K averaged over the azimuthal at 0.32 GeV (left), 1.0 GeV (middle) and 3.2 GeV (right).
cos θ=1 represents vertically downward going and cos θ=-1 represents vertically upward going neutrinos. Taken from (Honda et al., 2015).

Fig. 4.3 shows the direction-averaged neutrino flux as a function of neutrinoenergy at Super-K. The plots compare the flux predictions from Honda model withother models such as the Fluka model (Honda et al., 2015) and the Bartol model(Honda et al., 2015). The differences in the flux predictions between Honda model55

and other models are used to estimate the systematic uncertainties on the atmospheric neutrino flux.

Figure 4.3: The direction-averaged atmospheric neutrino flux at Super-K as afunction of neutrino energy. The left plot is absolute flux and the right plot is theflux ratios. Honda flux model is shown in red. Other model predictions are comparedwith Honda model to estimate the systematic uncertainties on the flux. Taken from
(Honda et al., 2015).

The MC events are generated using the predicted flux at Super-K assuming all
neutrinos reach the detector as their original flavor at the production in the atmosphere. Oscillations are applied by reweighting each MC event at the time of
performing an analysis according to its oscillation probability given the neutrinoenergy and direction.

4.2 Simulation of Neutrino Interactions
Atmospheric neutrinos are detected when they interact with water in the tank or
the rock surrounding the Super-K detector. Neutrino interactions at Super-K are
simulated by the neutrino event generator NEUT (Hayato, 2009). The software
simulates the neutrino interactions with nuclei in water and rock. The interactions
with atomic electrons are neglected due to the three orders of magnitude smaller56

cross section relative to that for nuclei at typical atmospheric neutrino energies. The
volumes simulated include the rock surrounding the Super-K detector, the OD andthe ID regions.
Neutrino interactions can classified into charged-current and neutral-current interactions based on which bosons take part in the interaction. A neutrino interacts
with other particles by exchanging a W ˘ boson in a CC interaction in which a charged
lepton of the same flavor as the neutrino is produced. The flavor of a neutrino canbe determined by detecting the charged lepton produced in its CC interactions. On
the other hand, a NC interaction is mediated by a Z 0 boson which does not leave
any signal of the neutrino flavor in the detector.
Five main modes are simulated:• quasi-elastic scattering: ν ` N Ñ l ` N 1• meson exchange current: ν ` N N ˚ Ñ l ` N 1 N 1˚• single meson production: ν ` N Ñ l ` N 1 ` m• coherent pion production: ν `16 O Ñ l `16 O ` π• deep inelastic scattering: ν ` N Ñ l ` N 1 ` hardons
where ν is a neutrino or anti-neutrino of any flavor, N /N ˚ and N 1 /N 1˚ are theincoming and outgoing nucleons respectively, and l is the outgoing lepton. All modes
except meson exchange current can occur in CC and NC types. The interactions aredescribed below.
4.2.1

Quasi-elastic Scattering

Charged-current quasi-elastic scattering (CCQE) is simulated for both free protons
and bound nucleons based on the models of (Llewellyn Smith, 1972) and (Smith
and Moniz, 1972) respectively. A relativistic Fermi gas model is used in the bound57

nucleon calculation. In both cases, the calculation uses vector and axial vector form
factors. The vector mass is set to 0.84 GeV/c2 , and the axial vector mass is set to 1.2
GeV/c2 . Both are tuned to match model prediction with measurement from K2K(Gran et al., 2006) and MiniBooNE (Aguilar-Arevalo et al., 2008). Fig. 4.4 shows

σ (10-38 cm2)

σ (10-38 cm2)

the calculated CCQE cross sections from NEUT and measurement results.
2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

ANL
BNL
SKAT

-

(a) νµn→µ p

0

2

4

–

6
8
10
Eν (GeV)

2

4

12

14

GGM (1977)
GGM (1979)
Serpukhov
SKAT+(b) νµp→µ n0

Serpukhov
GGM (1977)
GGM (1979)

6
8
10
Eν (GeV)12

14

Figure 4.4: CCQE cross sections of νµ (top) and ν̄µ from NEUT calculation, andoverlaid measurement data from ANL (Barish et al., 1977), Gargamelle (Bonetti
et al., 1977; Pohl et al., 1978), BNL (Ahrens et al., 1986), Serpukhov (Belikov et al.,1985) and SKAT (Brunner et al., 1990). Solid lines represent scattering off a free
proton, and dashed lines represent scattering off bound nucleons in 16 O. Taken from
(Hayato, 2009).

NC elastic cross sections are determined by scaling their charged-current quasi58

elastic counterparts (Hayato, 2009).

4.2.2

σpνp Ñ νpq “ 0.153 ˆ σpνn Ñ e´ pq,(4.1)

σpν̄p Ñ ν̄pq “ 0.218 ˆ σpν̄p Ñ e` nq,(4.2)

σpνn Ñ νnq “ 1.5 ˆ σpνp Ñ νpq,(4.3)

σpν̄n Ñ ν̄nq “ 1.0 ˆ σpν̄p Ñ ν̄pq

(4.4)

Meson Exchange Current

In the CCQE interactions on a bound nucleon described above, an assumption ismade that direct correlations between the nucleons in the nucleus are negligible suchthat a neutrino interacts with a single nucleon. However, the correlations among
the nucleons can result in neutrino interaction that involve multiple nucleons, which
is known as meson exchange current interactions (MEC). The MEC is conceived tobe responsible for the tension between the existing CCQE model and the data fromexperiments such as MiniBooNE (Aguilar-Arevalo et al., 2013). MEC interactions
simulation has been recently introduced to NEUT, based on the model by (Nieves
et al., 2011).
4.2.3

Single Meson Production

Single meson production is the dominant hadron production mechanism in the region
where the hadronic invariant mass is less than 2.0 GeV/c2 . Rein and Seghal’s method
is adopted to simulate the neutrino induced single pion production (Rein and Sehgal,
1981). In this method, the interaction is split into two steps as:ν ` N Ñ l ` N ˚, N ˚ Ñ π ` N 1,

(4.5)

where N and N 1 are nucleons, N ˚ is the baryon resonance. In total, 18 resonances
below 2 GeV/c2 are taken into account. Their amplitudes of production, multipliedby the decay probability of a resonance into one pion and one nucleon, are calculated59to obtain the cross sections. Interferences among the resonances are also taken
into account in the calculation. The pion angular distribution in the final stateis calculated from the Rein and Sehgal model for the ∆p1232q resonance, and isassumed to be isotropic for other resonances. The Pauli blocking effect in the decay ofresonance is taken into account by requiring the momentum of nucleon larger than the
Fermi surface momentum. Finally, absorption of pions is also considered, resulting
in about 20% of single pion production interactions having only an outgoing leptonand nucleon. Fig. 4.5 shows the calculated charged current single pion productioncm2)

νµ p → µ+ p π−

−38

1
0.9

σ(ν µ )(10σ(ν µ )(10

−38

cm2)

cross sections and experimental measurements.

0.8
0.7

cm2)
−38

0.7
0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0

0.1
1

1
0.9

νµ n → µ+ n π−

0.8

0.6

0.1

σ(ν µ )(10

1
0.9

10

+

νµ p → µ n π

0Eν (GeV)1

10Eν (GeV)0
BEBC Allen Nucl.Phys.B 264, 221 (1986)

0.8
0.7

BEBC Allasia Nucl.Phys.B 343, 285 (1990)

0.6
0.5

P. Schreiner, GGM, Proceedings of 1979 symposium at Fermilab

0.4
FNAL Barish Phys.Lett.B 91, 161 (1980)

0.3
0.2

BEBC Allisia, Z.Phys.C 20, 95 (1983)

0.1
0

1

10

Eν (GeV)Figure 4.5: Cross section for ν̄µ charged current single pion productions. Solidlines are NEUT calculations, points are experimental data. Taken from (Hayato,2009).

Based on the amplitudes of neutrino resonance productions in Rein and Seghal’smodel, it is also possible to calculate the cross sections of single photon, K, and η
production by changing the decay probabilities of the resonances. These interactions
are also included in the simulation package.
In coherent pion production, the incident neutrino interacts with the entire oxy60

gen nucleus producing a pion. Due to the low momentum transfer in the interaction,
the outgoing lepton and pion produced in this process are peaked in the forward
direction. The interaction is simulated with a modified Rein and Seghal model (Reinand Sehgal, 1983).
4.2.4

Deep Inelastic Scattering

In deep inelastic scattering (DIS), the incident neutrino interacts with a quark in thetarget nucleons and produces multiple hadrons. The process is included when the
hadronic invariant mass ą 1.3 GeV/c2 . The cross section is calculated by integrating
the following equation:
G2 MN Eνyd2 σ
“ F
pp1 ´ y ` y 2 ` C1 qF2 px, q 2 q ˘ ypp1 ´ ` C2 qqxF3 px, q 2 qq,
dxdy
π
2
C1 “

yMl2
xyMN
Ml2
´
´ 2 4Eν2 ´,
4MN Eν x
2Eν
Ml
2MN Eν x

C2 “ ´

Ml2,
4MN Eν x(4.6)

where MN and Ml are the nucleon and outgoing lepton mass, Eν is the neutrinoenergy, q 2 is the squared momentum transfer, and the variables x “ ´q 2 {p2MN pEν ´
El qq and y “ pEν ´ El q{Eν are the Bjorken scaling parameters. The nucleon form
factors, F2 and xF3 , are taken from the GRV98 parton distribution functions (Glucket al., 1998), with corrections from Bodek and Yang (Bodek and Yang, 2002) inthe low q 2 region. The NC cross sections are obtained by multiplying the CC crosssection by linear functions of neutrino energy (Hayato, 2009). The probability ofpion multiplicity, which is a function of W , provides the probability to generate morethan one pion in the region 1.3GeV/c2 ă W ă2.0 GeV/c2 . The average multiplicity
of charged pions is estimated from the Fermilab 15-foot hydrogen bubble chamber
experiment (Derrick et al., 1978) as
xnπ y “ 0.09 ` 1.83lnpW 2 q.
61

(4.7)

The number of pions in each event is determined with KNO scaling (Hayato, 2009).
The forward-backward asymmetry of pion multiplicity in the hadronic center of massis included by using the result of BEBC experiment (Barlag et al., 1982) as:nFπ0.35 ` 0.41lnW 2“.
nB
0.5 ` 0.09lnW 2
π

(4.8)

In the region of Wą2.0GeV/c2 , PYTHIA/JETSET (Sjostrand, 1995) is used todetermine the hadron kinematics.

4.3 Simulation of Charged-current Tau Neutrino Interactions
4.3.1

Simulation of Tau Neutrino Interactions in Water with NEUT

Although the original neutrino flux produced in the atmosphere only consists of electron and muon neutrinos, tau neutrinos are also detected in Super-K due to neutrino
oscillations. NEUT is used to simulate tau neutrino interactions as well. Because
NC simulation has included neutrinos of all three flavors, only charged-current tau
neutrino interactions are simulated and added to the final simulation sample. In
the NEUT simulation, the following charged-current tau neutrino interactions are
considered:• ντ N Ñ τ N 1 , charged-current quasi-elastic scatter• ντ N Ñ τ N 1 π, charged-current single π production• ντ N Ñ τ N 1 γ, charged-current single γ production• ντ N Ñ τ Λ1 π, charged-current single kaon production• ντ N Ñ τ N 1 η, charged-current single η production• ντ N Ñ τ N 1 hadrons, charged-current deep inelastic scattering62

Cross Section (10-38cm2)

Cross Section (10-38cm2)

101
total ντ xsecντ CC1πντ DISντ CCQE xsec

10-1

10-2

10

10

1

total ντ xsecντ CC1πντ DISντ CCQE xsec

10-1

10-2

2

10
Neutrino Energy (GeV)(a) The charged-current ντ neutrino
cross sections as a function of neutrino energy for different channels in
NEUT.

10

2

10
Neutrino Energy (GeV)(b) The charged-current ν̄τ neutrino
cross sections as a function of neutrino energy for different channels in
NEUT.

Figure 4.6: Charged-current tau neutrino cross sections in NEUT.
• ντ16 O Ñ τ πX, charged-current coherent π production
The cross sections for CC tau neutrino interactions are calculated with the samemodels for νµ and νe neutrinos as described in (Hayato, 2009). However, in order toaccount for the large mass of tau lepton, the kinematic effect of the tau lepton mass
is considered in the calculation. The kinematic effect of tau lepton mass requiresthat the minimum neutrino energy in the rest frame of the target nucleon isEmin “ mτ r1 `mτs “ 3.5 GeV
2mN

(4.9)

where mτ is the mass of tau lepton, and mN is the mass of the target nucleon.
Lepton mass corrections are employed in single and coherent pion production basedon models in (Berger and Sehgal, 2007; Rein and Sehgal, 2007). In the simulationof DIS interactions, the GRV98 model (Gluck et al., 1998) is used, and a correction
from (Bodek and Yang, 2002) is implemented to match the DIS cross section andresonance production cross section smoothly. Fig. 4.6 shows the charged-current tau
neutrino cross sections for the three most important channels and the total chargedcurrent cross sections for ντ and ν̄τ in NEUT. As shown in the figure, ντ has a larger
cross section than ν̄τ . CCQE dominates at the low energies (Eν ă 5 GeV), while
63

DIS has a much larger cross section than other channels at high energies (Eν ą10
GeV). At intermediate energies, resonance production is expected to dominate.
4.3.2

Modeling the Polarization of Tau Lepton

The tau leptons produced in the CC tau neutrino interactions are polarized, and thepolarization affects the decay particle distributions of the tau lepton. In the τ rest
frame, in which z-axis is taken along τ ’s momentum direction, a spin polarization
vector is defined as:PÑÝs “ psx , sy , sz q “ psin θP cos φP , sin θP sin φP , cos θP q
2

(4.10)

where θP and φP are the polar and azimuthal angle of the spin vector in the τ rest
frame respectively. The azimuthal angle is measured from the scattering plane whereÝÝ
φP “ π{2 is along the Ñpν ˆÑp τ direction in the lab frame. P represents the degreeof polarization, with P=1 for the fully polarized τ , and P=0 for the unpolarizedcase. The degree of polarization and the spin direction are functions of energy anddirection of outgoing τ . The polarization of tau leptons produced in the neutrino
interactions is modeled in the simulation with the predictions from (Hagiwara et al.,2003) for quasi-elastic, resonance production and deep inelastic respectively. Fig.
Ý
4.7 shows the polarization vector Ñs of τ ´ produced in ντ DIS interactions at Eν =10
GeV on the pτ cos θ-pτ sin θ plane, where pτ and θ are the momentum and scattering
angle of τ ´ in the lab frame.
4.3.3

Simulation of Tau Lepton Decay with Tauola

The tau lepton has a short lifetime of 2.9ˆ10´13 s, thus is not directly detected inSuper-K. The tau lepton decay has multiple channels, including leptonic decays andhadronic decays. The tau lepton decay channels with branching ratio larger than 1%
are presented in Table 4.1 from (Patrignani et al., 2016). Multiple charged particles
64Figure 4.7: The contour plot of the DIS τ ´ polarization in the plane of pτ cos θand pτ sin θ for ντ charged-current DIS process at Eν =10 GeV. The thick gray line
represents the kinematic boundary. The QE scattering contributed along the boundary, and the RES process contributes just inside of the boundary. The length of thearrows shows the degree of polarization, and the direction gives the τ ´ spin in theτ ´ rest frame. Taken from (Hagiwara et al., 2003).
ντν e, ν µ, uτ´W´e´ , µ´ , dFigure 4.8: Feynman diagram of a tau lepton decay.

are produced in most hadronic tau decays. Due to the multiple charged particles
and high-energy lepton, the leptonic decays of the tau look very similar to normal
CC DIS interactions from νµ and νe .
The decays of the tau leptons themselves are simulated by the tau lepton decay
library, TAUOLA (version 2.6) (Davidson et al., 2012). Secondary particles fromTAUOLA simulation are propagated to NEUT and the custom code for further
simulation is the same as that used for the simulation of other atmospheric neutrino

65

Table 4.1: Decay modes of tau leptons with branching ratio larger than 1%. Adaptedfrom (Patrignani et al., 2016).
Decay modeµ´ ν̄µ ντe´ ν̄e ντπ ´ ντπ ´ π 0 ντπ ´ 2π 0 ντπ ´ 3π 0 ντπ ´ π ` π ´ ντπ ´ π ` π ´ π 0 ντh´ ωντ

Branching ratio (%)
17.41 ˘ 0.04
17.83 ˘ 0.04
10.83 ˘ 0.06
25.52 ˘ 0.09
9.3 ˘ 0.11
1.05 ˘ 0.07
8.99 ˘ 0.06
8.99 ˘ 0.06
2.00 ˘ 0.08

interactions.

4.4 Final State Interaction
Hadrons produced inside of nuclei are known to interact with surrounding nucleons,
these interactions are known as final state interactions (FSI). Therefore, the FSI ofpions, kaons, etas and nucleons in the target are also simulated in NEUT. Pions are
the most common hadrons produced at Super-K, and it is important to accurately
simulate pion interactions. Pion FSI include four types of interactions: quasi-elastic
scattering in which a single pion of the same charge remains after the interaction,
charge exchange in which a charged pion is converted into a π 0 or vice versa, pion
absorption in which no pions remain after the interaction and pion production inwhich multiple pions are produced. Pion FSI are simulated by a NEUT cascade
model (de Perio, 2011). In the model, pions are propagated classically in finite
steps through a nuclear medium described by Woods-Saxon nuclear density profile(Woods and Saxon, 1954). At each step, the probabilities for various interactions
are calculated. The model has been tuned based on pion-nucleon and pion-nucleus
scattering data. Fig. 4.9 shows the tuned π ` ´12 C scattering cross sections simulated
in NEUT compared with experimental data.
66

σ (mb)Reactive (TunedFSI)
Reactive (OldFSI)
Quasi-elastic
Absorption
Single CX
π Production

450
400

π- 12
6 C

350
300
250
200
150
100
50
0
0

200

400

600

800

1000 1200 1400 1600
π- Initial Momentum (MeV/c)

Figure 4.9: π ` ´12 C scattering cross sections, where the cross sections for eachinteraction process are shown separately by the colors described in the legend. The
data pions are the measurement from π ` ´ C scattering experiments and the solidlines are the cross sections from the tuned NEUT pion cascade model. The dashedlines are for the previous NEUT pion FSI model. Taken from (de Perio, 2011).

4.5 Detector Simulation
After the simulations of the neutrino interactions and nuclear interactions in NEUT,
final state particles are propagated into a custom detector simulation package called
SKDETSIM. The simulation tracks the particles, models Cherenkov radiation andphoton propagation, and simulates the PMT response. The simulation is based on
GEANT3 (Brun et al., 1987).
Hadronic interactions are simulated with GCALOR (Zeitnitz and Gabriel, 1994)
except for pions below 500 MeV/c for which the NEUT cascade model is used.
Cherenkov radiation is modeled using inputs from the water calibration described in
Section 3.6.2. The scattering and absorption of photons in the water are modeled
based on the calibration result in Section 3.6.2. Reflection and absorption on different detector surfaces is modeled from calibration measurement. PMT response ismodeled using the measured quantum efficiency of each PMT. The charge and time
response of the PMTs and the electronics are also modeled based on the calibration
measurements.
67

5
Atmospheric Neutrino Event Reduction andReconstruction

Super-K records approximately 106 events per day. Most of the recorded events arefrom low background radioactivity surrounding the detector and cosmic ray muons,at a rate of 11 Hz and 3 Hz respectively. Therefore, dedicated reduction algorithms
are implemented to select neutrino events from large number of events in Super-K.
Atmospheric neutrino events in Super-K are classified into three types based on
the event topology: fully contained (FC), partially contained (PC), and upward-going
muon (UPMU), as shown in Fig. 5.1. FC events have the reconstructed neutrino
interaction vertex inside the ID, and the secondary particles deposit the energy in
ID. PC events have the reconstructed vertex in the ID, but the secondary particles
pass through the ID and produce Cherenkov radiation in both ID and OD. UPMU
events are muons produced by the neutrinos coming from the opposite side of the
Earth and interacting with the rock below the Super-K detector. UPMU events arefurther classified into upward stopping muons that lose energy and stop in the ID
and upward through-going muons that traverse the detector volume. These samples
are used in various Super-K atmospheric neutrino analyses. The analysis in this68

dissertation only uses FC events. The FC reduction and the updated PC reduction
are presented below.

UpwardUpward stopthrough-going
Fully contained
Partially containedping muon
muon
Figure 5.1: Schematic view of three types of neutrino events in Super-K.

5.1 Full Contained Sample
The FC reduction is divided into five stages, denoted as FC1-FC5, where events thatpass each step are further processed in order to optimize the processing time. The
reduction cuts were defined by comparing raw data with expectations from the SuperK atmospheric neutrino simulations, with a requirement of at least 99% efficiency forselecting true high-energy neutrino events at each step and rejection of backgroundat the level of less than 1%.
5.1.1

FC1

The first stage requires very basic cuts:
1a. PE300 ą200 (100 for SK-II), where PE300 is the total charge in the ID collected
in a 300 ns wide time window. 200 p.e.s is roughly the charge deposited bya 22 MeV/c electron, which is lower than the 30 MeV energy threshold for69

high energy neutrino analyses. This cut rejects most of the background fromlow-energy radioactivity.
1b. NHITA800 ă 50 (55 for SK-IV), where NHITA800 is the total number of OD
tubes hit within a time window (-500 ns, +300 ns) relative to the trigger time.
This cut removes most of the cosmic ray muons. The cut value is slightly larger
in SK-IV because of the replacement of the OD electronics and changes in the
PMT gains.
After the first reduction, the event rate is reduced from 106 per day to about3000 events per day. The efficiency for selecting atmospheric neutrino events withenergy more than 30 MeV is estimated by applying the FC1 cuts to the atmospheric
neutrino MC and is found to be 100.0%.
5.1.2

FC2

The cuts used in FC2 are defined as
2a. NHITA800 ď25 for events with PEtotal ă 100, 000 (50,000 for SK-II). Thisfurther removes cosmic ray muons, but allows extremely energetic events originating in the ID.
2b. PEmax /PEtotal ă0.5, where PEmax is the largest number of p.e’s observed in asingle PMT. This is to remove electrical noise events which have a single PMT
containing most of the total charge.
Fig. 5.2 illustrates the effects of two cuts in this step graphically. Following the FC2
reduction, the event rate is reduced to about 500 events per day, and the efficiency
of FC2 to select atmospheric neutrino events is estimated to be still 100%.

70

6 (a) DATA(a) DATA

10

5

10

4

10

3

10
5
10
4
10
3
10
2
10

rejected(b) MC10

number of events

number of events

4 (b) MC

10
3
10
2
10
10
1
10

rejected

3 (c) final sample
data

2

MC

10
1
0

5

rejected

10 15 20 25 30 35 40 45 50
FC2 NHITA800

4

10
3
10
2
10
10
1

rejected(c) final sample
3

data
10
MC
2
10
10
1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
FC2 PEmax/PE300

Figure 5.2: Distributions of the atmospheric MC and data before the FC2 cuts 2a(left) and 2b (right). MC events are high-energy FC neutrino events in FV, and data
represent the events that pass the FC1 reduction. Taken from (Lee, 2012).

5.1.3

FC3

The third reduction has a set of complicated cuts that are designed to remove remaining cosmic ray muons and low-energy radioactivity.
Through-going muon cut
Through-going muons are very energetic and deposit an extremely large amountof charge in the ID. In order to remove these events, a through-going muon fitter isapplied on events that have more than 1000 PMTs and collect more than 230 p.e.s.
The fitter finds the entering point of the muon in the ID by searching for the earliest
hit PMT with neighbor hit PMTs. The exiting point of the muon is defined as thecenter of the saturated ID PMTs. The goodness of the through-going muon fit isdefined as:
goodness “ ř

1

ÿ 1

1 ˆ
i σi2
i

σi2

e

´

pti ´Ti q2
2p1.5ˆσi q2

.

(5.1)

Events that meet the following criteria are classified as through-going muon, andtherefore are removed:
3a. goodnessą0.75
71

3b. N HIT Ain ě10 or N HIT Aout ě10, where N HIT Ain (N HIT Aout ) is the number of OD hit PMTs located within 8 m from the entrance (exit) point within
a 800 ns time window.

30(a) DATA

25

25

20

20
rejected

15

10

5

5

30

0

10

20

30

(c) final sample DATA

0
30

25

25

20

20

15

15

10

10

5

5

0

0

(b) MCrejected

15

10
0

NHITAout

30

10
20
30
FC3 NHITAin

0

0

10

20

30

(c) final sample MC

0

10

20

30

Figure 5.3: Two-dimensional scatter plot of NHITAin (x-axis) and NHITAout forthrough-going muons selection in FC3 reduction for SK-III. (a) shows the data after
the FC2 reduction. (b) shows MC events after the FC2 reduction. (c) shows the final
sample for data (bottom left) and Monte Carlo (bottom right) events. Cut values
are indicated by solid lines. Taken from (Lee, 2012).

Events that satisfy the above criteria are tagged as through-going muons andare rejected. Fig. 5.3 shows the scatter plot of NHITAin and NHITAout that havethrough-going muon goodness larger than 0.75 for data and MC after the 2nd reduction and the final FC sample.
72Additionally, a cut is applied to remove extremely energetic through-going muons
which are expected to deposit more energy via pair production and bremsstrahlung
than ionization.

The energy deposition per unit track length is higher than 2

MeV/cm. Events satisfying the following cut are removed:
3c. NHITA500 ą 400, where NHITA500 is the number of OD hits within 500 ns
sliding time window.
Stopping muon cut
The stopping muon selection follows a similar logic as for the through-going muon
selection, except there is no exit point. The direction of the muon is found by
maximizing the charge seen by the PMTs within a 42-degree opening angle fromthe entry point. A goodness parameter with the same definition as Equation 5.1 iscalculated. The criteria for selecting stopping muon events are defined as:
3d. NHITAin ě10 ( or NHITAin ě 5 if goodness of stopping muon fit ě 0.5 for
SK-I)
Cable hole cut
On top of the stainless steel structure where PMTs are mounted, there are twelve
cable holes for the signal and high voltage supply cables to pass through. There is
no OD PMT installed at the positions of these holes. Muons passing through theseholes may not be noticed by the OD. In order to compensate for the dead area in
the OD, veto counters (2 mˆ2.5 m plastic scintillation counters) were installed in
April, 1997 at four large cable holes. (For the other eight smaller holes without vetocounter, software background rejection will be applied in FC5.) Fig. 5.6 shows the
positions of the twelve cable holes and the 4 veto counters.
An event is removed if:
3e. cable hole veto counter is hit.
73

number of events(a) DATA

10

4

10

3

4
10
3
10
2

rejected(b) MCrejected

10101(c) final sample
3

10
2
10
10
1
0

data
MC

5

10 15 20 25 30 35 40 45 50
FC3 NHITAent

Figure 5.4: Number of hit OD PMTs around the entrance point of muons for eventsin FC3 stopping muon reduction in SK-IV. (a) shows the data after FC2 reduction.
(b) shows FC MC events after FC2 reduction. (c) shows the final samples for data
and MC after FC reduction. Taken from (Lee, 2012).

3f. Lveto ă4 m, where Lveto is the distance from the cable hole to the event vertex.
In SK-IV, a tighter selection is applied to reject stopping cable hole muons. Events
that satisfy the following criteria are removed:
3g. goodness ą 0.4;
3h. PEtotal ą 4000;
3i. cospθq ă-0.6, where θ is the muon’s zenith angle;
3j. The entry point of is within 2.5 m of the nearest cable hole.
Flasher event cut
A “flasher” is caused by repetitive discharge in the dynode of a PMT, which
may also affect nearby PMTs. Flasher events typically have broader hit timing
74

Figure 5.5: Schematic view of the positions of the twelve cable holes on the stainlesssteel support. Four holes in red circles are installed with veto counter. Taken from
(Fukuda et al., 2003).

Figure 5.6: Schematic view of a cable hole muon and a veto counter. Taken from
(Fukuda et al., 2003).

75

distributions than neutrino events, as shown in Fig. 5.7. Therefore, events thatsatisfy the following detector-specific criteria are rejected:
3k. (SK-I) NMIN100 ě 15, or NMIN100 ą10 if Nhit ď800, where NMIN100 is theminimum number of ID PMT hits in a 100 nsec time window from 200 ns to
700 ns after the trigger time, and Nhit is the total number of ID hits.
3l. (SK-II to IV) NMIN100 ě 20.

Figure 5.7: The timing distribution of a flasher (top) and an atmospheric neutrino
MC event (bottom). The arrows indicate the time window for counting NMIN100 .
Taken from (Lee, 2012).

Fig. 5.8 shows the NMIN100 distributions for data and FC MC after FC2 reduction and in the final samples.
Coincidence muons
Sometimes a cosmic ray muon enters in the trigger window of a low energy event.
The coincidence avoids the cuts designed to remove low-energy events because of theenergy deposition of the muon. Additionally, it avoids the high-energy cuts because
the OD activity of the muon is recorded later than the low-energy trigger. In orderto remove these event, the following cuts are applied:

76

10

5

10

4

10

5

10

10

4

3

10

3

10

2

10

2

(a) data

rejected

10

10

5

10

4

10

3

10

2

Number of events

Number of events

10(a) data

rejected(b) MC

10
1
10

4

10

3

(c) final sample

10

2

data
MC

10

10

5

10

4

10

3

10

2

(b) MC

10
1
10

4

10

3

10

2

(c) final sample
data
MC

10

1

1
0

5

10

15

20

25

30

35

40

45

50

NMIN100 (NHITID > 800)

0

5

10

15

20

25

30

35

40

45

50

NMIN100 (NHITID < 800)

Figure 5.8: NMIN100 distributions for events with more than 800 hit ID PMTs
(left) and less than 800 hit ID PMTs (right). (a) shows the data after the FC 2nd
reduction. (b) shows FC Monte Carlo events after the FC 2nd reduction. (c) shows
the final samples for data and the Monte Carlo events, in which the number of Monte
Carlo events is normalized to that of the data. Selection cuts are indicated by arrows.
Taken from (Lee, 2012)

3m. NHITAof f ě 20, where NHITA20 is the number of OD hits between 300 nsec
and 800 nsec after the trigger time;
3n. PEof f ą 5000 (2500 for SK-II), where PEof f is the total p.e.’s in the ID for thesame time window.
Low energy events
A specific cut is applied aiming to reject low-energy events. Events are removed
if NHIT50 ă 50 (25 for SK-II), where NHIT50 is the number of ID hits within asliding 50 nsec tie window. The time window uses the residual time that takes
into account the time-of-flight of the photons assumed to come from a single vertex.
NHIT50 corresponds to visible energy of 9 MeV, which is safely lower than the energy
threshold of 30 MeV for atmospheric neutrino analyses.
Following the FC3 reduction, the event rate is around 50 events per day. The
efficiency of FC3 to select true neutrino events is estimated to be 99.8%.
77

5.1.4

FC4

The fourth reduction is designed to remove flasher events that were not removed by
the above cuts. Flasher events usually emit light in a specific way, and create similar
hit patterns across several event timing windows. Therefore, a pattern-matchingalgorithm is applied to reject them. The ID tank wall is divided into 1450 roughly
equal 4 m2 patches, and the total charge (npe ) in each patch is computed. The
correlation between two events (A and B) is computed as:

r“

B
A
B
1 ÿ pQAi ´ xQ yq ˆ pQi ´ xQ yq,
N i
σA ˆ σB(5.2)

where N is the number of patches, Qi is the summed npe in each patch i, and σ isthe standard deviation of the Qi . If the distance between the PMTs with maximum
charge in the two events is smaller than 75 cm, an offset of 0.15 is added to r. Anothermatching parameter, dks is also computed by comparing the time distribution of thecharge by a Kolmogorov-Smirnov test. These two parameters are combined into alikelihood parameter “Prob”. The events with Prob exceeding 3.0 are removed.
The event rate after FC4 reduction is about 32 events per day, with a neutrino
efficiency of 99.6%.
5.1.5

FC5

The fifth reduction is designed to reject three specific categories of remaining background: invisible muons, coincidence muons, and long-tail flashers.
Invisible muons Cosmic ray muons with momenta below the Cherenkov threshold will pass the above reduction stages, but their decay electron can still be recorded.
These events are rejected by activity in the OD. Events satisfying the following criteria are removed:
5a. PEtotal ă 1000 (500 for SK-II);78

5b. NHITACearly + NHITAC500 ě10 if Distclust ă 500 cm; NHITACearly ą 9 otherwise, where NHITACearly is the maximum number of OD PMT hits in thecluster within in a 200 nsec sliding time window from -9000 to -200 nsec prior
to the primary trigger time, NHITAC500 is the hits in the time window of -200
to 300 nsec, and Distclust is the distance between the two OD clusters.
Coincidence muons
The remaining coincidence muon events are further removed by the following
cuts:
1. PE500 ă 300 p.e.s (150 p.e.s for SK-II), where PE500 is the total number of
p.e.s recorded in the ID in a fixed time 500 nsec time window from -100 nsec
to +400 nsec.
2. PElate ě 20 p.e.s, where PElate is the maximum number of OD hit PMTs in a
200 nsec sliding time window fro +400 nsec to +1600 nsec.
Sub-cable hole muons
Besides the cable hole muon cut in FC3, another cut is introduced in SK-IV toremove cosmic ray muons passing through eight sub-cable holes where veto counters
are not installed. The events satisfying the following criteria are rejected:
1. Goodness of stopping muon fit ě 0.4;
2. PEtot ą 4000;
3. downward direction (cos θ ă-0.6);
4. Lveto ă 2.5 m.
Long-tail flasher
A tighter flasher reduction than in FC3 reduction is applied, which removes events
that satisfy the following criteria:79

1. NMIN100 ě6 if goodness of Point-fită0.4, where NMIN100 is the minimum
number of ID hit PMTs in a sliding time window from +300 nsec to +800
nsec.
From SK-II, the following cuts at added
2. NMIN100 ď5 if goodness of Point-fită0.3.
Following FC5 reduction, FC event rate becomes „20 events per day, and the estimated neutrino efficiency is 99.9%.
5.1.6

Final Reduction

In the FC data set after FC1-5 reduction, neutrino events that satisfy the following
criteria are selected to be used in analyses:
1. The distance from the reconstructed neutrino event vertex to the ID wall is ą
2 m.
2. Number of OD PMTs in the OD cluster is ă 16 (10 for SK-II).
3. Evis ą 30 MeV. Evis is defined as the energy of an electron that would produce
the same amount of p.e.s as that in a event.
The cuts are expected to largely remove the remaining background events, suchas radioactivity. The event rate after the above cuts is approximately 8 events perday. The neutrino efficiency of these cuts is estimated to be 99.2%. The selection
efficiencies and the event rates in each Super-K periods are summarized in Table 5.1.
There is still a tiny portion of background events in the final sample, mostly
cosmic-ray muons and flasher events. The final sample is eye-scanned, and the background events are tagged; however they are not removed from the analysis, in order
to avoid a bias in the data and MC selection. The background contamination issummarized in Table. 5.2.
80

Table 5.1: Summary of selection efficiencies in each stages and event rates in the
final FC samples for each Super-K period.
FC
SK-I
SK-II
SK-III
SK-IV
reduction step
selection efficiency (%)
1st
100.0
99.97
100.0
100.0
2nd
100.0
99.92
99.98
99.99
3rd
99.93
99.78
99.81
99.82
4th
99.29
99.38
99.30
99.00
5th
99.26
99.30
99.24
98.95
5th (FV)
99.25
99.95
99.62
99.19
Event rate
event rate/day 8.21˘0.07 8.24˘0.10 8.41˘0.12 8.09˘0.09

Table 5.2: Estimated background contaminations in each Super-K period. Sub-GeV
means events with Evis ă 1.3 GeV and Multi-GeV means events with Evis ą 1.3
GeV.
BG
sub-GeVmulti-GeV

SK-I
cosmic ray flasher
0.2%
0.1%
0.3%
0.5%

SK-II
cosmic ray flasher
0.2%
0.4%
0.2%
0.2%

SK-III
cosmic ray flasher
0.3%
0.1%
0.3%
0.2%

SK-IV
cosmic ray flasher
0.1%
0.1%
0.1%
0.1%

5.2 Event reconstruction
Following the FC reduction, the remaining events are put through a reconstruction
algorithm in order to convert the charge and timing information of the PMTs into
physical variables. The reconstruction starts from finding the vertex position ofan event using PMT timing information. Following this, a ring-fitting algorithm
based on Hough transformations is applied to identify Cherenkov rings. Finally, theparticle identification (PID), and momentum, are calculated and decay electrons aresearched for. I describe the reconstruction of fully-contained events here following
the sequence of the reconstruction stages.

815.2.1Vertex Fitting

The event vertex is determined in three steps. In the first step, a simple fit called
Point-fit is applied. It assumes that the ID PMT hits are emitted from a point source,then uses ID PMT hit residual times relative to that point source to maximize agoodness of fit parameter:ˆ˙pti ´ t0 q2
1 ÿ
exp ´
,
GP “
N i
2p1.5 ˆ σq2(5.3)

where N is the total number of ID PMT hits, ti is the residual time of the ith hit PMT,
t0 is a free parameter to maximize the goodness GP , and σ is the time resolution
equal to 2.5 ns. The point that gives maximum goodness is taken as the event vertex.
Then, the overall direction of an event is estimated as the sum of the charge-weighted
directions from the reconstructed vertex to each hit PMT.
With the reconstructed vertex, the Cherenkov emission angle of the brightest ringis determined by searching for the edge of the ring. A goodness of fit parameter isdefined with the assumed angle of the ring as:şθedgeQpθedge q “0ˆ
˙
P Epθqdθ
pθedge ´ θexp q2
dP Epθq 2
q ˆ exp ´,ˆpsin θedge
dθ
2σθ2

(5.4)

where PE(θ) is the angular distribution of the observed charge as a function of
opening angle, θ, θexp and σθ are the Cherenkov opening angle expected from thecharge within the Cherenkov cone and its resolution respectively. The charge in
P Epθ) is corrected for attenuation and PMT angular acceptance. The ring direction
and opening angle θedge are varied to maximize Qpθedge ).
Finally, the vertex direction is determined by a more precise fitter, TDC-fit, thattakes into account scattering of light and the track length of the charged particle.
Instead of calculating the residual time relative to one assumed vertex, the residual82

time is calculated as:
1 Ñ
n Ñ
Ý
Ñ
Ý
Ý
Ñ
Ý
ti “t0i ´ X i ´ O i ´ P i ´ X i , for PMTs within the Cherenkov ring edge;
c
c
(5.5)
ti “t0i ´

n Ñ
Ý
Ñ
Ý
P i ´ X i , for PMTs outside the Cherenkov ring edge
c

(5.6)

Ñ
Ý
Ñ
Ý
where O i is the vertex position, X i is the emission point along the particle trackÑÝthat produce hit i, P i is the position of ith PMT, and n is the index of refraction.
The goodness parameter is again maximized to determine the vertex.
The vertex reconstruction is evaluated using MC by comparing the reconstructedvertex with the true MC vertex. The resolution for multi-ring samples ranges from50 to 100 cm.
5.2.2

Ring Counting

After the brightest Cherenkov ring and the vertex has been determined, other ringcandidates are searched for using a Hough transform method. Fig. 5.9 illustrates thebasic logic of the algorithm. By drawing charge-weighted virtual rings with 42-degree
angle around each hit PMT, the direction of the ring is found at the point wheremost virtual rings overlap.
Ring candidates are evaluated iteratively using a likelihood method. The hypothesis of (N+1) rings is tested against the hypothesis of N rings, up to a maximum of
five rings. During this process, previously found rings are kept fixed. The likelihood
is defined as:
˜
LN “

ÿ

log P pq obs ,

N
ÿ

¸
exp
αn qi,n
q ,

(5.7)

n“1

i

where αn is a weight for ring n that is varied to maximize the likelihood, and P is a
Gaussian probability:ˆ˙
pqiobs ´ qiexp q2
exp ´
.
2σ 2
2πσ 2

1
P pqiobs , qiexp q “ ?

83

(5.8)

42 deg. ringhit PMT(possible center)

centerCherenkov ring

(most probable)

Figure 5.9: The Hough transform method. Virtual rings with a 42-degree opening
angle are drawn around each hit PMT, and the ring center is found at the point at
which most virtual rings overlap.

After finding all rings, particle type and momenta are estimated for each of them.
Finally, the FC data are separated into single-ring and multi-ring events by cutting
the likelihood distribution at zero. The likelihood distributions for sub-GeV andmulti-GeV events are shown separated in Fig. 5.10.
Super Kamiokande I 1489.2 days
350

single ring

Super Kamiokande I 1489.2 days
400

multi ring

single ring

multi ring

350

300

300
250
250
200
200
150

150

100

100

50
0
-15

50
-10

-5

0

5

10

15

Ring counting likelihood multi-GeV(FC)

0
-15

-10

-5

0

5

10

15

Ring counting likelihood sub-GeV(FC)

Figure 5.10: The ring likelihood distribution for sub-GeV events (left) and multiGeV events (right) for SK-I. The data is shown in the black dots, and the MC withneutrino oscillations applied at ∆m2 =2.1ˆ10´3 eV2 and sin2 p2θq=1.0 is shown as thehistogram. A cut in the likelihood separates single-ring from multi-ring events.

845.2.3

Particle Identification

Each ring found in ring counting is classified in two broad categories: “e-like” or“µ”-like by particle identification (PID). “e-like” events that are produced by electromagnetic showers and have a diffused ring pattern. High-energy electrons produce
electromagnetic showers in the water via bremsstrahlung (e˘ Ñ e˘ γ), and photon
pair production (γ Ñ e´ e` ), and multiple scattering of low energy electrons. Highenergy photons also produce electromagnetic showers via pair production. Because
the cross section of bremsstrahlung is proportional to 1{M 2 , the more massive muons,charged mesons and protons rarely produce electromagnetic showers, especially inthe full-contained events.
At typical energies of atmospheric neutrinos, muons interact primarily via ionization, and produce rings with sharp edges. Moreover, Cherenkov rings from these
massive charged particles generally have smaller Cherenkov angles. Charged-mesons,particularly pions, also have sharp edges due to the lack of electromagnetic showers. However, the rings they form are usually thinner than those of muons because
hadronic interactions can suddenly alter their path. Fig. 5.11 shows examples of elike and µ-like rings. The particle identification algorithm (PID) works by comparing
the observed charge distribution with an expected charge distribution under e-like orµ-like hypothesis, together with the opening angle information. The expected charge
distribution for the ith PMT are defined as:R
1
qiexp peq “αe ˆ Qexp ppe , θi q ˆ p q1.5 ˆ ri {L ˆ f pΘq ` q scat ,rieqiexp pµq “pαµ ˆsin2 θxi
1` qiknock q ˆ ri {L ˆ f pΘq ` q scat .
dθeri psinθxi ` ri dx |x“xi q(5.9)(5.10)

where αe,µ are normalization factors, θi is the angle between the ith PMT and the ring
direction, f pΘq is the PMT acceptance as a function of incident angle Θi , ri is the
distance between the vertex and the ith PMT, L is the photon attenuation length,85

Super-Kamiokande IV
Run 999999 Sub 0 Event 1
17-01-23:21:39:55
Inner: 3589 hits, 8903 pe
Outer: 7 hits, 8 pe
Trigger: 0x03
D_wall: 1690.0 cm
Evis:

0.0 MeV

Charge(pe)>26.7
23.3-26.7
20.2-23.3
17.3-20.2
14.7-17.3
12.2-14.7
10.0-12.2
8.0-10.0
6.2- 8.0
4.7- 6.2
3.3- 4.7
2.2- 3.3
1.3- 2.2
0.7- 1.3
0.2- 0.7< 0.2
1700
0 mu-edecays

1360
1020
680
340
0
0

500

1000

1500

2000

Times (ns)Super-Kamiokande IV
Run 999999 Sub 0 Event 1
17-01-23:21:36:51
Inner: 2048 hits, 7555 pe
Outer: 8 hits, 10 pe
Trigger: 0x03
D_wall: 1690.0 cm
Evis:

0.0 MeV

Charge(pe)>26.7
23.3-26.7
20.2-23.3
17.3-20.2
14.7-17.3
12.2-14.7
10.0-12.2
8.0-10.0
6.2- 8.0
4.7- 6.2
3.3- 4.7
2.2- 3.3
1.3- 2.2
0.7- 1.3
0.2- 0.7< 0.2
950
1 mu-e
decay

760
570
380
190
0
0

500

1000

1500

2000

Times (ns)Figure 5.11: Simulated e´ and µ´ , both with momentum 1000 MeV/c. The ringpattern is sharper for µ´ and more diffuse for e´ . The color scale indicated theamount of charge deposited on the PMTs. The main display shows the ID, while
the small display shows on the top right shows the OD, in which there is negligible
activity.

86and qiscat is the expected charge due to scattering for the ith PMT. For the e-likehypothesis, MC simulations are used to estimate Qexp ppe , θi q at a vertex located at
R “ 16.9 m from the ID wall. The factor of pR{ri q1.5 accounts for the difference in
diffusion between the test track length and the actual track length. For the µ-likehypothesis, xi refers to the estimated point along the muon track from which photons
propagate to the ith PMT. The hypothesis also include contribution from “knock-on”
electrons (q knock ) that are ionized by muons.
A likelihood function is calculated with the observed and expected charge distribution for the nth ring as:
Ln pe, µq “

ź

ÿ

exp
P pqiobs , qi,npe, µq `

exp
qi,n
1 q,

(5.11)n1 ‰n

θi ă1.5ˆθC

where the product is over all the PMTs within the cone with an opening angle ofθc ˆ1.5, P is the Gaussian probability function defined in Equation 5.8, and thesum over n1 accounts for light from overlapping rings. The likelihood is maximizedby altering the direction and opening angle of each ring in order to optimize theexpected charge distributions. The charge distribution likelihood is then converted
to a χ2 parameter,χ2n pe, µq “ ´2 log Ln pe, µq ` const,(5.12)

which is then used to construct a probability parameter:
¸
2
2
2
2
pχ
pe,
µq
´
minpχ
peq,
χ
pµqqq
n
n
,
Pnpattern “ exp ´ n
2σχ2 2n˜

where σχ2n “(5.13)?
2N for N PMTs in the likelihood calculation.

For multi-ring events, Pnpattern pe, µq is the final PID probability. For single-ring
events, a second probability is calculated based on the Cherenkov angle:
P

angle

ˆ
˙
pθobs ´ θexp pe, µqq2
pe, µq “ exp ´
,
2δθ2
87

(5.14)

where θobs is the reconstructed Cherenkov angle, δθ is the angle fitting error, andθexp pe, µq is the expected angle based on the reconstructed e{µ-like momentum.
The final PID probability of single-ring events is the product of these two probability estimators:
P pe, µq “ P pattern pe, muq ˆ P angle pe, µq.

(5.15)

The final PID likelihood is defined as:
LP ID “aa
´ log P pµq ´ ´ log P peq.

(5.16)The distributions of the likelihood variable for single-ring events are shown in Fig
5.12. The FC data are separated into e-like and µ-like by cutting the PID variable
at zero.
Super Kamiokande I 1489.2 days

Super Kamiokande I 1489.2 days
160

300

e−like 3426

muon−like 3221

e−like 829

muon−like 694

140

250

120

200

100

CCQE muon

CCQE muon

80

150

60
100
40
50
0
−10 −8

20
−6

−4

−2

0

2

4

6

8

10

PID likelihood sub−GeV 1ring (FC)

0
−30

−20

−10

0

10

20

30

PID likelihood multi−GeV 1ring (FC)

Figure 5.12: The distribution of particle identification likelihood for sub-GeV (top)and multi-GeV (bottom) single-ring FC events. The hatched histogram shows the
MC contribution due to CC νµ interactions, the black dots show the data and theempty histogram shows the total MC with neutrino oscillations applied as in Fig.
5.10

5.2.4

Momentum Reconstruction

The momentum of a ring is determined by the observed p.e.s inside a Cherenkov
cone with an opening angle of 70˝ and the energy scale described in Chapter 3. The
88

observed p.e.s in hit ID PMTs are first separated into distributions from each ringbased on the expected p.e. distribution as:
exp
qi,n
obs
obs
qi,n “ qi ˆ ř exp ,n1 qi,n1(5.17)thwhere qobs
ring in the ith PMT, qobs
(qexpi,n is the fraction charge from the nii,n ) is theobserved (expected) p.e.s in the ith PMT. The sum of observed charge from a ringis used to determine its momentum. In order to avoid charge from decay electron,
the charge is restricted to a (-50, 250) ns time window relative to the time-of-flight
subtracted timing peak. After taking account of light attenuation, light scattering,
angular acceptance and PMT gain, the charge of a ring is calculated as:
RT OT “ÿÿGM C
cos Θipα ˆ pqiobs ˆ eri {L ˆq ´ Si q,
GDAT Af pΘi qii

(5.18)

MC
is the relative PMT gain for MC and data, α is a normalization factor,
where GGDATA

ri is the distance from the ith PMT to the event vertex, L is the water attenuation
length, f pΘq is the acceptance of the PMT for incident angle Θ, Si is the expectedcharge for the ith PMT due to scattered photons.

5.3 Decay Electron SearchMuons that stop in the detector decay into electrons that may be observed. Decay
electrons provide a relative independent information of event from the Cherenkov
rings produced by the primary particles. They are searched for in a 20 µs time
window after the event trigger. Decay electron candidates are divided into primaryor “in-gate” events or sub-events, depending on their time. In SK-I-III, the search
criteria for primary events are:• 0.1 µsătă 0.8 µs;
89• NHIT30 ě 60 (30 for SK-II).
where t is the time after the event trigger, and NHIT30 is the number of hits in a 30
ns time window. For sub-events, the requirements are defined similarly as:• 1.2 µsătă 20 µs;• NHIT30 ě 40 (20 for SK-II).
Additionally, both primary event and sub-event are required to have at least 50 hit
ID PMTs (25 for SK-II), and the vertex goodness (Equation 5.3) greater than 0.5.
The decay electron tagging efficiency is estimated with MC. The efficiency is 80%
for µ` and 63% for µ´ in SK-I-III. Due to the improvement of electronics in SK-IV,
the time between 0.8 µs and 1.2 µs is included in the search. Therefore, the efficiency
has increased to 96% for µ` and 80% for µ´ in SK-IV.

5.4 Summary of Reduction and Reconstruction
The data and MC are passed through the same reduction and reconstruction algorithms. The reduction selects the fully-contained atmospheric neutrino events in thefiducial volume with an efficiency of more than 99. The FC events in the fiducial volume are estimated to contain less than 1% background events. A measurement of thephysical quantities of an event including vertex, the number of Cherenkov rings, momentum, particle identification, and the number of decay electrons, is automatically
performed by the reconstruction algorithms with good quality. The fully-contained
events in fiducial volume is divided based on the observed energy to sub-GeV, below
1.33 GeV, and multi-GeV, above 1.33 GeV. The multi-GeV events are used in thisanalysis. The selection efficiencies for this cut of energy are 86% for the ντ CC signal
and 23% for the background respectively.

90

6
Atmospheric Tau Neutrino Appearance in Super-K

Super-K is credited for the discovery of neutrino oscillation by measuring a deficit
in atmospheric muon neutrinos. In three-flavor neutrino oscillations, the deficit isexplained by the change of muon neutrinos to tau neutrinos. An unambiguous measurement of tau appearance would be a direct evidence for three-flavor neutrino
oscillations. The details of atmospheric tau neutrino appearance and the technique
to measure it are presented in this chapter.

6.1 Atmospheric Tau Neutrino Appearance from Neutrino Oscillations
In three-flavor neutrino oscillations, tau neutrinos are expected to appear in theatmospheric neutrinos in the propagation of electron neutrinos and muon neutrinos.
The probabilities of tau neutrino appearance from muon neutrino or electron neutrino
in vacuum are calculated approximately as:
P pνµ Ñ ντ q “ cos2 θ13 sin2 p2θ23 q sin2 p1.27∆m223

Lq
E(6.1a)

Lq
E(6.1b)

P pνe Ñ ντ q “ sin2 p2θ13 q cos2 pθ23 q sin2 p1.27∆m223
91

Because neutrinos coming from below travel through the Earth, the oscillation
probabilities are also altered by the matter effect. Instead of using the tau appearance
probabilities in vacuum in Eq 4, a custom code by Wendell et al. (Wendell et al.,2012) is used to calculate the oscillation probabilties, which takes into account theeffect of neutrino types, path lengths, neutrino energies and matter effect. The
oscillation parameters used are sin2 2θ23 “ 1.0, ∆m223 “ 2.1 ˆ 10´3 eV 2 (Abe et al.,2013), ∆m221 “ 7.6 ˆ 10´5 eV 2 (Barger et al., 1980), sin2 2θ13 “ 0.099 (An et al.,2012; Ahn et al., 2012), δCP “ 0. The matter effect model from (Barger et al., 1980)
is included in the calculation of oscillation probabilities. Fig. 6.1b presents the tau
neutrino appearance probabilities from muon neutrinos and electron neutrinos inthree-flavor neutrino oscillation under the assumption of normal hierarchy. Givena neutrino energy and a neutrino path length, the probability of a muon neutrinochanging flavor to a tau neutrino is larger than that of electron neutrino. In addition,
muon neutrinos have a larger flux than electron neutrinos (shown in Fig. 4.3).
Therefore, tau neutrino appearance is dominantly from muon neutrino oscillation inatmospheric neutrinos.
By oscillating the electron neutrinos and muon neutrinos in the Honda flux prediction to tau neutrinos, we can calculate the expected atmospheric tau neutrino flux
from neutrino oscillations at Super-K asΦντ pEν q “

ż ÿ
e,µ
4π

Φνl pEν , θz , φA qOνl Ñντ pEν , θz , φA qdΩ

(6.2)

l

where Φνl pEν , θz , φA q is the differential flux as a function of neutrino direction (θz , φz )and
neutrino energy (Eν ) for electron neutrinos or muon neutrinos from Honda model,Oνl Ñντ pEν , θz , φA q is the probability of electron or muon neutrino changing flavor to
tau neutrino as a function of neutrino direction and neutrino energy from Prob3++(Wendell et al., 2012). Fig. 6.2 shows the expected atmospheric tau neutrino fluxes
92

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

0.5
0
-0.5
-1

1

Cosine Zenith Angle

Cosine Zenith Angle

1

10

1
0.6
0.5

0.5

0.4
0

0.3
0.2

-0.5

0.1
-1

1

0

10

Energy [GeV]

Energy [GeV](a) Probability of muon neutrinochanging flavor to tau neutrino asa function of neutrino energy andzenith angle.

(b) Probability of electron neutrinochanging flavor to tau neutrino asa function of neutrino energy andzenith angle.

Figure 6.1: Probability of tau neutrino appearance from muon neutrinos (left) andelectron neutrinos (right) as a function of neutrino energy and zenith angle.

Neutrino Flux (GeV-1/m2/s)

from neutrino oscillations at Super-K integrated over neutrino directions.

1ντ flux w osc

10-1ντ flux w osc

10-2
10-3
10-4
10-5
10-6
10-7

20

40

60

Neutrino Energy (GeV)

Figure 6.2: Flux of atmospheric ντ (red) and ν̄τ (blue) from neutrino oscillations
as a function of neutrino energy at Super-Kamiokande for the three-flavor oscillation
parameters given in the text.

93

q

ν̄τ

τ`

Z
W´

W`

ντ
ντ
q
ντ
τ´
Figure 6.3: Feynman diagrams of weak interactions that include a tau neutrino.

6.2 Charged-current Tau Neutrino Events in Super-Kamiokande
Because the neutral-current interaction is flavor blind, it is not possible to distinguish
the flavor of a neutrino in neutral-current interactions. However, a charged leptonof the same flavor as the neutrino is produced in the charged-current interactions.
The detection of the charged lepton is an unambiguous confirmation of the neutrino
flavor. Definite proof the appearance of tau neutrinos requires detection of tau
leptons produced in charged-current tau neutrino interactions.
1ντ ` N Ñ τ ´ ` N ` . . .

(6.3)

A charged-current tau neutrino interaction produces a tau lepton, which has amass of 1.77 GeV/c2 (Patrignani et al., 2016). The relatively large mass of thetau lepton greatly suppresses the cross section of charged-current tau neutrino crosssection at low energies and results in an energy threshold of 3.5 GeV. The chargedcurrent tau neutrino cross sections have been calculated by several groups (Jeongand Reno, 2010; Levy, 2009). In this study, we use the NEUT simulation package(Hayato, 2009) to model the neutrino cross section. Charged-current tau neutrino interactions are classified into three categories in NEUT: quasi-elastic charged-current
interactions, charged-current resonant interactions and deep inelastic interactions.
By folding the flux with the cross sections, we can calculate the expected tau neu94

Cross Section (10-38cm2)

10

1

10-1
NEUT CC ντ xsec

10-2
NEUT CC ντ xsec

10-30

20

40

60

Neutrino Energy (GeV)

Figure 6.4: Total charged-current cross section as a function of neutrino energy forντ (red) and ν̄τ (blue) for neutrino energy between 3.5 GeV and 70 GeV in NEUT.

trino events in Super-Kamiokande. The charged-current tau neutrino cross sections
suppressed by large tau lepton mass results in a low tau neutrino interaction rate,
especially at typical atmospheric neutrino energies. Moreover, the tau neutrino flux
from neutrino oscillations is smaller than that of electron neutrinos or muon neutrinos. Therefore, tau events are rare compared with the electron and muon neutrino
events, which will be background in this analysis.
Tau leptons are also difficult to detect because of the extremely short lifetimeof the tau lepton. The tau lepton decay has multiple channels, including leptonic
decays and hadronic decays. The tau lepton decay channels with branching ratio
larger than 1% are presented in Table 4.1 from (Patrignani et al., 2016). Leptonic
decays of tau look similar to normal CC DIS interactions from νµ and νe . Multiple
charged particles are produced in most hadronic tau decays. For this reason, theanalysis is developed to select hadronic tau decays.

6.3 A Neural Network Algorithm for Tau Neutrino Identification
As described in Sec. 6.2, tau leptons produced in CC ντ interactions decay quickly tosecondary particles, and the dominant hadronic decays produce extra pions. Becauseof the short lifetime of tau lepton, it is not possible to directly detect tau leptons in95ντν e, ν µ, uτ´W´e´ , µ´ , dFigure 6.5: Feynman diagram of a tau lepton decay.
Super-K. The extra pions from the hadronic decays of the tau allow the separationof the signal from νµ and νe and NC background statistically. Fig. 6.6 compares
the simulations of signal and background events in Super-K for example events. It
is noted that the signal event has multiple rings produced by charged particles inwater compared with a single ring in the background event. However, deep-inelastic
neutrino interactions can also produce multiple charged particles in the detector that
result in multi-ring patterns. Fig. 6.7 shows an example of background event that hasa deep-inelastic interaction. The existence of deep-inelastic νµ {νe interactions resultsin the difficulty of identifying tau signal events. The signal events cannot be identified
using a single variable in the standard Super-K reconstruction. The signatures ofhadronic tau decays are reflected in several different reconstructed variables. In
order to make use of multiple variables simultaneously in tau identification, a neuralnetwork is used in the analysis to look for extra pions from decays of heavy particles.
The neural network, configured with TMVA package (Hocker et al., 2007) has seven
input nodes, one hidden layer, and one output node.
6.3.1

Input Variables to the Neural Network

The seven input variables to the network are• The log base 10 of the total visible energy (Evis) of the event.
Evis is defined as the energy of a showering electron that produces the same96

Super-Kamiokande IV

Super-Kamiokande IV

Run 999999 Sub 1 Event 55

Run 999999 Sub 1 Event 27

16-03-10:18:48:25

16-04-13:05:43:20

Inner: 2701 hits, 18301 pe

Inner: 6059 hits, 29936 pe

Outer: 4 hits, 5 pe

Outer: 3 hits, 2 pe

Trigger: 0x07

Trigger: 0x07

D_wall: 588.8 cm

D_wall: 100.1 cm

Evis:

2.2 GeVEvis:

2.3 GeVmu-like, p = 1903.0 MeV/c

Charge(pe)Charge(pe)>26.7
23.3-26.7
20.2-23.3
17.3-20.2
14.7-17.3
12.2-14.7
10.0-12.2
8.0-10.0
6.2- 8.0
4.7- 6.2
3.3- 4.7
2.2- 3.3
1.3- 2.2
0.7- 1.3
0.2- 0.7
< 0.2

>26.7
23.3-26.7
20.2-23.3
17.3-20.2
14.7-17.3
12.2-14.7
10.0-12.2
8.0-10.0
6.2- 8.0
4.7- 6.2
3.3- 4.7
2.2- 3.3
1.3- 2.2
0.7- 1.3
0.2- 0.7
< 0.2
420

1100
1 mu-edecay

336

3 mu-edecays

880

252

660

168

440

84

220

0

0
0

500

1000

1500

2000

0

Times (ns)

500

1000

1500

2000

Times (ns)Figure 6.6: Simulation of a background (left) event with 2.2 GeV visible energy
and a CC ντ event (right) with 2.3 GeV visible energy in Super-K, which produces
multiple rings in the ID.

Figure 6.7: Simulation of a background event that has a deep-inelastic interaction,
which produces multiple rings in the detector.

amount of observed Cherenkov light. Due to the energy threshold of CC ντinteractions and the large mass of tau lepton, the signal events are expected
to have higher average visible energy than background events. The signal
distribution peaks around 4 GeV in Fig. 6.8. By contrast, the background
Evis spectrum falls with increasing Evis.
• The particle identification of the ring with maximum energy.
Tau leptons decay quickly to daughter particles after production through lep97

tonic decays and hadronic decays. Except for the leptonic decay to a muon,
most decay channels have at least one showering particle. A showering particle has a small value in the definition of particle identification likelihood,compared with a large value for a non-showering particle. The particle identification of the most energetic ring for signal has a distribution mostly in thenegative (showering) region, while the background has a broad distribution inboth negative and positive regions.
• The number of decay electron candidates in the event.
Naively, we expect more pion decay electrons in signal events from the tau
decays which are produced in hadronic tau decays. The variable does notdepend on ring reconstruction, so it is relatively independent of most othervariables.
• The maximum distance between the primary interaction point and any decay
electron from a pion or muon decay.
Energetic muons can travel a long distance in water. Therefore, CC νµ background involving a high energy neutrino is expected to have a large distance
between the primary interaction point and the decay electron from the muon.
In comparison, the pions from hadronic tau decay are expected to have smaller
momentum, resulting in a smaller value of the variable.
• The clustered sphericity of the event in the center of mass system.
Sphericity is a measure of how spherical an object is. A perfectly isotropic
object has sphericity 1, while a perfectly one-directional object has sphericity
0. By defining a sphericity tensor, sphericity can be deduced.

S

αβ

ř α β
i pi pi
“ ř
2
i pi

98

(6.4)

where α,β= 1, 2, 3 are the three cartesian directions of the momentum vectorsof the object, the sum is over the bins in geodesic space. Sphericity can beconstructed by finding the eigenvalues, λ1 ą λ2 ą λ3 , of the tensor.
3
S “ pλ2 ` λ3 q
2

(6.5)

The hadronic decay of heavy tau lepton is more isotropic than a typical νµor νe background. The spectrum of sphericity is centered near 1 for signal,
while the spectrum has an almost flat distribution between 0.1 and 0.8. The
background has a significant portion of events with large sphericity because of
the high percentage of DIS at high neutrino energies.
• The number of possible Cherenkov ring fragments.
A Hough transformation method is used to search for possible ring candidates
in Super-K (Abe et al., 2013). By Hough transformation, the observed photoelectrons are mapped to a circle centered on the hit PMT. By accumulating
the mapped circles, we can find peaks in the Hough space that are identified
as centers of possible ring candidates. We expect more rings in signal becauseof the multiple charged particles and pions in hadronic tau decay.
• The fraction of total number of photoelectrons in the event carried by the firstring.
This variable is the only one that has a different definition from that used in(Abe et al., 2013). A new variable is calculated with the charge in PMTs andthe reconstructed vertex and direction of the event as
ř
˝ qi
i ă48
rf rac “ θř
qi

(6.6)

where θi is the angle between the reconstructed direction of the first ring andthe direction of the reconstructed vertex to a PMT. The variable calculates
99the ratio of charge within 48˝ of the direction of first ring in the event. The
variable is expected to be smaller for the signal because the energy is carriedby multiple particles in the tau hadronic decay. The variable is useful in theselection of deep inelastic neutral-current background from the charged-current
background because energies are typically more equally divided into multiple
particles than the charged-current background.
Because the oscillation-induced tau events come from below, the downward sample is expected to have no tau events. Therefore, the data and MC are comparedin the downward-going sample to validate the agreement. Fig. 6.8 demonstrates theagreement between data and MC in the downward sample, along with the expected
tau signal, for the seven variables.
6.3.2

Training and Testing the Neural Network

500 years of MC simulations for signal and background are used in the analysis. The
MC samples are divided into three parts, one used to train the neural network, one
used to test the neural network performance and one used in the analysis of data. Tau
neutrinos from oscillations are expected to mostly come from below because of thepath lengths in the earth. Because of this, the oscillation weight described in Sec. 2
will encode the up-down asymmetry to the signal MC that only upward-going events
are selected. To select events based on event topology above and to avoid encoding
up-down biases into the neural network, the training sample is weighted by theenergy, as shown in Fig. 6.9, to cancel the up-down asymmetry in the oscillated MC.
In this way, the upward and downward samples are treated the same in the training
process, while the overall weight is still correct. The weights of the downward signal
simulation are not set to zero, thereby preserving signal MC statistics.
The neural network is trained to classify a tau-like event with a NN output close
to 1, and non-tau-like events have a NN output close to 0.
100

1000
3000

800
600

2000

400
1000
200
03

3.5

4

0
-300

4.5

Log of Visible Energy

-200

-100

0

ID of Maximum Energy Ring

2500

4000

2000

3000

1500

2000

1000
1000

500
00

2

4

6

00

8

Number of Decay Electrons

2000

4000

6000

Max Distance to Decay-e
1000

600

800
600

400

400
200
200
00

0.2

0.4

0.6

0.8

1

Clustered Sphericity

00

5

10

15

20

Ring Candidates

1500

1000

500

0

0.2

0.4

0.6

0.8

1

Fraction carried in First Ring

Figure 6.8: The input variables to the neural network for the downward MC (black
histogram) and data (black dots), along with tau MC (gray shaded histogram).

Table 6.1: MC sample sizes for each stage of the MLP. The same statistics are used
for SKI-IV.
SignalBackgroundtrainning„25,000
32559

testinganalysis
1500„6,600 (100 yr.)
1500 „82,000 (100 yr.)Fig. 6.10 shows the NN output for background and signal simulations. The spectrum of NN for signal is clearly separated from the background. The trained NN
is found to efficiently separate the CC tau signals from other background neutrino
events. When selecting tau-like events from the events passing preselection cuts
requiring the NN out greater than 0.5, 76% of signal events and only 28% of background remain. In the tau-like background, NC makes up 20% of the background. It
101

1
5
0.8

4

0.6

3

2

0.4

1

0.2

0
3

3.5

4

4.5

5

0
3

5.5
Log(evis)(a) Energy weights for backgroundin the training sample.

3.5

4

4.5

5

5.5
Log(evis)(b) Energy weights for tau signal in
the training sample.

Figure 6.9: Energy weights for the training sample.

5Signal (test sample)

Signal (training sample)

Background (test sample)

Background (training sample)

Kolmogorov-Smirnov test: signal (background) probability = 0.8 ( 0.82)

4
U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%(1/N) dN / dx

TMVA overtraining check for classifier: MLP

3

2

1

0

0

0.2

0.4

0.6

0.8

1

MLP response

Figure 6.10: Neural network output of training (filled histogram) and test (dots)
sample for background (red) and signal (blue) simulations.

is noted that multi-pion/DIS comprises 69% of the background. The neural network
is configured to identify hadronic tau decays which typically have multiple pions, so
the neural network also learns to idenfity multiple pion or DIS events. Table 6.2
summarizes the breakdown of the interaction modes in different samples.
By varying a cut on NN output, the efficiency of signal and background canbe changed. The background efficiency is anti-correlated with the signal efficiency,
102

Table 6.2: The breakdown of interaction modes of both background and expectedsignal for the SK-I period. By cutting the NN output at 0.5, each mode is separated
into tau-like and non-tau-like.
Interaction Mode
CC νe
CC νµ
CC ντ
NC

NN ă 0.5
818.1
1187.212.5
77.0

NNą0.5
428.4
216.6
39.4
161.6

All
1246.51403.851.9
238.5

which is shown in Fig. 6.11. The efficiency of signal and background by cutting NN
output is only plotted to assess the performance of the neural network in selecting
tau events. No actual cut is performed in the analysis itself.

Background rejection

Background rejection versus Signal efficiency
1
0.9
0.8
0.7
0.6
0.5
0.4
MVA Method:

0.3
0.2

MLP

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Signal efficiency

Figure 6.11: Efficiency of signal and background rejection by cutting the NN
output in SK-IV.

Atmospheric tau neutrinos are expected in Super-K due to the neutrino oscillations. However, the detection of tau neutrinos is challenging because charged-current
tau neutrino cross sections are suppressed by the relative large mass of the tau andthe short lifetime of tau. A neural network algorithm has been developed to identify
the events that have hadronic tau decays in Super-K, which shows good separationof charged-current tau neutrino events from other atmospheric neutrino events. The
algorithm will be used in the following analyses.
103

7
Search for Tau Neutrino Appearance

In this chapter, we describe the measurement of tau neutrino appearance. The data
are fitted against the PDFs built from MC simulations to measure a normalization
factor which characterizes the tau appearance signal.

7.1 Fit the Tau and Background Normalization
In order to extract maximum information from the sample, the analysis uses a 2dimensional unbinned likelihood fit. Using the neural network output and the reconstructed zenith angle of the events, two-dimensional particle distribution functions(PDFs) are built for background and tau signal. Fig. 7.1 is an example of the
2D distributions for oscillated signal on the left and background on the right. The
horizontal axis of the plots is the cosine of the reconstructed zenith angle that isdetermined by the energy-weighted sum of the ring directions in the event. The
vertical axis is the NN output, in which tau-like events have a value close to 1 andnon-tau-like events have a value close to 0. The signal events (left panel) are primarily tau-like and come from below (cospθq near -1), while the background eventsare more non-tau-like and come from all directions. The amount of signal and back104

ground events can be adjusted by varying the normalization of the distributions.

Figure 7.1: 2-dimensional probability distribution likelihood as a function of zenith
angle and neural network output for tau (left) and background (right) built with MC.

DIS background events typically produce multiple pions as in the hadronic tau
decays, thus are given high scores by the neural network tau identification. Therefore,
DIS background makes up a large component of the tau-like background. Fig. 7.2
shows the DIS and non-DIS component in the background simulation.
1500

Data
Signal Simulation
DIS Background
Non-DIS Background

1000

500

0

0

0.2

0.4

0.6

0.8

1

Neural Network Output

Figure 7.2: The breakdown of background simulation into DIS (orange) and NonDIS (red) as a function of neural network output.

The average neutrino energy in the DIS interactions in our sample is 14 GeV.
The DIS neutrino cross section in this region is not known to better than 10%.
105

In addition, the Bodek-Yang corrections (Bodek and Yang, 2002) implemented inour simulation of neutrino interactions tend to suppress DIS cross section at higher
energies by about 5%. Because of the large portion of DIS component and theuncertainties of the interaction, we simultaneously fit the tau normalization with
the DIS normalization with this uncertainty taken into account. The DIS error isintroduced into the fit as a 10% Gaussian error constraint. The fit of the data isdone to the sum of the signal and background PDFs varying the normalization,
Data “ α ˆ Signal ` β ˆ pBGN onDIS ` γ ˆ BGDIS q,(7.1)

where α, β and γ are the normalization of signal, non-DIS background and DIS
background respectively.
Separate PDFs are produced for each Super-K run period for both signal andbackground. The data sets are fit both individually for each run period and jointly.
Because of the low statistics of data events in the 2D distribution, we performed an
unbinned likelihood fit with RooFit (Verkerke and Kirkby, 2003). The likelihood in
the joint fit is constructed as

Lpα, β, γq “

I´IV
ź

ź

f pN N, cosθ; α, β, γq

SK run event i´pγ´1q2
1? e 2ˆ0.12
0.1 ˆ 2π

(7.2)

where f pN N, cosθ; α, β, γq is the normalized probability density function (PDF) constructed with the likelihood from simulation as shown in Fig. 7.1. The likelihood
is maximized by varying relative normalization of the three distributions, equivalent
to minimizing´ logrLpα, β, γqs “ ´ÿÿlogrf pN N, cosθ; α, β, γqs.

(7.3)

SKperiod event i

A sensitivity study for the joint fit of all Super-K run periods is performed with toy
MC generated with the signal and background PDFs, with α, β and γ set to be equal106to 1 in the study. 1,000 sets of toy MC were generated, and fit to extract the threeparameters. The distribution of fit tau normalization is shown in Fig. 7.3. The tau
normalization from the toy MC study has a Gaussian distribution centered at 1. The
tau normalizations in the fit of toy MC are fitted with a Gaussian function, and theσ of the Gaussian fit is 0.24.
Events / ( 0.118135 )A RooPlot of "Signal Fraction"
220
200
180
160
140
120
100
80
60
40
20
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2
2.2 2.4
Signal Fraction

Figure 7.3: Tau normalization fit results with 2D unbinned likelihood fit of toy
MC generated with α “ 1. The distribution has a σ of 0.24 with a Gaussian fit.

The data sets are fit individually and jointly for SK-I through SK-IV, and fit
results are summarized in Table 7.1, and are compared with the results from (Abeet al., 2013). The tau normalization in the joint fit is 1.37 with statistical uncertaintyof 0.23.
Table 7.1: Data fit results of tau normalization, background normalization and DIS
background normalization in each Super-K period.

SK-I
SK-II
SK-III
SK-IV
Total

α
1.25˘0.42
1.65˘0.54
1.99˘0.70
1.57˘0.35
1.37˘0.22

α in (Abe et al., 2013)
1.27˘0.49
1.47˘0.62
2.16˘0.78
1.42˘0.35

β
0.89
0.87
0.86
0.88
0.86

γ
1.11
1.18
1.09
1.11
1.18

The statistical uncertainty is 16% in the data fit, better than 24% of the fit
results of the toy MC study. The smaller uncertainty is due to the fitting of a tau
107normalization to a value larger than 1 in the data fit. Another set of toy MC isgenerated with α “ 1.37, β “ 1 and γ “ 1. The tau normalization is fit in the new
study, and the distribution of tau normalizations is shown in Fig. 7.4, which has astatistical error of 0.26. The statistical error in the new MC study is similar to thatin the data fit.
Events / ( 0.115963 )A RooPlot of "Signal Fraction"
200
180
160
140
120
100
80
60
40
20
0

0.5

1

1.5

2

2.5
Signal Fraction

Figure 7.4: Fit of tau normalization in toy MC study with tau normalizations set
to 1.37 in the toy MC generation. The distribution is centered at 1.37, and the fitof the distribution with a Gaussian function found a σ of 0.26.

To examine the results of the combined fit graphically, binned projections inzenith for both tau-like (NN output ą0.5) and non-tau-like (NN output ă 0.5)
events, and the projections in NN output for both upward-going (cospθq ą 0.1)
and downward-going (cospθq ă 0.1) events are shown in Fig. 7.5. The 2D PDFs ofsignal and background are projected into 1D and have been rescaled to the fitted
normalization values. Good agreements are seen in all distributions. As expected,
the tau neutrino signal shows up in the upward-going direction. The PDFs and datafrom all Super-K periods are combined in these plots.

108

Figure 7.5: Fit results showing the projections in NN output and zenith angle
distribution for tau-like (NN ą 0.5), upward-going (cospθq ą 0.1), non-tau-like (NNă 0.5) and downward-going (cospθq ă 0.1) events for both the two-dimensional PDFs
and data. The PDFs have been rescaled by the fit results. The PDFs and data sets
have been combined from SK-I through SK-IV. The fitted tau signal is shown ingray.

7.2 Systematic Errors
A subset of systematic uncertainties in atmospheric neutrinos is taken into account.
In total, there are 56 systematic uncertainties, 18 for the atmospheric neutrino flux,
18 for neutrino interactions, 15 for detector response and reconstruction and 5 oscillation parameter uncertainties. Out of the 56 uncertainties, 40 are common among
SK-I through SK-IV, while the other 16 are dependent on Super-K run period. A
brief description of the systematic uncertainties is presented below.
Atmospheric Neutrino Flux
1. Absolute normalization
The uncertainty of the absolute flux normalization is given by Honda flux
model (Honda et al., 2015). It is an energy-dependent error as shown in Fig.
109

7.6, where scaling factors for below and above 1 GeV independently scale theflux according to the function shown. Additionally, 5% normalization error on
multi-GeV events is considered in order to account for the differences above 10
GeV between the Honda flux model (Honda et al., 2015) and other flux models

Absolute normalization (%)

such as FLUKA (Battistoni et al., 2003) and Bartol (Gaisser, 2013).
35
30
25
20
15
10
5
0
10−1

1

10

102

103
104
Eν (GeV)Figure 7.6: Uncertainty on the absolute flux normalization as a function of neutrinoenergy.

2. Flux Ratios
The uncertainties on the various flux ratios are estimated by comparing the
Honda flux (Honda et al., 2015) to the FLUKA (Battistoni et al., 2003) andBartol (Gaisser, 2013) flux models. The flavor ratio (νµ ` ν̄µ )/(νe ` ν̄e ) in
FLUKA and Bartol model normalized by Honda flux is shown in Fig. 7.7. The
uncertainty is estimated to be 2% for Eν ă 1 GeV, 3% for 1 GeV ă Eν ă10
GeV, 5% for 10 GeVă Eν ă 30 GeV, then linearly increasing as a function oflog10 Eν to 30% at 1TeV.
The uncertainty on the anti-neutrino to neutrino ratio comes from π ` {π ´ andK ` {K ´ ratio in hadronic interaction of the flux calculation. For the ν̄e {νeratio, it is estimated to be 5% for Eν ă 10 GeV, 8% for 10 GeV ă Eν ă100
GeV, then linearly increasing as a function of log10 Eν to 30% at 1 TeV. For the110

(νµ+νµ)/(νe+νe) double ratio

1.4
1.3
1.2
1.1
1
0.9
0.8
0.7
0.6 −1101103Eν(GeV)

102

10

Figure 7.7: Flavor ratio (νµ ` ν̄µ )/(νe ` ν̄e ) in FLUKA model (solid) and Bartol
model (dashed) normalized by Honda flux model.

ν̄µ {νµ ratio, it is estimated to be 2% for Eν ă 1 GeV, 6% for 1 GeV ă Eν ă50
GeV, then linearly increasing as a function of log10 Eν to 60% at 1TeV. The
comparison of ν̄{ν ratios between Honda flux and other models is shown in

1.5
1.4
1.3
1.2
1.1
1
0.9
0.8
0.7
0.6
0.5 −110νµ/ νµ double ratioνe/ νe double ratio

Fig. 7.8.

1

10

102

103
Eν(GeV)

1.5
1.4
1.3
1.2
1.1
1
0.9
0.8
0.7
0.6
0.5 −1
10

1

10

102

103
Eν(GeV)

Figure 7.8: Anti-neutrino/neutrino ratio normalized by ratio in Honda flux modelfor ν e {νe (left) and ν µ {νµ (right) for FLUKA (solid) and Bartol (dashed) models.

The uncertainty on up/down ratio is due to the geomagnetic field. The uncertainty on horizontal/vertical ratio is mainly caused by the three-dimensional
calculation method in different flux models for neutrino energy below 3 GeV.
For neutrino energy above 3 GeV, it is caused by the difference in the calculated the K/π ratio in the hadronic interactions in the atmosphere. Bothuncertainties are estimated by comparing the zenith angle distributions using111

the three flux models.
3. K{π Ratio
Atmospheric neutrinos below 10 GeV are mainly produced by pion decays.
However, kaon contributions become important as the energy increases. The
SPY experiment measured the K/π production ratio with an accuracy of 3%
in proton interactions with beryllium for K and π momentum ranging from
7 GeV/c to 135 GeV/c (Ambrosini et al., 1998). Based on this measurement
and the correlation between hadron momentum and neutrino energy, the uncertainty on neutrino flux due to K/π ratio is estimated to be 5% for Eν below
100 GeV, linearly increasing to 20% at Eν “1 TeV, and 20% for Eν above 1
TeV.
4. Solar Activity
Primary cosmic rays are affected by the solar activity that periodically changes
in a 11-year cycle. The variation of neutrino flux caused by a ˘ 1 year variationin the solar cycle is assigned as an uncertainty on the flux, which is 20% for
SK-I, 50% for SK-II, 20% for SK-III, and 10% for SK-IV.
5. Neutrino Path Length
The production height of atmospheric neutrinos depends on the density of
the atmosphere, and the uncertainty affects the oscillation probability calculation. For upward-going events, the uncertainty produces a negligible effect.
For events near or above the horizon, it is treated as a systematic error. A
10% uncertainty on the atmosphere density structure is given by comparing
the US-standard 76 (NOAA and Force, 1976) and MSISE-90 (Hedin, 1991)
models for the atmosphere, and the resulting change in the production heights
of atmospheric neutrinos is propagated to the calculation of neutrino oscilla112

tions. Fig. 7.9 compares the neutrino path length as a function of zenith angle
calculated with the atmosphere density in Honda model and the atmosphere

path length (km)density compressed by 10%.
10

4

10

3

10

2

10

1

-0.2

0

0.2

0.4

0.6

0.8

1

Zenith angle (cosΘ)

Figure 7.9: Neutrino path length as a function of the zenith angle calculatedwith the atmosphere density in Honda model (solid) and the atmosphere density
compressed by 10% (dashed). Taken from (Lee, 2012)

6. Matter EffectAlthough the matter density structure of the Earth is well measured by (Dziewonski and Anderson, 1981), there is an uncertainty on the chemical compositionof the Earth core. The electron density that depends on the chemical composition of the Earth core affects the oscillation probability calculation. Therefore,
the uncertainty needs to be taken into account. The Earth core is assumed
to consist of mostly heavy elements such as iron, which has a electron density
about 6.8% less than that of light elements. The electron density in the Earthcore is estimated to be 6.8% smaller and is considered during the calculationof oscillation probabilities.
113

Neutrino Interactions
1. CCQE Cross Section
The cross sections of charged-current quasi-elastic scattering (CCQE) is calculated based on the Smith-Moniz CCQE model (Smith and Moniz, 1972) witha Relativistic Fermi Gas (RFG) nuclear model. The differences of the modelfrom Nieves model (Nieves et al., 2004) are assigned as the uncertainties of thetotal cross section as well as the cross section ratios ν̄{ν and pν̄µ ` νµ q{pν̄e ` νe q.
2. Axial Mass
The Super-K neutrino interaction MC has been generated at axial mass MA =1.2
GeV. A 10 % uncertainty on the axial mass is assumed in the analysis.
3. MECThe meson exchange current models (MEC) have large uncertainties, and theircross sections have not been explicitly measured. 100% uncertainty is assigned
to the total cross section of the MEC process.
4. Single Meson Production
The uncertainty on the single meson production cross section is estimated to
be 20% by comparing the model calculation with the experimental data ofνµ ` p Ñ µ´ ` p ` π ` measurement. Because of the lack of measurementof νµ ` p Ñ µ´ ` p ` π 0 , the uncertainty on its cross section is estimatedby comparing the base model in NEUT to the model from (Hernandez et al.,2008). The comparison of the theoretical cross sections is shown in the left
panel of Fig. 7.10. The double ratio of the ν̄{ν cross sections, shown in theright panel of Fig. 7.10, is also taken as an systematic uncertainty in the single
meson production. The ratio between the π 0 production to the charged pion
production is estimated to be 40% and is taken as the systematic error.
114

0CC ν/ ν double ratio

CC νµ cross section ratio

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0.5

1

1.5

2

2.5

3

3.5

4

4.5 5
Eν (GeV)

2.2
2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5 5
Eν (GeV)Figure 7.10: The Hernandez model normalized by the base model in NEUT forνµ p Ñ µ´ pπ ` (solid line), νµ n Ñ µ´ pπ 0 (dashed line) and νµ n Ñ µ´ pπ ` (dottedline) are shown in for the theoretical cross sections (left) and the double ratio forν̄{ν ratio (right). Figure is taken from (Lee, 2012).

The uncertainties on the parameters that characterize the Graczyk-Sobczyk
form factor (Graczyk and Sobczyk, 2008) were given by bubble chamber experiments data (Wilkinson et al., 2014).
5. Deep Inelastic Scattering Cross Section
The theoretical calculation of DIS cross section agrees with experimental measurements better than 5% for energies above 10 GeV, but the discrepancy
is larger in low energy region (ă10 GeV). A 10% uncertainty is assigned tothe cross section. Moreover, the application of Bodek-Yang correction tends
to suppress the DIS cross section by about 5%. An additional uncertainty istaken as the difference between the base model of GRV98 with the Bodek-Yang
correction with CKMT model (Capella et al., 1994).
Also, the uncertainty on the Q2 spectrum of the DIS interaction is given by
the difference between the GRV98 model with and without the Bodek-Yang
correction. Separate uncertainties are assigned to regions of W ă13 GeV/c2
and W ą13 GeV/c2 , where W is the invariant mass.
6. Coherent Pion Production
115

For CC νe and NC coherent pion production, 50% uncertainty is assignedby comparing the difference between the SciBooNE measurement (Nakajima,2011) and the Rein and Sehgal model. As SciBooNE had observed no evidenceof CC νµ coherent pion production, 100% uncertainty is assigned to its crosssection.
7. NC/CC Ratio
The uncertainty on NC/CC ratio is estimated to be 20%.
8. Nuclear Effect in 16 O Nucleus The uncertainty on pion final state interactions is estimated by comparing the NEUT simulation with external experiment
data (de Perio, 2011). The cross sections of these final state interactions arevaried with a weighing method while keeping the difference between the reweighted simulation from the experiment data less than 1σ, as shown in Fig.
7.11. The uncertainty on the nuclear effect in 16 O nucleus is given by the difference between the minimum weight and maximum weight on each simulated
event.
Detector Response and Reconstruction
1. FC Reduction
The uncertainty on FC reduction efficiency is estimated to be 0.2% for SK-I,
0.2% for SK-II, 0.8% for SK-III and 0.3% for SK-IV by comparing the distributions of each cut variable in data and MC. Background contamination in eachrun period is described in Section 5.1.
2. FC/PC separation
NHITAC, defined as the number of hits in the OD, is used to separate FC andPC events, and the uncertainty on the separation is estimated by the differenceof NHITAC distribution between the data and MC.
116

σAbsorption (mb)σQuasi-elastic (mb)300

qe1.0_ab1.1_cx1.0
qe0.6_ab0.7_cx0.5
qe0.6_ab0.7_cx1.6
qe0.7_ab1.6_cx0.4
qe0.7_ab1.6_cx1.6
qe1.3_ab0.6_cx0.6
qe1.3_ab0.7_cx1.6
qe1.5_ab1.5_cx0.4
qe1.6_ab1.6_cx1.6

200

200
150
100

100
50
0
0

200
400
600
π+ Initial Momentum (MeV/c)σSingle CX (mb)σReactive (mb)

0
0

400

200
400
600
π+ Initial Momentum (MeV/c)

100

50
200

0
0

200
400
600
π+ Initial Momentum (MeV/c)

0
0

200
400
600
π+ Initial Momentum (MeV/c)

Figure 7.11: Re-weighted simulation of π ` scattering on carbon target with NEUT
and experimental data in: quasi-elastic channel (top left), absorption channel (topright), total reaction channel (bottom left) and single charge exchange channel (bottom right). Taken from (Lee, 2012).

3. Fiducial Volume
The uncertainty on the fiducial volume is estimated to be 2%, corresponding
to the change in fiducial volume with a inward/outward shift of the vertex
position by „ 10 cm near the fiducial boundary.
4. Ring Counting
The ring separation error is considered by examining the differences in the
distributions of ring-counting likelihood between the data and MC.
5. Particle Identification
Particle identification errors are estimated similarly by comparing the distributions of e{µ likelihood ratio between the data and MC. The errors are estimatedseparately for single-ring and multi-ring events.

117

105

105
FC

PC

104

104

103

103

102

102

10

10

1
0

20

40

60

1
0

80

105

FC

20

FC

60

80

103

102

102

10

10

20

40

60

80

PC

104

PC

103

1
0

40

105
FC

104

PC

40

60

1
0

80

20

Figure 7.12: NHITAC distributions for FC and PC final samples for SK-I (topleft), SK-II (top right), SK-III (bottom left) and SK-IV (bottom right). Data areshown by dots with error bar, MC events assuming no oscillation are shown withdashed lines and MC with νµ Ñ ντ 2-flavor oscillation with (sin2 2θ, ∆m2 ) = (1.00,
2.5 ˆ 10´3 eV2 ) are shown with solid lines.

6. Energy Calibration
Energy scale errors are described in Section 3.6. The error on the absoluteenergy scale is estimated to be 1.1% for SK-I, 1.7% for SK-II, 2.7% for SKIII, and 2.3% for SK-IV. The up/down asymmetry of the energy scale is also
evaluated from the momentum of decay electrons with different directions and
0.3% uncertainty is assigned.
7. Decay Electron Tagging
The uncertainty on decay electron detection efficiency is estimated to be 1.5%
for SK-I to SK-III and 0.8% for SK-IV. In addition, π ` decay probability
uncertainty is taken into account, which ranges between 3% and 13% depending118

on the π ` momentum.
Oscillation Parameters
1. sin2 pθ23 q and ∆m223 are varied within 1σ limits of a combined SK-I+SK-II+SKIII atmospheric oscillation result assuming normal mass hierarchy (Abe et al.,2013). The ∆m232 is varied between 1.92ˆ10´3 and 2.22ˆ10´3 eV2 , sin2 2θ23 isvaried between 0.93 and 1.0.
2. The sin2 pθ13 q values are varied within our combined Daya Bay (An et al., 2012)
and RENO (Ahn et al., 2012) results of sin2 2θ13 “ 0.099 ˘ 0.014.
7.2.1

Individually Estimate Systematic Uncertainties

The systematic errors in the atmospheric neutrino analysis (Wendell et al., 2010) are
reevaluated, and 28 systematic errors are included in this method. The systematic
errors are updated regularly. The values of systematic errors that are common in four
Super-K run periods are summarized in Table 7.2. Table 7.3 shows the systematic
errors that are dependent on Super-K run periods.
Table 7.2: Systematic errors used in tau neutrino appearance search that are common
to all Super-K run periods. 1σ size of the systematic error is shown as a percentage.
Systematic error
Flux normalizationν̄e {νe ratioEν ą 1GeV

σp%q
7

1 GeV ă Eν ă 10GeVEν ą 10GeV

35

1 GeV ă Eν ă 10GeVEν ą 10GeV

58

1 GeV ă Eν ă 10GeVEν ą 10GeV

66ν̄µ {νµ ratioν¯µ {νe ratio

Up/down ratio
Multi-GeV

119

e-likeµ-like0.7
0.2

Multi-ring Sub-GeV
Multi-ring Multi-GeVµ-likee-likeµ-like

0.2
2.8
1.5

Multi-GeV

e-likeµ-likeµ-likee-likeµ-like

0.7
0.2
0.2
2.8
1.5
10
10
5
10
1.0
20
1.0
5
100
20
30
10
10
40
10
10

horizontal/vertical ratio

Multi-ring Sub-GeV
Multi-ring Multi-GeV
K{π ratio in flux
Neutrino path length
sample by sample normalization
MA in QE and single π
CCQE cross section
Single meson production cross section
DIS model difference
DIS cross section
Coherent π production
NC/CC ratio
Nuclear effect in pion spectrum
CCQE ν{ν̄ ratio for νe
CCQE cross section µ{e ratio
π 0 {π ˘ ratio
single pion production ν̄{νratio
Hadron simulation

FC Multi-GeV

The effect of an individual systematic error on expected signal and backgroundis calculated by reweighting every event in simulations after shifting the systematic
error value from the central value by +1σ/-1σ. The neural network is re-trained withthe reweighted MC, and the data are re-classified with the new network. The PDFs
are updated with distributions of the reweighted and reclassified MC. The unbinned
likelihood fit is repeated with the reprocessed data and new PDFs. The differenceof the refitted tau normalization with systematic error shift from the tau normalization with systematic errors at the central value is calculated as the uncertainty
from this systematic error term. By iterating the process for each systematic error
one at a time, the systematic uncertainties of tau normalization are estimated for120

Table 7.3: Systematic errors used in tau neutrino appearance search that are independent between Super-K run periods. 1σ size of the systematic error is shown as apercentage.
Systematic error
solar activity
FC reduction
Multi-GeV e-like BG
Multi-GeV mu-like BG
PID for single-ring e-like
PID for single-ring mu-like
PID for multi-ring e-like
PID for multi-ring mu-like
Energy calibration
up/down symmetry for energy calibration
fiducial volume

SK-I
20
0.2
0.1
0.1
0.2
0.2
6.5
2.9
1.1
0.6
2

SK-II
50
0.2
0.1
0.1
0.2
0.3
9.7
3.8
1.7
0.6
2

SK-III
20
0.8
0.4
0.2
0.3
0.3
4.5
2.7
2.7
1.3
2

SK-IV
10
0.3
0.1
0.8
0.4
0.4
3.3
1.6
2.3
0.3
2

each systematic error for both plus 1σ and minus 1σ. The central value of the tau
normalization is taken as the fit result without changing systematic errors, which
does not change in the estimation of systematic errors. The changes in fitted tau
normalization after shifting the systematic errors are summarized in Table 7.4.
Table 7.4: Change in fitted tau normalization after shifting systematic error by`1σ{ ´ 1σ.

Systematic Error
ν{ν̄ ratio for Eν ą 1GeVνe {ν̄e ratio for Eν ą 1GeV
νµ {ν̄µ ratio for Eν ą 1GeV
up/down ratio
horizontal/vertical ratio
K/π ratio
neutrino path length
absolute neutrino Eą1GeV
relative norm
axial mass
CCQE cross section
single meson
DIS model difference
coherent pi
NC/CC ratio121+1σ (%)
2.5
0.2
1.2
3.3
-4.8
1.9
0.1
7.7
4.7
4.8
-1.9
3.3
-7.1
1.0
2.5

-1σ (%)
-2.6
-0.1
-1.1
-3.3
4.8
-1.9
-0.1
-9.0
-5.1
-5.7
1.8
-3.8
6.2
-1.0
-2.9

nuclear pi spectrumCC ντ cross section
FC reduction
hadron simulationnon νe backgroundnon ν background
PID 1 ring
PID multi-ringenergy calibration
up/down energy calibration
solar activity
CCQE ν{ν̄ ratio
CCQE νµ {νe ratio
π{π ˘ ratioν̄{ν 1 pi ratio
ν{ν̄ ration 1GeV ă E ă 10GeVνe {ν̄e ration 1GeV ă E ă 10GeV

0
-25.0
0.3
0.01
0.01
0.1
0.1
2.4
2.0
2.3
1.0
0.8
-0.3
-2.8
-0.4
0.6
0.0

0
25.0
-0.3
-0.01
-0.01
-0.1
-0.1
-2.7
-1.9
-2.2
0.8
-0.8
0.3
2.7
0.5
-0.7
0.3

The systematic errors are divided into two sets. The first set includes the systematic errors used to quote the error on the expected number of tau events, which
does not affect the significance of the measured tau normalization. The tau crosssection error, which is 25% error, is the largest in this set. The error on the tau crosssection was estimated by comparing the charged-current tau neutrino cross sectionsin NEUT with several other theoretical models, especially the difference between
NEUT and cross section model by (Hagiwara et al., 2003). Detector biases on selection and fitting are also included in this set, but are small compared to the tau crosssection uncertainty.
The second set of errors would affect the observed signal but not the background,or the significance of the tau normalization in the fit. The five errors in this set
are all expressed as ratios: upward to downward neutrino flux, horizontal to verticalneutrino flux, kaon to pion originated neutrino flux, NC to CC cross section andthe upward to downward detector energy scale difference. The dominant error on
the signal is the horizontal/vertical ratio changing the fitted number of tau events of122

about ˘5%.
The systematic errors of the oscillation parameters are also included in the errors
which can change the measured results and significance. In this study, they arevaried within the 1σ limit of a combined SK-I+SK-II+SK-III atmospheric oscillation
analysis result assuming the normal hierarchy (Abe et al., 2013). The ∆m232 is varied
between 1.92 ˆ 10´3 and 2.22 ˆ 10´3 eV 2 , sin2 2θ23 is varied between 0.93 and 1.0.
The θ13 value is varied within the world average result of sin2 p2θ13 q “ 0.099 ˘ 0.14
in (Patrignani et al., 2016). Due to the small uncertainty of θ13 , the variation in θ13
around this central value results in less than 1% change in the fit result. However, theuse of non-zero θ13 value results in a 14% reduction of the fitted tau normalization
because three-flavor oscillations produce high energy upward-going electron neutrino
in the expected background thus reducing the needed number of tau neutrinos toexplain the signal region. Varying the value of δcp results in less than 1% change
in the fitted number of tau events. Therefore, the value of δcp is set to zero in thisanalysis, and the uncertainty is neglected.
Table 7.5: Summary of errors of oscillation parameters.
Systematic uncertainties of ντ events∆m232
sin2 p2θ23 q
sin2 p2θ13 q

+%
2.4
0.5

-%
-4.5
-4.6
-0.7

The systematic errors are summarized in Table 7.6. Combining the systematic errors on the observed tau events, the measured value of tau normalizationis 1.37˘0.23`0.13´0.11 . The expected number of tau events in the fiducial volume is 230.9with combined systematic errors of (`28.7 ´ 29.1%). After rescaling the MC by thefitting factors and corrections for efficiency, the observed number of tau events fit for`30.0
the entire running period is calculated to be 316.3˘53.1´25.4.

The observed number of tau events is converted to the significance level for reject123

Table 7.6: Summary of systematic errors for both the expected and measured numberof ντ events. The errors for each category have been added in quadrature.
Systematic uncertainties of ντ events
Super-K atmospheric ν oscillation errors
28 error term (expected events)
5 error terms (observed events)
Tau neutrino cross section (expected events)
Oscillation parameters (observed events)+%-%

14.0 14.9
7.0 7.1
25.0 25.0
6.5 2.4

ing the non-tau-appearance hypothesis. The measured signal normalization with anormalization of zero corresponds to the hypothesis of no-tau-appearance. A normalized asymmetric Gaussian function centered at 1.37 is constructed with the combined
error of measured tau normalization as shown in Fig. 7.13. A p value calculated as
the integral of the function below 0 is calculated to be 1.1 ˆ 10´8 , corresponding toa significance level of 5.1σ.

Figure 7.13: An asymmetric Gaussian function is constructed with the fitted tau
normalization and it uncertainty. The significance of tau appearance is calculated
by integrating the negative part of the function.

If the inverted hierarchy is chosen in calculating the oscillation probabilities instead of the normal hierarchy, the expected number of θ13 -induced upward-going
electron neutrino background is reduced, resulting in a higher fitted value (1.56) andcorresponding to a higher significance.
1247.2.2

Simultaneous Fitting of the Tau Normalization and Systematic Errors

In the above method, the tau normalization is fit with the systematic errors at thecentral value, and the systematic errors are estimated individually and added inquadrature. The correlation between signal and systematic errors and correlation
within systematic errors are thus neglected in the fit. To address the problem, a
new fit method is developed to simultaneously fit the tau normalization and thesystematic errors. This new method is more accurate and conservative. The fit
method is also based on RooFit. The systematic errors in Super-K atmospheric
neutrino analysis (Abe et al., 2013) are taken into account in the fit. PDFs for
signal and background are built as described in the above subsection for each run
period, but the DIS component is not built separately. PDFs for systematic errors
are built out of two-dimensional histograms prepared from the MC simulation andthe three-flavor oscillation software (Osc3++) (Wendell et al., 2010). Osc3++ has afunction to calculate the change of event distribution for a 1σ change in an individual
systematic error. Using the function, we calculate the change of the two-dimensional
event distribution of zenith angle and neural network output for each systematic
error. Assuming the systematic coefficients have a linear effect on the bin content in
the two-dimensional distribution, the PDF of each systematic error is calculated as
the difference of distributions between ˘σ.

P DFi “pP DFBG p`σi q ` P DFτ p`σi qq ´ pP DFBG p´σi q ` P DFτ p´σi qq
,
2

(7.4)

where P DFBGpτ q p˘σi q is the two-dimensional event distribution for signal (background) after shifting the systematic error by `1σ{ ´ 1σ. An example of twodimensional distribution for 1σ change is shown in Fig. 7.14. The data are fitted tothe summation of background PDF, signal PDF and systematic PDFs varying the125normalizations.
Data “ P DGBG ` α ˆ P DFtau `ÿi ˆ P DFi ,(7.5)

Neural Network Output

Instead of only fitting one signal normalization and two background normalizations,

1
0.8
0.6
0.4
0.2
0
-1

-0.5

0

0.5

1

Cos of the Zenith Angle

Figure 7.14: An example histogram of the PDF of K/π ratio uncertainty for SK-I
background. The vertical axis is the output of the NN; the horizontal axis is the
cosine of the event zenith distribution.

the new method includes coefficients of all systematic errors and normalization of
signal distribution into the fit. An unbinned likelihood fit of data is done to the sumof the signal, background and systematic uncertainty PDFs by varying the overall
normalizations. The likelihood is a function of tau normalizations and magnitudes
of systematic errors where the magnitude of every systematic error has a Gaussian
error constraint.

Lpα, 0 , 1 , . . . , m q “

I´IV
ź

ź

f pN N, cosθ; α, 0 , 1 , . . . , m q

SK run event imź12? e´i {2 , (7.6)
2π0where 0 , 1 , . . . , m are the magnitudes of the m systematic errors and P DFi is the
PDF for each systematic error built with the MC simulations.
All the systematic errors from atmospheric neutrino oscillations, except the uncertainty of CC tau neutrino cross section, are taken into account in the fit. The
tau cross section uncertainty is excluded because it affects the overall normalization126

of tau events, which acts the same as α in the fit. The uncertainties of oscillation
parameters can also change the measured results and significance, which are fitted
together. sin2 p2θ23 q and ∆m223 are varied between 1σ limits of Super-K result (Abeet al., 2013), and sin2 p2θ13 q is varied between 1σ limit of reactor neutrino measurements (Beringer et al., 2012).
Toy MC is generated with the PDFs, and the fit of tau normalization in the
toy MC is shown in Fig. 7.15. The fraction of toy MC with negative fitted tau

Events per Bin

normalization is calculated to be 3.6ˆ10´4 .
1.4

×103

1.2
1
0.8
0.6
0.4
0.2
0

0

0.5

1

1.5

2

Tau Normalization
Figure 7.15: Tau normalization in toy MC generated and fitted with all atmospheric neutrino analysis systematic errors, which is centered at 1.0 and has a pvalue of 3.6ˆ10´4 .

However, some of the systematic errors cause tiny or even no change in the event
distributions after shifting the systematic error, meaning that these systematic errors have a small effect in this fit. An example of two-dimensional distribution forsuch a systematic error is shown in Fig. 7.16, in which the fractional change in thetwo-dimensional event distribution is less than 0.1%. To exclude these negligible
systematic errors in the fit, the systematic errors for this analysis are selected by
cutting on the maximum fractional change in any bin of the two-dimensional histograms. A series of toy MC’s are generated by varying the cut value, and the tau
127normalization is fit in the toy MC’s. By requiring the significance level for rejecting
no-tau-appearance to be the same as the result without cutting on the systematic
errors in the toy MC, we can select the maximum cut value to remove unrelated systematic errors. Fig. 7.17 demonstrates the distribution of fitted tau normalizationsin the toy MC after selecting the systematic errors with the maximum cut value of
2.5%. The p value in the distribution is 4.4ˆ10´4 , corresponding to a significancelevel of 3.3σ. The significance level of rejecting no-tau-appearance in the toy MCchanges to 3.8σ (p=6.0ˆ10´5 ) when selecting the systematic error with a cut value
of 3%. Therefore, the systematic errors used in the analysis are selected by requiringthat the maximum fractional change in the bins of two-dimensional histograms not
be less than 2.5%. The fitted tau normalization in the toy MC generated and fittedwith selected systematic errors and uncertainties of oscillation parameters is shown

Neural Network Output

in Fig. 7.17.
×10-3
1

10

0.8

8

0.6

6

0.4
4
0.2
2
0
-1

-0.5

0

0.5

1

0

Cos of Zenith Angle
Figure 7.16: Histograms of the PDFs for uncertainty on absolute normalization ofneutrino flux below 1 GeV that causes small change in event distributions in SK-I.
The color demonstrates the change of events number in the bin after shifting thesystematic error by 1σ.

The fit is performed jointly with all data periods being fit at the same time.
Relative to the expectation of unity, the tau normalization is fitted to be 1.47˘ 0.32
(stat+syst) in the joint fit, corresponding to a 4.6σ rejection of the no-tau-appearance
128Events per Bin

800
600
400
200
0

0

1

2

Tau Normalization
Figure 7.17: Tau normalization in toy MC generated and fitted with selected
atmospheric neutrino analysis systematic errors, which is centered at 1.0 and has ap value of 4.4ˆ10´4 .

hypothesis. This result is less significant than that of the quadratic summation ofsystematic errors because of the more conservative method of simultaneous fitting
tau normalization and systematic errors. The fitted values and sizes of systematic
errors in the analysis are summarized in Table 7.7 for the systematic errors common
in Super-K run periods and Table 7.8 for systematic errors independent between
Super-K runs. Fig. 7.18a illustrates the distribution of the values of fitted systematic
errors. In order to evaluate the fitted systematic errors, we use a variable defined asEcor “Ef it ´ Einput,
sf it ´ sinput(7.7)

where Ecor is the value of systematic error after correction. Ef it (Einput ) is the value ofthe systematic error in the fit (input) and sf it (sinput ) is the 1σ size of the systematic
error in the fit (input). The corrected errors are shown in Fig. 7.18b.
The results of the final combined fit are examined graphically by plotting thebinned projections of the fitted results. Fig. 7.19 demonstrates the projection inzenith angle for both tau-like (NN output ą 0.5) and non-tau-like (NN output ă0.5)
events, along with the projections in NN output for both upward-going (cospθq ą 0.1)
129

20
15
10
5
0
-2

-1

0

1

2

Fit Error(a) Distribution of the fitted values of systematic error
from the fit.

12
10
8
6
4
2
0
-3

χ2 / ndf
12.59 / 8Prob0.1269Constant
7.054 ± 1.999
Mean-0.05919 ± 0.10178
Sigma
0.4883 ± 0.1140

-2

-1

0

1

2

3

Corrected Fit Error(b) Distribution of the variable defined in Equation
7.7. The distribution is expected to have a Gaussian
distribution with a width of 1.

Figure 7.18: Fitted values of systematic error in simultaneous fit of tau normalization and systematic errors.

130

Table 7.7: Fit results for common systematic errors in all Super-K run periods insimultaneous fit of tau normalization and systematic errors.
Systematic error
fit mean
CCQE ν{ν̄ ratio
6.8ˆ10´4
CCQE cross section ratio
0.10
DIS model difference-0.51
DIS cross section
0.20
K{π ratio
0.01
NC/CC ratio
0.22
Flux normalization (Eν ą 1GeV)-0.04
coherent π cross section-0.03
DIS q 2 dependence for low W
1.03
2
-0.27
DIS q dependence for high W
1π CAaa5 error-0.41
1π axial coupling-0.68
1π background scale factor-0.23
horizontal/vertical ratio
0.01
2∆m23 error-0.14
matter effect0.003
meson exchange on/off
0.05
ν{ν̄ ratio (Eν ą 10GeV )0.03
ν̄e {νe ratio (Eν ą 10GeV )0.008
neutrino path length-3.8ˆ10´4single-pion ν̄{ν ratio
0.002
ν̄µ {νµ ratio (Eν ą 10 GeV)
0.09
0˘
single pion π {π ratio
0.15
sample by sample normalization FC multi-GeV-0.07
sin2 pθ23 q
0.19
single meson cross section-0.17
sin2 p2θ13 q
0.02

fitted uncertainty
0.98
0.48
0.37
0.47
0.66
0.44
0.48
0.59
0.49
0.45
0.43
0.35
0.59
0.57
0.59
0.70
0.44
0.7
0.7
1.36
0.04
0.69
0.43
0.4
0.66
0.68
0.98

Table 7.8: Fit result of systematic errors that are independent between Super-K run
periods in the simultaneous fit tau normalization and systematic errors.
Systematic error

Mean
1-ring e-like background0.01
multi-ring e-like background 0.32
ring separation
0.01
multi-ring PID
0.18

SK-I
SK-II
SK-III
SK-IV
uncertainty Mean uncertainty Mean uncertainty Mean uncertainty
1.35
0.01
1.30
0.31
1.29
-0.04
1.34
0.83
0.66
0.88
0.33
1.07
0.26
0.70
1.19
-0.41
1.24
-0.08
1.35
-0.36
1.30
0.81
0.52
0.67
-0.26
1.1
0.1
0.98

131

and downward-going (cospθq ă0.1) events. In these plots, the signal PDFs have been
rescaled to the fitted normalization values, and PDFs of systematic errors for signal
and background have been rescaled by the fitted magnitudes of systematic errors
and added to signal and background respectively. The fitted tau signal is shown ingray. All distributions have good agreement between data and MC simulations. In

300

EventsEventsthese plots, the PDFs and data from all of the run periods are combined.
BG after fit

1000

Tau after fit
Data

800

200
600
400

100

200

Tau-like
0
-1

-0.5

0

0.5

0

1

Upward
0

0.2

0.4

0.6

0.8

1

0.8

1

NN output
EventsEvents

CosΘ

1200

600

1000
800

400

600
400

200

200

Non tau-like
0
-1

-0.5

0

0.5

0

1

CosΘ

Downward
0

0.2

0.4

0.6

NN output

Figure 7.19: Fit results showing projections in the NN output and zenith angle
distribution for tau-like (NNą0.5), upward-going [cospθq ą0.1], non-tau-like (NNă0.5) and downward-going [cospθq ă0.1] events for both the two-dimensional PDFs
and data. The PDFs and data sets have been combined from SK-I through SK-IV.
The fitted tau signal is shown in gray.

A significance level of 3.3σ is expected for the nominal expected signal, which issmaller than the measured result. The larger significance of the measurement is dueto the fact that more signal was measured than expected. Another set of toy MC isgenerated with 1.47 times the expected signal, and the fitted tau normalization in the
toy MC is shown in Fig. 7.20. The significance level has increased with more signal
132

than nominally expected. Because of the enormous amount of computing required
to generate toy MC for the p value calculation, the significance level is determinedby fitting the distribution of fitted tau normalizations in the toy MC by a Gaussianfunction. The fitted Gaussian function in centered at 1.47 with a width of 0.30,
corresponding to a significance level of 4.9σ, which is in relative good agreement

Events per Bin

with the measured result.
×103
1
0.8
0.6
0.4
0.2
0

0.5

1

1.5

2

2.5

Tau Normalization

Figure 7.20: Distribution of the fitted tau normalization in toy MC generated with
1.47 times of the nominal expected tau events.

The fit is repeated with inverted hierarchy when calculating the oscillation probabilities, resulting in a reduction in the expected number of θ13 -induced upward-going
electron neutrino. Under the assumption of inverted hierarchy, the fit results in ahigher fitted value (1.57˘0.31) and a correspondingly higher significance (5.0σ).

1338Charged-current Tau Neutrino Cross Section
Measurement

In this chapter, we describe the measurement of charged-current tau neutrino crosssections using atmospheric tau neutrino events in Super-K. The result is comparedto the measurement of the DONUT experiment.

8.1 Charged-current Tau Neutrino Cross Section Measurement
In the tau neutrino search, we developed a neural network algorithm to identify tau
events and measure the number of tau events in Super-K. By scaling the theoreticalcross section in the MC simulations to match the data, we can measure the inclusive
charged-current tau neutrino cross section in water.
σmeasured “ Sτ ˆ xσtheory y,(8.1)

where Sτ is the factor that is used to scale the theoretical cross section to match simulations and data. xσtheory y is the theoretical cross sections used in MC simulations.
We have already fitted a normalization factor to match data and MC simulation witha two-dimensional unbinned likelihood fit in the last section, which can be used to134scale the theoretical cross section. The scale factor can be also obtained by counting
the number of tau events in the data to be compared with the expected number oftau events in the MC simulations, which is basically an alternative one-dimensional
fit of data against MC. The two methods will be discussed and compared in thefollowing subsections.
8.1.1

Charged-current Tau Neutrino Cross Section Measurement Based on Twodimensional Fit of Tau Normalization

In the tau neutrino search, the normalization of tau events is fitted to be 1.47˘0.32
with systematic errors taken into account. This normalization matches the data
and MC simulations by scaling up the overall normalization of two-dimensional tau
distributions. Therefore, the normalization factor can also be converted to a chargecurrent tau neutrino cross section by scaling the theoretical cross sections in the MC
simulation with the factor,
σmeasured “ p1.47 ˘ 0.32q ˆ xσtheory y,(8.2)

Here, xσtheory y is a flux-averaged theoretical charged-current tau neutrino crosssections included in NEUT. To calculate it, we need to convolute the differential tau
neutrino cross sections as a function of neutrino energy with the energy spectrum of
atmospheric tau neutrinos from neutrino oscillations, which have been discussed inSec. 3 and Sec. 2 respectively. xσtheory y is calculated as follows:ř
xσtheory y “ş dΦpEν q

σpEν qdEν
dEν,
ş dΦpEν qdEνντ ,ν̄τdEνντ ,ν̄τ

ř

(8.3)

νq
where dΦpE
is the differential flux of tau neutrinos as a function of neutrino energy,dEν

σpEν q is the differential charged-current tau neutrino cross sections used in NEUT
νqsimulation. dΦpE
has been obtained from integrating the Honda flux over directionsdEν

in Sec. 2, and σpEν q is extracted from NEUT simulation. To calculate xσtheory y, we135

need to choose the energy range of the integral in Eq. 8.3, which is the neutrinoenergy range of interest in our measurement. Fig. 8.1a shows the energy spectrumof tau events in MC simulation in SK-I through SK-IV after weighting by oscillation probabilities and live time. In the MC simulations, tau events have an energy
range spanning from 3.5 GeV, the threshold energy of charged-current tau neutrino
interactions, to energies more than 60 GeV. However, it is noted that tau eventswith energy more than 60 GeV are rare. The number of tau events is integrated
for events with energy more than a specific value in simulation to determine themaximum energy of expected tau events. Fig. 8.1b shows the integral of the tau
neutrino events with energies more than the energy value in the horizontal axis. The
integral of events becomes less than 1 at the energy value of 67 GeV, meaning thatless than one tau event is expected with more than 67 GeV energy after selection.
The neutrino energy region of interest in the analysis is between 3.5 GeV and 67

200

Integral of Events

Events per Bin

GeV. The upper edge of the energy region is rounded to 70 GeV.

12
10

150

8

100

6
4

50

2
0
0

20

40

60

0
0

80

Neutrino Energy (GeV)(a) Energy spectrum of expected tau
events for 5326 days of live time.

20

40

60

80

Neutrino Energy (GeV)(b) Integral of events shown in theleft plot as a function of minimum
energy.

Figure 8.1: Energy spectrum of tau neutrino events after selection in SK-I through
SK-IV simulation after weighting by oscillation probabilities and live time (left), andintegral of number of events in the left plot with energy above the energy value inthe bin (right).

Given the theoretical cross section and tau neutrino energy spectrum by oscillating the electron and muon neutrinos to tau neutrinos, the flux-averaged theoretical136

charged-current tau neutrino cross section is calculated between 3.5 GeV and 70 GeV
to be 0.64ˆ10´38 cm2 . Fig. 8.2 shows the theoretical charged-current tau neutrino
cross sections in NEUT and tau neutrino fluxes as a function of neutrino energy by
oscillation electron and muon neutrinos in Honda model prediction to tau neutrinos.
The flux-averaged charged-current tau neutrino cross section is shown by the gray

10

10

1

1

10-1

NEUT CC ν τ xsec

Flux (GeV-1/m2/s)

Cross Section (10-38cm2)

dashed line.

10-1

NEUT CC ν τ xsec

10-2

Flux averaged NEUT xsec
ν τ fluxν τ flux

10-3

10-2
10-3

10-4
10-50

20

40

60

Neutrino Energy (GeV)

Figure 8.2: Theoretical charged-current tau neutrino cross sections (solid lines) as
a function of neutrino energy for ντ (blue) and ν̄τ (red), atmospheric tau neutrino
flux (dashed lines) from neutrino oscillations by oscillating the electron and muons
neutrinos in HONDA model prediction to tau neutrinos for ντ (blue) and ν̄τ (red)and flux-averaged tau neutrino cross section (gray dashed line).

Based on Eq. 8.2, the flux-averaged charged current tau neutrino cross section
is measured to be p0.94 ˘ 0.20q ˆ 10´38 cm2 . The measured cross section is shown
together with the theoretical cross sections and the MC simulations in Fig. 8.3. The
center and range of the energy axis (10`25´5 GeV) is calculated from the 90% range ofthe MC simulations after selection in Sec. 4, and the uncertainty on the cross section
axis (0.94˘0.20ˆ10´38 cm2 ) is from the uncertainty in the fit of tau normalization.
8.1.2

Comparing the Number of Measured Tau Events with Expectation

Instead of using the tau normalization in the two-dimensional fit, Sτ in Eq. 8.1 canalso be calculated from comparing the observed tau events with the expected tau
137

10

1

1

Events per Bin

Cross Section (10-38cm2)

10

10-1

10-1

10-2

10-2

10-30

20

40

60

10-3

Neutrino Energy (GeV)

Figure 8.3: Measured flux-averaged charged-current tau neutrino cross section(black), together with theoretical differential cross sections (ντ in blue and ν̄τ in red),flux-averaged theoretical cross section (dashed gray) and tau events after selectionin MC simulations (gray histogram).

events as follows
Sτ “

exp
obs
´ NBG
Ntot
Nτexp

(8.4)

exp
obs
where Ntot
is the total number of events in the sample after selection, NBG
is thenumber of background events in the MC simulations after selection, and Nτexp is thenumber of tau events in the MC simulations after selection.
The samples used in the tau appearance search have pre-selections, and the tau
neutrino concentration is only 1.8% in the MC after pre-selections. The statisticaland systematic uncertainties arising from the background are overwhelmingly larger
than Nτexp , resulting in large statistical and systematic errors in the measured Sτ .
Therefore, multiple cuts are applied on the sample used in the tau neutrino search
to enrich the tau signal.
• Cosine of zenith angle is less than -0.2. Tau neutrinos are expected mostly inthe upward-going sample, while background events show up in all directions.
The cut removes more than half of the background, but keeps most of the tau
signal events.

138• Neural network output larger than 0.54 (0.54 is the bin edge closest to 0.5 in
the two-dimensional PDF histograms). As shown in Fig. 6.10, tau events tend
to have neural network value close to 1 with background events having values
close to 0. The cut further increases the tau purity in the sample.
Table 8.1 summarizes the number of events and signal and background efficiency,
tau purity and number of events passing the cuts in the data. The cuts applied
are efficient in rejecting background, while retaining high signal efficiency. After thecuts, the expected tau concentration increases to 13.6%.
Table 8.1: Number of background, tau signal and data passing the event selection ineach Super-K period.
Tau

Tau
BG
BG
efficiency
efficiency
Multi-GeV, FCFV 193.5
10320.0
and CosThetaă-0.2 167.2
86.4%
3305.932.0%
and NNą0.54
118.5
61.2%
752.1
7.3%

Tau Data
purity1.8% 9884
4.8% 3262
13.6% 937

In order to obtain the information about change in expected events due to thesystematic errors, the cuts are applied on the two-dimensional event distributions forexp
the systematic errors built in Sec 6.3. The change in NBG
and Nτexp after shifting

one systematic error can be determined by counting the number of events in the
distributions after the cuts, which are summarized in Table 8.2 and 8.3.
Since the fitted results of the systematic errors are mostly different from thecentral value in Table 7.7 and 7.8, the expected number of signal and backgroundare calculated with the changes taken into account as follows
Nτexp “ Nτ0 `ÿ

∆Nτi ˆ σif it ,(8.5a)i∆NBG
ˆ σif it ,(8.5b)

i

exp
0
NBG
“ NBG
`

ÿ
i

139

exp
Table 8.2: Change in Nτexp and NBG
after shifting one systematic error by 1σ forcommon systematic errors in all Super-K run periods after selection cuts.

Systematic error∆Nτexp
CCQE ν{ν̄ ratio
0.13 (0.0)
CCQE cross section ratio
1.8 (1.2)
DIS model difference
105.0 (27.7)
DIS cross section
32.6 (15.5)
K{π ratio
6.4 (3.7)
NC/CC ratio
47.3 (20.3)
Flux normalization (Eν ą 1GeV)
57.4 (29.2)
coherent π cross section
2.6 (1.4)
2
DIS q dependence for low W
35.6 (22.9)
DIS q 2 dependence for high W
7.7 (2.8)
1π CAaa5 error
6.6 (3.2)
1π axial coupling
12.7 (5.5)
1π background scale factor
4.3 (1.9)
horizontal/vertical ratio
2.2 (0.8)
2.2 (2.9)∆m223 error
matter effect
0.5 (0.1)
meson exchange on/off
3.3 (2.5)
ν{ν̄ ratio (Eν ą 10GeV )
2.6 (1.3)ν̄e {νe ratio (Eν ą 10GeV )
2.7 (2.5)
neutrino path length
0.12 (0.0)
single-pion ν̄{ν ratio
1.9 (0.5)ν̄µ {νµ ratio (Eν ą 10 GeV)
3.5 (2.9)
single pion π 0 {π ˘ ratio
4.1 (2.7)
sample by sample normalization FC multi-GeV 37.5 (17.8)
sin2 pθ23 q
13.0 (5.6)
single meson cross section
3.7 (0.8)
2
sin p2θ13 q
1.7 (1.2)

exp∆NBG
0.12 (0.0)
0 (0.0)1.9 (0.1)
3.8 (2.4)
1.5 (0.9)
0.0 (0.0)
8.7 (5.1)
0.02 (0.0)
9.7 (7.0)
0.7 (0.1)
1.3 (0.5)
3.1 (1.1)
0.4 (0.10
0.9 (0.6)
1.7 (2.0)
0.1 (0.0)
0.0 (0.0)
2.4 (1.7)0 (0)0.1 (0.0)
0.2 (0.0)
1.4 (1.0)
1.3 (0.5)
5.3 (3.0)
0.1 (0.0)
0.6 (0.3)
0.03 (0.1)

exp
Table 8.3: Change in Nτexp and NBG
after shifting one systematic error by 1σ forsystematic errors that are independent between Super-K run periods.

SK-I
SK-II
SK-III
SK-IV
exp
exp
exp
exp
∆Nτexp∆NBG∆Nτexp∆NBG∆Nτexp∆NBG∆Nτexp∆NBG
1-ring e-like BG
0.4 (0.3) 0.0 (0.0) 0.5 (0.3) 0.0 (0.0) 0.3 (0.3) 0.0 (0.0) 0.8 (0.7) 0.0 (0.0)
multi-ring e-like BG 9.3 (3.8) 0.0 (0.0) 4.8 (3.5) 0.0 (0.0) 3.1 (1.3) 0.0 (0.0) 15.8 (6.6) 0.0 (0.0)
ring separation
3.1 (0.6) 0.5 (0.3) 1.1 (0.5) 0.1 (0.1) 0.4 (0.1) 0.1 (0.0) 1.2 (0.4) 0.3 (0.2)
multi-ring PID
9.3 (3.9) 1.5 (0.8) 7.5 (3.5) 1.0 (0.5) 2.4 (1.0) 0.4 (0.2) 8.6 (3.7) 1.4 (0.8)
Systematic error

140

exp
where Nτexp (NBG
) is the expected signal (background), Nτ0 is the signal (background)i)before adding the effect of systematic errors, as shown in Table 8.1, ∆Nτi (∆NBG

is the change of signal (background) after shifting the i’th systematic error by 1σ
and σif it is the fitted value of the i’th systematic error shown in Table 7.7 and 7.8.
The expected signal is shown in Fig. 8.4, along with data after subtracting expectedbackground from it. The number of expected tau events is 126.0 (81.3), while thenumber of background events is 748.0 after the cuts in Table 8.1 and adding the

Events per bin

effect of systematic errors.

Data-BG with stat&sys errors
Data-BG with stat error

80

Tau MC with stat&sys errors

60

Tau MC with stat error

40
20
0
5

10

15

20

25

30

Evis (GeV)
Figure 8.4: Visible energy of events passing the cuts in Table 8.1. Expectedbackground is subtracted from data, shown in dots, with statistical uncertainties fromdata (orange) and systematic uncertainties from background simulations. Expected
tau events are shown in the blue histogram, with statistical (blue) and systematic(shaded green) uncertainties.

A set of toy MC experiments is produced to estimate the statistical and systematic uncertainties of the measured Sτ . We assigned statistical and systematic
uncertainties to signal and background events separately from Table 8.2 and 8.3.
The fitted sizes of the systematic errors from Table 7.7 and 7.8 are introduced toconstrain the systematic errors in the toy MC experiments.
i “ Gausp0, σif it q,
141(8.6a)

Nτ1 “ Nτexp `

ÿ

exp
1
NBG
“ NBG
`

i ˆ ∆Nτi ,

i

Sτ “

ÿ

i
i ˆ ∆NBG,(8.6b)

i

GauspN obs ,a1
GauspN obs q ´ NBG,
Nτ1(8.6c)

where i is the value of systematic error randomly generated with a Gaussianfunction centered at 0 with a width of the fitted error size of that systematic error,
1
and Nτ1 (NBG
) is the number of expected signal (background) including the changes

caused by systematic errors. The number of events in the observation is randomized
with a Gaussian function. By forcing σi to be 0, we can estimate the statisticaluncertainty of Sτ . The distribution of Sτ in the toy MC experiments is shown in Fig.
8.5.
×103
MC study for stat error

30
MC study for total error

20

10

0
-1

0

1

2

3

4

xsec normalization
Figure 8.5: Sτ distribution in toy MC experiments generated with expected signal,
background and data after selections. The solid black line shows the central valueof Sτ as 1.5. The dashed line denotes the boundaries of 68% C.L. region, which is0.91-2.15.

By multiplying Sτ with the flux-averaged theoretical charged current tau neutrino`0.42
cross section, we obtain the measured cross section as 0.96´0.38
ˆ 10´38 cm2 , to becompared with the previous result of 0.94˘0.20ˆ10´38 cm2 . The different uncertainties of two methods are basically due to the difference between the one-dimensional
fit and two-dimensional fit.
142

10

1

1

Events per Bin

Cross Section (10-38cm2)

10

10-1

10-1

10-2

10-2

10-30

20

40

60

10-3

Neutrino Energy (GeV)

Figure 8.6: Flux-averaged tau neutrino cross section between 3 GeV and 70 GeV.

We expect that multidimensional fits have better sensitivity because more dataand information are used. The tau normalization is fitted with a single dimensional
fit, with only the zenith angle information used in the fit. The distribution of fitted
tau normalization in a set of toy MC is shown in Fig. 8.7, which has a width of
0.46 in a Gaussian fit. It is noted that the width is much larger than 0.24 in thetwo-dimensional fit.
Events / ( 0.14611 )A RooPlot of "Signal Fraction"
160
140
120
100
80
60
40
20
0
0

0.5

1

1.5

2

2.5
Signal Fraction

Figure 8.7: Tau normalization in toy MC generated and fitted with one-dimensional
PDFs.

Due to the better sensitivity of the full two-dimensional fit, we report 0.94˘0.20ˆ10´38 cm2
143

as the final result of the cross section study.

8.2 Comparisons with Previously Reported Results
8.2.1

Flux-averaged Tau Neutrino Cross Section

Because of the difficulties in tau neutrino production and detection, charged-current
tau neutrino cross sections have not been well measured. DONUT is the only experiment that previously measured the charged-current DIS tau neutrino cross section.
They used 9 observed charged-current tau neutrino events with an estimated background of 1.5 events. 800 GeV protons from Fermilab Tevatron were used to produce
neutrino beam by colliding with a beam dump, and tau neutrinos were produced via
decays of charmed mesons. The mean energies of the detected tau neutrino interactions were calculated to be 111 GeV (Kodama et al., 2001), in which deep inelastic
interactions are dominant. Assuming that the DIS charged-current tau neutrino
cross section has a linear dependence on neutrino energy, DONUT measured theslope of the energy dependence of the cross section including the kinematic effect of
tau lepton mass from the standard model calculation.
σpEq “ σconst ¨ E ¨ KpEq

(8.7)

where σpEq is the charged-current deep inelastic cross section as a function of neutrino energy, σconst is the slope that is to be measured, and KpEq is the kinematiceffect of the tau lepton mass from standard model calculation (shown in Fig. 8.8).
σconst is measured to be p0.39 ˘ 0.13 ˘ 0.13q ˆ 10´38 cm2 GeV´1 (Kodama et al., 2001)
in the final DONUT result. DONUT are incapable of distinguishing the charge ofthe τ lepton. Therefore, the measurement result is an average of ντ and ν̄τ crosssections assuming equal number of ντ and ν̄τ in the neutrino flux (Kodama et al.,2001). The ντ flux is larger than the ν̄τ flux, which are used to reweight the DONUT
cross sections.
144

Figure 8.8: Suppression of tau neutrino cross section relative to electron and muon
neutrinos due to the tau lepton mass. Taken from (Kodama et al., 2001).

In order to compare the Super-K measurement with the DONUT measurement,
the DONUT cross section result is reweighted by the Super-K tau neutrino flux. σpEq
is calculated for neutrino energy between 3.5 GeV and 70 GeV from σconst , KpEq andthe ratio of cross sections between νµ and ν̄µ in (Patrignani et al., 2016). σpEq and
xσy are shown in Fig. 8.9. The flux-weighted charged-current tau neutrino DIS crosssection is p0.37 ˘ 0.18q ˆ 10´38 cm2 , to be compared with p0.94 ˘ 0.20q ˆ 10´38 cm2 in
the Super-K measurement. Our result is larger than the reweighted DONUT crosssection because of the difference in neutrino energies. DONUT measured the DIS
cross section with a neutrino beam that has a much higher average energy than thatof tau neutrinos at Super-K. CCQE and CC coherent production is negligible in the
DONUT measurement. However, tau neutrinos at Super-K have a large componentof neutrinos below 10 GeV, where CCQE and CC coherent production cross sections
make a significant contribution to the final result.

145

Cross Section (10-38cm2)

10

1

10-1

NEUT CC ν τ xsec

NEUT CC ν τ xsec

DONUT CC ν τ xsec

DONUT CC ν τ xsec

Flux averaged NEUT xsec

SK measuread CC tau xsec

Flux weighted DONUT xsec

10-20

20

40

60

Neutrino Energy (GeV)

Figure 8.9: Comparison of the Super-K measured flux-averaged CC tau neutrino
cross section with DONUT result reweighted by the Super-K tau neutrino flux.
DONUT cross section digitized from (Kodama et al., 2008), and weighted by thetau neutrino from at Super-Kamiokande, ă σ ą“ p0.37 ˘ 0.18q ˆ 10´38 cm2 .

8.2.2

Comparison of σ{E with DONUT Result

Instead of measuring a flux-average cross section, we can also measure σ{E, as forthe DONUT result. The theoretical value of σ{E for ντ and ν̄τ can be calculated with
the NEUT model, which is a function of energy. Using the same method of scaling
MC to match data, we can calculate the value of σ{E at energies between 3.5 GeV
and 70 GeV in our measurement. Because ντ and ν̄τ have different cross sections(σντ {E “ 0.455ˆ10´38 cm2 GeV´1 , σντ̄ {E=0.226ˆ10´38 cm2 GeV´1 ) and fluxes in our
measurement, the measured σ{E has to be weighted by their fluxes. The measured
σ{E at 70 GeV is 0.51˘0.11ˆ10´38 cm2 GeV´1 .
In order to compare our result with DONUT result, the DONUT cross section as afunction of energy is calculated based on Equation 8.7. The kinematic function usedin the DONUT cross section is digitized, and the value at 70 GeV is taken to calculate
σ{E (E=70 GeV) to be 0.364ˆ10´38 cm2 GeV´1 for ντ and 0.164ˆ10´38 cm2 GeV´1 forν̄τ . The σ{E weighted by Super-K tau neutrino fluxes is 0.27˘0.13ˆ10´38 cm2 GeV´1 .
Our σ{E result is more consistent with the DONUT result than the flux-averaged146cross section because CC DIS is dominant at 70 GeV.

147

9
Three-flavor Oscillation Analysis with Tau
Identification

In spite of the dominant effect of νµ Ñ ντ oscillation in atmospheric neutrinos, subdominant effects are expected in atmospheric neutrino oscillations. The non-zero θ13from reactor neutrino experiment measurements implies that νµ oscillates into νe .
The oscillation has a considerable effect on atmospheric neutrino oscillations whenmatter effects are involved. A full three-flavor oscillation analysis could extract information about the neutrino mass hierarchy from Super-K atmospheric neutrino databy measuring the resonance of electron neutrinos in the multi-GeV region. Unfortunately, the tau neutrino appearance shows up in the same upward-going multi-GeV
sample, which looks similar to electron neutrino in Super-K reconstruction. The CC
tau events are an important background that is not well constrained due to the large
uncertainty on its cross section. The tau identification tools developed for the tau
neutrino appearance search have the capability to separate tau events from otherevents statistically, thus reducing this background. The tool has been incorporated
into the full three-flavor neutrino oscillation analysis, and is capable of improving
the sensitivity to the neutrino mass hierarchy by constraining the tau background.
1489.1 Three-flavor Atmospheric Neutrino Oscillation
9.1.1

Electron Neutrino Resonance Induced by Matter Effectνe oscillation is a sub-dominant effect in atmospheric neutrino oscillations. It includes
νe disappearance and νe appearance from νµ . The oscillation probabilities in vacuum
can be expressed as:
P pνe Ñ νe q “1 ´ sin2 p2θ13 q sin2 p

1.27∆m223 LqE

P pνµ Ø νe q “ sin2 pθ23 q sin2 p2θ13 q sin2 p

1.27∆m223 L
q
E

(9.1)

As discussed in Section 1.4.4, the matter effect plays an important role in νeoscillations. Given a non-zero θ13 , the matter effect results in a resonance in the νe(ν̄e ) appearance in the normal hierarchy (inverted hierarchy). Fig. 9.1 illustrates
the probabilities of νµ Ñ νe and ν̄µ Ñ ν̄e under the assumptions of the normal
hierarchy or inverted hierarchy. If the mass hierarchy turns out to be normal, thereis expected to be a resonance of νe appearance in the upward-going multi-GeV region
neutrinos but not in antineutrinos, and vice versa. Therefore, a definite measurementof the difference between electron neutrinos and antineutrinos in this region could
determine the mass hierarchy.
9.1.2

Oscillation Analysis Framework

In order to measure the oscillation parameters including the mass hierarchy, SuperK performs a full three-flavor oscillation analysis (Wendell et al., 2010; Abe et al.,2011b). The analysis framework is designed to compare the data against the MCfor different oscillation models. A particular oscillation model specifies a value foreach oscillation parameter, and the MC is oscillated based on that model. A model’s
compatibility with the data is defined as a goodness of fit statistic between the two.
A set of oscillation models forms a grid in the parameter space, and the analysis
149

0.6
0.5

0.5

Cosine Zenith Angle

Cosine Zenith Angle

1

1
0.5
0.5

0.4

0

0.3

0.4
0

0.3

0.2

0.2

-0.5

-0.5

0.1

0.1
-1

1

0

10

-1

1

0.6
0.5
0.4
0

Cosine Zenith Angle

Cosine Zenith Angle

Energy [GeV]

1
0.5

0.3
0.2

-0.5

0

10

Energy [GeV]

1
0.5
0.5

0.4
0.3

0

0.2
-0.5
0.1

0.1
-1

1

0

10

-1

Energy [GeV]

1

0

10
Energy [GeV]

Figure 9.1: Probabilities of νµ Ñ νe (left) and ν̄µ Ñ νe under the assumptionof the normal hierarchy (top) and inverted hierarchy (bottom). It is noted that
resonance in the upward multi-GeV region shows up in different neutrino types under
the assumptions of the normal or inverted hierarchy. Oscillation probabilities arecalculated assuming the oscillation parameters: ∆m212 “ 7.6 ˆ 10´5 eV 2 , ∆m232 “
2.4 ˆ 10´3 eV 2 , sin2 θ12 “ 0.3, sin2 θ23 “ 0.5, sin2 θ13 “ 0.0238 and δCP “ 0.

performs a search over the grid for the best agreement with the data.
All data samples including FC, PC and UPMU are used in the analysis. In order
to best extract the information from the data, the MC and data are divided into
subsamples. The details of the event separation is described in (Lee, 2012). The
subsamples are further binned in energy and cosine of reconstructed zenith angle ofthe event, as shown in Fig. 9.2.
Super-K performs a full three-flavor neutrino oscillation analysis with the atmospheric neutrino sample, in which the sub-leading effects are taken into account. In
150

Upmu
shower
FC
FC
FC
FC
FC
(1x10 bins)
1ring
Multi-ring Multi-ring Multi-ring
1ring
1ring
PC
PC
Upmu
µ-likeνe-likeanti-νe-likeνe-likeanti-νe-likeµ-likestop
throughstop(4x10 bins) (4x10 bins) (2x10 bins) (3x10 bins) (3x10 bins) (4x10 bins)(2x10 bins) (4x10 bins) (3x10 bins)

Log10 P (MeV/c)

5.0

Upmu
non-shower(1x10 bins)

4.0
3.7
3.4
3.0

3.2
3.0
2.8
2.6
2.4

FC
Multi-ring
Unclassified(4x10 bins)

FC Sub-GeV
1-ring e-like

FC Sub-GeV
1-ring π0-like

0-decay
1-decay(5x10 bins) (5x1 bins)(5x1 bins)

FC Sub-GeV
2-ring π0-like

FC Sub-GeV 1-ring µ-like
0-decay
1-decay
2-decay(5x10 bins) (5x10 bins) (5x1 bins)
3.2
3.0
2.8
2.6
2.4
2.2

2.0

3.2

2-decay(5x1 bins)

2.8
2.6
2.4
2.2
2.0

Figure 9.2: Definition of the binning scheme in the three-flavor oscillation analysis.
410 bins in FC, 60 bins in PC, 50 bins in UPMU. White boxes indicate that momentum bin is uniformly divided into 10 zenith angle bins between -1 and 1 for FC andPC, between -1 and 0 for UPMU. Shaded boxes indicate that the momentum bin is
not further divided into zenith angle bins.
order to examine the compatibility between the data and the oscillation parameters,
data are compared the MC expectation using a pull χ2 method based on Poisson
probability distribution. The binned likelihood is defined as

LpN

exp

,N

obs

q“

obs
exp
n
ź
e´Ni N exp Ni

i

Niobs !

i“1

,

(9.2)

where Niobs (Niexp ) is the number of observed (expected) events in the ith bin. χ2 isdefined as the log likelihood ratio
χ2 “ ´2 lnnÿLpN exp , N obs q
Niobs
exp
obs
obs
“
2
pN
´
N
`
N
ln
exp q.
ii
i
LpN obs , N obs q
N
i
i“1

(9.3)

The effect of systematic errors are included into the bin with the pull method (Wendell et al., 2010) by replacing Niexp with
Niexp Ñ Niexp p1 `

m
ÿ
j

151

fij j q,(9.4)

where j is a variation corresponding to the j th systematic uncertainty, fij is thefractional change of the event rate in the ith bin due to a change of 1σ in the jthsystematic uncertainty, the sum is over all the systematic uncertainties considered in
the oscillation analysis (Wendell et al., 2010). For each systematic uncertainty, thequantity j {σj is assumed to follow Gaussian distribution, where σj is the 1σ valueof the jth systematic uncertainty. Therefore, an additional term is added to the χ2
to constrain the size of systematic uncertainties in the fit. The final form of χ2 isdefined as
2

χ “2

n
ÿ

rpNiexp p1 `i“1

ÿfij j q ´ Niobs q ` Niobs ln

j

ÿ j
Niobs
ř
s
`
p q2 . (9.5)
Niexp p1 ` j fijj q
σjj

The fit varies the systematic uncertainties to minimize χ2 for a given MC expectation.
In the analysis, θ23 , ∆m223 , and δcp are fitted simultaneously, while θ13 is constrained to the recent reactor neutrino measurements. The best-fit information forthis fit is summarized in Table. 1. One-dimensional ∆χ2 distributions for theoscillation parameters are shown in Fig. 5. Normal hierarchy (NH) is preferred,χ2N H ´ χ2IH “ ´4.3, by Super-K atmospheric neutrino data. Toy Monte Carlo samples were generated with the assumption of either normal or inverted hierarchy totest the significance of the preference. The probability of having ∆χ2 smaller than
-4.3 is 0.45 in the NH assumption, and is between 0.007 (sin2 pθ23 q=0.4) and 0.031(sin2 pθ23 q=0.6) in the IH assumption. The best fit information of this fit is also
summarized in Table. 1.
Table 9.1: Summary of best-fit information for fits assuming the normal or inverted
hierarchy with the Super-K neutrino data only (517 d.o.f.).
Fit
SK NH
SK IH

∆2min δCP sin2 θ23
571.7 4.19
0.59
576.0 4.19
0.58

152

∆m223,13 reV 2 s
2.5ˆ10´3
2.5ˆ10´3

20

15

15

15

10

∆ χ2

20

∆ χ2

∆ χ2

20

10

99%

10

99%

5

99%

5

5

95%

95%

95%

90%

90%

90%

68%

68%

0
0.001

0.002

0.003

0.004

| ∆ m232 | , | ∆ m213 | eV2

0.005

0
0.2

68%

0.4

0.6

0.8

sin2 θ23

0
0

2

4

6

δcp

Figure 9.3: Constraints on neutrino oscillation parameters from the SKI-IV atmospheric neutrino data. The orange line denotes the result in the IH assumption,
which has been offset from the NH result, shown in blue, by the difference of their∆χ2 .

9.2 Constraining the Tau Background with Tau Identification
In the standard Super-K reconstruction and analysis, tau events are not separated
from other event categories. Due to the dominant hadronic decay channels, the tau
events are mostly categorized into e-like samples in the reconstruction and analysis
samples. Moreover, the tau appearance shows up in the same region as the electron
neutrino appearance induced by θ13 and matter effect, which becomes an importantbackground in the three-flavor neutrino oscillation analysis. Fig. 9.4 shows the
MC expectations under the assumption of the normal or inverted hierarchy andthe tau events in four e-like samples in the three-flavor oscillation analysis. It isnoted that the tau appearance shows up in the same region as the electron neutrino
appearance resonance, and the two have comparable magnitude. In addition, thecharged-current tau neutrino cross sections are known not better than 25%, leading
to a large uncertainty in the tau background in the analysis.
We have shown the capability of identifying the tau events with a neural network
algorithm in Section 6.3. It is possible to use the tau identification tools in the threeflavor oscillation analysis. In order to use tau identification, the multi-GeV eventsin the MC simulations and data are processed with the trained neural network, in
153

Multi-GeV e-like νeEventsEventsMulti-GeV e-like νeNormalInverted

100

CC ντ

NormalInverted
CC ντ

300

200
50
100

0
-1

-0.5

0

0.5

0
-1

1

-0.5

0

0.5

cos( θ )Multi-Ring e-like νeEventsEventsMulti-Ring e-like νeNormalInverted150

1
cos( θ )150

NormalInverted

CC ντ

CC ντ

100
100

50
50

0
-1

-0.5

0

0.5

0
-1

1
cos( θ )

-0.5

0

0.5

1
cos( θ )EventsMulti-Ring Other
300

NormalInverted
CC ντ

200

100

0
-1

-0.5

0

0.5

1
cos( θ )Figure 9.4: Zenith angle distribution of five samples that include electron appearance signal in the Super-K three-flavor neutrino oscillation analysis. The MC
expected under the assumption of the normal (inverted) hierarchy is shown in blue
(orange) histogram, without taking tau events into account. The tau events areshown in the green filled histogram.

154

which each event is given the neural network output of the tau appearance search.
In the modified three-flavor oscillation analysis, the existing bins for the samples
that are expected to have electron neutrino appearance, such as multi-GeV/multiring e-like νe {ν̄e samples, and the multi-ring unclassified sample, are further divided
into two bins based on the tau neutral network output. As shown in Fig. 9.5, one
bin in the multi-GeV e-like νe sample is divided into two bins at the neural network
output value of 0.25. The bin with NN output less than 0.25 has almost no tau

Number of Events

events, while the other bin has almost all the tau events.

8
6
4
2
0
0

0.2

0.4

0.6

0.8

1

NN Output
Figure 9.5: Distribution of neural network output for events with 3.0 ă log Evis ă
3.4 and ă ´1 cos θ ă-0.8 in the multi-GeV e-like νe sample, corresponding to one
bin the standard three-flavor oscillation analysis. The simulations are scaled to thelive time of 5326 days in SK-I through SK-IV. The νe and νµ events are shown in theblack histogram, while the tau events are shown in red filled histogram on top of theother events. The single bin is divided into two bins at the neural network equal to
0.25.

A sensitivity study was performed to study the effect of using the tau identification in the oscillation analysis. The study was carried out with SK-I, SK-II, SK-III,and SK-IV MC scaled to the data live time. As in the standard three-flavor oscillation analysis, ∆m232 , sin2 θ23 and δCP are free parameters and fitted in the study.
Fig. 9.6 shows the ∆χ2 as a function of sin2 θ23 between the normal and inverted
hierarchies under the assumption of true normal or inverted hierarchy. The incor155

poration of tau identification in the analysis on average improves the sensitivity tomass hierarchy by „ 10%.

Figure 9.6: ∆χ2 distribution as a function of sin2 θ23 for standard binning (red) andwith additional binning in neural network output (yellow). Under the assumption oftrue (inverted) normal hierarchy, the ∆χ2 between normal and inverted hierarchies
is shown in solid (dashed) lines.

As shown in this section, we have developed a framework to use the tau identification to constrain tau background in the oscillation analysis. The use of tau
identification improves the sensitivity to mass hierarchy by „ 10%. The method is
new being incorporated into the Super-K three-flavor oscillation analysis.

156

10
Conclusion

Super-Kamiokande has collected atmospheric neutrino data for more than 20 years.
This dissertation describes the result of a search for atmospheric tau neutrino appearance using 5,326 days of atmospheric neutrino data. Under the assumption of
the normal mass hierarchy, Super-K measured the tau normalization at 1.47 ˘ 0.32
(stat+sys), relative to the expectation of unity for nominal oscillation parameters
and assumed charged-current tau neutrino cross section. This result excludes the
hypothesis of no tau appearance with a significance of 4.6σ. The result is direct
proof of three-flavor neutrino oscillation.
Super-K measured 338.1˘72.7 charged-current tau neutrino interactions fully
contained in Super-Kamiokande. This charged-current tau neutrino interaction sample is larger than that of the DONUT experiment and the OPERA experiment, whichmeasured 9 and 5 charged-current tau neutrino interactions respectively. Using thislarge sample, the charged-current tau neutrino cross section is measured by scaling the theoretical cross section in the Super-K simulation to match the measured
charged-current tau neutrino interactions. A flux-averaged cross section is measured
to be p0.94 ˘ 0.20q ˆ 10´38 cm2 for neutrino energy between 3.5 GeV and 70 GeV,
157

compared with the flux-averaged theoretical cross section of 0.64 ˆ 10´38 cm2 . The
measured cross section is consistent with the Standard Model prediction within 1.5σ, and is only the second reported measurement of this cross section.
An event-by-event tau identification is achieved using tools developed for the
search for atmospheric tau neutrino appearance. The tau identification has been
implemented in the Super-K three-flavor oscillation analysis. A Monte Carlo study
has shown that the sensitivity to the mass hierarchy is improved by „10% with thetau identification implemented in the analysis.

158

Appendix A
Ratio of Charge in the Most Energetic Ring

In Ref. (Abe et al., 2013), the fraction of the charge carried by the most energetic
ring is defined as
rf rac˚ “

RT OT0
Evis(A.1)

where RT OT0 is the momentum of the most energetic ring in the reconstruction.
However, the data and MC do not agree well in this variable in SK-II. The MC on
average has a larger value than data, as shown in Fig. A.1. The ratio of the meanvalue of the variable in data and MC is 0.95. To correct the discrepancy, (Abe et al.,2013) introduced a normalization factor of 0.95 for SK-II MC to reduce the variableof each event.
In this analysis, a newly-defined variable is used to replace the problematic variable. The new definition is shown in Eq. 6.6. Instead of using the event reconstruction result to calculate the ratio, the PMT charges are taken as input to thecalculation. The new variable has good agreement between data and MC in all fourSK run periods, as shown in Fig. A.2.

159

Super Kamiokande I 1489.2 days

Super Kamiokande II 798.6 days

400

220

350

200
180

300

160
140

250

120

200

100
150

80

100

60
40

50
0
0

20
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

001

rfrac (FC multi-GeV)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

rfrac (FC multi-GeV)

Figure A.1: The fraction of energy carried by the most energetic ring for data
(dots) and MC (histogram) in SK-I (left) and SK-II (right). SK-III and SK-IV have
similar data and MC agreement as in SK-I.

Figure A.2: Ratio of charge carried by most energetic ring calculated with PMT
charge in each SK period.

160

Appendix B
Neural Networks and Their Application in Particle
Physics

An artificial neural network (ANN) is an information processing method inspiredby biological systems. It is composed of a large number of highly inter-connected
processing elements (neurons) working together to solve a problem. An ANN is typically organized in layers of inter-connected neurons. The external input is presented
to the network via input layer, which communicates to one or more hidden layers
where the processing is done via a system of weighted connections. The hidden layer
is then linked to an output layer. The ANN is capable of deriving meaning fromcomplicated or imprecise data, which is advantageous over conventional computing.
The weights of the connections in the network are adjusted in the learning processaccording to the external input examples, a process known as training. Neural networks are widely used in particle physics. In this analysis, a multilayer perceptron ischosen to identify tau events. The neural network used is configured with the TMVA
package (Hocker et al., 2007). Fig. B.1 demonstrates the structure of a multi-layer
perceptron with one hidden layer.

161

Input Layer

Hidden Layer

Output Layer

x1
1
x2
x3
x4

...

Outputx5x6
10x7
Figure B.1: Structure of a multilayer perceptron with seven input nodes in the
input layer, one output node in the output layer, and ten nodes in the hidden layer.

162

Appendix C
Partially Contained Reduction

C.1 Partially Contained Reduction
Partially contained events are different from fully contained events in the OD activities. A dedicated background reduction for the PC events is applied. During thereconstruction of the detector before SK-III, a opaque segmentation was installedin the OD to divide it into three regions - top, side barrel and bottom. The PC
reduction is substantially adjusted for SK-III onwards. There are also five stages for
PC reduction.
Here PC reduction in SK-III and SK-IV is summarized. Description of PC reduction for SK-I and SK-II could be found in Ref. (Lee, 2012).
C.1.1

PC1

Events that satisfy the following criteria are selected:
1. PEtot ě 1000 p.e.s (500 p.e.s for SK-II). PC events used in the atmospheric
neutrino analysis are required to have vertex at least 2 m from the ID wall,corresponding to a minimum momentum of 500 MeV/c. The cut is applied to163conservatively reject events with energy lower than approximately 310 MeV/c.
2. NHITAtop ă 10 or NHITAbottom ă 10, where NHITAtop (NHITAbottom ) is thenumber of OD hits at the top (bottom) of the OD. The rejects through-going
muons from top to bottom.
3. NHITAendcap ă 25 or NHITAside ă 70, where NHITAside is the number of OD
hits on the side barrel of the OD, and “endcap” refers to either top or bottom.
This rejects corner clipper muons that produce light in both the endcap andside barrel of the OD.
4. ODRmean ă 2100 cm, if OD hitsă20 in a 500 nsec time window. where
ODRmean is the average distance between all hit pairs that is defined as:

ODRmean “

1

N
´1
ÿ

Npair i“1ÝÝ|Ñxi ´ Ñ
xj |.

(C.1)

Following the first reduction, the event rate is around 2ˆ104 events per day.
C.1.2

PC2

PC2 is targeted to remove muon backgrounds based on a cluster algorithm. Events
satisfying the following criteria are selected:
1. NHITAendcap ă 20 or NHITAendcap ă MAX(NHITAside ), where MAX(NHITAside )is defined as epα´βˆN HIT Aside q with α=5.8, β=0.023 if NHITAside ă 75, α=4.675,
β=0.008 otherwise.
2. Nouter2ď10, where Nouter2 is the number of hits in the second highest charge
OD cluster.
After PC2, the event rate is reduced to around 8,000 events per day.

164

C.1.3PC3PC3 is is designed to reject flasher events following the same logic as that in the FC3
reduction. Events satisfying the following criteria are removed:
1 NHIT100 ą14 (9 if the number of ID hits ă 800).
C.1.4

PC4

PC4 reduction is based on a muon fitter “muboy” which classifies the remaining
cosmic muon background events into stopping muon, through-going muon, multiple muon or corner clipper muon. The six selection criteria based on the muboy
information are defined below:
1. Muboy angle ă 90˝ , where Muboy angle is the angle between the muboy fitted
direction and the vector between the pf it vertex and the center of the highest
charge OD cluster. The cut is targeted to remove stopping muon events.
2. Muboy dotprodą-0.8, where muboy dotprod is similar to muboy angle but instead of using the center of the highest charge OD cluster, the earliest saturated
ID PMT is used. This is designed to remove both stopping and through-goingcosmic ray muons.
3. Muboy track lengthă1750 cm, where muboy length is the length of the muontrack fitted by muboy.
4. Muboy goodnessă0.52. This is the goodness of fit given by muboy. A good fit
indicates a cosmic muon.
5. Muboy cornerě300 cm, where muboy corner is the distance between the muboy
entry point and the corners of the tank.

165

6. Muboy ehit8mě300 cm, where ehit8m is the number of OD hits located within
8 m of the entrance point from muboy fit in a fixed 5 nsec time window.
If an event is classified as a through-going muon by muboy, then it has to passfour of the first five cuts. If an event is classified as a stopping muon by muboy,
then it has to pass four out of the first five cuts including the muboy dotprod cut.
Additionally, it has to pass either the muboy ehit8m cut or have a muboy goodness
ă0.5. Otherwise, the event only needs to pass two of the first five cuts.
Finally, a cut on the total ID charge is applied to remove low-energy events:• PEtot ă 2900 p.e.
After the PC4 reduction, the event rate is „300 events/day.
C.1.5

PC5

PC5 consists of a series of cuts dedicated for different types of background sources.
The cuts are classified into hard cuts and soft cuts. Hard cuts are cuts that an eventis required to pass in order to be kept, including a through-going muon cut, stopping
muon cut, cable hole muon cut and corner clipper muon cut. The events are allowed
to fail the soft cuts once before being rejected. Soft cuts include through-going muoncut, stopping muon cut, corner clipper muon cut and a decay muon cut. A detailed
description of the cuts can be found in Ref. (Lee, 2012).
A corner clipper based on clusters in OD segmented regions is applied on eventswith visible energy less than 1 GeV, aiming to remove corner clipper muons. Insteadof clustering the OD PMTs together, the cluster algorithm is applied on the PMTs
on top, bottom and side walls separately. The largest charged cluster on endcap
and side wall are selected. In addition, the ID PMTs are clustered and the largest
charged cluster is used. The maximum distance between either of the three clusters (Distsegclust ) and the smaller charge of the two OD clusters (NHITAsegclust ) are
166

calculated. The events satisfying the following criterion are rejected:• Distsigclust ă1500 and NHITAsegclust ą30.
Fig. C.1 shows the scatter plot of Distsegclust and NHITAsegclust in MC and data.
4 ×10

500

3

400
300

2
200
1

100

0
0

4 ×10

3

20

40

60

80

clusters distance (cm)clusters distance (cm)30Number of hits in smaller cluster

7
6

3

5
4

2
3
2

1

1
0
0

20

40

60

80

0

Number of hits in smaller cluster

Figure C.1: Scatter plot of Distsegclust and NHITAsegclust in MC (left) anddata (right) before the last corner clipper cut in PC5. The data events withNHITAsegclust ą30 and Evis ă 1500 MeV are mostly background events tagged in
the eye-scan.

Following the PC5 reduction, the event rate is „1.2 events/day.
C.1.6

PC Reduction Summary

The final dataset of PC events are selected by the following cuts:
1. Reconstructed vertex is located inside the fiducial volume (at least 2 m fromthe ID wall),
2. NHITACě16 (10 for SK-I),
3.0 Visible energy ą 350 MeV.
The final sample is also eye-scanned to tag the background events (but not to remove
them from the final sample). Table C.1 shows the estimated background contamination in different run periods.

167

Table C.1: Estimated contamination by cosmic ray muons for PC.
SK-I SK-II SK-III SK-IV
cosmic ray BG 0.2% 0.7%
1.8%
0.6%

168

Bibliography
Aad, G. et al. (2012), “Observation of a new particle in the search for the Standard
Model Higgs boson with the ATLAS detector at the LHC,” Phys. Lett., B716,
1–29.
Aartsen, M. G. et al. (2016), “Neutrino oscillation studies with IceCube-DeepCore,”
Nucl. Phys., B908, 161–177.
Abdurashitov, J. N. et al. (1999), “Measurement of the solar neutrino capture ratewith gallium metal,” Phys. Rev., C60, 055801.
Abe, K. et al. (2006), “A Measurement of atmospheric neutrino flux consistent with
tau neutrino appearance,” Phys. Rev. Lett., 97, 171801.
Abe, K. et al. (2011a), “Indication of Electron Neutrino Appearance from an
Accelerator-produced Off-axis Muon Neutrino Beam,” Phys. Rev. Lett., 107,041801.
Abe, K. et al. (2011b), “Search for Differences in Oscillation Parameters for Atmospheric Neutrinos and Antineutrinos at Super-Kamiokande,” Phys. Rev. Lett., 107,241801.
Abe, K. et al. (2012a), “First Muon-Neutrino Disappearance Study with an Off-Axis
Beam,” Phys. Rev., D85, 031103.
Abe, K. et al. (2013), “Evidence for the Appearance of Atmospheric Tau Neutrinos
in Super-Kamiokande,” Phys. Rev. Lett., 110, 181802.
Abe, K. et al. (2014), “Calibration of the Super-Kamiokande Detector,” Nucl. Instrum. Meth., A737, 253–272.
Abe, K. et al. (2017), “First combined analysis of neutrino and antineutrino oscillations at T2K,” .
Abe, Y. et al. (2012b), “Indication for the disappearance of reactor electron antineutrinos in the Double Chooz experiment,” Phys. Rev. Lett., 108, 131801.
Achard, P. et al. (2004), “Measurement of the atmospheric muon spectrum from
20-GeV to 3000-GeV,” Phys. Lett., B598, 15–32.
169

Adamson, P. et al. (2011), “Improved search for muon-neutrino to electron-neutrino
oscillations in MINOS,” Phys. Rev. Lett., 107, 181802.
Adamson, P. et al. (2016), “First measurement of muon-neutrino disappearance in
NOvA,” Phys. Rev., D93, 051104.
Agafonova, N. et al. (2015), “Discovery of τ Neutrino Appearance in the CNGS
Neutrino Beam with the OPERA Experiment,” Phys. Rev. Lett., 115, 121802.
Aguilar-Arevalo, A., Bazarko, A., Brice, S., Brown, B., Bugel, L., Cao, J., Coney, L.,
Conrad, J., Cox, D., Curioni, A., et al. (2008), “Measurement of muon neutrino
quasielastic scattering on carbon,” Physical review letters, 100, 032301.
Aguilar-Arevalo, A. A. et al. (2013), “First measurement of the muon antineutrino
double-differential charged-current quasielastic cross section,” Phys. Rev., D88,
032001.
Ahmad, Q. R. et al. (2002), “Direct evidence for neutrino flavor transformation fromneutral current interactions in the Sudbury Neutrino Observatory,” Phys. Rev.
Lett., 89, 011301.
Ahn, J. K. et al. (2012), “Observation of Reactor Electron Antineutrino Disappearance in the RENO Experiment,” Phys. Rev. Lett., 108, 191802.
Ahn, M. H. et al. (2003), “Indications of neutrino oscillation in a 250 km long baseline
experiment,” Phys. Rev. Lett., 90, 041801.
Ahrens, L. A. et al. (1986), “Precise Determination of Sin2 θw From Measurementsof the Differential Cross-sections for νµ p Ñ νµ p and ν̄µ p Ñ ν̄µ p,” Phys. Rev. Lett.,
56, 1107, [Erratum: Phys. Rev. Lett.56,1883(1986)].
Alcaraz, J. et al. (2000), “Cosmic protons,” Phys. Lett., B490, 27–35.
Allison, W. W. M. et al. (1999), “The Atmospheric neutrino flavor ratio from a 3.9
fiducial kiloton year exposure of Soudan-2,” Phys. Lett., B449, 137–144.
Ambrosini, G. et al. (1998), “K / pi production ratios from 450-GeV/c protons onberyllium,” Phys. Lett., B420, 225.
An, F. P. et al. (2012), “Observation of electron-antineutrino disappearance at Daya
Bay,” Phys. Rev. Lett., 108, 171803.
Araki, T. et al. (2005), “Measurement of neutrino oscillation with KamLAND: Evidence of spectral distortion,” Phys. Rev. Lett., 94, 081801.
Arnison, G. et al. (1983), “Experimental Observation ?
of Isolated Large Transverse
Energy Electrons with Associated Missing Energy at s = 540-GeV,” Phys. Lett.,
122B, 103–116, [,611(1983)].
170

Bahcall, J. N., Calaprice, F., McDonald, A. B., and Totsuka, Y. (1996), “Solar
neutrino experiments: The next generation,” Physics Today, 49, 30–36.
Banner, M. et al. (1983), “Observation of Single Isolated Electrons of High Transverse Momentum in Events with Missing Transverse Energy at the CERN anti-p
p Collider,” Phys. Lett., 122B, 476–485.
Barger, V. D., Whisnant, K., Pakvasa, S., and Phillips, R. J. N. (1980), “Matter
Effects on Three-Neutrino Oscillations,” Phys. Rev., D22, 2718.
Barish, S. J. et al. (1977), “Study of Neutrino Interactions in Hydrogen and Deuterium. 1. Description of the Experiment and Study of the Reaction Neutrinoµ`dÑ µ´ +p+ps ,” Phys. Rev., D16, 3103.
Barlag, S. et al. (1982), “Charged Hadron Multiplicities in High-energy Anti-muon
Neutrino n and Anti-muon Neutrino p Interactions,” Z. Phys., C11, 283, [Erratum:
Z. Phys.C14,281(1982)].
Battistoni, G., Ferrari, A., Montaruli, T., and Sala, P. R. (2003), “The FLUKA
atmospheric neutrino flux calculation,” Astropart. Phys., 19, 269–290, [Erratum:
Astropart. Phys.19,291(2003)].
Becker-Szendy, R. et al. (1992), “The Electron-neutrino and muon-neutrino content
of the atmospheric flux,” Phys. Rev., D46, 3720–3724.
Becker-Szendy, R. et al. (1993), “IMB-3: A Large water Cherenkov detector for
nucleon decay and neutrino interactions,” Nucl. Instrum. Meth., A324, 363–382.
Belikov, S. V. et al. (1985), “Quasielastic Neutrino and Anti-neutrinos Scattering:
Total Cross-sections, Axial Vector Form-factor,” Z. Phys., A320, 625.
Benvenuti, A. et al. (1974), “Observation of Muonless Neutrino Induced Inelastic
Interactions,” Phys. Rev. Lett., 32, 800–803.
Berger, C. and Sehgal, L. M. (2007), “Lepton mass effects in single pion production
by neutrinos,” Phys. Rev., D76, 113004.
Beringer, J. et al. (2012), “Review of Particle Physics (RPP),” Phys. Rev., D86,
010001.
Bienlein, J. K. et al. (1964), “Spark chamber study of high-energy neutrino interactions,” Phys. Lett., 13, 80–86.
Bodek, A. and Yang, U. K. (2002), “Modeling deep inelastic cross-sections in the
few GeV region,” Nucl. Phys. Proc. Suppl., 112, 70–76, [,70(2002)].

171

Bonetti, S., Carnesecchi, G., Cavalli, D., Negri, P., Pullia, A., Rollier, M., Romano,F., and Schira, R. (1977), “Study of Quasielastic Reactions of Neutrino and antineutrino in Gargamelle,” Nuovo Cim., A38, 260–270.
Brun, R., Bruyant, F., Maire, M., McPherson, A. C., and Zanarini, P. (1987),“GEANT3,” .
Brunner, J. et al. (1990), “Quasielastic Nucleon and Hyperon Production by Neutrinos and Anti-neutrinos With Energies Below 30-GeV,” Z. Phys., C45, 551.
Capella, A., Kaidalov, A., Merino, C., and Tran Thanh Van, J. (1994), “Structure
functions and low x physics,” Phys. Lett., B337, 358–366.
Chatrchyan, S. et al. (2012), “Observation of a new boson at a mass of 125 GeV with
the CMS experiment at the LHC,” Phys. Lett., B716, 30–61.
Cowan, C. L., Reines, F., Harrison, F. B., Kruse, H. W., and McGuire, A. D. (1956),“Detection of the Free Neutrino: a Confirmation,” Science, 124, 103–104.
Danby, G., Gaillard, J.-M., Goulianos, K., Lederman, L. M., Mistry, N., Schwartz,
M., and Steinberger, J. (1962), “Observation of High-Energy Neutrino Reactions
and the Existence of Two Kinds of Neutrinos,” Phys. Rev. Lett., 9, 36–44.
Davidson, N., Nanava, G., Przedzinski, T., Richter-Was, E., and Was, Z. (2012),“Universal Interface of TAUOLA Technical and Physics Documentation,” Comput.
Phys. Commun., 183, 821–843.
Davis, R. (1994), “A review of the Homestake solar neutrino experiment,” Prog.
Part. Nucl. Phys., 32, 13–32.
Davis, Jr., R. (1955), “Attempt to detect the antineutrinos from a nuclear reactor
by the Cl37(anti-nu,e-) A37 reaction,” Phys. Rev., 97, 766–769.
de Perio, P. (2011), “NEUT pion FSI,” AIP Conf. Proc., 1405, 223–228.
Derrick, M. et al. (1978), “Properties of the Hadronic System Resulting from antiMuon-neutrino p Interactions,” Phys. Rev., D17, 1.
Dimopoulos, S. and Susskind, L. (1978), “On the Baryon Number of the Universe,”
Phys. Rev., D18, 4500–4509.
Dziewonski, A. M. and Anderson, D. L. (1981), “Preliminary reference Earth model,”Physics of the earth and planetary interiors, 25, 297–356.
Fermi, E. (1934), “Trends to a Theory of beta Radiation. (In Italian),” Nuovo Cim.,
11, 1–19, [,535(1934)].

172

Feynman, R. P. and Gell-Mann, M. (1958), “Theory of Fermi interaction,” Phys.
Rev., 109, 193–198.
Fukuda, S. et al. (2001), “Solar B-8 and hep neutrino measurements from 1258 daysof Super-Kamiokande data,” Phys. Rev. Lett., 86, 5651–5655.
Fukuda, Y. et al. (1994), “Atmospheric muon-neutrino / electron-neutrino ratio in
the multiGeV energy range,” Phys. Lett., B335, 237–245.
Fukuda, Y. et al. (1996), “Solar neutrino data covering solar cycle 22,” Phys. Rev.
Lett., 77, 1683–1686.
Fukuda, Y. et al. (1998), “Evidence for oscillation of atmospheric neutrinos,” Phys.
Rev. Lett., 81, 1562–1567.
Fukuda, Y. et al. (2003), “The Super-Kamiokande detector,” Nucl. Instrum. Meth.,
A501, 418–462.
Gaisser, T. K. (2013), “Atmospheric leptons,” EPJ Web Conf., 52, 09004.
Gamow, G. and Teller, E. (1937), “Some Generalizations of the beta Transformation
Theory,” Phys. Rev., 51, 289–289.
Glashow, S. L. (1961), “Partial Symmetries of Weak Interactions,” Nucl. Phys., 22,579–588.
Gluck, M., Reya, E., and Vogt, A. (1998), “Dynamical parton distributions revisited,” Eur. Phys. J., C5, 461–470.
Goldhaber, M., Grodzins, L., and Sunyar, A. W. (1958), “Helicity of Neutrinos,”
Phys. Rev., 109, 1015–1017.
Graczyk, K. M. and Sobczyk, J. T. (2008), “Form Factors in the Quark Resonance
Model,” Phys. Rev., D77, 053001, [Erratum: Phys. Rev.D79,079903(2009)].
Gran, R., Jeon, E., Aliu, E., Andringa, S., Aoki, S., Argyriades, J., Asakura, K.,
Ashie, R., Berghaus, F., Berns, H., et al. (2006), “Measurement of the quasielastic
axial vector mass in neutrino interactions on oxygen,” Physical Review D, 74,052002.
Hagiwara, K., Mawatari, K., and Yokoya, H. (2003), “Tau polarization in tau
neutrino nucleon scattering,” Nucl. Phys., B668, 364–384, [Erratum: Nucl.
Phys.B701,405(2004)].
Haino, S. et al. (2004), “Measurements of primary and atmospheric cosmic - rayspectra with the BESS-TeV spectrometer,” Phys. Lett., B594, 35–46.

173Hampel, W. et al. (1999), “GALLEX solar neutrino observations: Results for
GALLEX IV,” Phys. Lett., B447, 127–133.
Hasert, F. J. et al. (1973), “Observation of Neutrino Like Interactions Without Muon
Or Electron in the Gargamelle Neutrino Experiment,” Phys. Lett., 46B, 138–140.
Hayato, Y. (2009), “A neutrino interaction simulation program library NEUT,” Acta
Phys. Polon., B40, 2477–2489.
Hedin, A. E. (1991), “Extension of the MSIS thermosphere model into the middle
and lower atmosphere,” Journal of Geophysical Research: Space Physics, 96, 1159–
1172.
Hernandez, E., Nieves, J., and Valverde, M. (2008), “Neutrino induced weak pion
production off the nucleon,” Mod. Phys. Lett., A23, 2317–2320.
Hocker, A. et al. (2007), “TMVA - Toolkit for Multivariate Data Analysis,” PoS,ACAT, 040.
Honda, M., Kajita, T., Kasahara, K., and Midorikawa, S. (2004), “A New calculationof the atmospheric neutrino flux in a 3-dimensional scheme,” Phys. Rev., D70,
043008.
Honda, M., Sajjad Athar, M., Kajita, T., Kasahara, K., and Midorikawa, S. (2015),“Atmospheric neutrino flux calculation using the NRLMSISE-00 atmosphericmodel,” Phys. Rev., D92, 023004.
Ishihara, C. (2010), “Full three flavor oscillation analysis of atmospheric neutrino
data observed in Super-Kamiokande,” Ph.D. thesis, Tokyo U.
Jeong, Y. S. and Reno, M. H. (2010), “Tau neutrino and antineutrino cross sections,”
Phys. Rev., D82, 033010.
Kajita, T. (2006), “Discovery of neutrino oscillations,” Reports on Progress inPhysics, 69, 1607.
King, S. F. and Luhn, C. (2013), “Neutrino Mass and Mixing with Discrete Symmetry,” Rept. Prog. Phys., 76, 056201.
Kodama, K. et al. (2001), “Observation of tau neutrino interactions,” Phys. Lett.,
B504, 218–224.
Kodama, K. et al. (2008), “Final tau-neutrino results from the DONuT experiment,”
Phys. Rev., D78, 052002.
Kume, H., Sawaki, S., Ito, M., Arisaka, K., Kajita, T., Nishimura, A., and Suzuki,
A. (1983), “20-INCH DIAMETER PHOTOMULTIPLIER,” Nucl. Instrum. Meth.,
205, 443–449.
174

Lee, K.-P. (2012), “Study of the neutrino mass hierarchy with the atmospheric neutrino data observed in Super-Kamiokande,” Ph.D. thesis, University of Tokyo.
Lee, T. D. and Yang, C.-N. (1956), “Question of Parity Conservation in Weak Interactions,” Phys. Rev., 104, 254–258.
Levy, J.-M. (2009), “Cross-section and polarization of neutrino-produced tau’s madesimple,” J. Phys., G36, 055002.
Llewellyn Smith, C. H. (1972), “Neutrino Reactions at Accelerator Energies,” Phys.
Rept., 3, 261–379.
Mead, C. A. (1958), “Quantum Theory of the Refractive Index,” Phys. Rev., 110,
359–369.
Michael, D. G. et al. (2006), “Observation of muon neutrino disappearance with the
MINOS detectors and the NuMI neutrino beam,” Phys. Rev. Lett., 97, 191801.
Mikheev, S. P. and Smirnov, A. Yu. (1985), “Resonance Amplification of Oscillationsin Matter and Spectroscopy of Solar Neutrinos,” Sov. J. Nucl. Phys., 42, 913–917,[Yad. Fiz.42,1441(1985)].
Motoki, M. et al. (2003), “Precise measurements of atmospheric muon fluxes with
the BESS spectrometer,” Astropart. Phys., 19, 113–126.
Nakajima, Y. (2011), “A Measurement of Neutrino Charged Current Interactions anda Search for Muon Neutrino Disappearance with the Fermilab Booster Neutrino
Beam,” Ph.D. thesis, Kyoto U.
Nieves, J., Amaro, J. E., and Valverde, M. (2004), “Inclusive quasielastic chargedcurrent neutrino-nucleus reactions,” Physical Review C, 70, 055503.
Nieves, J., Ruiz Simo, I., and Vicente Vacas, M. J. (2011), “Inclusive Charged–
Current Neutrino–Nucleus Reactions,” Phys. Rev., C83, 045501.
Niita, K., Sato, T., Nakashima, H., Iwase, H., Nose, H., and Sihver, L. (2006),“PHITS: A particle and heavy ion transport code system,” Radiat. Meas., 41,
1080–1090.
Nishino, H. (2009), “Search for Nucleon Decay into Charged Antilepton plus Meson
in Super-Kamiokande,” Ph.D. thesis, Tokyo U.
Nishino, H., Awai, K., Hayato, Y., Nakayama, S., Okumura, K., Shiozawa, M.,
Takeda, A., Ishikawa, K., Minegishi, A., and Arai, Y. (2009), “High-speed chargeto-time converter ASIC for the Super-Kamiokande detector,” Nucl. Instrum.
Meth., A610, 710–717.

175

NOAA, U. and Force, U. A. (1976), “US standard atmosphere, 1976,” Tech. rep.,
NOAA-S/T.
Olive, K. A. et al. (2014), “Review of Particle Physics,” Chin. Phys., C38, 090001.
Patrignani, C. et al. (2016), “Review of Particle Physics,” Chin. Phys., C40, 100001.
Pauli, W. (1930), “Pauli’s Letter-Dear Radioactive Ladies and Gentlemen,”Tübingen, December, 4.
Pohl, M. et al. (1978), “Study of Neutrino Proton Elastic Scattering in the
GARGAMELLE Freon Experiment,” Phys. Lett., B72, 489–492.
Rein, D. and Sehgal, L. M. (1981), “Neutrino Excitation of Baryon Resonances andSingle Pion Production,” Annals Phys., 133, 79–153.
Rein, D. and Sehgal, L. M. (1983), “Coherent pi0 Production in Neutrino Reactions,”
Nucl. Phys., B223, 29–44.
Rein, D. and Sehgal, L. M. (2007), “PCAC and the Deficit of Forward Muons in pi+
Production by Neutrinos,” Phys. Lett., B657, 207–209.
Roesler, S., Engel, R., and Ranft, J. (2000), “The Monte Carlo event generator
DPMJET-III,” in Advanced Monte Carlo for radiation physics, particle transport
simulation and applications. Proceedings, Conference, MC2000, Lisbon, Portugal,
October 23-26, 2000, pp. 1033–1038.
Salam, A. (1968), “Elementary Particle Physics: Relativistic Groups and Analyticity(Nobel Symp. 8), ed,” N. Svartholm, p. 367.
Salam, A. and Ward, J. C. (1964), “Electromagnetic and weak interactions,” Phys.
Lett., 13, 168–171.
Sanuki, T. et al. (2000), “Precise measurement of cosmic ray proton and helium
spectra with the BESS spectrometer,” Astrophys. J., 545, 1135.
Sanuki, T. et al. (2002), “Measurements of atmospheric muon spectra at mountain
altitude,” Phys. Lett., B541, 234–242, [Erratum: Phys. Lett.B581,272(2004)].
Schael, S. et al. (2006), “Precision electroweak measurements on the Z resonance,”
Phys. Rept., 427, 257–454.
Schwinger, J. S. (1957), “A Theory of the Fundamental Interactions,” Annals Phys.,
2, 407–434.
Sekiya, H. (2013), “Solar neutrino analysis of Super-Kamiokande,” in Proceedings,
33rd International Cosmic Ray Conference (ICRC2013): Rio de Janeiro, Brazil,
July 2-9, 2013, p. 0048.
176

Sjostrand, T. (1995), “PYTHIA 5.7 and JETSET 7.4: Physics and manual,” .
Smith, R. A. and Moniz, E. J. (1972), “NEUTRINO REACTIONS ON NUCLEAR
TARGETS,” Nucl. Phys., B43, 605, [Erratum: Nucl. Phys.B101,547(1975)].
Sudarshan, E. C. G. and Marshak, R. E. (1994), “Origin of the Universal V?Atheory,” AIP Conf. Proc., 300, 110–124, [,0001(1984)].
Suzuki, A., Mori, M., Kaneyuki, K., Tanimori, T., Takeuchi, J., Kyushima, H.,
and Ohashi, Y. (1993), “Improvement of 20-inch diameter photomultiplier tubes,”
Nucl. Instrum. Meth., A329, 299–313.
Thébault, E., Finlay, C. C., Beggan, C. D., Alken, P., Aubert, J., Barrois, O.,
Bertrand, F., Bondar, T., Boness, A., Brocco, L., et al. (2015), “International
geomagnetic reference field: the 12th generation,” Earth, Planets and Space, 67,
1–19.
Verkerke, W. and Kirkby, D. P. (2003), “The RooFit toolkit for data modeling,”
eConf, C0303241, MOLT007, [,186(2003)].
Weinberg, S. (1967), “A Model of Leptons,” Phys. Rev. Lett., 19, 1264–1266.
Wendell, R. et al. (2010), “Atmospheric neutrino oscillation analysis with sub-leading
effects in Super-Kamiokande I, II, and III,” Phys. Rev., D81, 092004.
Wendell, R. et al. (2012), “Prob3++,” http://www.phy.duke.edu/~raw22/public/
Prob3++/, Accessed: 2017-01-10.
Wilkinson, C., Rodrigues, P., Cartwright, S., Thompson, L., and McFarland, K.
(2014), “Reanalysis of bubble chamber measurements of muon-neutrino induced
single pion production,” Phys. Rev., D90, 112017.
Wolfenstein, L. (1978), “Neutrino Oscillations in Matter,” Phys. Rev., D17, 2369–
2374.
Woods, R. D. and Saxon, D. S. (1954), “Diffuse Surface Optical Model for NucleonNuclei Scattering,” Phys. Rev., 95, 577–578.
Wu, C. S., Ambler, E., Hayward, R. W., Hoppes, D. D., and Hudson, R. P. (1957),“Experimental Test of Parity Conservation in Beta Decay,” Phys. Rev., 105, 1413–
1414.
Xing, Z.-z. and Zhao, Z.-h. (2016), “A review of µ ´ τ flavor symmetry in neutrino
physics,” Rept. Prog. Phys., 79, 076201.
Yamada, S. et al. (2010), “Commissioning of the new electronics and online system
for the Super-Kamiokande experiment,” IEEE Trans. Nucl. Sci., 57, 428–432.
177

Zeitnitz, C. and Gabriel, T. A. (1994), “The GEANT - CALOR interface and benchmark calculations of ZEUS test calorimeters,” Nucl. Instrum. Meth., A349, 106–
111.

178

Biography
Zepeng Li was born February 20, 1989, in Yantai, China. He attended University of
Science and Technology of China between 2007 and 2011. He received his Bachelor
degree in Physics in 2011. In the same year, he started graduate studies in physics at
Duke University. Under the direction of Professor Chris Walter and Professor Kate
Scholberg, he worked on the Long-Baseline Neutrino Experiment and the SuperKamiokande experiment.

179