Master Thesis

Development of the Java-based Simulator
for IceCube

02UM1102

Rie Ishibashi

A thesis submitted in partial fulfillmentof the requirements of the degree of Master of ScienceatGraduate School of Science and Technology
Chiba University

Abstract
The IceCube is an underground neutrino observatory now under construction at
the south pole. The advantages of the observation of the universe by neutrinos arethat they can arrive on the earth directly from the sources without deflection bymagnetic fields, and that they penetrate stars, gas and galaxies because they rarelyinteract with matters and radiations during their journey. Therefore the neutrinos
can propagate without losing the information of their sources.
In order to determine the energy and direction of extremely high energy (EHE) neutrinos with energies up to ∼ 1020 eV, the energy deposit as Cherenkov light emission
from electromagnetic and hadronic cascades in ice must be measured. However, theenergy deposit behavior has the stochastic nature in ultra high energy range which
can not be properly calculated by analytic calculations, because the interactions anddecay give rise to not a continuous energy loss but a stochastic loss.
Therefore we developed the Monte-Carlo/Numerical simulator named JULIeT –Java-based Ultra high energy Leptons IntEgral Transportor. It provides
two features to simulate the EHE particle propagation. The numerical calculations
are used to study the propagation of neutrinos and charged leptons in Earth overlong courses. Secondary produced EHE cosmic neutrinos and charged leptons arealso followed by our numerical calculation. The Monte-Carlo simulations are used toobtain the distributions of energy deposit in the deep-ice detector. By convoluting
the propagation in the earth and the energy deposit in the detector, we can evaluate
the IceCube sensitivity on EHE neutrino detection.
The 90% C.L. IceCube sensitivity for primary EHE neutrino flux has been foundto be placed at E 2 dF/dE ' 2.1 × 10−8 [GeV cm−2 sec −1 sr −1 ] for the primary νµand ντ with 109 GeV in absence of signals with energy deposit of 10 PeV or greater.

Acknowledgments
I would like to express sincere gratitude to my supervisors, Prof. S. Yoshida andProf. H. Kawai, for their support and advice in this work.
I would like to thank all of the members of Chiba IceCube group, Dr. K.Hoshina,
H. Miyamoto, K. Fujimoto, T. Noda for their suggestions and help.
I would also like to thank H. Nakayama, T. Ooba, Y. Shiino, Y. Unno, N. Orita,
T. Tsuda, H. Uchiyama, M. Konishi, T. Nunomura, Y. Ono, M. Tabata, N. Inadama,
S. Saitoh. They always supported and helped me, and I enjoyed my study and daylife with them.
I am very grateful to my parents for their support and encouragement.

Contents
1 Introduction
1.1 Introduction of EHE Neutrino Astrophysics . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Detection Method : Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 IceCube : The Underground Neutrino Observatory at the South Pole . . . . . . . . . . . . .

1
1
1
2

2 Propagation of Neutrinos and Charged Leptons
2.1 Interactions of Neutrinos and Charged Leptons . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Neutrino Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 Charged Lepton Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Decays of Charged Leptons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Issues and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5
5
5
6
11
12

3 Development of Java-based MonteCarlo/Numerical Simulator
14
3.1 Design Policy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2 Structure of JULIeT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2.1 Points Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2.2 Particles Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2.3 Interactions Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2.4 Decay Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2.5 Propagation Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2.6 Event Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.3 Simulated Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4 IceCube Sensitivity for EHE Neutrinos
30
4.1 Flux of Neutrinos and Charged Leptons at IceCube Depth . . . . . . . . . . . . . . . . . . . 30
4.1.1 GZK Neutrinos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.1.2 Z-burst . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.1.3 Top Down Scenario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.3 Flux at the IceCube Depth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.4 Evaluation of IceCube Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
5 Summary and Outlook
38
5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38i

List of Figures
1.1
1.2
1.3
1.4

IceCube:Underground Neutrino Observatory . . . . . . . . . . . . . . . . . . . . . . . . . . .
Schematic View of Digital Optical Module . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Top View of IceCube Observatory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
DOM Arrangement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
4
4

2.1

Cross Sections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

3.1 Package Structure of the JULIeT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Class Structure of the JULIeT –propagation . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Class Structure of the JULIeT –event . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Definition of geometric parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 Distribution of Energy Deposit in Ice for Muon . . . . . . . . . . . . . . . . . . . . . . . . .
3.6 Distribution of Energy Deposit in Ice for Muon . . . . . . . . . . . . . . . . . . . . . . . . .
3.7 Distribution of Energy Deposit in Ice for Tau . . . . . . . . . . . . . . . . . . . . . . . . . .
3.8 Distribution of Energy Deposit in Ice for Tau . . . . . . . . . . . . . . . . . . . . . . . . . .
3.9 Dependence of Energy Deposit on the Primary Energy for Muon . . . . . . . . . . . . . . .
3.10 Dependence of Energy Deposit on the Primary Energy for Tau . . . . . . . . . . . . . . . .
3.11 Dependence of the Range on the Primary Energy for µ and τ . . . . . . . . . . . . . . . . .
3.12 Distribution of Energy Deposit of Each Interaction/Decay in Ice for ν µ . . . . . . . . . . . .

15
17
18
19
25
25
26
26
27
27
28
29

4.1
4.2
4.3
4.4

Angular Distribution of Charged Lepton flux . . . . . . . . . . . . . . . . . . . . . . . . . . 33
Angular Distribution of Energy Deposit in the Detector . . . . . . . . . . . . . . . . . . . . 33
IceCube Sensitivities on EHE Neutrino Fluxes for νµ and for ντ . . . . . . . . . . . . . . . . 35
IceCube Sensitivities on EHE Neutrino Fluxes . . . . . . . . . . . . . . . . . . . . . . . . . 36

ii

List of Tables
2.1

Parameters in ALLM97 parametrization . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

3.1
3.2

Channels of Infinitesimal Propagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Channels of Interaction/Decay for Monte-Carlo Simulation . . . . . . . . . . . . . . . . . .

15
21

4.1
4.2

EHE Neutrino Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
90% C.L. Upper Limit for EHE Neutrino Fluxes . . . . . . . . . . . . . . . . . . . . . . . .

31
37

iii

Chapter 1

Introduction
1.1

Introduction of EHE Neutrino Astrophysics

The observation of high energy neutrinos provides a new window to search for high energy phenomena
in universe. Neutrinos can arrive on the earth directly from the sources, because they have no charge
and they are not deflected by magnetic fields. Additionally, neutrinos interact only weakly, so they rarelyinteract in their journey.
The extremely high energy (EHE) neutrinos with energies up to 1020 eV may be produced by mechanisms such as the Greisen-Zatsepin-Kuzmin (GZK) mechanism [1]. The GZK neutrinos via decays ofcharged pions produced by photo-pion production of EHE cosmic rays (EHECRs) colliding cosmic microwave background (CMB) photons. In this production mechanism no ντ s are generated because thepion-muon decay chain,pγ→ π ± + anything→ µ± + νµ (ν̄µ )→ e± + νe (ν̄e ) + ν̄µ (νµ ),does not have any channel to produce ντ . By the interaction of EHECRs and CMBs, EHECRs lose theirown energy during their propagation. Therefore the spectrum of EHECRs should have a cut-off feature
at around 5 × 1019 eV. This “GZK cut-off” has not clearly seen in the observed data. Whether the GZKmechanism, the most solid model of EHE neutrino radiation, takes place in the present universe still
remains a mystery.
There are a lot of other hypotheses to predict EHE neutrinos. The observation of EHE neutrinos andthe search of their sources lead to new understanding of high energy universe.

1.2

Detection Method : Overview

Because the cross sections of the neutrino interactions in this energy range are enhanced, they interact
with the nuclei of medium in the earth and produce neutrinos or charged leptons via neutral current
interaction (NC) or charged current interaction (CC);νl (ν̄l ) + N
νl (ν̄l ) + N→ νl (ν̄l ) + X(N C),→ l +X

(CC),

∓

1

where N is a nucleon of the medium, l is charged lepton (e, µ, τ ), X is produced hadrons.
Hadrons produced via NC or CC induce a hadronic cascade, and an electron produced via CC initiates
an electromagnetic cascade. Therefore EHE particles arriving at an underground detector are not justprimary neutrinos but secondary produced ν’s, µ’s, and τ ’s.
When a charged particle propagates in transparent medium with the speed faster than that of lightin the medium, the Cherenkov light is emitted. The number of radiated photons is proportional to path
length, and the emission angle depends on the refractive index of the medium. By seeing emitted Cherenkov
light and its angle, we can detect the charged particle and reconstruct the event track.
When an EHE particle propagate in a detector, the induced cascades contain many charged particles
which emit a bulk of the Cherenkov photons. Amount of such Cherenkov light packets is an estimator ofenergy deposit of an EHE particle.

1.3

IceCube : The Underground Neutrino Observatory at theSouth PoleThe IceCube is the underground neutrino observatory now under construction at the south pole. It
is the telescope to explore the high energy universe with cosmic neutrinos emitted from energetic distantsources. Since only neutrinos can penetrate underground, the earth is an effective shield to eliminate highenergy primary particles except neutrinos (See Figure 1.1). The main backgrounds are atmospheric muonsand neutrinos which will be described in details in section 4.2.
The flux of EHE neutrinos that arrive at the earth is very small. Therefore we need a detector withlarge target volume. The detector of the IceCube is made of polar ice with the volume of 1km 3 . It is madeof natural glacial ice with the depth from 1400m to 2400m under the ice surface. The glacial polar ice issuitable for the neutrino telescope like the IceCube [2] as demonstrated by the AMANDA experiment.
Neutrinos produced via NC, or charged leptons produced via CC can produce electromagnetic and/or
hadronic cascades by interactions. When a cascade is produced in the detector, the cascade particles
induce the Cherenkov light. The IceCube observes the neutrinos by detecting the Cherenkov light fromthese cascades using many photomultiplier tubes (PMTs) instrumented in ice.
The IceCube detector will instrument the 4800 PMTs which have 10 inch diameter. These PMTs will
be buried in vertical holes of ice drilled by hot water, so they have to endure very high pressure and low
temperature. Therefore each PMT is put in the pressure glass sphere comprise a digital optical module(DOM). The IceCube DOMs contain a signal processing electronics board, a LED flasher for calibration,
and a high voltage supplier. Figure 1.2 shows the schematic view of the DOM. The PMT in the glass
sphere is supported by gel as a shock absorber.
These DOMs will be arrayed on 80 strings, which contains 60 DOMs spaced by 17m. Each string arespaced 125m ∼ (Figure 1.3) and instrumented in the vertical holes of ice such as Figure 1.4.
We can obtain the profile of the EHE ν-induced cascades by detecting the Cherenkov light hit at many
DOMs arranged three dimensionally.

2

Figure 1.1: IceCube is the underground neutrino observatory. It uses the earth as a shield to cut thebackgrounds.

Figure 1.2: Schematic view of the Digital Optical Module (DOM) for the IceCube. It contains a PMT,
electronics board, LED flasher for calibration and a high voltage supplier. They are protected by pressure
glass sphere.

3

Figure 1.3: Top view of the IceCube observatory. The 60 DOMs will be arrayed on 80 strings which arespaced by 125m ∼. Each string is instrumented in the vertical holes of ice.

Figure 1.4: DOM arrangement in the IceCube detector. The each string has instrumented in the verticalholes of ice had 60 DOMs. They will be arranged three dimensionally.

4

Chapter 2

Propagation of Neutrinos andCharged Leptons
2.1

Interactions of Neutrinos and Charged Leptons

The EHE neutrinos that arrive at the earth interact with the medium in the earth as described in chapter 1.
The produced neutrinos and charged leptons via NC or CC can propagate in the earth interacting with
nucleons of matter around them.
The relevant interactions for EHE neutrinos are• neutral current interaction• charged current interaction
and that for EHE charged leptons are• pair creation• bremsstrahlung• knock-on electrons• photo-nuclear interaction
Figure 2.1 shows differential cross sections of interactions described in this section. The details of
calculations are described below.

2.1.1

Neutrino Interaction

Neutrinos interact with nucleon of the propagation medium, and produced neutrinos or charged leptons
via neutral current interaction (NC) or charged current interaction (CC);νl (ν̄l ) + N
νl (ν̄l ) + N→ νl (ν̄l ) + X(N C),→ l +X(CC),∓

where N is a nucleon of the medium, l is charged lepton (e, µ, τ ), X is produced hadrons.
For EHE neutrinos, the cross sections of those interactions depend on the behavior of parton distributions at very small values of the momentum fraction, however, where there are no direct experimental
measurement. So we used the cross sections calculated by CTEQ [3, 4].
52(log) Cross Section [cm]2(log) Cross Section [cm]−25
Mu to e Pair Creation
Bremsstrahlung
Photo−Nuclear
Mu to mu Pair Creation−30
CC

Knock−on ElectronNC−25
Tau to e Pair Creation

Photo−Nuclear
Tau to mu Pair Creation
Bremsstrahlung−30

Mu to tau Pair Creation
6

8

10

CC

Knock−on Electron

NC
6

12

8

Tau to tau Pair Creation
1012(log) Primary Energy [GeV](log) Primary Energy [GeV]

Figure 2.1: The left panel shows the cross sections in ice for νµ and µ. The right panel shows that for ντand τ . The most dominant interaction is pair creation producing e± pair. The photo-nuclear cross section
has large uncertainty from its parameterization.

2.1.2

Charged Lepton Interaction

Pair Creation
When high energy charged leptons (µ and τ ) interact with an nucleon of the propagation medium, a pairof the e± or µ± or τ ± is produced. The pair creation of an e± pair is dominant.
We define the inelasticity y as
y = 1−

Epro
Erec
=
,
Ein
Ein

(2.1)

where Erec is the energy of recoiling charged lepton, Ein is that of incoming charged leptons and Epro isthat of a produced pair (e± or µ± or τ ± ).
The differential cross section is calculated as Ref. [6];Z  m 2dσpairC
me 2 1 − ypro(Ein , y) = α2 Re2 Z (Z + ζ(Ein , Z))
φpro +
φrec dρ(2.2)
dy
mpro 3π yminρwhere α is the structure constant, Re is the classical electron radius and Z is the charge of propagation
medium. The ζ(Ein , Z) is the correction factor for screening effect on an atomic electron. It is parametrizedas
0
E ≥ 35min








Ein /min(2.3)
ζ(Ein , Z) =0.073 ln
2/3 E /m

1
+
γ
Z
1
in
in


otherwise
Ein /min
 0.058 ln
1 + γ2 Z 1/3 Ein /min

6

γ1

γ2

=


4.4 × 10−5=
4.8 × 10−5

for Hydrogen1.95 × 10−5

otherwise5.30 × 10−5

otherwise.

for Hydrogen

The ratio of electron mass and produced particle’s mass, me /mpro is the scaling factor for the flavorof produced particles.
The last term integraled over ρ is an asymmetry term, it is defined by!r
6m2in
E+ − E−
4mpro, |ρ| ≤ 1 − 2
,
(2.4)
ρ =
1−
E+ + E−
Ein (1 − y)Ein ywhere E + , E − is the energy of produced l + , l− , respectively.
The functions φpro and φrec is defined as






1 − ρ2 − β1+φpro =− (3 + ρ2 ) Lpro ,(2 + ρ2 )(1 + β) + ξ(3 + ρ3 ) ln 1 +
ξ
1+ξ
φrec

=

Lpro

Lrec



(1 − ρ2 )(1 + 2β)3ln(1 + ξ)(1 + ρ2 )(1 + β) −2ξξ(1 − ρ2 − β)+ (1 − ρ2 )(1 + 2β) Lrec ,
+
1+ξ



(1 + ξ)(1 + Ypro )RZ√
2mrec eRZ −1/3 (1 + ξ)(1 + Ypro ) 1+
Ein y(1 − ρ2 )"#
2
3mpro 1/3
1Z
(1 + ξ)(1 + Ypro ) ,+ ln 1 +2
2min
3mpro 1/3
Z


2min
,
= ln 
√−1/3
2mrec eRZ
(1 + ξ)(1 + Ypro ) 1+
Ein y(1 − ρ2 )

Yrec

=

Ypro

=

5 − ρ2 + 4β(1 + ρ2 ),2(1 + 3β) ln(3 + 1/ξ) − ρ2 − 2β(2 − ρ2 )4 + ρ2 + 3β(1 + ρ2 )(1 + ρ2 )(3/2 + 2β) ln(3 + ξ) + 1 − 3ρ2 /2y2,
β=2(1 − y)

(2.6)



p
−1/3


= ln 


(2.5)

ξ=


7

min y
2mpro

2

1 − ρ2
,
1−y

(2.7)

(2.8)

(2.9)

,

(2.10)

where R is the radiation constant, and e is the base of natural logarithm. The integration limits for paircreation is given by3 √ min 1/3
4mpro,
ymax = 1 −Z .
(2.11)
ymin =eEin
4Ein

Bremsstrahlung
When charged leptons interact with the electromagnetic field produced by atoms of the propagation
medium, an emitted high energy photon induced electromagnetic cascade. The differential cross section ofbremsstrahlung is parametrized as Ref. [10];
2 


1
min
dσbremss2
2 2
2(Ein , y) = 4αRe Z(2 − 2y + y )Ψ1 (qmin , Z) − (1 − y)Ψ2 (qmin , Z) , (2.12)
dy
mpro
y
3
where the functions Ψ1 (qmin , Z) and Ψ2 (qmin , Z) is defined byΨ1,2 (qmin , Z) = Ψ01,2 (qmin , Z) − ∆1,2 (qmin , Z),
=

1
2



m2 a 2
1 + ln in 12
1 + x1

Ψ02 (qmin , Z) =

1
2



2
m2 a 2+ ln in 12
3
1 + x1

Ψ01 (qmin , Z)

+

1
Z



(2.13)

1
1+− x1 arctan
x1
Z



 


1
m2in a22
1
1 + ln− x2 arctan(,2.14)2
1 + x22
x2



1
3
x21
+ 2x21 1 − x1 arctan+ ln
x1
4 1 + x21

 



m2 a 23x22
1 2
1+ ln in 22 + 2x22 1 − x2 arctan+ ln, (2.15)2 3
1 + x2x2
4 1 + x220(Z = 0)mζζ+1
ln in + ln
(Z 6= 0),qc
2 ζ −1 0∆2 (qmin , Z) =minζ
ζ + 1 m2in+ (3 − ζ 2 ) ln
+ 2
ln
qc
4
ζ −1
qc

∆1 (qmin , Z) =qmin =

m2in y,2Ein (1 − y)xi = ai qmin (i = 1, 2),

ζ=

s

1+

4m2in
,
qc2

a1 =

qc =

111.7
Z 1/3 mpro

(Z = 0)(Z 6= 0),

,

a2 =

724.2
Z 2/3 m,pro

1.9min.
Z 1/3

The integration limits for bremsstrahlung is given byymin = 0,ymax = 1 −8

3 √ min 1/3
Z .
e
4
Ein

(2.16)

Knock-On Electrons
The charged leptons also can “knock on” the atomic electrons of propagation medium. The differential
cross section of knock-on electrons is parametrized as Ref. [10];mpro 11 β2
1dσknock(1 + ∆(Ein , y)) ,(2.17)(Ein , y) = 2πRe2 Z 2−+dy
β E y2y ymax2∆(Ein , y) =  2




α
2yEin
2yEin
4Ein (1 − y)lnln 1 +− ln 1 +
2π
mpro
m2in
mpro

(2.18)

where β = p/Ein , p is the four momentum of the incoming charged lepton, and ymax is the upper limit ofintegration. The integration limits for knock-on electrons is given byymin =I¯,
Einymax =

2mpro β 2 Ein.
m2pro + m2in + 2mpro Ein(2.19)

where I¯ is the mean of ionization potential of the propagation medium.

Photo-nuclear Interaction
The photo-nuclear is the interaction of charged lepton with nuclei via virtual photon. The photo-nuclear
cross section is calculated using the deep-inelastic scattering formalism with the ALLM parametrization
for the nucleon structure function [13, 14].
The photo-nuclear interaction is described asl(kin )N (p) → l(krec )X,
where l is a charged lepton, and kin , krec is the four momentum of incoming and recoiling charged lepton,
respectively. The N is nucleon, and X is hadron changed via interaction. The X induces the hadronic
cascade.
The variables used in the calculations are defined asq2= (kin − krec )2 = −Q2 ,(2.20)x =Q2
2p · q(2.21)yp·q.
p · kin(2.22)=The differential cross section of photo-nuclear interaction is given by"!#
4Mp2 x2Q2dσphoto4πα22min y 2
1−y−1+
FN (x, Q2 )= 2
+ 1−
dQ2 dx
Q x
4Ein
Q2
2
Q2

9

(2.23)

where Mp is the proton mass, min is the incoming particle mass and Ein is the incoming particle energy.
The nucleon structure function is given by Ref. [13];FN (x, Q2 ) = a(A, x)FP (x, Q2 )(Z + (A − Z)P (x)),
P (x)= 1 − 1.85x + 2.45x2 + 2.35x3 + x4 .

(2.24)(2.25)The FP (x, Q2 ) is the proton structure function parametrized by ALLM97 [14]. It is given by
Q2

FP (x, Q2 )=F P (x, Q2 )= cP (t)xPP

F R (x, Q2 )= cR (t)xRR

Q2 + m20

F P (x, Q2 ) + F R (x, Q2 ) ,(2.26)a (t)(1 − x)bP (t) ,(2.27)a (t)(1 − x)bR (t) .

(2.28)

where xP and xR are defined by
xP

xR

=

W 2 − Mp2
1+ 2
Q + m2P!−1,(2.29)=

W 2 − Mp2
1+ 2
Q + m2R!−1.

(2.30)

where W2 is the center of mass energy (GeV2 ), m2P and m2R are given by Table 2.1.
The parameters in (2.27) and (2.28) is defined as functions of t;ln
t = ln Q2 + Q20

Λ2

Q20
ln 2
Λ

(2.31)

where Q0 and Λ are given by Table 2.1. Some of the parameters which is the function of t; aR (t), bR (t),cR (t) and bP (t) is given by
f (t) = f1 + f2 tf3 ,(2.32)

while aP (t) and cP (t) are given byg(t) = g1 + (g1 − g2 )




1
−1 .
1 + t g3

(2.33)The parameters in (2.32) and (2.33) (cP 1 , cP 2 etc.) are also given by Table 2.1.
The a(A, x) in (2.24) is the ratio of the nuclear structure function to that of nucleon. It is given bya(A, x)


−0.1


A
=
A0.069 log10 x+0.097



1

10

x ≤ 0.0014
0.0014 < x ≤ 0.04
0.04 < x.

(2.34)

m20 GeV2
m2P GeV2
m2R GeV2
Q20 GeV2
Λ2 GeV2
cP 1
cP 2
cP 3
aP 1
aP 2
aP 3

0.3198549.4570.150520.460170.065270.28067
0.22291
2.1979−0.0808−0.44812
1.1709

bP 1
bP 2
bP 3
cR1
cR2
cR3
aR1
aR2
aR3
bR1
bR2
bR3

0.60243
1.3754
1.84390.801070.973073.49420.58400
0.37888
2.60630.107111.9386
0.49338

Table 2.1: The parameters for proton structure function by ALLM97 [14].

The differential cross section dσ/dy is given by integration of dσ/dQ2 dx over dQ2 using valiable transformation from x to y. The valiable x and y are related by
x =

Q2
,
2Mp Ein y(2.35)

dx =

Q2
dy.
2Mp Ein y 2(2.36)The differential cross section is given by
dσ
dy

2.2

=

Z

dQ2

dσ dx
.
dxdQ2 dy

(2.37)

Decays of Charged Leptons

The µ decay and τ to lepton decay is calculated using the approximation that produced lepton mass is
negligible compared to parent lepton mass [11]. The differential decay rate is given bydn= 2 − 6yν + 4yν3 − (−2 + 12yν − 18yν2 + 8yν3 ),dyν

(2.38)

where yν = Eνe /Eνµ (for µ decay), Eνe,µ /Eντ (for τ decay), and
dn
5
4
1
8
= − 3zl2 + zl3 − ( − 3zl2 + zl3 ),
dzl
3
3
3
3

(2.39)

where zl = Eνl /El , l is charged lepton (µ or τ ).
The distribution of τ to hadron decay is given by
1
dn X
=
Bi
,
dz
1
−
ri
i

11

ri =

m2i
m2τ

(2.40)

where i indicates produced hadrons such as π, ρ, a1 , and X(other hadrons). The Bi is the branching ratioof the each decay mode which is normalized byZXdn
Bi = 1.
(2.41)
dz =
dz
i

2.3

Issues and Motivation

The accurate understanding of the EHE particles propagation in the earth is essential to search for EHE
cosmic neutrinos by an underground neutrino telescope. Let me explain why.
In the case of EHE particles, because the chain reactions of the interactions and decays occur during thedecayCC
P airCreation
propagation in the earth, e.g., νµ −−→ µ −−−−−−−−−→ τ −−−→ µ · · · , the behavior of particle propagationin EHE range is rather complex, and full calculations to follow all the interactions/decays and secondary
produced particles is indispensable to estimate the detection sensitivities. Full Monte-Carlo simulator tofollow all these particles may take too long time in practice, however, so we should develop the effective
calculation method for saving CPU time. It should be noted that any method for particle propagation in
EHE region demands a capability to treat µ’s and τ ’s on the same footing to handle the chain reactions.
Furthermore, as the particle intensity is expected to depend strongly on the nadir angle, one has to calculate
the particle flux in all phase space of solid angles because the IceCube is three dimensional array with 4π
detection sensitivity. The most practical way to satisfy these requirements is to build the proper transport
equations and numerically resolve them.
The transport equations which describe the propagation of neutrinos and charged leptons in the earth
are given by Ref. [15], using differential cross sections;for neutrino propagation:
Z
dNν
dNν
ml
1 dnl dNl
= −NA σCC+N C (El )
+
dEl
dxdEν
dEν
cτl ρ
El dEν dEl
Z
Z
dσ
dNl
dσ
dNν
+ NA dEl0 lN ·CC+ NA dEν0 νN ·N C
0dEν dEldEν dEν0for charged lepton propagation:
dNl
ml 1 dNl
dNl= −NA σCC (El )
−
dxdEl
dEl
cτl ρ El dEl
Z
1 dnl dNlml
dEl0 0+cτl ρEl dEl dEl0
Z
Z
dσ
dNν
dσ dNl+ NA dEν0 νN ·CC
dEl0 lN+NA
dEl dEν0dEl dEl0(2.42)(2.43)

where NA is the Avogadro’s number, c is the speed of light, ml is the mass of the charged lepton l, τl is
the lifetime of the l, x is the column density, unit in [g/cm2 ], which is defined byx≡

Z R

ρ(r)dr.

0

The ρ(r) is the density of the propagation medium, which depends on particle location.

12(2.44)

The first term of eq.(2.42) indicates the neutrino reduction due to the CC and NC. The second term
indicates the neutrino appearance due to the lepton decay. The third and fourth terms indicate theneutrino appearance due to the CC and NC, respectively. The dσlN ·CC /dEν term represents that theenergy interactions of the charged lepton l (due to the pair creation, bremsstrahlung, etc.) should betaken account.
The first and second terms of eq.(2.43) represent the charged lepton reduction due to the CC anddecay, respectively. The third term indicates the appearance of charged lepton due to the lepton decay.
The fourth term represents the neutrino appearance due to the CC. The last term accounts the interactions
such as pair creation, bremsstrahlung etc.
We need the numerical simulator which can solve the transport equations (eq.(2.42) and (2.43)) for thelong distance and for the all solid angles with saving CPU time.
The IceCube does not directly measure the energy of the EHE particle, but the energy deposit in thedetector. The energy deposit of an EHE particle are given by
−

dE
dx

= α + βE

(2.45)

' βE
Z

dσ
=
y dy E,dy(2.46)(2.47)

where the α is the ionization continuous loss term and β is the stochastic energy loss term by otherinteractions such as pair creation, bremsstrahlung, knock-on electrons, and photo-nuclear interactions.
The β depends on weakly the particle’s energy, E. The β term is dominant for EHE particles. In order tounderstand the stochastic energy deposit in the detector, we should handle it as probability distribution byusing the Monte-Carlo method. The Monte-Carlo method is the technique to sample a variable following
a given probability distribution using the random number. The relation between the random number r
and the probability distribution f (x) is given by
Z x
Z
r = F (x) =f (t)dt,
f (x)dx = 1.
−∞

By sorting the F (x), the uniform random number r samples a value of x.
In order to understand both the EHE particles propagation in the earth and in the detector, we developed
the simulator which can calculate propagation in the earth numerically, and the energy deposit in thedetector by the Monte-Carlo method.

13

Chapter 3

Development of Java-basedMonteCarlo/Numerical Simulator
3.1

Design Policy

We developed Java-based Monte-Carlo/Numerical simulator named JULIeT ∼ Java-based Ultrahigh energy Leptons IntEgral Transportor. In order to solve the issues described in section 2.3, the
JULIeT must meet a couple of requirements.
The first requirement is that the JULIeT should be able to trace neutrinos and charged leptons propagating in the earth by solving the transport equations (eq.(2.42) and eq.(2.43)) numerically. Therefore
the JULIeT is designed to solve the transport equations and produce the matrices describing the particle
propagation in the earth in order to obtain the fluxes of the particles arrived at the detector depth withall the solid angle. The second requirement is that the JULIeT should be also able to trace particles inthe detector i.e. in ice by Monte-Carlo method in order to obtain the distribution of the energy deposit
as electromagnetic and hadronic cascades produced in the detector.
The JULIeT is designed to have capability of calculating the IceCube sensitivity on EHE neutrino
detection with reasonable CPU time (chapter 3).
We introduced Java technology, one of the most modern programing languages, because• the Java is the object-oriented programming language which makes it easy to build large source codes
and update them over long years,• amount of source codes can be shorter than that in C++,• it provides an OS-independent platform.

3.2

Structure of JULIeTThe JULIeT is based on a package system. Each package contains classes under similar concepts. Fig.3.1
shows the structure of the packages. The details of the roles of each packages and their classes are described
below.
Table.3.1 shows the first order interactions/decays channels considered in the JULIeT. The details ofeach interactions/decays channels are described in section 2.1 and 2.2. The particles handled in the JULIeT

14

Table 3.1: Channels of infinitesimal propagation.
Rows are primary and columns are produced/survival particles.
νeNCa /Sh

νµ

ντ

νe
νµ
NC/S
ντ
NC/S
c
µ
D
CC/D
τ
D
D
CC/D
a
Neutral Current interactionbCharged Current interaction
c
Decay
d
Pair CreationeBremsstrahlungfKnock-on Electrong
Photo-nuclear interactionhSurvive

e/γ
CCb

µ

τ

CC
d

e

f

PC /B /K /D
PC/B/K/D

PC/S
PC/D

CC
PC
PC/S

hadron
NC/CC
NC/CC
NC/CC
NC/CC/PNg
NC/CC/Dare νe , νµ , ντ , µ and τ . Allowed region of primary energy of particles is from 106 [GeV] to 1012 [GeV], andthat of produced energy of particles is from 102 [GeV] to 1012 [GeV].

JULIET
Event

Particle
Interactions

Decay
Point

Propagation

Figure 3.1: Package Structure of the JULIeT.

The representative packages among them are propagation and event both of which has its own method
to run the simulation. The propagation package contains the classes to calculate the particle propagation
by numerically solving the relevant transport equations (eq.(2.42) and eq.(2.43)). The event package is agroup of the classes to calculate the particle propagation with Monte-Carlo methods. It has the interactive
15

mode to run the simulation. The propagation and event packages rely on classes in all the other packages,
such as the interaction package which is responsible for handling the particle interactions. The relation
between the classes are shown in Figure 3.2 for the propagation package, and Figure 3.3 for the event
package.

3.2.1

Points Package

This package contains the classes concerning the particle location and the propagation medium. The
class in this package, ParticlePoint class, provides the point vector used by calculations on propagationof EHE particles in the earth. Figure 3.4 shows the definition of geometric parameters. It also provides
the parameters of the medium such as the density which are used for the calculations of differential crosssections.

3.2.2

Particles Package

This package contains the Particle class to define particles and their property such as their mass, name,flavor, doublet, and lifetime. It also provides their energy.

3.2.3

Interactions Package

This package contains the classes to define interactions such as PairCreation class, Bremsstrahlung classetc. They are classes inherit from Interactions class. It also has the classes (MakeBremsstrahlungMtx
class etc.) to make dumped matrices which elements are differential cross sections corresponding to logenergy of primary and produced particles. The dumped matrices defined in InteractionsMatrix class are“serialized” objects, which provide “get” methods to return the differential/total cross sections, particle’s
flavor, doublet, energy and so on. The differential cross sections are calculated as described in section 2.1.
The matrices made in this package are used by calculations of propagation in the earth.
Additionally, this package has InteractionsBase class and ElectronBase class which are extended from its
super class, MonteCarloBase class in event package (See 3.2.6). These “Base” classes exist in order thatclasses in event package which calculate the particle propagation by Monte-Carlo method can handle all of
the interactions/decays under the same footing. These classes essentially provide the probability tables tosample the pathlength of interaction and produced energy with Monte-Carlo methods (See 3.2.6). It also
supplies the method to return the pathlength and produced energy for a given random number.

3.2.4

Decay Package

The classes in this package handle the lepton decay processes such as τ → µ →e, and make decay matrices
which elements are differential decay rate.
The Tau(Mu)DecayYMatrix class are similar with the InteractionsMatrix class in the interaction package.
These classes make the matrices of the energy transfer dN/dLogy for the decays. The matrix elements
are calculated by the methods supplied by the Decay class. The decay rates are calculated as described in
section 2.2. The biggest difference between Tau(Mu)DecayYMatrix class and InteractionsMatrix class is, thedecay matrix element is calculated whenever the object is generated while interaction matrices are usually
pre-calculated and stored in form of files to be read out for the Event/PropagationMatrix class. This isbecause it takes much less time to calculate the decay matrices than the interaction matrices.
There are also “Base” classes, MuDecayBase and TauDecayBase, for calculations of the particle propagation by Monte-Carlo methods. Similar with the InteractionsBase class, they calculate the Monte-Carlo16

propagation package

particle package

MainRunParticle

RunPropagationMatrix
ParticlePoint
PropagationMatrix

points package
interaction package

Dumped Data

Interactions(InteractionsMatrix objects)
MakeNeutrinoChargeMtx

NeutrinoCharge

MakeNeutrinoNeutralMtx

NeutrinoNeutral

MakeBremsstrahlungMtx

Bremsstrahlung

MakeKnockOnElectronsMtx

KnockOnElectrons

MakePairCreationFitMtx

PairCreation

MakePhotoNuclearFitMtx

PhotoNuclear

PairCreationFit

MakePairCreationTable

PhotoNuclearFit

MakePhotoNuclearTable

MuDecayMatrix
Decay
TauDecayMatrix

decay package
Figure 3.2: Class structure of the JULIeT from aspects of the Propagation package. MainRun class has amain program for tracing neutrinos and charged leptons in the earth.

17

event package

particle package

RunManager

Particle

Event
ParticlePoint
MonteCarloBase

points package

IntearctionsBase
Interactions
Dumped Data
(InteractionsMatrix
objects)

ElectronBase

interaction package

MakeNeutrinoChargeMtx

NeutrinoCharge

MakeNeutrinoNeutralMtx

NeutrinoNeutral

MakeBremsstrahlungMtx

Bremsstrahlung

MakeKnockOnElectronsMtx

KnockOnElectrons

MakePairCreationFitMtx

PairCreation

MakePhotoNuclearFitMtx

PhotoNuclear

PairCreationFit

MakePairCreationTable

PhotoNuclearFit

MakePhotoNuclearTable

MuDecayBase

MuDecayYMatrix

TauDecayBase

TauDecayYMatrix

Decay

decay package
Figure 3.3: Class structure of the JULIeT from aspects of the event package. RunManager class has a main
program for tracing neutrinos and charged leptons in the IceCube detection volume.

18

AxisLengthθαREarth

laxis

α

R

particle

Earth

Figure 3.4: Definition of geometric parameters. The α is nadia angle, θ is zenith angle.

probability table using the decay matrices calculated by the MuDecayYMatix class and TauDecayYMatrix
class.
The τ decay has some decay modes such as τ → µ, e, and hadrons. The branching ratio of each modeare supplied in the Decay class. The hadrons produced via τ decay may be π, ρ and so on. But we regard
them as a hypothetical ’hadron’ which is made of superpositions of π, ρ and so on.
The mode of τ decay is determined by the method supplied in TauDecayBase class. In the method, acumulative table of branching ratio is produced, and one of the decay mode is chosen by random number.

3.2.5

Propagation Package

This package contains the classes for numerically solving the transport equations (eq.(2.42) and eq.(2.43)),
which provides a faster and more accurate evaluation of the particle propagation in long distances.
The propagation matrices calculated in this package describe the energy and particle species transfer by
the interactions/decays. They are built by pre-calculated InteractionsMatrix and MuDecayYMatix/TauDecayYMatrix
objects. The elements of the matrices indicate the energy distribution and the flux of particles after propagation of an infinitesimal distance in the earth. In this package, there is the MainRun class which has amain method to calculate the propagation matrices.
The procedure of making propagation matrices is described below. Figure 3.2 shows the relation between the classes in this package.
Calculation of the propagation matrices:
1) In order to avoid sizable error due to the rock density profile, the trajectory length is divided to
5 sections.
2) The following calculations 3) ∼7) is repeated 5 times (for each section).
3) The elementary infinitesimal propagation matrices over infinitesimal distance (dX[g/cm 2 ]) are built
using differential cross sections calculated in the interactions package.
4) The particles are propagated over X = 10 × dX[g/cm2 ], i.e. infinitesimal propagation matrices
19

produced in 3) is multiplied 10 times.
5) The particles are propagated over X × 2k [g/cm2 ], where k is the number of steps needed to propagate
the section, i.e. finite propagation matrix produced in 5) is multiplied 2k times.
6) During 5), the matrices are stored when particles reaches 3 % of the total path in their journey.
7) After 5), the resultant matrices propagated over the section are stored.
8) For a rest of trajectory length, particles are propagated by multiple the matrix stored in 5) until
the particles reach the final point.

3.2.6

Event Package

This package contains the classes related to ’occurrence of event’. These classes define the behavior ofan event for propagating particles by Monte Carlo methods.
There is the MonteCarloBase class in this package. It is a super class of InteractionsBase and ElectronBase
which are classes in the interactions package and MuDecayBase and TauDecayBase which are classes in the
decay package. The Event class and RunManager class handle the interactions/decays in form of “list”.
In order to register interactions and decays in the list, all the interactions/decays should have a common
super class. The MonteCarloBase class is an abstract class defines the methods for both interactions anddecay that determines the pathlength and produced energy with the Monte-Carlo methods. The actual
implements of the methods are made by each “Base” class. It also defines methods to get the informationof interactions/decays so that the Event class can know that information.
The Event class actually makes an event. It is the fundamental class for propagating particles by theMonte-Carlo methods. Table 3.2 shows channels of the interaction/decay considered in the Monte-Carlo
methods. The CC interaction from charged lepton to neutrino is negligible because propagation distance(detector size) is short enough. The neutrinos produced via charged lepton decay are also neglected by
the same reason. In order to save CPU time, the differential cross sections of these interactions are
pre-calculated by MakeBremsstrahlungMtx class etc. in Interactions package. An interaction/decay which
returns the shortest pathlength is chosen, and an “event” has occurred then particles propagate. So theinformation of an interaction/decay which has occurred, such as transfered energy, produced particle,
pathlength etc. are supplied the methods of Event class.
The procedure of choosing an interaction/decay which returns the shortest pathlength, and choosing
the produced energy via the interaction/decay by Monte-Carlo method are described below. Figure 3.3
shows the relation between the classes in this package.
Calculation of the pathlength by the Monte-Carlo method:a) In the case of the changed lepton interactions
The number of particles when they propagate X[g/cm2 ] is given bydNdX= −σ(E)N,(3.1)

N (E, X) = N (E, X0 ) exp[−(X − X0 )NA σ(E)],(3.2)

where NA is the Avogadro’s number. The probability that interaction has occurred while a particle
propagates X[g/cm2 ] is given by
P (E, X) = 1 − exp[−XNA σ(E)].
20

(3.3)

Table 3.2: Channels of interaction/decay for Monte-Carlo simulation. Rows are primary and columns areproduced particles. The CC interaction from charged lepton to neutrino is negligible for a propagation distance (detector size) considered here. The neutrinos produced via charged lepton decay are also neglectedby the same reason.
νe
NCa

νµ

ντ

e/γ
CCb

νe
νµ
NC
ντ
NC
PCd /Be /Kf /Dµτ
PC/B/K/Da
Neutral Current interactionbCharged Current interaction
c
Decay
d
Pair CreationeBremsstrahlungfKnock-on Electrong
Photo-nuclear interactionhSurvive

µ

τ

PC
PC/Dc

PC
PC

hadron
NC/CC
NC/CC
NC/C
PNg
D

The P(E,X) is determined by a random number(0∼1). So the pathlength isX−1ln(1 − r),
NA σ(E)=(3.4)

where r is a random number.
b) In the case of the neutrino interactions
The probability that interaction has occurred while a neutrino propagates X[g/cm 2 ] is given by

P (E, X) = 1 − exp[−XNA σ(E)](3.5)' XNA σ(E),(3.6)

where under an approximation that cross section of neutrino interaction is very small.
The P(E,X) is determined by random number(0∼1). So the pathlength is
X

=

r
,
NA σ(E)

where r is random number.
c) In the case of decay

21

(3.7)

The number of particles when they propagate while t[sec] is given bydN
dt

1
N (t),
τ (E)
h
t i
,
N (E, t) = N (t0 ) exp −
τ (E)
= −

(3.8)
(3.9)

where τ is the lifetime of the particle. The probability that decay has occurred while a particle
propagates t[sec] is given byhP (E, t) = 1 − exp −

t i
.
τ (E)

(3.10)

From X = cτ ρ (c is light speed and ρ is the density of the medium[g/cm3 ]),h
P (E, X) = 1 − exp −X i.
cτ (E)ρ

(3.11)The P(E,X) is determined by random number(0∼1). So the pathlength is
X = −cτ (E)ρ ln(1 − r),(3.12)

where r is random number.

Calculation of the Produced Energy by the Monte-Carlo method:
The produced particle’s energy is determined by Monte-Carlo method. Energy deposit anddimensionless energy y = Epro /Ein are related by eq.(2.47). A y is sampled using a random number
by
r

=

1
σ

Z

dσ
dy
dy

(3.13)

where r is a random number.
The RunManager class controls the run, so it has a main method which runs the events and records
some types of data files. It calculates the distribution of the energy deposit includes their fluctuation. It
can run some types of runs for recording some types of data e.g., full-text-based data which records all
the profile of the propagation for debugging, binary data which records distributions of energy depositin the detector for each registered interaction/decay, and so on. For example, one of the procedure ofcontroling the run which records binary data of the distributions of the energy deposit to electromagnetic
and hadronic cascades, respectively, with given monocromatic primary energy, is described below;
Procedure of the run:1) The constructor of the RunManager class generates a list of MonteCarloBase object whose elements
are “Base” classes for each interaction/decay registered by an user interactively.
2) A Particle object is generated as a primary particle (See 3.2.2).
3) A ParticlePoint object is generated with given material(ice or rock) as a starting point (See 3.2.1).
4) An Event object is generated (See above).
5) The pathlength of interactions/decays in the list generated in the constructor is sampled and22compared with each other, then an interaction/decay which returned the shortest pathlength is determined
to take place.
6) The particle is propagated over the pathlength determined in 5).
7) If the propagating particle is still in the detector, calculations is continued.
8) The energy of the produced particle is sampled by Monte-Carlo method (See above).
9) The procedure 5) ∼ 9) is repeated until the particle stops or changes into an electron(photon)/hadronor passes through the detector.
10) The distributions of energy deposit as electromagnetic/hadronic cascades and total energy deposit
are recorded during the procedure of 2) ∼ 10)

The simulated results calculated in this RunManager class by Monte-Carlo methods will be shown insection 3.3.

3.3

Simulated Results

As described in section 3.2, the RunManager class in the event package can calculate the particle propagation with the Monte-Carlo methods.
We ran the particles 1km in the ice, and calculated the distribution of energy deposit of them. The
produced particles via interactions/decays induce electromagnetic or hadronic cascades. The resultant
distributions are compared to analytic calculations.
Figure 3.5 and 3.6 shows the distribution of energy deposit for muon propagating in 1km ice (thedensity of the ice is ρ = 0.917). The energy of primary muon is 10EeV (= 10 10 GeV). The analytic meanvalue of energy deposit in Figure 3.5 is given by∆E= 1 − eβx ,−Ewhere x is the path length and β is the inelasticity defined as
Z
X
dσi
βi ,
βi = yβ=dy,dyi(3.14)

where the i is the index of relevant interactions. Figure 3.5 shows the distribution of energy loss toelectromagnetic and hadronic cascades, respectively. The mean value of these energy deposit is consistent
with analytic calculation, but you can see the fluctuation of the energy deposit shown in Figure 3.5 is not
negligible.
In order to understand causes of the fluctuation, the contribution of the e± pair creation, bremsstrahlung
and photo-nuclear interaction to the energy deposit are separately shown in Figure 3.6. While the energy
deposit via e± pair creation has a peak, those via bremsstrahlung and photo-nuclear interaction exhibits
the divergent structure. Therefore these stochastic energy loss contribute to make the distribution wide.
While the following figures show the same as above, but for τ . Figure 3.7 shows the same as Figure 3.5,
but for the case of τ propagation, the contribution from decay is larger than the case of µ, because
the lifetime of τ is far short compared with that of µ. The contribution of decay can be shown in
Figure 3.8. The cross section of photo-nuclear interaction has large uncertainty (See Figure 2.1), because
the parameterization of cross section depends on the nuclei structure function. Therefore we may need theother parameterization for photo-nuclear interaction in this energy range for systematic investigation.
23

The relation between the energy deposit obtained by the Monte-Carlo simulation and analytical calculation are shown in Figure 3.9 for µ and Figure 3.10 for τ as a function of primary energy. The energy
deposit of µ is larger than that of τ above 108 GeV region. But, below 108 GeV region, the contributionof τ decay increases, which make a difference between the simulated result and analytical calculation.
We can also shows the contribution of interactions which induce stochastic energy loss as follows. The
average distance which particles propagated until its energy becomes 1/e times of the primary energy, i.e.,“range” is determined by Monte-Carlo simulation. Figure 3.11 shows the range for µ and τ . The blue
and magenta lines indicate expected range obtained by analytic calculation (blue for the τ , magenta forthe µ) 1/β by eq.(3.14). The instability of the photo-nuclear cross section adversely affects both analytic
and simulated values. The green line in indicates the contribution of the τ decay. Below 10 9 GeV, thecontribution of τ decay plays a dominant role on the τ propagation. You can see that the range of τ islonger than that of µ. It suggests that the detection volume of τ is larger than that of µ, i.e. the τ flux
arrive on the detector is larger than µ flux above ∼ 108 GeV. However, since the contribution of τ decay
is dominant, the range of τ becomes short below ∼ 108 GeV.
We have discussed the case when incoming particles are muons and taus. But the IceCube detector
will occasionally see events induced by incoming neutrinos events in EHE range. We, therefore, need toinvestigate the case when neutrinos are incoming into the 1km3 volume, and producing muons and tausinternally. Figure 3.12 shows the energy deposit distributions when νµ s are incoming into the detector.
The νµ produces µ via CC interaction, but the location where CC interaction has occurred is random.
The long tail shown in Figure 3.12 reflects the difference of the distance which produced µ propagated.
In this way, we obtained the distributions of energy deposit in the detector for neutrinos and charged
leptons propagation. These results are used to evaluate the IceCube sensitivity on EHE cosmic neutrinos
given by various models.

24

−1

(log) dF/dLogE

Muon 10EeV
Total Cascade
Electromagnetic Cascade
Hadronic Cascade

Mean (MC)Mean (Analytic)

−2

−3

−4
6

8

10

12

(log) Energy Deposit [GeV]
Figure 3.5: Distribution of energy deposit in 1km ice for muon. The green line shows the energy deposit byelectromagnetic cascade, the red line for hadronic cascades, and the blue line for the total energy deposit
in 1km ice. The unit of dF/dLogE is arbitrary.

−1

Muon 10EeV(log) dF/dLogEe Pair Creation−2
Photo−nuclear
Bremsstrahlung

−3

−4
6

8

10

12

(log) Energy Deposit [GeV]
Figure 3.6: Distribution of energy deposit in 1km ice for muon. The contributions from pair creation,bremsstrahlung and photo-nuclear interaction are shown respectively.
25

−1Mean (Analytic)

Total Cascade
Electromagnetic Cascade
Hadronic Cascade−2(log) dF/dLogE

Tau 10EeV

Mean (MC)

−3

−4

−5

−6
6

8

10

12

(log) Energy Deposit [GeV]
Figure 3.7: Distribution of energy deposit in 1km ice for tau. The green line shows the energy depositby electromagnetic cascade, the red line shows that by hadronic cascades. The blue line shows the totalenergy deposit in 1km ice. The unit of dF/dLogE is arbitrary.

−1

Tau 10EeV(log) dF/dLogE−2

Tau to e Pair Creation
Photo−nucear−3

Mu to e Pair Creation−4

Tau to Hadron Decay−5

Bremsstrahlung
Tau to e Decay

−6
6

8

10

12

(log) Energy Deposit [GeV]
Figure 3.8: Distribution of energy deposit in 1km ice for tau. The line of photo-nuclear interaction is very
instable because the cross section of photo-nuclear interaction has large uncertainty.
26

(log) Energy Deposit [GeV]

Mu
10

5

5

10

(log) Primary Energy [GeV](log) Energy Deposit [GeV]

Figure 3.9: Dependence of energy deposit on the primary energy for muon. The red points shows theenergy deposit obtained by Monte-Carlo simulation, and blue line shows the analytical value.

Tau
10

5

5

10

(log) Primary Energy [GeV]
Figure 3.10: Dependence of energy deposit on the primary energy for tau. The red points shows the energy
deposit obtained by Monte-Carlo simulation, and blue line shows the analytical value. The departure fromthe analytical estimation which is visible below 108 GeV is due to the contributions by the τ decay.

27

40

Range [km]

30

20

10

0
6

8

10

12

(log) Primary Energy [GeV]
Figure 3.11: Dependence of the range on the primary energy for µ (magenta) and for τ (blue). The point
data shows the energy deposit obtained by Monte-Carlo simulation, and line data shows the analyticalvalue.

28

Muon Neutrino 10EeV(log) dF/dLogE−2e Pair Creation−4

CC

BremsstrahlungNC−6
Photo−nuclear

−8
6

8

10

12

(log) Energy Deposit [GeV]
Figure 3.12: Distribution of energy deposit of each interaction in 1km ice for ν µ . The NC interaction
produce hadronic cascades. The CC interaction also produce hadronic cascades, and produced µ induces
the interaction such as pair creation, bremsstrahlung and photo-nuclear interaction. The long tail shownin pair creation indicates the difference of the location where the µ was propagated.

29

Chapter 4

IceCube Sensitivity for EHE
Neutrinos
4.1

Flux of Neutrinos and Charged Leptons at IceCube Depth

As described in chapter 3, the JULIeT has a capability to simulate propagation of neutrinos and charged
leptons in the earth. Here we calculate the fluxes of neutrinos and charged leptons by propagating the
primary EHE neutrinos at the detector depth using the transfer matrices which are generated a priori by
PropagationMatrix class in the propagation package (See 3.2.5). Full mixing of ν oscillation is assumed andmade proper adjustment in the primary flux of EHE cosmic neutrinos.
Table 4.1 summarized the EHE neutrino models and their parameters, which we used for estimationof the EHE cosmic neutrino flux at the surface.
They are input spectra to our propagation calculation, then we discuss the IceCube sensitivity for
EHE neutrino fluxes taking into account their propagation effects in the earth and the energy deposit
characteristics of the secondary leptons.

4.1.1

GZK Neutrinos

The models 1∼6 [16] and 10, 11 [20] in Table 4.1 are the GZK neutrino models. The GZK neutrinos are
generated via decays of charged pions which are produced by photo-pion production of EHECRs colliding
CMB photons. This production chain,pγ→ π ± + anything→ µ± + νµ (ν̄µ )→ e± + νe (ν̄e ) + ν̄µ (νµ ),does not have any channel to produce ντ . There is no ντ s are primarily produced within the GZK neutrino
production.
The GZK model 1∼6 from Ref. [16, 17] has a few parameters concerning the cosmic ray source distribution. m is the evolution parameter which relevant to the luminosity of EHECR as a function of red
shift (m = 0 corresponds to the no evolution). Zmax is the “turn-on time”, which means the time when
EHECR were created. γ indicates the energy spectrum of primary EHECRs following E −γ .

30

1
2
3
4
5
6
7
8
9
10
11

model
The GZK Neutrinos
The GZK Neutrinos
The GZK Neutrinos
The GZK Neutrinos
The GZK Neutrinos
The GZK Neutrinos
The Z-burst
The Top Down (SUSY)The Top Down (QCD)
The GZK Neutrinos (by Sigl)
The GZK Neutrinos (by Sigl)

m
0
2
2
4
4
7

Zmax
2
2
4
4
5
5

γ

1.5
1.5

3
5

2
2

1
2

Emax[ZeV ]100
10

Table 4.1: EHE neutrino models. Model-1∼6 are GZK Neutrinos [16], Model-7 is Z-burst model [18],Model-8,9 are Top Down models [19], Model-10,11 are GZK Neutrinos [20]. The m, Z max , γ, andEmax[ZeV ] are parameters which individualize the models. Note that there is no ντ in GZK model.
The another GZK model 10 and 11 from Ref. [20] has one more parameter, Emax is the maximum
energy of primary EHECRs.

4.1.2

Z-burstThe Z-burst model [18] is about interactions of EHECRs with cosmic background neutrinos which are arelic of the Big Bang. The cross section of the interaction has a resonance when the collision energy is themass of Z0 boson. The produced charged leptons or hadrons induce the neutrino cascade. This cascadecan produce the EHE particles above GZK threshold.

4.1.3

Top Down Scenario

The model 8, 9 in the Table. 4.1 is the top down scenario [19]. This scenario supposes that the EHE
cosmic rays and neutrinos are produced from massive particles via an interaction/decay. The top down
scenario may be able to explain the existence of EHE cosmic rays beyond GZK cut-off.

4.2

BackgroundThe backgrounds in detecting EHE cosmic neutrinos are atmospheric muons and neutrinos. They are
produced via the decay of mesons produced by the interaction of cosmic rays with atmosphere. The
atmospheric muons and neutrinos can arrive at the detector depth.
Since the energy spectrum of atmospheric neutrinos falls steeply like dNν /dEν ∼ E −3.7 , the EHE
cosmic neutrino energies extend to higher values. The dependence of the flux on the zenith angle can becharacterized by


dF
dE


= kE −(γ+1) 


1

1
0.054
+
1 + 1.1(E/E1 ) cos θ 1 + 1.1(E/E2 ) cos θ
31


,


(4.1)

where the cosmic ray flux is proportional to E −γ , E1 =115GeV and E2 =850GeV are the factors due tothe dependence of path length on the zenith angle. The atmospheric muons should contribute to only thedown-going events. The energy spectrum of atmospheric muons may be steeply enough to be distinguishedfrom signal events by setting a propaer threshold energy in filtering.

4.3

Flux at the IceCube Depth

Figure 4.1 shows the angular distributions of the relevant fluxes above 10 PeV of µ and τ at the detector
after propagation in the earth. The fluxes of µ and τ are calculated by propagating the primary GZK
neutrinos (with m = 4, Zmax = 4 in Ref. [16]) in the earth. It should be noted that the full mixing ofν oscillation of the primary cosmic neutrino is assumed i.e., νe : νµ : ντ = 1 : 1 : 1. The intensity of
the fluxes strongly depends on zenith angle. The up-going particles are strongly attenuated by the earth,because they have to propagate a distance longer than down-going particles The main contribution to thefluxes is the down-going events (cos θ ≥ 0, θ is the zenith angle of incident particle). The flux of τ is larger
than that of µ in up-going and horizontal region. It indicates that the detecting τ events with up-goingand horizontal angle is advantageous. The reason the τ has larger flux is that the τ propagate in th earth
with losing its energy slowly compared with µ regeneration process like ντ → τ → ντ also contributes theflux.
There are two lines of the fluxes of atmospheric muons are also shown in Figure 4.1. The one iscalculated by the extrapolation from the flux of 5 TeV with assuming that cosmic ray spectrum has E −2.7
form and the another is given by Corsika [22] with assuming E −3 . The contribution of atmospheric muonsis, however, also large where cos θ ≥ 0. Figure 4.2 shows the angular distributions of fluxes which deposit
their energy above 10 PeV. The contribution of atmospheric muon is small because the energy of them isrelatively lower than that of µ and τ .
It also shows that the horizontal events (cos θ ∼ 0) are the mostly responsible for an observed bulk of
EHE events.

32

Upgoing

Downgoing

3

Flux [1019 cm−2 sec −1 sr −1]

GZK Neutrino

Atm MuonE > 10PeV

τ

2

θ

µ
1

0
−1

−0.5

0

0.5

1

Cos θ
Figure 4.1: The angular distributions of fluxes of µ and τ after propagation in the earth. The fluxes above 10
PeV are plotted. θ is zenith angle of incident neutrinos defined as the right-side figure. In the case of cos θ ≥ 0, it
indicates down-going events. The primary neutrinos are given by GZK model in Ref. [16] with m = 4, Z max = 4.
The fluxes of atmospheric µ are also plotted, solid line is given by the extrapolation from 5 TeV, dashed line isgiven by Corsika estimation [22].

Flux [1019 cm−2 sec−1 sr−1]

3

GZK Neutrino
E Dep. > 10PeV
2µ
1

Atm Muon

τ
0
−1

−0.5

0

0.5

1

Cos θ
Figure 4.2: The angular distributions of energy deposit in the detector for µ and τ . The fluxes which deposit theirenergy above 10 PeV are plotted. θ is zenith angle of incident neutrinos. The primary neutrinos are given by GZKmodel in Ref. [16] with m = 4, Zmax = 4. The fluxes of atmospheric µ are also plotted, solid line is given by the
extrapolation from 5 TeV, dashed line is given by Corsika estimation [22].
33

4.4

Evaluation of IceCube Sensitivity

The IceCube sensitivity is evaluated by the event rate per energy decade dN/dLogE. In order toobtain the fluxes of neutrinos and charged leptons incident to the detector as a function of zenith angle,
the numerical transportation in the earth is carried out for a given flux of primary neutrinos. Then
convoluting the results of Monte-Carlo simulation, we obtained the IceCube sensitivity.
Figure 4.3 shows the resultant sensitivity with 1km2 detection area for νµ and for ντ , respectively.
The various model predictions are also plotted together for comparison. The event flux which lose theenergy above threshold is calculated by Monte-Carlo simulation described in section 3.3. The obtainedenergy rate is integrated over zenith angle. The 90% C.L. upper limit, i.e., (2.3 events)/(10 years)/(1km 2
detection area)/(energy decade), for the threshold of the energy deposit in the detector for 1 PeV and 10
PeV are shown.
The ντ sensitivity is better than that of νµ above 109 GeV region because the detection volume of τis larger than that of µ. While, below 109 GeV region, the νµ sensitivity is better than ντ sensitivity,
because the energy deposit of µ is basically larger than the case of τ (See Figure 3.10 and 3.9). But, theνµ sensitivity has a bump structure around 108 GeV region. It is due to the contribution of the τ decay
which induces large energy deposit.
Additionally, Figure 4.4 shows the sensitivity which include both ντ and νµ contributions. This figure
indicates that the model above the lines of threshold can be excluded in absence of signals with 90% C.L..
The solid blue lines show the limits obtained by the case of the incident particle is charged leptons. Whilethe dashed lines show the limits in the case of the charged leptons produced via CC inside the detector.
Below 108 GeV region, charged leptons must be produced near the detector volume to remain their energies
above the threshold, which leads to reduction of the detector volume, while energy deposit at the CC/NC
vertex contributes total energy loss in significant way in case of neutrino primary events. The resultant
sensitivity around the threshold energy will be, consequently, enhanced by looking at not just muon/tau
tracks but neutrinos inducing cascades inside the detector volume. From Figure 4.4, the IceCube has thesensitivity enough to exclude the GZK model which is the most solid model for the EHE cosmic neutrinos.
The Table 4.2 shows 90% C.L. upper limit of primary EHE neutrino fluxes, which are placed by
E 2 dF/dE [GeV cm−2 sec−1 sr−1 ] form.

34

1PeV
E Dep.
thres.

−6

−2
2

νµ

10PeV E dep. thres.

−1−1(log) j(E)E [GeV cm sec sr ]

−5

TD

−7

GZK

−8
Z−burst
GZK

−9

−10
6

8

10

12

(log) j(E)E2 [GeV cm−2 sec−1 sr−1 ](log) Primary Energy [GeV]−5
1PeV
E Dep.
thres.

−6ντ

10PeV E Dep. thres.

TD
−7
GZK

−8

GZK

Z−burst
−9

−10
6

8

10

12

(log) Primary Energy [GeV]

Figure 4.3: The IceCube sensitivities on the EHE neutrino fluxes. The upper figure shows the ν µ sensitivity,
and the lower figure shows the ντ sensitivity. The 90% C.L. upper limits of the primary neutrino flux by a
1km2 detection area with 10 years operation for the threshold energies with 1 PeV and 10 PeV are plotted.
The solid blue lines show the limits obtained in the case of the incident particle is charged leptons. Whilethe dashed lines show the limits by the case of the charged leptons are produced via CC inside the detector.
The other lines indicate the various model predictions given by, GZK ([16] with (m, Zmax )=(4, 4) for the
red line, [20] with (m, Zmax , γ, Emax )=(5, 2, 2, 10 ZeV) for the magenta line), TD [19] and Z-burst [18].
35

−1

−1

(log) j(E)E [GeV cm sec sr ]−5
1PeV
E Dep.
thres.

2−2−6

νµ + ντ

10PeV E Dep. thres.

TD

−7
GZK
−8

Z−burst
GZK

−9

−10
6

8

10

12

(log) Primary Energy [GeV]
Figure 4.4: The IceCube sensitivities on the EHE neutrino fluxes. The blue lines show the sensitivities
which include both ντ and νµ contributions. The solid thick lines show the limits obtained by the case of
the incident particle is charged leptons. While the dashed lines show the limits by the case of the charged
leptons are produced via CC inside the detector. Same as Figure 4.3, the various models predictions arealso plotted.

36

1 PeV thres.
10 PeV thres.

νµ3.3 × 10−8
3.7 × 10−8ντ
3.0 × 10−8
4.6 × 10−8

νµ + ν τ
3.2 × 10−8
4.1 × 10−8

νµ + ντ /flavor
1.6 × 10−8
2.1 × 10−8

Table 4.2: 90% C.L. upper limit for EHE neutrino fluxes, which are given by E 2 dF/dE [GeV cm−2 sec−1
sr−1 ] form are shown with the primary EHE neutrinos’ energies of 109 GeV. Rows are threshold energy
deposit and columns are primary neutrinos. The threshold energy with 1 PeV and 10 PeV are shown.

37

Chapter 5

Summary and Outlook
5.1

Summary

We developed the MonteCarlo/Numerical simulator “JULIeT” for transporting the neutrinos and charged
leptons for the IceCube.
The JULIeT was developed in order to understand both the EHE particles propagation in the earth
and in the detector. As fundamental checks of JULIeT, we calculated the distributions of energy deposit
in 1km ice for µ and τ . They shows the fluctuation which are affected by stochastic interactions such as
bremsstrahlung, photo-nuclear and so on. In the case of τ , especially, the contribution of decay is not
negligible in lower energy range. The dependence of energy deposit on the primary particle’s energy also
shows that the contribution of τ decay is sizable.
The JULIeT still has a problem, arisen from the uncertainty of the cross section of the photo-nuclear
interaction. The instability of this cross section adversely affect on all the results we obtained.
We calculated the propagation in the earth by numerically solving transport equations. We obtained
the incident fluxes of these particle to the detector by transporting the primary EHE neutrinos given byvarious models. While, by using the Monte-Carlo methods, we calculated the energy deposit distributionof the particles in the detector made of ice.
By combining those results, we obtained the energy deposit distribution of the neutrinos and charged
leptons produced by EHE cosmic neutrinos. The event rate for given EHE neutrino fluxes was obtained
to evaluate the IceCube sensitivity. The 90% C.L. upper limits for the primary EHE neutrino fluxes, i.e.,(2.3 events)/(10 years)/(1km2 detectionarea)/(energy decade), for the threshold energies with 1 PeV and
10 PeV were estimated as a function of primary neutrino energy. The resultant upper limits for 10 PeV
threshold would be placed at E 2 dF/dE ' 3.7 × 10−8 [GeV cm−2 sec −1 sr −1 ] for νµ , and 4.6 × 10−8 [GeV
cm−2 sec −1 sr −1 ] for ντ with 109 GeV. We also calculated the 90% C.L. upper limits for the sum of νµand ντ origin events; 2.1 × 10−8 [GeV cm−2 sec −1 sr −1 flavor−1] with 109 GeV threshold.

5.2

Future Work

The highest-priority issue for the progress of the JULIeT is the improvement of the parameterization ofdifferential cross sections of photo-nuclear interaction. The instability of the photo-nuclear cross section
adversely affect on all the results we obtained.
The JULIeT in the current version is not able to simulate cascade development. Instead, it sim-38

ply records primary energy of the electron(or photon)/hadron (i.e. primary energy of the electromagnetic/hadronic cascades) and its generated location, i.e. it sees the cascades as point like. In order tounderstand the Cherenkov light from cascades with accuracy, the JULIeT should be able to simulate thecascade development with Monte-Carlo methods.
The detector defined in the current version of the JULIeT is just an ice, in which the DOMs have notinstrumented yet. It needs a class to define the geometry of DOMs in the points package of the JULIeT.

39

Bibliography[1] K. Greisen, Phys. Rev. Lett. 16, 748, 1966.
[2] IceCube Collaboration, astro-ph/0305196, (2003).
[3] R. Gandhi, C. Quigg, M. H. Reno and I. Sarcevic, Phys. Rev. D 58, 093009, 1998.
[4] H. L. Lai et al., Phys. Rev. D 55, 1280, 1997.
[5] P. Miočinović, Ph.D. thesis, (2001).
[6] R. P. Kokoulin and A. A. Petrukhin, Proc. of the 12th ICRC, Vol.6 p.A2436, (1971).
[7] S. R. Kelner, R. P. Kokoulin and A. A. Petrukhin, Phys. At. Nucl. 62, 1894, (1999).
[8] Yu. M. Andreev, L. B. Bezrukov and E. V. Bugaev, Phys. At. Nucl. 57, 2006, (1994).
[9] S. R. Kelner, R. P. Kokoulin and A. A. Petrukhin, Phys. At. Nucl. 60, 576, (1997).
[10] I. A. Sokalski, E. V. Bugaev and S. I. Klimushin, Phys. Rev. D 64, 074015, (2001).
[11] S. Yoshida, Ultra-High Energy Particle Astrophysics, (Nova Science Publishers, Inc., New York, U.S.,
2003).
[12] S. I. Dutta, M. H. Reno, I. Sarcevic, Phys. Rev. D 62, 123001, (2000).
[13] S. I. Dutta, M. H. Reno, I. Sarcevic and D. Seckel, Phys. Rev. D 63, 094020, (2001).
[14] H. Abramowicz and A. Levy, hep-ph/9712415, (1997).
[15] S. Yoshida, R. Ishibashi and H. Miyamoto, astro-ph/0312078, accepted for publication in Phys. Rev.
D 69, (2004).
[16] S.Yoshida et al., The Astrophysical Journal, 479, 547-559, (1997).
[17] S. Yoshida and M. Teshima, Prog. Theor. Phys., 89, 833, (1992).
[18] S.Yoshida, G. Sigl and S. Lee, Phys. Rev. Lett., 81, 5505, (1998).
[19] G.Sigl et al., Phys. Rev. Lett., D 59, 043504, (1998).
[20] O. E. Kalahsev, G. Sigl et al., Phys. Rev. Lett., D 66, 063004, (2002).
[21] V. Agrawal et al., Phys. Rev. Rev., D 55, 1314, (1996).
[22] D. Heck et al., Report FZKA 6019, (1998).

40