Search for Ultra High Energy Photons
with the Hybrid Detectorof the Telescope Array Experiment(テレスコープアレイ実験ハイブリッド検出器による極高エネルギーガンマ線探索)

理学研究科
数物系専攻

平成 26 年度
Katsuya Yamazaki(山﨑 勝也)

Search for Ultra High Energy Photons with theHybrid Detector of the Telescope Array
ExperimentKatsuya Yamazaki

Abstract
After the discovery of cosmic rays many of the features of cosmic rays, such as theenergy spectrum, primary mass composition, and anisotropies are studied in thewide energy range from 109 eV to beyond 1020 eV. In particular, recent progressesof cosmic ray studies in the energy range above 1018 eV, called ultra high energy
cosmic rays (UHECRs), are remarkable. However, the origins of UHECRs are still
unknown. This is one of the largest problem in the area of UHECR study.
Currently, there are two major experiments for observing and studying UHECRs,
that is the Telescope Array (TA) experiment in the northern and Pierre Auger
Observatory (PAO) in the southern hemisphere. These two observatories publish theresults of the energy spectrum, primary mass composition, and anisotropy studies.
However, the nature of UHECRs has not been fully revealed yet, especially origins
of UHECRs have not been discovered.
In this thesis I will report the studies for the nature of extensive air showers
induced by primary photons and the observational results of searching for UHE
photons as another approach to understand the nature of UHECRs. The flux of UHE
photons depends on the major composition of UHECRs due to their generation andpropagation mechanisms. Furthermore, they arrive at the Earth with no magnetic
deflection. Thus, the detection of UHE photons and determination of their flux canbe critical clues to reveal the nature and origins UHECRs.
TA is a hybrid detector consisting of an array of scintillator detectors and fluorescence detector stations. This hybrid detector measures incoming cosmic ray induced
air showers with about 700 km2 eﬀective detection area in the central Utah. We
searched for UHE photons using the TA detectors in the hybrid mode. The hybrid
data set used in this thesis were accumulated in the period from March 2008 to July2013.
The depth of shower maximum (Xmax ) is used to select photon-like events since
photon induced air showers are expected to have significantly larger Xmax in theatmosphere than hadroni induced showers. In order to determine a selection criterion
a full Monte Carlo (MC) simulation including air shower simulations and hybrid

iii
detector simulations are used.
As a result of UHE photon search we found 24, 11, 3, and 0 photon-like eventsin the energy range greater than 2, 3, 5, and 10 EeV, respectively. The number ofphoton-like events found in the hybrid data set are consistent with the number of
expected false selected events, which are caused by primary protons. Therefore, thephoton fraction upper limits with 95% confidence level are derived from the numbersof photon-like events. In order to calculate the upper limits the detection bias ofprimary photons and primary protons and the selection eﬃciency are calculatedwith the MC simulations. The systematic uncertainty on the photon fractions whichcomes from the uncertainty on the energy determination with hybrid analysis are also
considered. Finally, the upper limits of photon fractions are derived as 9.4, 9.0, 9.8,
and 26.6% in the energy range greater than 2, 3, 5, and 10 EeV, respectively. These
upper limits are the first result which measured by the hybrid detection techniquein the northern hemisphere.
As a result of the upper limits the predicted photon fraction with top-down
models which considered in this thesis are not constrained, but in consideration of
combination with the TA surface detector result, super heavy dark matter models
and topological defect model of the UHECR generation are constrained with 95%
confidence level, and Z-Burst model are survived with these upper limits in thenorther hemisphere. In addition, the result ensure the uncertainties caused by photon contaminations on other analyses, such as composition analysis using average ofthe Xmax , are reasonably small.

Acknowledgements
I wish to express my sincere appreciation to my supervisor Dr. Shoichi Ogio1 , Associate Professor, for the opportunity to be part of the Telescope Array collaboration,and for his many advices and supports for my thesis. I learned many things fromhim. I also would like to oﬀer my great appreciation to Dr. Yoshio Hayashi1 , Professor, and Dr. Koji Kojima1 , Visiting Professor, they has given me many helpfuladvices and discussions.
I would like oﬀer my appreciation to Dr. Masaki Fukushima2 , Professor, and Dr.
Hiroyuki Sagawa2 , Associate Professor for the great management of the Telescope
Array experiment and for the helpful discussions. I also express my appreciation to
Dr. Gordon B. Thomson3 , Professor, Dr. John N. Matthews3 , Research Professor,
for their kindly support at the US.
I wish to oﬀer my deep gratitude to Dr. Toshihiro Fujii4 , Dr. Daisuke Ikeda2 ,and Dr. Yuichiro Tameda2 . They have been done the pioneer work on developing
the analysis software package for fluorescence detectors, and also they have given
me great helps and advices. I also would like to express my deep gratitude to Dr.
Benjamin Stokes3 for helpful technical advices and helps.
Telescope Array is supported by many staﬀ members that help with issues bothtechnical and administrative work. Without their support neither this work nor any
other would be possible. These people include Mr. Gary McDonough, Mr. Frank
Misak, Mr. Jeremy Smith, Dr. Robert Cady, Mr. Stan Thomas, Mr. Dave Barr,and Ms. Cindy Barr.
I would like to thank the members of the Telescope Array collaboration for making the chance to participate in the grate collaboration, especially I would like tothank Dr. Kazumasa Kawata2 , Project Assistant Professor, Dr. Eiji Kido2 , Project
Assistant Professor, Dr. Toshiyuki Nonaka2 , Assistant Professor, Dr. Takeshi
Okuda5 Project Assistant Professor,, Dr. Nobuyuki Sakurai2 , Project Assistant Professor, Dr. Tatsunobu Shibata2 , Project Assistant Professor, Dr. Masahiro Takeda2 ,Assistant Professor, Dr. Hisao Tokuno2 , Project Assistant Professor, Dr. Takayuki
Tomida6 , Assistant Professor, Dr. Yoshiki Tsunesada7 , Assistant Professor, and Dr.
Shigeharu Udo8 , Assistant Professor, for many helps and comments.
iv

v

I would like to thank the members of our laboratory including graduated students, Ms. Mayuko Minamino, Mr. Yu Yamashita, Mr. Daisuke Matsumiya, Mr.
Kazuyuki Kuramoto, Mr. Ryo Yamane, Mr. Yasuhisa Yoneda, Mr. Hyouei Nanpei,
Mr. Takashi Goto, Mr. Shogo Kobayashi, Ms. Ayaka Mori, Mr. Junya Kobayashi,
Mr. Atsushi Terano, Mr. Yasuhiro Yamaguchi, Mr. Shogo Konishi, Mr. Shingo
Wachi, Mr. Yoshiki Nishimoto, Ms. Mika Yamamoto, Mr. Shoichi Kishigami,
Mr. Ryota Onogi, Mr. Yuichi Takahashi for their pleasant communications andencouragements.
Finally, I am deeply grateful to my family, especially Yasuo Yamazaki and Miki
Yamazki for numerous supports and encouragements, and also I would like to thankNaho Sasaki for kindly encouragements.
1 Graduate School of Science, Osaka City University
2 Institute for Cosmic Ray Research, University of Tokyo
3 Physics and Astronomy, University of Utah
4 Kavli Institute for Cosmological Physics, University of Chicago
5 Department of Physical Sciences, Ritsumeikan University
6 Faculty of Engineering, Department of Computer Science&Engineering, Shinshu Univer-sity
7 Graduate School of Science and Engineering, Tokyo Institute of Technology
8 Faculty of Engineering, Kanagawa University

The Telescope Array experiment is supported by the Japan Society for the Promotionof Science through Grants-in-Aid for Scientific Research on Specially Promoted Research(21000002) “Extreme Phenomena in the Universe Explored by Highest Energy Cosmic
Rays”, Grant-in-Aid for Scientific Research(B) (20340057) “Calibrations of atmosphericfluorescence detectors of the Telescope Array experiment using a portable UV laser”,
and the Inter-University Research Program of the Institute for Cosmic Ray Research; by
the U.S. National Science Foundation awards PHY-0307098, PHY-0601915, PHY-0703893,
PHY-0758342, PHY-0848320, PHY-1069280, and PHY-1069286 (Utah) and PHY-0649681(Rutgers); by the National Research Foundation of Korea (2006-0050031, 2007-0056005,
2007-0093860, 2010-0011378, 2010-0028071, R32-10130); by the Russian Academy of Sciences, RFBR grants 10-02-01406a and 11-02-01528a (INR), IISN project No. 4.4509.10
and Belgian Science Policy under IUAP VI/11 (ULB). The foundations of Dr. Ezekiel R.
and Edna Wattis Dumke, Willard L. Eccles and the George S. and Dolores Dore Eccles
all helped with generous donations. The State of Utah supported the project through its
Economic Development Board, and the University of Utah through the Oﬃce of the Vice
President for Research. The experimental site became available through the cooperationof the Utah School and Institutional Trust Lands Administration (SITLA), U.S. Bureauof Land Management and the U.S. Air Force. We also wish to thank the people and theoﬃcials of Millard County, Utah, for their steadfast and warm support. We gratefully

vi

acknowledge the contributions from the technical staﬀs of our home institutions as well as
the University of Utah Center for High Performance Computing (CHPC) and computing
system at Institute for Cosmic Ray Research (ICRR) of University of Tokyo.
This work is supported by the Japan Society for the Promotion of Science (JSPS)
through Grant-in-Aid for JSPS Fellows and by the Japan Student Services Organization(JASSO) through Student Exchange Support Program (Scholarship for Short-Term Study
Abroad).

Contents
Abstractii

Acknowledgements

iv

1 Introduction

1

2 Cosmic rays
2.1 Overview of cosmic ray physics . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Energy spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 Chemical composition . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.3 Anisotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.4 Acceleration mechanisms . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Ultra high energy cosmic rays . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Propagations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Origins of UHECRs . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Extensive air showers . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.4 Detection techniques . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.5 Review of the recent results . . . . . . . . . . . . . . . . . . . . . . .

3
3
3
6
6
8
10
10
14
15
19
20

3 Ultra high energy photons
3.1 Production mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Propagations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Flux predictions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Photon-induced air showers . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1 Pre-shower eﬀect . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 Landau-Pomeranchuk-Migdal eﬀect . . . . . . . . . . . . . . . . . . .
3.4.3 Features of photon-induce EASs . . . . . . . . . . . . . . . . . . . .

29
29
30
32
32
33
33
34

4 The Telescope Array experiment
4.1 Surface detectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Detector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.2 Trigger and data acquisition system . . . . . . . . . . . . . . . . . .
4.2 Fluorescence detectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38
38
38
40
42

vii

CONTENTS

4.2.1
4.2.2
4.2.3
4.2.4

viii

Detector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Trigger and data acquisition system . . . . . . . . . . . . . . . . . .
Atmospheric monitoring . . . . . . . . . . . . . . . . . . . . . . . . .
Cloud monitoring . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42
44
47
47

5 Monte Carlo simulation and shower reconstruction procedures of thehybrid analysis
50
5.1 Monte Carlo simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.1.1 Extensive air shower simulations . . . . . . . . . . . . . . . . . . . . 50
5.1.2 Detector simulations for surface detectors . . . . . . . . . . . . . . . 52
5.1.3 Detector simulations for fluorescence detectors . . . . . . . . . . . . 55
5.2 Event reconstructions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2.1 PMT selections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2.2 Geometrical reconstruction . . . . . . . . . . . . . . . . . . . . . . . 57
5.2.3 Longitudinal profile reconstruction . . . . . . . . . . . . . . . . . . . 59
5.2.4 Quality cuts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.3 Monte Carlo study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.3.1 Monte Carlo data set . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.3.2 Resolutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.4 Data and Monte Carlo comparisons . . . . . . . . . . . . . . . . . . . . . . . 62
5.4.1 Observed data set . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.4.2 Comparisons of observable parameters . . . . . . . . . . . . . . . . . 65
6 Search for ultra high energy photons
6.1 Acceptance bias and data reduction . . . . . . . . . . . . . . . . . . . . . .
6.2 Selection criterion for photon primary events . . . . . . . . . . . . . . . . .
6.3 Data analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4 Upper limits of the photon fraction . . . . . . . . . . . . . . . . . . . . . . .
6.5 Systematic error from energy uncertainty . . . . . . . . . . . . . . . . . . .

70
70
71
72
75
76

7 Discussion

80

8 Summary and conclusion

82

References

84List of Figures
2.1

The all particle cosmic ray energy spectrum in the broad energy range [1]. .

4

2.2

The chemical abundance of cosmic rays measured at the Earth comparedwith the solar system abundances [8]. Vertical axis is the relative abundance
based on that silicon equal 100. Filled circle shows the chemical abundanceof cosmic rays measured in low energy region (70–280 MeV/nucleon), andopen circle shows the chemical abundance of cosmic rays in high energy
region (1000–2000 MeV/nucleon). Diamond shows the chemical abundanceof the solar system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.3

Acceleration at supernova shock front [11]. . . . . . . . . . . . . . . . . . . . 12

2.4

The attenuation length of UHECR protons in the CMBR photon field [15].
This calculation including the photo-pion production and electron-positronpair production. First dip around 1019 eV is caused by the pair production,and the second dip is created by the photo-pion production. . . . . . . . . . 12

2.5

Energy loss lengths for diﬀerent nuclei [18]. . . . . . . . . . . . . . . . . . . 13

2.6

Projected view of 20 trajectories of proton primaries injecting from a point
source for several energies. The trajectories are plotted until they reach thephysical distance from the source of 40Mpc [19]. . . . . . . . . . . . . . . . 22

2.7

The Hillas plot, which shows the relation between the magnetic field strengthand the typical size of astronomical objects [21]. Thick line and dashed line
indicate the lower limits which can accelerate particles up to 1020 eV witheach shock condition β (written in the figure), i.e., objects drawn in upperside of the line is feasible as candidates of UHECR origins. . . . . . . . . . 23

2.8

An schematic diagram of the basic processes in an EAS. It consists of threemajor components; the electromagnetic cascade, the hadronic component,
and the muonic component . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.9

The Heitler model for the cascading of the electromagnetic component EAS.
Every steps after the propagation of radiation length, λ, the number ofparticles in the electromagnetic cascade doubles, and the amount of energy
per particle decreases in half. . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.10 The relative fluorescence spectrum measured by the FLASH experiment [31]. 25ix

LIST OF FIGURESx

2.11 The energy spectra obtained by AGASA (blue crosses), HiRes-1/HiRes2 (open pink squares and pink circles), Auger (open blue triangles), TA
SD (filled red circles), TA MD monocular (filled green triangles), and TA
BR/LR hybrid (filled black squares) [71]. These spectra multiplied by E 3 .
The systematic uncertainty of the flux scaled by E 3 due to the uncertaintyof the energy of 21% is indicated by arrow . . . . . . . . . . . . . . . . . . . 26
2.12 Upper figure shows the averaged Xmax observed by the TA MD hybridmeasurement [40], and lower figure shows that of the PAO hybrid measurement [41]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.13 The correlation between UHECR arrival directions observed by PAO andthe astronomical objects [42], and this figure is in the galactic coordinate.
The circles are arrival directions of UHECRs with energies above 57 EeV,
and the red points indicate positions of astronomical objects included inthe Veron-Cetty and Veron catalogue. Blue contour shows the exposure of
the PAO surface detectors, and dashed line shows the super galactic plane.

28

2.14 The significance of the intermediate scale anisotropy observed by the TA
SD [45] in the galactic coordinate. Arrival directions of UHECRs observedby the TA SD with energies above 57 EeV are oversampled within 20◦ circle,
and then the significances of the deviation from isotropic distribution arecalculated and plotted. The solid gray line is the super galactic plane. . . . 28

3.1

Typical energy loss lengths of UHE photons, protons, and iron nuclei inthe universe [47]. There are three lines indicating the typical energy loss
lengths for each primary, the primary of each line is indicated beside theeach line. Interactions with infrared (IR), cosmic microwave backgroundradiation (CMBR), and universal radio background (URB) photons are
considered in this calculation. The dashed horizontal line indicates theenergy loss length via adiabatic energy losses with redshift. . . . . . . . . . 31

3.2

Predicted ratios of photons to integrated all particle flux above each energy
as a function of primary energy for AGASA (left) and HiRes (right) [48]energy spectra in percentage, and upper limits from AGASA (A), AGASA
and Yakutsk collaboration (AY), Haverah Park (H), and PAO. The shaded
regions show the range of GZK-photon fractions. The upper bands andlower bands are the minimum and maximum fraction derived from diﬀerent assumptions, respectively. The lines indicate minimal photon fractions
calculated by top-down scenarios, SHDM (blue), Z-burst (ZB, pink), andTD (green) models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

LIST OF FIGURES

xi

3.3

Averaged Xmax obtained from simulations for diﬀerent primaries and fromexperimental data [47]. The lines show predicted averaged Xmax with different high energy hadronic interaction models. The elongation rate forthe photon primaries changes at around 1019 eV due to the LPM eﬀect.
The splitting of photon line shows diﬀerent trajectories of primary photons,
because interactions related UHE photons depend on specific trajectories
through the geomagnetic field. . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.4

Average of Xmax for photons, protons, and irons induced EASs [54]. Low
and high energy interaction models used in this calculation are AIRES andQGSJET. There are five lines for photons, the dots curve indicates the caseof no geomagnetic field. The thick solid lines and dashed lines are calculated
zenith angle of 54.0◦ and 61.6◦ , respectively . . . . . . . . . . . . . . . . . . 37

4.1

A map of the TA experiment showing the location of the FD station andSDs indicated by the green squares and the black squares, respectively. The
blue cross shows the location of the Central Laser Facility, and the orange
circles are that of the communication towers. The arrows show azimuthal
extent of the field of view of each FD station. . . . . . . . . . . . . . . . . . 39

4.2

(Left) the schematic view of inside of the scintillation counter. (Right) theside view of the two layers plastic scintillators. . . . . . . . . . . . . . . . . 41

4.3

Example of FADC count distribution of a SD calculated from level-0 triggered waveform. The hatched histogram is a pedestal distribution. This
distribution is collected for every PMT every 10 minutes. . . . . . . . . . . 41

4.4

The possible hit patterns to generate a level-2 trigger. . . . . . . . . . . . . 43

4.5

(Left) an appearance photograph of the Black Rock Mesa FD station building. (Right) the photograph of one of the pair of upper and lower ringtelescopes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.6

(Left) the typical transmittance of the acrylic plate, UV00, the data points
are the median value for three samples, and errors indicate the diﬀerences
between the median and the otter two samples. (Right) the typical transmittance of the UV filter, BG3 filter, and the data points are the medians
for all sampled filters with the bars corresponding to one standard deviation. 45

4.7

A block diagram of the triggering electronics and data acquisition system
of FD [60]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4.8

A schematic diagram of the complete track finding process [60]. . . . . . . . 46

4.9

A schematic diagram of the complete track finding process to find a partial
track pattern across boundaries of adjoining two cameras [60]. . . . . . . . . 46

4.10 The appearance of the IR camera (left) and the twelve pictures taken by
the IR camera at two diﬀerent times (right). The numbers shown in each
picture indicate the existence of clouds or not by 1/0 for in four separated
regions in each picture. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

LIST OF FIGURES

xii

5.1

Schematic view for geometry of a dethinned particle [64]. Vertex position
is chosen on the trajectory of the weighted particle as described in text. . . 53

5.2

Comparison of secondary electron spectra with and without dethinning fora thinned simulation of a proton-induced EAS with primary energy E =1019.0 eV and primary zenith angle θ = 45◦ [64]. Horizontal axis is the
kinetic energy, and vertical axis is the number of particles. Spectra forthinned simulation are showed by gray, and dethinned spectra are showedby black. The secondary electrons within a region from −30◦ to 30◦ inazimuthal angles with respect to the azimuthal angle of the primary particle
direction, and within a region from lateral distances from 500 m to 1000 m
are tabulated in these figures. The angle showed in each figure is incident
angle with respect to the ground. . . . . . . . . . . . . . . . . . . . . . . . 53

5.3

Overview of the simulated SD configuration in the GEANT4 simulation. It
includes many of the components for realistic simulation, such as scintillators, mounting parts, a stainless-steel box, a battery, a solar panel, wirelessantenna, and so on. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

5.4

The cosmic ray energy spectra for some primary species measured by AMS[66, 67]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

5.5

A digital model of a FD station building in the FD detector simulation. . . 56

5.6

Schematic view of the parameters related to the geometrical reconstruction[71]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5.7

An example of the result of the longitudinal profile reconstruction with the
IMC technique [71]. The horizontal axis is the slant depth along the shower
axis, and the vertical axis is the number of photo-electrons observed by the
FD. The black points show the observed data and the red and the blue
histograms indicate the fluorescence and the scattered photon components
optimally estimated by the IMC. Shower parameters of this event are estimated as follows: zenith angle is 4.0 degrees, azimuthal angle is 313.1
degrees, impact parameter (distance from the station to the shower axis) is
17.7 km, primary energy is 8.1 × 1019 eV, and the Xmax is 756.4 g/cm2 . . . 60

5.8

The histogram of diﬀerence between simulated and reconstructed core position. The horizontal axis is the diﬀerence, and the vertical axis is normalized
entry. The resolution is calculated with region of 68%, and it is 111 m. . . . 63

5.9

The histogram of diﬀerence between simulated and reconstructed open angleof arrival direction. The horizontal axis is the diﬀerence, and the verticalaxis is normalized entry. The resolution is calculated with region of 68%,
and it is 0.7◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

5.10 (Left) the histogram of energy ratio of simulated and reconstructed energy.
The systematic shift of reconstructed energy is about +2% and the resolution is about 9.5%. (Right) the energy dependence of the energy ratio.
Reconstructed energies have small energy dependence. . . . . . . . . . . . 64

LIST OF FIGURES

5.11 (Left) the histogram of diﬀerence between simulated and reconstructed
Xmax . The systematic shift of reconstructed Xmax is about −14 g/cm2 andthe resolution is about 29 g/cm2 . (Right) the energy dependence of the
Xmax . Almost no energy dependence can be found in Xmax reconstruction.

xiii

64

5.12 The Distributions of the number of event. Each MC events are weighted
to fit the energy spectrum to the obtained energy spectrum from the SD
analysis [6]. The proton spectrum is in good agreement with that of real
data without energy of 1020 eV bin. The discrepancy is come from thechoice of the bin size. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.13 The Distributions of reconstructed Xmax . There is systematic shift between
proton MC and real data distribution. The amount of shift is about 20
g/cm2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.14 The distributions of core positions along with the x-axis from west to east
direction centered at the CLF position. The MC data set distribution is ingood agreement with that of the real data set. . . . . . . . . . . . . . . . . . 67
5.15 The distributions of core positions along with the y-axis from south to northdirection centered at the CLF position. The MC data set distribution is ingood agreement with that of the real data set. . . . . . . . . . . . . . . . . . 67
5.16 The distributions of the zenith angle of EASs. The MC data set distribution
is in good agreement with that of the real data set with in statistical errors. 68
5.17 The distributions of the azimuthal angle of EASs. The MC data set distribution is in good agreement with that of the real data set. . . . . . . . . . . 68
5.18 The distributions of the impact parameter of EASs. The MC data set
distribution is in good agreement with that of the real data set. . . . . . . . 69
5.19 The distributions of the ψ angle, which is the angle between the shower axis
and the direction from the FD station to the shower core position on the
SDP (see Figure 5.6). The MC data set distribution is in good agreementwith that of the real data set below 120◦ . Above 120◦ , there is systematic
diﬀerence between the real data set and MC data set. . . . . . . . . . . . . 69
6.1

The survival fraction R. Red and blue points indicate R for protons and thatfor photons, respectively. The lines show the fit results with an empiricalequation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

6.2

The blue dots are Xmax of each photon-induced EAS. The blue open circles
show the averaged Xmax at each energy bin, the size of each bin is 0.1 inlogarithm of energy. The lines show the fitting results of photon averaged
Xmax . There is a bend at 1019.6 eV for averaged Xmax of photons due tothe preshower eﬀect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

LIST OF FIGURES

6.3

6.4

6.5

6.6

6.7

6.8

7.1

The interaction model dependence for Xmax of primary photon showers.
The filled circles indicates the averaged Xmax for each energy range withQGSJETII–03, and the solid line is fitted to the filled circles. The open
circles and dashed line are calculated with EPOS–LHC. These values arecalculated with CORSIKA. . . . . . . . . . . . . . . . . . . . . . . . . . . .
The red dots are reconstructed Xmax for MC simulated events of primary
protons, and blue line is the selection criterion determined by Figure 6.2.
Open circles show the averaged Xmax of protons at each energy bin, the sizeof each bin is 0.1 in logarithm of energy. We can find the contaminationfrom protons to the photon selection criterion. . . . . . . . . . . . . . . . .
Xmax histograms for several energy ranges in 1018.2 eV, 1018.5 eV, 1018.8 eV,
1019.1 eV, 1019.4 eV, and 1019.7 eV from upper to lower figures, respectively.
The red histograms and the blue histograms indicate protons and photons,
respectively. The blue arrow in each figure shows the photon selection
criterion at each energy range. . . . . . . . . . . . . . . . . . . . . . . . . .
Reconstructed Xmax of observed events and the selection criterion. The
circles show data points and the line is the selection criterion described in
Section 6.2. The observed data are passed through the quality cut and thephoton enhance cut. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The photon fraction upper limits. Black arrows show the results of upper
limits from this work. The other arrows show other experimental results:
gray is the TA mono analysis result, red is the TA SD analysis result [74],light green is the Auger hybrid result presented in 2009 [75], green is theAuger hybrid result presented in 2011 [76], light blue is the Auger SD result[77], blue is the Haverah Park [78, 79], yellow is the Yakutsk [80, 81], andpink is the AGASA result [82, 83]. The lines show the predictions fromtop-down models (see [84–86]) . . . . . . . . . . . . . . . . . . . . . . . . . .
The photon fraction upper limits derived from the worst case of upper limits
with shifting energy scale calculations. Notation in this figure is the sameas Figure 6.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xiv

73

73

74

75

77

79

The photon fraction upper limits with enough (10 times) statistics. Notation in this figure is the same as Figure 6.7. Fraction upper limits at lower
energies are higher than that of 10 EeV, because separation of average Xmax
between protons and photons is worse at lower energy ranges. . . . . . . . 81

List of Tables
5.1
5.2
5.3

6.1
6.2
6.3

A list of most important parameters used for the CORSIKA simulations inthis work. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
Parameters for making the MC data set. . . . . . . . . . . . . . . . . . . . . 61
The statistics of the MC data for protons and photons. The parent CORSIKA showers are repeatedly used for the simulations (see Section 5.1.2). . 62
The results of the photon fraction upper limits with 95% confidence level.
See text for details. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
The systematic uncertainties on energy determination [71]. . . . . . . . . . 78
The photon fraction upper limits calculated with shifting energy scales. . . 78

xv

Chapter 1
Introduction
Cosmic rays are relativistic energy particles traveling through the universe, some of theseparticles reach the Earth. They were firstly discovered by V. F. Hess in 1912. After the
discovery, cosmic rays have been measured in the broad energy range which is form 108 eV
to beyond 1020 eV.
The energy spectrum of cosmic rays roughly follow a power low with the index of −3,
but some fine structures, such as bends and dips, can be found in the spectrum [1]. These
structures are changes of the power low index, and they are called “knee” (at around 1015.5
eV), “second knee” (at around 1017.5 eV), “ankle” (at around 1018.5 eV) and “cutoﬀ” (ataround 1019.6 eV), respectively. In order to understand the nature of these structures, thechemical composition of cosmic rays in each energy range is also important. Furthermore,
the anisotropies in the distribution of cosmic ray arrival directions is also important keys
to discover origins of cosmic rays and reveal the nature of the cosmic ray physics.
We have not understood all of the nature of cosmic rays, but it is widely believed low
energy (≤ 1015 eV) cosmic rays are accelerated by SNR shocks and propagate throughthe galactic halo. On the other hand, the nature of cosmic rays with energies greater
then 1018 eV, which are called Ultra High Energy Cosmic Rays (UHECRs), has not been
understood yet, because the flux of UHECRs are quite low.
Currently, some experiments continue observing UHECRs. One of them is the Telescope Array (TA) experiment, which is the largest observatory in the northern hemisphere,
and other one is the Pierre Auger Observatory (PAO), which is the largest one in the southern hemisphere. These two observatories have published the results of energy spectrum,
primary composition and anisotropy studies. The reported anisotropies results are almost
consistent, and the energy spectra have a small diﬀerence at the cutoﬀ energy region, but
it is consistent within systematic uncertainties. In contrast with the spectra, the interpretation on the observed results of the primary composition have large discrepancy between
the observatories. TA reported that the primary composition of UHECRs is consistent
with the purely proton at E > 1018.2 eV, on the other hand, the PAO claimed that thedominant component is gradually changing from protons to heavier nuclei at E > 1018.2
eV. This discrepancy is critical to understand the nature of UHECRs. To overcome this1

Chapter 1. Introduction

2

problem, we need to check and compare both of the results carefully. Today, they aretrying that exchange each detectors and compare the detector biases.
In addition to understanding the discrepancy, we should attempt other approaches tounderstand the nature of UHECRs, such as using neutral particles. It is expected that
there are neutral particles, such as photons and neutrinos, in UHECRs. Charged particle
are deflected by the magnetic field in the universe in their propagation. In contrast, neutral
particles are not deflected by the magnetic filed
Observation of neutral particles is one of the approaches for understanding the natureof UHECRs. Fluxes of these particles can depend on the primary composition of the
UHECRs because of their generation and propagation mechanism. Furthermore, they
arrive at the earth with no deflection by the magnetic fields in the universe. Especially for
the UHE photon, the arrival directions can be a direct key to search for UHECR sources.
Arrival directions of these photons directory point to their sources, and source candidatesof these photons are limited to nearby sources, because the mean free path of such highenergy photon is constrained by the interaction with the cosmic microwave backgroundradiation (CMBR). In addition to that, if UHE photons are generated with transientphenomena, such as GRB, simultaneous detection with multi-wavelength experiments is
expected. That can be a evidence of discovery of UHECR sources. Thus, UHE photons
can be a “smoking-gun” of UHECR sources.
In this thesis, Chapter 2 gives a overview of cosmic rays, UHECRs physics and review
the recent results. In Chapter 4, the TA detectors are reviewed. Monte Carlo simulationand the reconstruction procedures are given in Chapter 5. Chapter 6, 7 and 8 give analyses
results, discussion and conclusion.

Chapter 2
Cosmic rays
After the discovery of cosmic rays, they have been measured and studied over a long
period of time. Today, a great deal is known about cosmic ray physics. These knowledges
obtained from many studies are reviewed in this chapter.
This chapter consists of two sections. First, in Section 2.1, overview of cosmic ray
physics in the broad energy range is presented. Second, in Section 2.2, I will review thephysics related to UHECRs, which is important topics in the field of cosmic ray physics
and is the main topic of this thesis.

2.1

Overview of cosmic ray physics

There are three major approaches to cosmic ray studies, such as measurement of energy
spectrum, chemical composition and anisotropy in cosmic ray arrival directions. In thissection, I describe an overview of these three topics, and additionally, I introduce cosmic
ray acceleration models and theoretical discussion about origins of cosmic rays.

2.1.1

Energy spectrum

A cosmic ray energy spectrum for a wide energy range is shown in Figure 2.1. The energy
spectrum with energies below about 1010 eV is flatten than that of higher energy range
due to the shielding eﬀect with the magnetic field of the heliosphere. In contrast to that,above this energy range, the spectrum follows a power low with an index of approximately−3, but there are fine structures called “knee”, “2nd knee”, “ankle” and “cutoﬀ”.
The knee is at around energy of 1015.5 eV. Below this energy, the power low index
is −2.7, and above the energy it is −3. One of possible explanations for this change ofthe spectral index is that this is caused by the change of dominant origins or acceleration
mechanisms. KASCADE and other experiments reported that the dominant chemicalcomposition changes from light to heavy component with increasing energy around theknee region [3]. In common electromagnetic acceleration mechanisms, the maximum accelerated energy depends on the charge of particle. Therefore, the experimental result32.1. Overview of cosmic ray physics

4

Figure 2.1: The all particle cosmic ray energy spectrum in the broad energy range [1].

2.1. Overview of cosmic ray physics

5

supports the model in which a dominant mechanism reaches its acceleration limit withincreasing energy at the knee. On the other hand, there is another explanation for theknee, which is caused with the energy and charge dependent leakage of cosmic rays fromthe confinement volume of our galaxy. The length scale of the motion of the relativistic
charged particle with energy of E and with charge of q = Ze (e : the elementary charge)
moving in the magnetic field B is equivalent to an order of the its Larmor radius [2],rL =E
E [PeV]= (1.08 [pc]).
qcB
Z (B [µG])(2.1.1)

Then, for a particle with energy of 1015 eV in the galactic magnetic field assumed tobe 3 µG, the Larmor radius is 0.36 pc. This is enough smaller than the thickness ofthe galactic disk (about 100 pc), but above this energy the Larmor radius becomes morelarger. Therefore confinement of cosmic rays in the galaxy becomes less eﬀective withincreasing energy, and the spectrum appears to be steep.
This is also a natural explanation of the 2nd knee which is small steepening of thespectral shape at energy around 1017.5 eV. The maximum energy of the acceleration isdepend on an atomic number Z. Thus the acceleration limit of nuclei with charge of Ze
is expected to be
Emax,Z = ZEmax,p ,(2.1.2)where Emax,p is the acceleration limit of protons. Therefore, it can be understood that 2nd
knee is caused by the acceleration limit of heavy nuclei. However, the chemical compositionat the 2nd knee region is poorly probed, and at this time this model is not confirmed. Near
future, the low energy extensions of PAO, which are called HEAT and AMIGA, and thelow energy extension of TA, which is called TALE, will provide clues to the solution of
the problem.
The ankle is at around 1018.5 eV, which is hardening of the energy spectrum. The index
changes to −2.7 above this energy. There are also some models to explain this feature.
One of the models predicts that the ankle is caused by the transition from galactic toextragalactic origins of cosmic rays with increasing energy. In the model, we expect atransition from heavy to light compositions at the ankle with increasing energy. Anothermodel is that the ankle is caused by the flux decreasing of extragalactic protons with theenergy loss through electron-positron pair production with cosmic microwave backgroundradiation (CMBR) photons, i.e., p + γ → p + e+ + e− . This model predicts that thedominant composition is proton in lower and higher energy ranges than the ankle.
If the dominant composition remains proton up to 1020 eV, a cutoﬀ is expected at1019.6 eV, and it is produced through photo-pion productions with CMBR photons,
p + γ → ∆+1232 →{
π 0 + p,π + + n.

(2.1.3)This mechanism had been proposed by Greisen, Zatpin and Kuz’min in 1966 for the first2.1. Overview of cosmic ray physics

6time, and then so called “GZK eﬀect” [4] [5]. The cutoﬀ has been observed by boththe TA [6] and PAO [7] experiments. In spite of the fact these two experiments propose
diﬀerent scenarios to explain the cutoﬀ. TA has mentioned that GZK eﬀect is favored to
explain the cutoﬀ, with their result of the proton dominant chemical composition. On theother hand, PAO has claimed that the cutoﬀ is caused by an acceleration limit for heavy
nuclei from extragalactic origins. This conflict has not been resolved yet. The details of
the results of these two experiments are reviewed in Section 2.2.5.

2.1.2

Chemical composition

As discussed in previous section, the chemical composition is a key quantity to examine theoretical theoretical acceleration and propagation models for all the energy ranges.
Below 1014 eV, the flux of primary cosmic rays is so suﬃciently large that we expect toobtain suﬃcient statistics of events from direct observations of primaries with balloonor satellite-borne detectors. Thus the chemical abundance can be measured from direct
observations of primary cosmic rays. Figure 2.2 shows the chemical abundance of cosmic
rays and that of solar system [8]. The cosmic ray composition shows the overabundancesat some elements, such as Li–B and Sc–Mn. This is believed to be caused by the spallation
process of cosmic ray nuclei with the interstellar matter on the way of travel to Earth.
On the other hand, the flux of cosmic rays above 1014 eV is very small, hence wecannot expect to obtain enough statistics with direct measurements. Therefore, such highenergy cosmic rays are observed by measuring secondary cosmic rays. i.e. an Extensive
Air Shower (EAS), which is a complex of electromagnetic and hadronic cascades (detailed
in section 2.2.3). Instead of the direct measurements of the mass and the charge ofprimaries, we use some observable quantities of EAS in order to determine the particle
type of primaries. One of such quantities is the depth of the maximum development of
EAS, which is called Xmax . This quantity reflects the atomic mass number of primaryparticles. The derivation of the dependence of Xmax on atomic mass numbers for UHECRs
is particularly discussed in Section 2.2.3.

2.1.3

Anisotropy

The anisotropy of cosmic ray arrival directions is studied by various methods in order todirectly identify sources of cosmic rays or in order to reveal the distribution of sources.
One of the simple way to examine the anisotropy is the auto-correlation analysis, which issearching for event concentrations on an arrival direction map. Another approach such as
investigation of correlations between the arrival direction of cosmic rays and the positionof known astronomical objects is one of the methods to search for sources. The otherapproach is a comparison between the arrival direction distribution of observed eventsand simulated one, based on the assumption of some source distributions, matter andphoton distributions, and interactions with matter, photons and the magnetic field during
their propagations.

2.1. Overview of cosmic ray physics

7

Figure 2.2: The chemical abundance of cosmic rays measured at the Earth comparedwith the solar system abundances [8]. Vertical axis is the relative abundance basedon that silicon equal 100. Filled circle shows the chemical abundance of cosmic rays
measured in low energy region (70–280 MeV/nucleon), and open circle shows thechemical abundance of cosmic rays in high energy region (1000–2000 MeV/nucleon).
Diamond shows the chemical abundance of the solar system.

2.1. Overview of cosmic ray physics

8

Unfortunately, most of the cosmic rays have a non-zero charge, and there is the magnetic field in the universe. The arrival direction of a cosmic ray is expected to be refracted
and deviated from its source position, and the distribution of cosmic ray arrival directions is smeared. Therefore, we need a number of statistics to reveal the anisotropy ofcosmic ray arrival directions from smeared distribution. However, the magnetic deflectionof UHECRs is smaller than that of lower energy cosmic rays, because it depends on therigidity of particles. Thus, study of the anisotropy for UHECRs is an eﬀective method to
identify their sources.

2.1.4

Acceleration mechanisms

Acceleration mechanisms of cosmic rays are not fully understood, and there are sometheoretical predictions. The most basic mechanism is a stochastic acceleration model ofcharged particles, which was originally proposed by Fermi in 1949 [9], and it is called“Fermi mechanism”.

Fermi mechanism
The original picture of the Fermi acceleration is transferring macroscopic kinetic energy ofa moving magnetized plasma to charged particles. If such energy transfer is repeated manytimes, then the particle is stochastically accelerated to nonthermal energies. We consider
a process that a charged particle interact with magnetized plasma cloud. The particle isaccelerated to proportional to its energy with each interaction. Then, the particle obtain
the energy ∆E = ξE per an encounter, and after n times interactions the energy of theparticle reaches En
En = E0 (1 + ξ)n ,(2.1.4)
where E0 is the initial energy when it is injected into the accelerator. The number ofinteractions needed to reach an energy E is( )E
n = ln/ln(1 + ξ).
(2.1.5)
E0
When the escape probability from the acceleration region is Pesc per interaction, the energy
spectrum of particles accelerated to energies greater than E is
N (≥ E) ∝∞∑(1 − Pesc )m =m=n

and hence
N (≥ E) ∝

1
Pesc

(

)

where
γ ≡ ln1
1 − Pesc(EE0(1 − Pesc )n,Pesc(2.1.6))−γ,/ln(1 + ξ) ≈(2.1.7)Pesc.
ξ

(2.1.8)

2.1. Overview of cosmic ray physics

9

Consequently, Fermi mechanism leads to a power low energy spectrum as the observed
cosmic ray energy spectrum.

Acceleration at supernova
Supernova (SN) explosions eject materials to the interstellar medium (ISM), and they
drive a shock wave with magnetic field where charged particles can be accelerated. The
finite lifetime of the SN shock wave limits the maximum energy per particle that can beachieved with the acceleration at SN. The acceleration rate is
dE
ξE
=
,
dt
Tcycle(2.1.9)where ξ is the fractional energy gain per single acceleration cycle, and Tcycle is the characteristic time for one cycle. In order to estimate the maximum energy Emax by integrating(2.1.9), we need to know Tcycle .
Firstly, let us consider the upstream region, which is the outer region of a SN shock
sphere (Figure 2.3). The particle flow is given by
J = −D∇N + uN,(2.1.10)where D is a diﬀusion coeﬃcient and N is the number density of particles. In the upstream
region, the fluid velocity u1 is negative in the rest frame of the shock wave. There is no
net flow in equilibrium, therebydN
D1= −u1 N,
(2.1.11)dz
and so in the upstream region()
zu1
N (z) = ρCR exp −,(2.1.12)D1
where ρCR is the number density of cosmic rays at the shock front. The total number ofparticles per unit area in the upstream region is ρCR D1 /u1 . The rate per unit area atwhich relativistic cosmic rays cross a shock plane is cρCR /4. Thus, the mean time thata particle is remained in the upstream region is 4D1 /(u1 c). The downstream region issomewhat more complicated to analyse because it is necessary to average the residence
time only over those particles that do not escape. The analysis is shown explicitly by
Druri [12]. According to that, the form of averaged residence time is identical to that in
the upstream region. Therefore,()
4 D1 D2
Tcycle =+.
(2.1.13)c u1
u2
To proceed the calculation, we need to estimate the diﬀusion coeﬃcients. Lagage andCesarsky argued that the diﬀusion length of charged particles cannot be smaller than its

2.2. Ultra high energy cosmic rays

10

Larmor radius rL = E/(ZeB) [13], where Z is the charge of the particle, E is the totalenergy and B is the magnetic field strength. Hence, the minimum diﬀusion coeﬃcient,which gives the maximum acceleration rate is calculated with rL as
Dmin =

rL c
1 Ec
∼
.
3
3 ZeB

(2.1.14)

Therefore, when we substitute D1 = D2 = Dmin and u2 = u1 /4 with a strong shock
assumption,
Tcycle ≥ 20E/(3u1 ZeB),(2.1.15)
Finally, the maximum energy of charged particles in the SN shock wave acceleration isderived from eq. (2.1.9) and (2.1.15) as
Emax ≤3 u1
ZeB(u1 TA ),20 c(2.1.16)

where TA is lifetime of the SN shock wave. For example, the SN which ejects 10 M⊙ at
5 × 106 m/s during TA ∼ 1000 years, into the ISM with the density of 1 cm−3 and withthe magnetic field the order of µG. In such situation, the maximum energy is derived fromeq (2.1.16) as
Emax ≤ Z × 1014 eV.
(2.1.17)
As a result of the rough estimation, SN shock waves can only accelerate cosmic ray protons
up to 1014 eV.

2.2

Ultra high energy cosmic rays

The cosmic ray energy spectrum continues to above 1018 eV and such high energy cosmic
rays are called ultra high energy cosmic rays (UHECRs). The flux of UHECRs is extremely
small reflecting the stochastic nature of unknown acceleration mechanisms to such ultrahigh energies. Moreover, UHECRs are aﬀected by the magnetic fields and the cosmic
microwave background radiation (CMBR) photons in their propagation to the Earth. In
this section, a general overview and recent results of UHECR studies are described.

2.2.1

Propagations

UHECRs interact with CMBR photons and with intergalactic and galactic magnetic
fields in their propagations. The interaction with CMBR photons causes important processes, which are photo-pion production, electron positron pair-production and photodisintegration. These processes are important to understand origins of UHECRs.

2.2. Ultra high energy cosmic rays

11

Photo-pion production and electron-positron pair production
Greisen, Zatsepin and Kuzmin pointed out that the universe is opaque for cosmic rays
which have suﬃciently high energies. This is called GZK eﬀect. When cosmic ray protons
have energies around 1020 eV, the energy of CMBR photons approximately corresponds
to the pion rest mass in the rest frame of protons. In this case, the interaction of an UHE
proton with a CMBR photon creates pions, and as a result the proton loses its energy.
The most important process which contributes the total cross-section of the photo-pion
production is the ∆+ (1232) resonance production, such as
p + γ → ∆+1232 →{
π0 + pπ+ + n.

(2.2.18)

The energy threshold of photons for this interaction in the rest farme of the proton isapproximately equivalent to the rest mass of π. It correspond to the primary protonenergy of 6.79 × 1019 eV for 2.7 K photons.
UHECR protons also lose their energies by the electron-positron pair production,p + γ → p + e+ + e − .

(2.2.19)

Figure 2.4 shows the attenuation length of UHECRs in the CMBR photon field. The
attenuation length for UHECRs energies above 1020 eV is less than 100 Mpc, therefore,
these extremely high energy cosmic rays can only reach us from relatively nearby sources.

Photo-disintegration
UHE nuclei also interact with CMBR photons, and they are broken in fragments, i.e., into
lighter nuclei. This process is called photo-disintegration. This is an important energy
loss mechanism for UHECR nuclei. UHECR nuclei are strongly aﬀected by this process
in the energy range from 10 MeV to 30 MeV in the rest frame of UHECR nuclei.
The attenuation lengths for four species are shown in Figure 2.5 [18]. For the nuclei,
this process is dominant energy loss mechanism in the energy range above 1018.0 eV. In
addition, this mechanism also predicts the cutoﬀ at approximately same energy as the
GZK cutoﬀ. Thus, the mass composition of cosmic rays and the dominant energy loss
mechanism are important to discuss the propagation of UHECRs.

Propagation in the magnetic fields
Most of cosmic rays are charged particles, therefore they are aﬀected and deflected bymagnetic fields in the universe. The Larmor radius of cosmic ray protons with energy of1019 eV is about 3 kpc in the 3 µG interstellar magnetic field. This radius is larger than
the thickness of the galactic disk, thus cosmic rays which have energies above 1019 eV
cannot be confined in our galaxy. Therefore, the hypothesis that origins of UHECRs are2.2. Ultra high energy cosmic rays

12

Figure 2.3: Acceleration at supernova shock front [11].

Figure 2.4: The attenuation length of UHECR protons in the CMBR photon field[15]. This calculation including the photo-pion production and electron-positronpair production. First dip around 1019 eV is caused by the pair production, and thesecond dip is created by the photo-pion production.

2.2. Ultra high energy cosmic rays

Figure 2.5: Energy loss lengths for diﬀerent nuclei [18].

13

2.2. Ultra high energy cosmic rays

14

extragalactic is favored. The strength of the extragalactic magnetic field is measured by
the Faraday rotation techniques, but there are large uncertainties the measurements, andthe order of magnitude is believe to be 1 nG. Figure 2.6 shows simulated trajectories of
UHE protons in the 1 nG magnetic field. Trajectories of 1020 eV is almost straight, andthey can be traced back to origins of cosmic rays.

2.2.2

Origins of UHECRs

Origins of UHECRs have not been resolved yet, and there are many theoretical models for
UHECR origins. They are divided two types of models. One is the “Bottom-up” models,
and the other is the “Top-down” models.

Bottom-up scenarios
The bottom-up models predict that UHECRs are generated through the accelerations of
low energy cosmic rays. In this model, the acceleration of lower energy particles is believe
to be basically performed with the Fermi mechanism by shock waves, thus cadidates of
UHECR origins are limited to a certain type of astronomical objects. For UHECR origins,
their magnetic field strength and the size of the objects are required to be large enoughto accelerate particles up to 1020 eV. The relation between the maximum energy Emax ofan accelerated particle and parameters of an acceleration site is given by
Emax ∝ ZeBLβ,(2.2.20)

where Ze is the charge of the particle, B is the strength of the magnetic field at theacceleration site, L is the size of the acceleration site, and βc is the velocity of the shock
wave at the site. Figure 2.7 shows the relationship between the magnetic field strengthand the typical size of the acceleration region of astronomical objects. This figure shows
candidates of UHECR origins which can accelerate cosmic rays up to 1020 eV. If the sizeof the acceleration site is small, then the magnetic field strength has to be large enoughto accelerate particles to ultra high energies, because the Larmor radius must be smaller
than the acceleration region size in order to eﬃciently confine UHE particles in the region.
Furthermore, these candidates have to be found within the distance of 100 Mpc because
UHE particles lose their energies by the GZK eﬀect or the photo-disintegration process.

Top-down scenarios
In contrast with the bottom-up models, top-down scenarios predict that UHECRs are
generated through the decay or the interaction of energetic particles, for example, theinteraction of the topological defects [22, 23], the decay or annihilation of super heavy
dark matters [25], Z-burst model, which is caused by the resonant production of Z bosonby the interaction between relic neutrinos and UHE cosmic neutrinos [24]. Most of the topdown scenarios predict relatively large abundance of the UHE photons and UHE neutrinos.

2.2. Ultra high energy cosmic rays

15

Therefore, it can be proved when we measure significant fraction of photons in UHECRs
and a large flux of UHE neutrinos.
In addition, there are other possibilities to verify these models. One of the ways isthat probing the anisotropy of cosmic ray arrival directions. If UHECRs are produced bysome sort of dark matter, we expect an excess of the UHECR flux from the galactic centerregion.

2.2.3

Extensive air showers

When primary cosmic rays with suﬃcient energy enter the atmosphere, they interact
with molecules of the atmosphere and produce fluxes of secondary particles. All theseparticles together create cascades known as extensive air showers (EASs). This phenomena
is discussed in detail here. Since the flux of UHECRs is quite small, it is diﬃcult toobtain suﬃcient statistics of events by direct observations with ballon- or satellite-borne
detectors, because it is not realistic to construct such a huge detector which can obtain
suﬃcient statistics of UHECR events at the top of the atmosphere. Therefore, the indirect
measurement for EAS is the main method for observations of UHECRs.
In the early stage of the EAS development, the number of secondary particles increases
exponentially. The number of secondary particles in this cascade finally reaches its maximum, and then the number of secondary particles decreases. The EAS cascade consists ofthree major components; the electromagnetic cascade, the hadronic component and themuonic component. Figure 2.8 shows a schematic diagram of these basic processes.
For the dominant component of UHECR primaries, the initial interaction between
primaries and the atmosphere is the hadronic interaction, and mostly produces mesons,
where the most abundant secondary is pion. Pions with three diﬀerent charge, i.e., π + , π 0 ,and π − , are produced equal probability. The charged pion, which has relatively long life
time, continues producing more hadronic collisions, or it decays into muon and neutrino,as π ± → µ± + ν. The neutral pion which has relatively short life time decays into two
photons, as π 0 → 2γ, and they originate electromagnetic cascades.

Electromagnetic cascade
Electromagnetic cascades are induced by interactions of high energy photons, electrons,
and positrons with atmospheric molecules. One high energy photon generates a pairof electron and positron with the pair production, γ + N → e+ + e− + N, where N is a
nucleus. And a high energy electron and positron emits photons via bremsstrahlung,e± + N → e± + N + γ. As a result, the number of particles increases as long as energy ofparticles are suﬃciently for continuing the cascading process. Figure 2.9 gives the idea
which is proposed by Heitler [26] how the bremsstrahlung and the pair production processes
work in the electromagnetic cascade. It should be noted that the radiation length forbremsstrahlung of electrons and positrons is 36.5g/cm2 in the air, and it is approximately
equal to the attenuation length of high energy photons with the pair production. As the2.2. Ultra high energy cosmic rays

16

air shower which particles have suﬃciently high energy travels one radiation length, thenumber of particles is doubled through these two processes, and the average energy per
particle is halved.
The critical energy is determined from the balance of the two energy loss processes ofelectrons and positrons, i.e., bremsstrahlung, which drives the growth of the number ofparticles, and ionization. The probabilities of these two processes become equivalent atenergy of 84 MeV in the air, and that is called the critical energy Ec . Along with the growthof the number of secondary particles by cascading processes as the shower development,
the primary energy is divided up into the secondary particles, and finally the averagedenergy per secondary particle reaches the critical energy. At this point, the number of thesecondary particles in the cascade reaches a maximum, called Nmax , and the slant depthof the point along with the shower axis is called the depth of shower maximum, Xmax .
As mentioned in Section 2.1.2, Xmax is one of the most important quantities sensitive tothe atomic mass number of primary cosmic rays. In addition, Nmax is a quantity which
sensitive to the energy of primary cosmic rays because it linearly depends on the energiesof the primaries. After the shower maximum, the number of secondary particles decreases
along with the shower development.

Hadronic component
The energy of the hadronic component in EAS is constantly converted into the electromagnetic cascade via the decay of π 0 which are produced in the hadronic interactions.
π ± and π 0 are produced with an equal probability of π ± and π 0 in each hadronic collision. Because of this energy transfer process, the electromagnetic portion for the numberof secondary particles becomes dominant after only a few generations of collisions. The
hadronic component constitute about 1% of the total number of secondary particles in
one EAS.

Muonic component
Most of low energy charged pions from hadronic collisions at the later stage of the air
shower development decay into muons and neutrinos (π ± → µ± + ν). The mean lifetimeof the charged pions is relatively long (2.6 × 10−8 s) because the decay is via the weakinteraction. Most of the muonic component in the shower is generated through these decays. Higher energy muons are generated in the upper atmosphere, where the atmosphericdensity is suﬃciently low such that the charged pions have more of a chance to decay before colliding, and lower energy muons are generated in the lower atmosphere. Secondary
muons and neutrinos undergo hard collisions only very rarely in the atmosphere. Thus,
most of the muons in EASs deeply penetrate the atmosphere, and muons which are notcollide with the atmospheric molecules can arrive at the ground.

2.2. Ultra high energy cosmic rays

17

Longitudinal development

As mentioned above, secondary particles in EAS are dominated by the electromagnetic
component. In this section, let us consider a very simple toy model of the longitudinaldevelopment of a purely electromagnetic cascade caused by a photon with very high energyof E0 . In this model, a suﬃciently high energy photon produces an electron-positron pairhaving the half of the photon energy after the photon has traveled a distance of λ, whichis the radiation length for the pair production process in the air. Moreover, a suﬃciently
high energy electrons or positrons emit one photon by bremsstrahlung with giving a halfof their kinetic energy on average to the photon after traveling a distance of λ. With thesetwo processes, the primary photon creates an electron and a positron in the first step, andin the second step, the number of particles increases to four, i.e., two photons, an electron
and a positron. Each particle has the energy of E0 /4. It can be denoted that after i steps
the number of particles is Ni = 2i , and each particle has the energy of Ei = E0 /2i . Using
the grammage distance, x, instead of i, N and E are written in the form
N (x) = 2x/λ ,(2.2.21)and,E0.
(2.2.22)2x/λThis process continues until the average energy per particle reaches the critical energy, Ec .
At the Xmax , the number of particles reaches the maximum and N (Xmax ) can be writtenas
N (Xmax ) = Nmax = E0 /Ec .
(2.2.23)E(x) =

Using the relation (2.2.21) we obtain Xmax as,
Xmax = λln(E0 /Ec ).
ln2(2.2.24)

Therefore, Xmax has the logarithmic dependence on E0 , and dXmax /d(ln E0 ) is called
called the elongation rate.
We discussed about the longitudinal development of air showers induced by primary
photons. Now, let us revisit the longitudinal development for primary protons and heavy
nuclei. Primary protons initially create the hadronic component of showers, and most of
their energies are transported to the electromagnetic component. For primary nuclei, it
can be seen as a superposition of low energy nucleons, each of which has the average energyof E0 /A, where E0 is the energy and A is the mass number of the primary nucleus. Thus,
the air shower induced by the nucleus of A is considered to consist of A subshowers withthe primary energy of E0 /A. Therefore, in the toy model of the longitudinal development
described above the primary energy is E0 /A, except E0 , for the nucleus A. Then, the2.2. Ultra high energy cosmic rays

18

mass number dependence for Xmax is expressed as,
Xmax ∝ ln

E0.
A(2.2.25)

Thus, Xmax can be used for determining the composition of primary cosmic rays.
In addition, the shape of the distribution of Xmax also depends on the composition ofprimary cosmic rays. As described above, the EAS induced by nuclei of A is considered toconsist of A subshowers. The bunch of these subshowers results in an averaging eﬀect on
the single EAS. Therefore, EASs which is induced by heavy nuclei have smaller fluctuations
than the EASs induced by protons.
Though Heitler model gives the Xmax dependency for the chemical composition ofprimary cosmic rays, there is another treatment of longitudinal development. The GaisserHillas function is generally used for fitting the development of charged particles in EAS[11, 27],(
N (x) = Nmax

x − X0
Xmax − X0)(Xmax −X0 )/λ(exp

Xmax − xλ),(2.2.26)

where, X0 is the first interaction point of primary particles, and λ is the attenuation
length, typically λ = 70 g/cm2 . The variable x used above is the grammage distance that
an EAS has traveled through the atmosphere. That is called slant depth which unit is[g/cm2 ].

Lateral distribution
Particles in an EAS have a lateral spread. Lateral distributions of electromagnetic cascades
are parametrized by Nishimura, Kamata and Greisen called NKG function [28],Ne f (r, s),2
rM)s−4.5( )s−2 (rr
f (r, s) = C(s)
+1
,
rM
rM
ρ(r, s) =(2.2.27)(2.2.28)

where r is the distance from the shower axis of EAS, rM is the Moliére unit, s is the shower
age, C(s) is a normalization factor, and Ne is the number of particles at any point in theshower.
In order to represent lateral distributions for hadron-induced EASs, the TA experimentuses the empirical function established by the AGASA experiment [29]:
(

)
(
)
[
( r )2 ]−0.6
r −1.2r −(η(θ)−1.2)
ρ = A
,
1+
1+
rM
rM
1000 mη(θ) = 3.97 − 1.79 [sec(θ) − 1] ,(2.2.29)(2.2.30)

where A is a normalization factor, and θ is the zenith angle of the shower axis. This

2.2. Ultra high energy cosmic rays

19

empirical equation includes the zenith angle dependency.

2.2.4

Detection techniques

As mentioned above, in order to study the nature of UHECRs we detect the EAS particles
initiated by incident primaries instead of direct observations of the primaries. In thissection, the techniques for detecting EASs are briefly reviewed.

Particle detection techniques on the ground
The most commonly-used method for the EAS detections is based on the coincident detections of the remnants of EAS particles with a particle detector array at the groundlevel. In general, this type of detectors are called surface detector (SD) arrays. Each SD
records the number of particles and the arrival timings at the ground level. Then, thenumber of particle distribution as a function of the distance from the EAS core and thearrival timing distribution are used to reconstruct the energy and the arrival direction of
EASs. Since the typical size of lateral spread of EAS particles at the ground level is abouta few kilometers in ultra high energies, SDs should be deployed with about one kilometerspacing.
Energies of primary particles are estimated from the comparison between the measured
and the simulated particle distributions. In the simulation we need to assume the chemicalcomposition of primaries. However, the composition is not known exactly in the indirect
measurement energy range. Moreover, hadronic interaction models used in the simulations
are extrapolations of experimentally proved results with accelerator experiments, because
the interaction energies at current particle colliders are significantly lower than that of air–
UHECR interactions. Thus, the primary energy estimation with the SD array detections
has systematic uncertainty by these ambiguity and extrapolations. However, SD arrays
are generally stable observational equipments and are able to achieve large exposure.

Fluorescence detection techniques
In the fluorescence method, we uses the nitrogen fluorescence emission to measure thelongitudinal development of EASs as they move through the atmosphere. In other words,
the atmosphere plays a role as a huge scintillator. Since the number of secondary particles
in a EAS induced by UHECRs is very large, fluorescence photons can be detected byordinary photodetectors, such as using PMTs and mirrors. This is one of the detection
techniques for UHECRs measurement used as fluorescence detectors (FDs) by TA, PAO
and other experiments.
Charged particles with suﬃcient energy excite atmospheric molecules, such as nitrogens
or oxygens, and these excited molecules emit fluorescence photons, in ultra-violet andoptical wavelength band, mainly from 300 nm to 400 nm by nitrogens. Fluorescenceemissions are isotropic, and so the emitted photons can be measured from the side of
the EAS axis, and, in particular, showers can be detected from a long distance when the2.2. Ultra high energy cosmic rays

20

shower size is suﬃciently large. There is no need for EASs to point toward the detector
or even to strike it directly as it is required for SD arrays.
Fluorescence yield is the measurement, which provides the number of fluorescencephotons per energy deposit in the atmosphere. The fluorescence yield depends mostly onenergy deposit, but has a weak dependence on air pressure, temperature and humidity.
The fluorescence yield spectrum was measured by multiple experiments. In this thesis,
the FLASH spectrum [31] and the Kakimoto absolute fluorescence yield [30] are used for
simulations and reconstructions of EAS profiles. The measured yield spectrum is shownin Figure 2.10.

2.2.5

Review of the recent results

UHECRs have been observed by various observatories, firstly it was reported in 1963 [32].
After this report various experiments, such as Haverah Park [33], Yakutsk [34], Fly’s
Eye [35], AGASA [36], HiRes [37, 38], TA, and PAO have also observed UHECRs until
now. Recently, progress of UHECR studies is remarkable, however, the nature of UHECR
origins is not revealed yet. In this subsection, a brief review of these recent progress
focused on TA and PAO results related to three topics, the energy spectrum, primary
mass composition, and anisotropies, are described.

Energy spectrum
Figure 2.11 shows the energy spectra measured by several experiments. Although theabsolute values of the energy spectra are diﬀerent within the systematic uncertainties ofeach experiment, these spectra show similar shape without the highest energy edge of the
AGASA spectrum.
According to the recent results of TA and PAO it is established that the existenceof a dip and a cutoﬀ in the energy spectrum around the energy of 1018.7 , 1019.5 eV,
respectively. The energy spectra of these two experiments are consistent within theirsystematic uncertainties, however, there is a systematic diﬀerence in energy scales between
these two experiments. In order to reveal the nature of UHECRs accurate determinationof absolute energies at the spectral bending points is necessary. Therefore, TA and PAO
have started the collaborative work to understand the diﬀerence of their energy scales [39]Mass composition
In the study of mass composition for UHECRs Xmax is used widely, because of the sensitivity in the mass number. Figure 2.12 shows results of mass composition analysis using
Xmax . TA and PAO have diﬀerent statement about the interpretation of UHECR composition. TA has reported the primary mass composition is dominated by protons above
the energy of 1018.2 eV. In contrast, PAO has mentioned that the dominant composition
is changing to heavier component along with increasing energy above the energy of 1018.3
eV. The discrepancy makes it diﬃcult to understand the nature of UHECRs. There are2.2. Ultra high energy cosmic rays

21

two possible explanation, one is that the discrepancy arise from errors of one experimentor both experiments, the other possibility is that the discrepancy arise from the diﬀerenceof UHECR sources which can be observed in the norther or the southern hemisphere. This
has not been concluded yet, and now these two experiments are working to resolve thisproblem.

Anisotropies
In the energy of E > 57 EeV, PAO found correlations of the cosmic ray directions within
a 3.1◦ radius circle centered at nearby AGNs (within 75 Mpc) in the southern sky (Figure2.13) [42]. The correlation becomes weak thereafter in the updated measurement [43].
For the TA experiment, this type of correlation analysis is searched, and TA does notfind significant correlation between UHECR arrival directions observed by the TA SD
with energies above 40 EeV and positions of astronomical object listed in six diﬀerent
catalogs [44].
On the other hand, TA recently reported the intermediate scale anisotropy, hotspot, byoversampling using 20◦ radius circles in the norther hemisphere [45]. Figure 2.14 shows theresult of the intermediate scale anisotropy. This hotspot can be a key to reveal UHECR
origins.
In the anisotropy study, there is also collaborative work between TA and PAO. Then,
now we can search the whole sky anisotropies of UHECR arrival directions [46].

2.2. Ultra high energy cosmic rays

22

Figure 2.6: Projected view of 20 trajectories of proton primaries injecting from apoint source for several energies. The trajectories are plotted until they reach thephysical distance from the source of 40Mpc [19].

2.2. Ultra high energy cosmic rays

23

Figure 2.7: The Hillas plot, which shows the relation between the magnetic fieldstrength and the typical size of astronomical objects [21]. Thick line and dashedline indicate the lower limits which can accelerate particles up to 1020 eV with eachshock condition β (written in the figure), i.e., objects drawn in upper side of theline is feasible as candidates of UHECR origins.

2.2. Ultra high energy cosmic rays

24

Figure 2.8: An schematic diagram of the basic processes in an EAS. It consists ofthree major components; the electromagnetic cascade, the hadronic component, andthe muonic component

Figure 2.9: The Heitler model for the cascading of the electromagnetic component
EAS. Every steps after the propagation of radiation length, λ, the number of particles in the electromagnetic cascade doubles, and the amount of energy per particle
decreases in half.

2.2. Ultra high energy cosmic rays

25

Figure 2.10: The relative fluorescence spectrum measured by the FLASH experiment[31].

2.2. Ultra high energy cosmic rays

26

Figure 2.11: The energy spectra obtained by AGASA (blue crosses), HiRes-1/HiRes2 (open pink squares and pink circles), Auger (open blue triangles), TA SD (filled
red circles), TA MD monocular (filled green triangles), and TA BR/LR hybrid (filled
black squares) [71]. These spectra multiplied by E 3 . The systematic uncertainty of
the flux scaled by E 3 due to the uncertainty of the energy of 21% is indicated byarrow2.2. Ultra high energy cosmic rays

27

Figure 2.12: Upper figure shows the averaged Xmax observed by the TA MD hybridmeasurement [40], and lower figure shows that of the PAO hybrid measurement [41].

2.2. Ultra high energy cosmic rays

28

Figure 2.13: The correlation between UHECR arrival directions observed by PAO
and the astronomical objects [42], and this figure is in the galactic coordinate. The
circles are arrival directions of UHECRs with energies above 57 EeV, and the redpoints indicate positions of astronomical objects included in the Veron-Cetty andVeron catalogue. Blue contour shows the exposure of the PAO surface detectors,
and dashed line shows the super galactic plane.

Figure 2.14: The significance of the intermediate scale anisotropy observed by the
TA SD [45] in the galactic coordinate. Arrival directions of UHECRs observed by the
TA SD with energies above 57 EeV are oversampled within 20◦ circle, and then thesignificances of the deviation from isotropic distribution are calculated and plotted.
The solid gray line is the super galactic plane.

Chapter 3
Ultra high energy photons
Although many features of UHECRs have been found as discussed above, the origins of
UHECRs are still unknown.
Most of UHECRs are believed to be charged particles, so that they are deflected by
the magnetic field. In order to discover sources of UHECRs in this situation, quite large
statistics of UHECR events are necessary. For neutral particles, however, the situation
is diﬀerent. Neutral particles are not deflected by the magnetic field, and they travel ina straight line from the origin to the Earth. Thus, the arrival directions point to their
originated points. In particular for UHE photons, the mean free path suppressed by theinteraction with CMBR photons as discussed in the later section. It is much larger than
the size of our galaxy, but it is less than the cosmological distance scale. The number ofsource candidates of UHE photons is not so large because of their short mean free path, andalso we can obtain position and timing information of the source candidates from opticaland radio observation results. If UHE photons are generated by transient astronomical
phenomena, such as GRBs and SNe, simultaneous detection with these observations are
possibly expected, and it can be a evidence of their sources. Furthermore, the detailed
study of the UHE photon flux is an important key to reveal the origin and the propagationof UHECRs. In particular, top-down scenarios predict higher photon fluxes than that ofbottom-up scenarios.
In this chapter, physics connected to UHE photons are reviewed.

3.1

Production mechanisms

There are many models that predict UHE photon production.
In top-down scenarios, photons are produced from decay or annihilation of exotic
particles such as topological defects (TD) [22, 23] and super heavy dark matter (SHDM)[25] models, and also the Z-burst model [24] predicts UHE photon generation. The TD,
such as cosmic strings, produces extremely heavy particles, and they typically decay intoquarks and leptons. The quarks become hadrons and some leptons are decay. As a result29

3.2. Propagations

30of the interactions, a large amount of photons, neutrinos, leptons, and also nucleons are
generated. In the SHDM model, SHDM particles, such as named as cryptons, wimpzillas,and so on, may decay or annihilate, and then UHE photons are created. The Z-burst
model is based on the annihilation of neutrinos with Z resonance, which occur between
UHE neutrinos and relic background neutrinos. Then, the Z boson decay into protons,
neutrinos, and photons. These top-down scenarios predict the relatively higher dominationof photon component than bottom-up scenarios in the energy region of UHE. Therefore,
the measurement of the UHE photon flux can be a crucial test for these models.
On the other hand, in bottom-up scenarios, UHE photons are mainly produced by
the decay of neutral pions which are generated by GZK mechanism in the propagation of
UHECRs, called GZK-photons. The flux of UHE photons predicted by bottom-up models
is not so large as the top-down scenarios, but the amount of the flux depend on theenergy spectra and environment, such as strength of the magnetic fields, of their sources.
Thus, the GZK-photon flux can be a complementary clue to understand the acceleration
mechanism at UHECR sources.

3.2

Propagations

UHE photons lose their energies through the interaction with background photons as is thecase with UHECRs. In this interaction, UHE photons induce electromagnetic cascades.
There is significant uncertainty for the background radiation spectrum in the universe,
thus the energy loss length of UHE photons has the large uncertainty. Figure 3.1 shows
the calculated energy loss length of photons [47]. Typical energy loss length is the orderof 10 Mpc at the energy of 1019 eV. Therefore, UHE photons can arrive at the Earth fromour galaxy or relatively nearby our galaxy.
It is important that photons are neutral particles, and they are not deflected by themagnetic fields in the universe. Thus, their arrival directions can be directly traced backto their origins. Because of the no deflection, their arrival timing would be correlatedwith other optical observations. If UHE photons are emitted by transient sources, suchas gamma ray bursts and super novae, with other wavelength photons, origins of UHE
photons are possibly identified with consideration for other wavelength observations.
In the case of SHDM models, SHDM particles may be concentrated at the galactic
center region, and UHE photons as decay or annihilation products from SHDM particles
would be observed as anisotropically from the galactic center region. Other top-downparents which are distributed relatively far from the Earth, thus the arrival directions ofobserved photons would have the relatively isotropic distribution. In the case of bottom-up
scenarios, photons are mainly generated by GZK mechanism. Thus, their arrival directions
may not point directly to origins of UHECRs but the distribution of the directions could
show up the informations of sources.

3.2. Propagations

31

Figure 3.1: Typical energy loss lengths of UHE photons, protons, and iron nuclei
in the universe [47]. There are three lines indicating the typical energy loss lengths
for each primary, the primary of each line is indicated beside the each line. Interactions with infrared (IR), cosmic microwave background radiation (CMBR), anduniversal radio background (URB) photons are considered in this calculation. The
dashed horizontal line indicates the energy loss length via adiabatic energy losses
with redshift.

3.3. Flux predictions

3.3

32

Flux predictions

In order to estimate expected photon flux at the Earth, we need to include the energy loss
process in the propagation of UHE photons. For these calculations we need to assume theenergy spectrum of background photons, and this uncertainty causes a systematic error on
the predicted photon flux. Moreover, we also need to assume source information for flux
estimations. For example, in the estimation based on a top-down model, the density andthe lifetime of parent SHDM particles are assumed. On the other hand, in the calculation
based on the bottom-up model the source distribution, the UHECR spectrum, and thecomposition are assumed. These uncertainties can be reduced by assuming and fitting
observed cosmic ray energy spectra.
In the paper by Gelmini et al. [48, 49], various scenarios were used to fit the observed
UHECR energy spectra from AGASA and HiRes. Figure 3.2 shows the predicted photon
fractions for top-down models and GZK-photons. In this calculation, the energy spectra
reported by AGASA and HiRes are used to fit the the simulated energy spectra of UHECRs. There are large diﬀerence between the photon fractions with using AGASA and
HiRes energy spectra. This diﬀerence is caused by the diﬀerence between AGASA and
HiRes energy spectra.
In order to fit the spectra 1019 eV, low energy component from galactic cosmic rays
are also assumed.
In order to fit these energy spectra, diﬀerent parameters are assumed in the study,
the source spectrum, the maximum energy of UHECRs, and the minimal source distance,
and also low energy components from galactic cosmic rays are assumed to fit the spectra
below 1019 eV. These assumed parameters make the major diﬀerence of photon fraction.
In addition, the abundance of radio background and the strength of the extragalacticmagnetic field, which are alter the probability of the GZK interaction and mean free pathof the UHE photons, are assumed to calculate GZK-photon fluxes, and it represent theband width of predictions of GZK-photon fraction.
In the case of top-down scenarios, these parameters are also assumed to fit the AGASA
and HiRes energy spectra, but the lines showed in the figure are indicate minimal case ofphoton fraction with in the assumed parameter spaces in the study.

3.4

Photon-induced air showers

Photons which have suﬃcient energies initiate an almost purely electromagnetic cascadein the atmosphere. In addition to the ordinary cascading processes, it is necessary toconsider additional processes for UHE photons, such as the pre-shower eﬀect [50] andLandau-Pomeranchuk-Migdal (LPM) eﬀect [51]. These processes are reviewed in thissection. Furthermore, some specific and observable features of photon-induced EAS are
discussed later in this section.

3.4. Photon-induced air showers

3.4.1

33

Pre-shower eﬀect

In general, high energy photons interact with a magnetic field, and then they produce
an electron and a positron through the pair production process. Erber showed single
parameter, Γ, which is used for the determination of the threshold condition of this process
[50],
(
)(
)
E
H⊥
m2e c3(3.4.1)
Γ≡,H≡≃ 4.414 × 1013 [G]
crit
me c2
Hcrit
eℏ
where E is the energy of the projectile electron, H⊥ is the strength of the magnetic field
perpendicular to the momentum direction of the primary photon, me c2 is the rest massenergy of electrons, and Hcrit is called critical field strength.
This interaction can be occur even in weak magnetic fields, when the energy of primary
photons is suﬃciently high, such as UHE photons. The probability of this interaction
become non-negligible at Γ > 0.5, and it corresponds to the primary photon energy of∼ 1019 eV. Thus, a significant proportion of photons with energies above 1019 eV converts
to electron-positron pairs in the geomagnetic field, and projectile electron components
emit synchrotron radiation photons. Then, the averaged first interaction point of air
showers reaches far high above the atmosphere, ∼ 1000 km above sea level for the primaryenergy of 1020 eV, as comparing with the mean free path of 100 g/cm2 for the proton–
air inelastic hadronic collations. This phenomenon is called pre-shower eﬀect. Since theErber’s parameter Γ depends on H⊥ , the probability of the pre-shower eﬀect depends
strongly on the trajectories of primary photons in the geomagnetic filed and observationsites.
The typical height of pre-shower interaction is about 1000 km above sea level, andthen subsequent electrons, positrons, and photons are induced in the atmosphere, andeach of them create sub-showers in the atmosphere. The primary energy is divided into
these sub-showers, and each sub-shower develops as lower energy shower than the primary
photon. This bundle of sub-showers is observed as a single shower, because the lateral
spread of the pre-shower particles and the diﬀerence of arrival timings are suﬃciently small
for the current UHECR detectors. As a result, Xmax of the showers, which undergo thepre-shower eﬀect, become small.

3.4.2

Landau-Pomeranchuk-Migdal eﬀect

In a medium, the cross-section of the pair production process is suppressed due to aquantum mechanical interference in very high energies [51–53]. This is called LPM eﬀect.
The reduced cross-section σLPM is approximately derived as√
Eγ ELPMσLPM = σBH,
Ee (Eγ − Ee )(3.4.2)
ELPM ≈ (7.7 [TeV/cm]) × X0 ,

3.4. Photon-induced air showers

34

where σBH is Bethe-Heitler cross-section (≈ 0.51 b in air), Eγ is primary photon energy,
Ee and (Eγ − Ee ) are electron and positron energies emitted by the pair production, X0and ρ are the radiation length and the density of media, respectively. For example, in thereference [52], ELPM ∼ 2.8 × 1017 eV at 300 m a.s.l. and ELPM ∼ 1019 eV at the upper
atmosphere.
It causes delay in the shower development at the initial phase of electromagnetic cascades initiated by UHE photons, because primary photons can penetrate deeper into theatmosphere than if the cross-section is not reduced by the eﬀect. It should be noticed
that the fluctuation of the shower development becomes large, since σLPM decreases withincreasing the atmospheric depth.

3.4.3

Features of photon-induce EASs

In contrast to nucleus-induced EASs, photon-induced EASs suﬀer the additional eﬀects
as reviewed above in this section. Therefore, photon-induced EASs show diﬀerent developments from nucleus-induced EASs. The diﬀerence on the shower developments canbe seen in the diﬀerence on Xmax . Figure 3.3 shows averaged Xmax values for diﬀerent
primaries [47]. Photons have a large slope of the averaged Xmax , called elongation rate
(dXmax /d ln E), and thus their averaged Xmax are significantly larger than that of nucleus
primaries at energies above 1016 eV. At energies around 1019 eV the elongation rate gradually increases due to the LPM eﬀect, and also the Xmax rail of photons has a sharp cutoﬀ
at 1020 eV due to the pre-shower eﬀect. The energy threshold of the pre-shower depend on
the trajectories of EASs as discussed above. Figure 3.4 shows average Xmax dependency
on shower geometries at the TA site [54]. If photon-induced EASs come from the northdirection, the pre-shower eﬀect begins lower energy than that of the south direction.

3.4. Photon-induced air showers

35

Figure 3.2: Predicted ratios of photons to integrated all particle flux above each energy as a function of primary energy for AGASA (left) and HiRes (right) [48] energy
spectra in percentage, and upper limits from AGASA (A), AGASA and Yakutsk collaboration (AY), Haverah Park (H), and PAO. The shaded regions show the rangeof GZK-photon fractions. The upper bands and lower bands are the minimum andmaximum fraction derived from diﬀerent assumptions, respectively. The lines indicate minimal photon fractions calculated by top-down scenarios, SHDM (blue),Z-burst (ZB, pink), and TD (green) models.

3.4. Photon-induced air showers

36

Figure 3.3: Averaged Xmax obtained from simulations for diﬀerent primaries andfrom experimental data [47]. The lines show predicted averaged Xmax with diﬀerent high energy hadronic interaction models. The elongation rate for the photon
primaries changes at around 1019 eV due to the LPM eﬀect. The splitting of photon line shows diﬀerent trajectories of primary photons, because interactions related
UHE photons depend on specific trajectories through the geomagnetic field.

3.4. Photon-induced air showers

37

Figure 3.4: Average of Xmax for photons, protons, and irons induced EASs [54].
Low and high energy interaction models used in this calculation are AIRES andQGSJET. There are five lines for photons, the dots curve indicates the case of no
geomagnetic field. The thick solid lines and dashed lines are calculated zenith angleof 54.0◦ and 61.6◦ , respectively

Chapter 4
The Telescope Array experiment
The Telescope Array (TA) experiment is the largest observatory for UHECR in the northern hemisphere and it is located at 39.3◦ N, 112.9◦ W, and about 1400 m above sea
level, Utah, USA. The experiment is operated by the international collaboration of Japan,
USA, Korea, Russia, and Belgium, since 2003. TA consists of three fluorescence detectors(FDs) [55] and an array of 507 surface detectors (SDs) [56]. Figure 4.1 shows a map of
the TA experiment. The area covered by the SD array is about 680 km2 , and the FDs arelocated at around the SD array viewing over the SD array. The FDs and SDs sample cosmic ray air showers by diﬀerent methods, and each sampled data are used complementary
to reconstruct the air showers.

4.1

Surface detectors

4.1.1

Detector

The air shower array consists of 507 surface detectors (SDs) deployed on a 1.2 km spacing
square grid. The SD for the array is autonomic detector, and each SD consists of two layersof 3 m2 plastic scintillation counters, a data acquisition and control electronics, a GPS
receiver, a wireless LAN interface. Each SD is powered by one solar panel and one deep
cycle buttery. The charging power of the solar panel is 125 W, and the buttery capacity is
100 Ah. Figure 4.2 shows a schematic view of inside of the scintillation counter. Thicknessof the plastic scintillator is 1.2 cm, and a stainless plate is inserted between the two layers
to avoid cross-talk photons. Scintillation photons induced by passages of charged particles
are sampled by wavelength shifting fibers laid on the surface of the scintillator and are led
into a photo multiplier tube (PMT) for each layer. Output of each PMT are successively
digitized by a 12 bit FADC with sampling time of 20 ns.

38

4.1. Surface detectors

39

Figure 4.1: A map of the TA experiment showing the location of the FD stationand SDs indicated by the green squares and the black squares, respectively. The
blue cross shows the location of the Central Laser Facility, and the orange circles
are that of the communication towers. The arrows show azimuthal extent of the
field of view of each FD station.

4.1. Surface detectors

4.1.2

40

Trigger and data acquisition system

There are three communication towers at the edges of the SD array, that provide wirelesscommunication to each SD. This network is used to triggering the SDs and to store thedata. The SD trigger routine consists of two parts. The first part is the local trigger oneach detector, and it makes for an SD electronics to store waveforms to its buﬀer memory.
The second is the array trigger, and it makes for all SDs to send the buﬀered waveforms
the communication towers.
The local trigger has two levels labeled level-0 and level-1 trigger. The level-0 trigger
is generated when the pulse height signals reaches the height equivalent to 0.3 minimum
ionizing particle (MIP) passage, and the level-0 triggered SD stores the waveform including
the triggered signal with 2.56 µs span into the buﬀer memory. The FADC counts equivalent
to the 1 MIP is monitored with this waveform by integrating the bins between −4 binsfrom trigger timing and +8 bins after trigger timing. Figure 4.3 shows an example of
1 MIP histogram. This histogram is collected every 10 minutes for monitoring 1 MIP
variation.
When the signals from the upper and the lower layers coincidently exceed three MIPs
height within 160 ns, the level-1 trigger is generated, and the trigger timing informationis stored. Timings of level-1 triggers are accumulated for one second interval, and theset of the information, which is called trigger table, is sent to the host electronics at thecommunication tower in the following interval.
The accumulated trigger tables are checked by the trigger-decision logic on the hostelectronics. An event trigger, called level-2 trigger, is generated by the host electronics
when three or more adjacent SDs have recorded level-1 triggers within the timing diﬀerenceof 8 µs. Figure 4.4 shows possible hit patterns to generate a level-2 trigger. Then, the hostelectronics distributes the event trigger information to all the SDs. Each SD electronic
searches for the coincident waveform with the level-2 trigger in the buﬀer memory. If
it is found, the SD electronics send it to the host electronics together with its timing
information.
The average trigger rates are about 750 Hz and 30 Hz for level-0 and level-1 trigger,
respectively. Level-2 trigger rate is about 0.01 Hz, and corresponding trigger eﬃciency of
the the SD array is about 50% at the energy of 1018.5 eV.
In addition to the level-1 trigger information, SDs send monitor data to the host electronics every second. The monitor data include the information about the gains of PMTs,
environmental conditions, GPS reception conditions, power generation and consumption
status, etc. Details of the monitor data can be found in the reference [56].

4.1. Surface detectors

41

Figure 4.2: (Left) the schematic view of inside of the scintillation counter. (Right)
the side view of the two layers plastic scintillators.

Figure 4.3: Example of FADC count distribution of a SD calculated from level-0
triggered waveform. The hatched histogram is a pedestal distribution. This distribution is collected for every PMT every 10 minutes.

4.2. Fluorescence detectors

4.2

Fluorescence detectors

4.2.1

Detector

42

The TA experiment has three fluorescence detectors (FDs). One of them consists ofrefurbished HiRes-I mirrors, PMTs, and re-engineered electronics. This FD is located in
the northern part of the TA site which called Middle Drum (MD), station named after
the neighboring mountain (see Figure 4.1). The others are newly designed and produced
detectors. These two FDs are located in the south-east and south-west of the TA site, andthey are called Black Rock Mesa (BR) and Long Ridge (LR) stations, respectively. Dataobtained with BR and LR stations are used for the analyses in this thesis, thus details ofthese two FDs are mainly described in the followings in this section.
BR and LR stations have almost the same detector configuration. Each station has
twelve telescopes, and they are grouped into two groups, called ring, such as the lower
viewing ring and the upper viewing ring. Each station covers the field of view (FOV) of 3◦
to 33◦ in elevation angle and 108◦ in azimuthal angle. The field of view covers the wholearea of the SD array. Figure 4.5 shows the appearances of the BR station buildings andtelescopes in the station.
Each telescope consists of a large spherical mirror, which consists of 18 segments.
The diameter of the mirror is 3.3 m and the curvature is 6.067 m. The imaging camera is
installed on the prime focus of each mirror, and It consists of 256 PMTs. The front windowof the camera box is made of UV transparent acrylic plate, paraglas UV00 by Kuraray co.,ltd., and a UV transparent bandpass optical filter, BG3 by Schott, is attached on every
photo cathode window. The night sky background photons are reduced by the filter,
and the typical background level is ∼ 30 photoelectrons in 100 ns. The transmittancespectrum of the filter has a necessary and appropriate width to transmit air fluorescencephotons. Figure 4.6 shows the typical transmittance of the UV00 and BG3. The PMTs
are manufactured by Hamamatsu Photonics K.K., and R9508 based on R6234-01 with aprinted circuit board include a bleeder and a preamplifier, which are used for the PMTs at
BR and LR FD stations. The PMT has a hexagonal bialkali photo-cathode and borosilicate
glass window, and the size of the window is a distance of 60 mm between the parallelsides. A printed circuit board which includes a bleeder and a preamplifier is installed on
the bottom of each PMT. The PMT has eight dynodes, and the gain is ∼ 8 × 104 with 800
V. The FOV of each PMT is approximately 1◦ . The PMT gains are relatively monitored
and adjusted to absolutely calibrated standard PMTs [57]. The gains of standard PMTs
are monitored by small light pulser of YAP (YAlO3 :Ce-Am light pulser) [58], which ismounted on the photo cathode surface of each standard PMT. The gains of other PMTs
are monitored relative to the standard PMTs with a Xe flasher installed at the center of
each mirror [59].

4.2. Fluorescence detectors

43

Figure 4.4: The possible hit patterns to generate a level-2 trigger.

Figure 4.5: (Left) an appearance photograph of the Black Rock Mesa FD stationbuilding. (Right) the photograph of one of the pair of upper and lower ring telescopes.

4.2. Fluorescence detectors

4.2.2

44

Trigger and data acquisition system

FDs have three types of VME electronics modules, signal digitizer/finder (SDF), trackfinder (TF), and central trigger distributor (CTD). Figure 4.7 shows a block diagram ofthese modules. Output signals from each PMT pass through a DC-coupled pre-amplifier,and it is digitized by a FADC on SDF with 12 bit 40 MHz sampling rate. Then, four
consecutive bins are combined to one bin in the SDF, i.e., it works as 14 bit and 10 MHz
sampling rate digitizer. One SDF module processes 16 input channels and 16 SDF modules
are installed for each camera. The output from 16 SDFs are processed on the TF module
to find a track in the signals recorded in the camera. The TF module is installed for eachcamera, and output from 12 TF modules are processed on CTD to trigger the all cameras.

The triggering logic of the FDs consists of three levels [60]. The level-1 trigger onindividual PMTs is generated in the SDFs. The SDF slices the data into 25.6 µs frames
containing 256 bins of FADC data samples, which has an overlap of 12.8 µs between
adjacent frames. The SDF calculates moving averages for four time windows (1.6, 3.2, 6.4
and 12.8 µs) and compares them with the background level. The background average andstandard deviation is calculated from past 1.6 ms, and they are updated every 0.4 ms.
When the SDF finds an excess greater than 6 sigma above the background level, it sends
an level-1 trigger to the TF module of the camera. The average trigger rate is about 3 Hz
with 6 sigma threshold level for the level-1 trigger. The level-2 triggers are generated by
the TF that performs pattern recognition for level-1 trigger maps of each camera.
There are two modes to search for air shower tracks on hit maps by TF. One is searching
for well-contained tracks in the camera’s FOV. TF scans over hit maps of cameras witha 5 × 5 pixel subarray window (Figure 4.8). When a hit pattern of the subarray matches
one of the patterns on the look-up table, TF generates a level-2 trigger and sends it withtagging complete track code to CTD.
The other mode is searching for shower tracks across the edges of cameras. In thismode, TF scans hit maps of cameras along the camera edges with 4 × 4 pixel subarray
window (Figure 4.9). When a hit pattern matches the look-up table for this mode, TF
send the level-2 trigger with the partial track code to CTD.
CTD judges a final trigger decision as a level-3 trigger. CTD generates a level-3 trigger,
when CTD receives at least one level-2 trigger with the complete track code, or when it
receives more than two level-2 triggers with the partial track code from adjacent cameras.
Then, CTD provides a level-3 trigger to all the cameras in a station. The level-3 trigger
makes all the PCs which control individual cameras to start storing waveforms from all
the PMTs in each camera. This trigger rate is about 1 ∼ 3 Hz in a stable observationnight.

4.2. Fluorescence detectors

45

Figure 4.6: (Left) the typical transmittance of the acrylic plate, UV00, the data
points are the median value for three samples, and errors indicate the diﬀerences
between the median and the otter two samples. (Right) the typical transmittance of
the UV filter, BG3 filter, and the data points are the medians for all sampled filterswith the bars corresponding to one standard deviation.

Figure 4.7: A block diagram of the triggering electronics and data acquisition system
of FD [60].

4.2. Fluorescence detectors

46

Figure 4.8: A schematic diagram of the complete track finding process [60].

Figure 4.9: A schematic diagram of the complete track finding process to find apartial track pattern across boundaries of adjoining two cameras [60].

4.2. Fluorescence detectors

4.2.3

47

Atmospheric monitoring

In order to reduce systematic errors and to keep the high event reconstruction accuracy
the TA experiment has various monitoring systems for the atmospheric condition.

Central laser facility
The central laser facility (CLF) is located in the center of the TA site, 20 km away fromall three FDs. CLF has a 355 nm Nd:YAG laser, Ultra CFR by Quantel. In usual air
shower observations with FDs, CLF shoots the laser into the FOV of FDs for 300 times
every half hour with 10 Hz repetition rate, which are observed by all the FDs, and FDs
observe side scattered photons of the laser. The total amount of side scattered photons
induced by a 5mJ laser shot is equivalent to the total fluorescence photons induced by the
air shower with primary energy of 1020 eV.
CLF is an important tool for understanding the atmosphere at the TA site. The
atmospheric transmittance and the vertical aerosol optical depth (VAOD) can be estimatedby comparing observed CLF events to that of simulations.

LIDAR
Light detection and ranging (LIDAR) is widely used technology for monitoring aerosols. In
the TA experiment, the LIDAR system was installed at 100 m apart from the BR station.
The LIDAR system consists of 355 nm Nd:YAG laser and telescope. The LIDAR systemobserve backscattered laser light by atmospheric molecules via Rayleigh scattering and byaerosols via Mie scattering. The time profile of backscattered light relates to the density
profiles of atmospheric molecules and aerosols along with the laser beam. Therefore, theatmospheric parameters, such as VAOD, can be estimated by LIDAR observation [61].
The typical attenuation length of Mie scattering on the ground and scale hight of aerosol
distribution are 29.4 km and 1.0 km, respectively, that corresponds the VAOD at the 3.5
km is 0.033. These values are determined by this LIDAR measurements, and they areused in the analysis of this thesis.

4.2.4

Cloud monitoring

The existence of clouds in the FOV causes unexpected photon scattering, and then it causes
the systematic errors on EAS reconstructions. Thus, we need to continue to monitor clouds
during every FD observation night. The TA experiment has two methods to monitor toevaluate the cloudiness.

IR camera
The TA experiment has an infrared (IR) camera to take IR pictures, which correspondto temperature measurement of the sky, of the FOV of the BR station. If it is cloudy,
the temperature in the cloud region is relatively higher than the clear sky region. In

4.2. Fluorescence detectors

48

the usual FD operations, the IR camera take twelve pictures which correspond to thetelescopes’ FOV every 30 minutes. Figure 4.10 shows the appearance of the IR camera andphotographs taken by the camera. In the analysis of the IR camera pictures, each picture
is divided into four separated regions in the horizontal direction. Then, the existence ofcloud is determined by the temperature in each region [62]. However, the IR camera is
not operated all the observation night, there is another cloud monitoring method in our
experiment.

WEAT code
In the MD station, an operator watches the existence of clouds every hour by own eyes.
The result for each time is recorded as a numerical code called WEAT code. The code
indicates the existence of clouds in the north, east, south, west, and zenith directions.
Additionally, the code includes the cloud coverage parameter in the whole sky and thetransparency of the horizontal direction. Since WEAT code is recorded in all of observationnights when MD FD is operated, the code can be used all of the observation terms. The
consistency between the WEAT code and the IR score is confirmed [62].

4.2. Fluorescence detectors

49

Figure 4.10: The appearance of the IR camera (left) and the twelve pictures takenby the IR camera at two diﬀerent times (right). The numbers shown in each picture
indicate the existence of clouds or not by 1/0 for in four separated regions in each
picture.

Chapter 5Monte Carlo simulation andshower reconstruction proceduresof the hybrid analysis
In order to evaluate reconstruction accuracies and systematic errors for observed resultswith a complex detector we use Monte Carlo (MC) simulations. In this chapter, the showerreconstruction and the simulation methods of hybrid technique are presented.

5.1

Monte Carlo simulations

For the hybrid event analysis in the TA experiment, there are three steps for EAS event
simulation. First, EASs are simulated with de facto standard MC simulation package called
CORSIKA, which is described in Section 5.1.1. Second, detector responses of SDs arecalculated using GEANT-4, which is a standard simulation package developed by CERN,
and it is discussed in Section 5.1.2. Third, detector responses of FDs are calculated using
the software which is developed by the TA collaboration, and the detailed descriptions arepresented in Section 5.1.3.

5.1.1

Extensive air shower simulations

In the study described in this thesis, EASs are simulated with the Cosmic Ray Simulations
for KASCADE (CORSIKA) package [63]. This package was firstly developed for the
KASCADE experiment, and now, is generally used as a de facto standard in the field ofcosmic ray physics. CORSIKA allows us to simulate EASs with various primary particle
types and primary energies. Moreover, we can set various parameters to control EAS
simulations. Table 5.1 shows a list of the most important parameters, options and set
values which are used in this work. Ecut in the table gives the lower energy limit ofparticles. When the particle’s energy reaches this limit, CORSIKA stops tracking the50

5.1. Monte Carlo simulations

51

Table 5.1: A list of most important parameters used for the CORSIKA simulationsin this work.
Ecut (hadrons)
Ecut (µ± )Ecut (e± , γ)Thinning factor
High energy interaction models
Primary particles
Energy range

50 MeV
50 MeV250 keV
10−6
QGSJET-II-03 & EPOS-LHC
Proton & Photon
from 1018.0 to 1020.5

particle and deposit its all remaining energy to the atmosphere.
For primary cosmic rays above 1018 eV, the number of secondary particles becomes toolarge to trace all of the particles. Therefore, for high energy simulations with CORSIKA,
we use the thinning option to reduce the number of secondary particles to trace, and it
can be set a weight to a particle. The weighted particle represents many particles of thesame species with the weight to ensure that the energy and the momentum are conserved.
When the thinning option is applied to CORSIKA, secondary particles are combined atany stage in the EAS development when their energies reach lower than the energy of,
Eparticle < ϵth E0 ,(5.1.1)

where E0 is the energy of the primary cosmic ray, and ϵth is the thinning factor. The
weighted particle is chosen randomly from all of the low energy particles. The thinning
factor of 10−6 is used in this work.
The thinning causes large fluctuation in the lateral distribution of the particle density.
The fluctuation in the lateral particle distribution does not aﬀect the production of fluorescence photons, and therefore, thinned showers work well for FD simulations, because
production of fluorescence photons depends on the overall deposited energy in each stepof the shower development. On the other hand, the fluctuations in the lateral density distribution at the ground level cause systematic eﬀects on SD simulations. Then, thinned
air showers can not be used directly for practical SD simulations. In order to respond thissituation, a method has been developed by the TA collaboration to recover the information
lost by thinning and to successfully use thinned showers for practical simulations. This
process is called dethinning, and the overview of this process is described in the followings
in this subsection.

Dethinning CORSIKA air showers
The overview of the dethinning process is described in this subsection, and the detailsof the dethinning process is described in the reference [64]. In the dethinning process,5.1. Monte Carlo simulations

52

weighted particles are spread out before the ground to restore the realistic information,
which is lost by the particle thinning in CORSIKA simulations [64]. In the dethinning
process, first, the vertex point of the weighted particle is chosen in the trajectory of theweighted particle with requirement that no particle can have an arrival times that precedes
the arrival timing of the shower front. An early arrival timing occurs when the total timeof-flight from the first interaction point, x0 , to the vertex point and then to the position on
the ground of the dethinned particle is less than the time-of-flight directory from x0 to thearrival position of the dethinned particle. In fact, there is a maximum distance between
the vertex point and the arrival position of the weighted particle. When the time andposition of first interaction are t0 , and x0 , and the arrival time and position of weighted
particle on the ground are ti and xi , respectively. The maximum distance, Dmax , along
the unit vector of the weighted particle trajectory, p̂i , is calculated as
Dmax =c2 (ti − t0 )2 − |xi − x0 |2,2 {c(ti − t0 ) − (xi − x0 ) · p̂i }(5.1.2)where c is a speed of light. The Dmax is the maximum distance, and any shorter distancecan be chosen for the vertex point.
Second, trajectories of the dethinned particles are chosen by a two-dimensional Gaussian with the σ of a few degrees. The σ is set to βd, where d is the lateral distance fromthe shower core to the weighted particle and β = 3◦ /km for electromagnetic particles
and β = 1◦ /km for muons and hadrons. These values for β is empirically determined bycomparison between unthinned showers and dethinned showers.
Third, crossing points between the chosen trajectories and the ground level are calculated as arrival points of the dethinned particles, and arrival timings and energies are
assigned to the dethinned particles. In fact, the timing is calculated by the trajectory ofeach particle, and the energy is determined by the Gaussian distribution centered on theenergy of the weighted particle with ±10% energy of sigma. Figure 5.1 shows an example
of a geometry determination for a dethinned particle.
To validate the dethinning process, approximately 100 unthinned air showers withprimary energies greater than 1018 eV were generated and compared with the air shower
data produced by the dethinning process. As a result, the lateral distributions producedby the dethinning procedure agree with the unthinned lateral distributions [65]. Figure 5.2
shows the comparison of secondary electron spectra with and without dethinning [64]. For
each histogram, good agreement is observed between thinned and dethinned simulations,and dethinned spectra are smoothed by the dethinning procedure.

5.1.2

Detector simulations for surface detectors

The dethinned CORSIKA simulations are used for the SD simulations. Each CORSIKA
shower is used repeatedly with random diﬀerent core positions to reduce the calculation
time. For SD simulations with proton primaries, each proton shower is also used repeatedly
with random diﬀerent azimuthal angle. In contrast, for photon primaries, the shower

5.1. Monte Carlo simulations

53

Figure 5.1: Schematic view for geometry of a dethinned particle [64]. Vertex position
is chosen on the trajectory of the weighted particle as described in text.

Figure 5.2: Comparison of secondary electron spectra with and without dethinning
for a thinned simulation of a proton-induced EAS with primary energy E = 1019.0
eV and primary zenith angle θ = 45◦ [64]. Horizontal axis is the kinetic energy, andvertical axis is the number of particles. Spectra for thinned simulation are showed
by gray, and dethinned spectra are showed by black. The secondary electrons withina region from −30◦ to 30◦ in azimuthal angles with respect to the azimuthal angleof the primary particle direction, and within a region from lateral distances from500 m to 1000 m are tabulated in these figures. The angle showed in each figure isincident angle with respect to the ground.

5.1. Monte Carlo simulations

54

Figure 5.3: Overview of the simulated SD configuration in the GEANT4 simulation.
It includes many of the components for realistic simulation, such as scintillators,
mounting parts, a stainless-steel box, a battery, a solar panel, wireless antenna, andso on.

development depends on their azimuthal angles because of preshower eﬀect as discussed
in Section 3.4.1.
The energy deposit processes on each SD are simulated by GEANT4 simulation package. The SD construction is simulated accurately in the GEANT4 simulation. The TA SD
is accurately constructed as a digital model using GEANT4 library as shown in Figure 5.3.
Deposited energies for passing particles with various diﬀerent incident angles, momenta,
and species are calculated.
To save conputing time of the SD simulations, energy depositions in the scintillators
have been simulated thousands of times and the results are combined to a look-up table.
The table is read out with the index made from the parameters of injected particles, which
are outputs of dethinned CORSIKA simulations.
Then, energies deposited in the scintillators and time dependent SD calibration information are combined for simulation of digital output waveforms by the SD electronics. The
triggering process described in Section 4.1.2 is also included in the simulation calculations.
Background signals induced by secondary cosmic rays are also simulated to makerealistic simulation data. The simulated background signals are generated from COSMOS
simulated secondary cosmic rays based on the energy spectra of primary cosmic rays for
all the primary species measured by AMS [66, 67]. Figure 5.4 shows the energy spectra of5.1. Monte Carlo simulations

55

Figure 5.4: The cosmic ray energy spectra for some primary species measured by
AMS [66, 67].

the considered primary cosmic rays. The simulated background signals agree well with themeasured background signals. Details of this analysis were described in the references [68]and [69]
Finally, the waveforms are packaged in the same format as the real data.

5.1.3

Detector simulations for fluorescence detectors

Simulation tools for FDs has been developed in Java with the TA collaboration. The
longitudinal development of energy depositions simulated with CORSIKA are used in
the FD simulations, because productions of fluorescence photons are induced by energies
deposited in the atmosphere. To simulate EASs as hybrid events, it is necessary to use
the same EASs in both the SD and the FD simulations. Thus, the shower developmentof simulated showers are packed in a look-up shower library, and it also enable to reduce
computing time.
Energy depositions are converted to fluorescence photons taking into account the wavelength spectrum described in Section 2.2.4, and parameterized atmospheric conditions.

5.1. Monte Carlo simulations

56

Figure 5.5: A digital model of a FD station building in the FD detector simulation.

The atmospheric condition parameters are obtained from results of the LIDAR operations
and radiosonde data. The radiosonde data are allowed universal access to all in at the website [70]. We take into account not only fluorescence photons, but also Cherenkov photons. After photon production processes, the number of photons is reduced with suﬀering
atmospheric attenuations and the solid angle of FDs.
Remaining photons trace to PMTs with considering the obscurations by the structuresof the FD station. Figure 5.5 shows a digital model of a FD station building in the FD detector simulation. The photons which hit a mirror are reflected with additionally including
a random fluctuation following the Rayleigh distribution with approximately σ = 0.06 degrees (diﬀerent measured values are used for each mirror). The parameter, σ, corresponds
to the averaged adjustment accuracy of segment mirrors. The optical parameters and thearrangement of FDs are checked with starlight calibrations with following two methods.
One is a way to use the time variation of pedestals of PMT outputs caused by momentsof stars in the FOV. Starlight can be observed as high pedestal value with DC-coupling
electronics of our FDs. Since positions of stars are well known, the pointing direction of
the optical axis of the PMTs and the curvatures of composite mirrors can be calibrated
by comparing measured and expected pedestal variations.
The other is a way to use photographs of star images on the focus plane, i.e., thecamera surface. The method is that taking pictures of star images on a fluorescent screen
on a camera surface from the center of a mirror. Then, the pointing directions of thecameras and curvature of the mirrors are calibrated by comparing the photographs withsimulated images on the camera surface. The results from these two methods are ingood agreement. Therefore, we use the averaged geometry of these two methods, and thesystematic uncertainty is less than 0.2 degrees.
The number of injected photons are reduced taking into account the transmissivities of5.2. Event reconstructions

57

the camera window and the UV filter in front of PMTs and converted to photo-electrons
taking into account the quantum eﬃciency, the correction eﬃciency, and the uniformityof PMTs. Then, signals measured at FDs are finally reproduced.
In addition, background light is important for FD simulations, especially for evaluating
the detection eﬃciency. The background light is caused by night sky background, stars
and artificial light. The mean and fluctuation of the background are recorded in every 10
minutes, which corresponds to the typical passing time of a star through the field of viewof a PMT.

5.2

Event reconstructions

This section explains an overview of the hybrid event reconstruction procedures [71, 72]
used in this study. The hybrid reconstruction consists of four steps: PMT selection, showergeometry reconstruction, longitudinal shower profile reconstruction, and quality cuts.
In the hybrid analysis, both FD and SD data are used. The timing information in
SD data is used for geometrical reconstructions. It allows significant improvement in theaccuracy of the determination of shower geometries comparing with the FD monocularanalysis. After the geometrical reconstruction, the longitudinal profile of the EAS isreconstructed using the FD data with precise geometry.

5.2.1

PMT selections

Before the reconstruction process, PMTs which truly detect shower signals must be identified and noise only pixels must be discarded. There are four steps of this PMT selection.
The first step of the PMT selection is done based on the strength of PMT signals.
When the signal of PMTs in a triggered camera is greater than 3σ above the backgroundlevel, the PMTs are selected.
Then, the PMTs selected with the first level PMT selection are used to draw a shower
track on the cameras. The PMTs’ positions and signal timing are spatially and temporallycompared with the track, and PMTs isolated from the track are discarded with the second
and third level PMT selection. At this step, the shower detector plane (SDP) is defined,
which is the plane include the shower axis and the FD.
Only high significant signal PMTs are selected before this selection, however in the
fourth level selection, PMTs which have lower significance of signal are also considered.
Signal timings of PMTs are fitted and distances between each PMT and the SDP arecalculated. Then, isolated PMTs are discarded, and PMTs which are well fitted and close
to the SDP are added to the selected PMT list.

5.2.2

Geometrical reconstruction

The energy threshold of the analysis in this thesis is lower than the analysis which uses
only SD data. For such low energy EASs, the number of triggered SDs is too small to5.2. Event reconstructions

58

Figure 5.6: Schematic view of the parameters related to the geometrical reconstruction [71].

estimate the core position of the EAS by using only SD information. Therefore, in thishybrid analysis, the timing information of one triggered SD and FD information are used
to reconstruct the shower geometry.
In the geometrical reconstruction of the monocular FD analysis, the geometry is determined taking into account both the pointing direction and the timing of hit PMTs fit
to the following equation,
Texp,i = Tcore +sin ψ − sin αi
Rcore ,
c sin(ψ + αi )(5.2.3)where Texp,i and αi are the expected timing and the elevation angle on the SDP for thei-th PMT, ψ is the angle between shower axis and the ground, and Tcore is the timing
when the shower hit the ground, and Rcore is the distance between the core position of
the EAS and the FD (see also Figure 5.6).
Tcore can be expressed with the information of the SD which is the nearest to the EAS
core as follows,1′
Tcore = TSD+ (Rcore − RSD ) cos ψ,c(5.2.4)1 ⃗′′TSD = TSD − {(P SD − P⃗SD ) · P⃗ },cwhere P⃗SD is the position of the SD, P⃗ ′ SD is the projected SD position on the SDP, P⃗ is the′ is the corrected timing
direction of the shower axis, TSD is the timing of the SD signal, TSD
at the projected SD position, and RSD is the distance between the FD and P⃗ ′ SD . Then,5.2. Event reconstructions

59

to determine the shower geometry, signal timings are fitted to minimize the following χ2 ,
χ2 =∑ (Texp,i − Ti )2
i

2
σT,i

,

(5.2.5)where σT,i is the fluctuation of the signal timing.
SDs which are used in the geometrical reconstruction are selected from the SDs withincloser than 1.5 km from the core position, which is determined by weighted average pointof SD signal distribution in the first calculation, and closer than 1.2 km from the lineof intersection of the SDP and the ground to reduce the fake triggered SDs, which are
triggered with chance incident particles. When the core position is determined by theprocedures described above, then these procedures are iterated by using the calculated
core position as initial value of it to get the minimum χ2 .

5.2.3

Longitudinal profile reconstruction

For the longitudinal profile reconstruction we use only FD data. Since the shower axis
for each event is already known as the result of the geometrical reconstruction, the development of the EAS along the the shower axis is evaluated with the profile reconstruction
process. For the profile reconstruction, the amount of fluorescence photons is important, but there are contamination from Cherenkov photons in detected signals. There are
four types of photon components in detected signals, that is fluorescence photons, direct
Cherenkov photons, scattered photons by atmospheric molecules and aerosols. On orderto reconstruct the longitudinal profile analyzing those various photon components we use
a method called inverse Monte Carlo (IMC) technique.
The idea of the IMC technique is that searching the best fit EAS profile for observed
air showers through comparisons between the observed data and MC simulations. By
using the IMC technique, processes which can not be reproduced from observed signals,
such as non-uniformity of the the PMT sensitivity, can be realistically considered in thereconstruction process.
For the shower longitudinal profile calculations in the IMC we use numerical calculations based on Gaisser-Hillas function [27] instead of using full MC simulations, such as
CORSIKA. In the numerical calculations the first interaction point, X0 , and the interaction length, λ, is fixed at 0 g/cm2 and 70 g/cm2 , respectively, because these parameters
are less sensitive for the energy determination. Then, the adjustable parameters in the
IMC are Xmax and Nmax .
At the first step of the IMC process, the Nmax is fixed to 1.0, and the Xmax is theadjustable parameter to be optimized. With altering values substituted into Xmax wesimulate the expected number of photo-electrons for all the PMTs taking into account thedetector response as the same routine described in Section 5.1.3. The optimized value for5.2. Event reconstructions

60

Figure 5.7: An example of the result of the longitudinal profile reconstruction with
the IMC technique [71]. The horizontal axis is the slant depth along the shower
axis, and the vertical axis is the number of photo-electrons observed by the FD. The
black points show the observed data and the red and the blue histograms indicate
the fluorescence and the scattered photon components optimally estimated by the
IMC. Shower parameters of this event are estimated as follows: zenith angle is
4.0 degrees, azimuthal angle is 313.1 degrees, impact parameter (distance from thestation to the shower axis) is 17.7 km, primary energy is 8.1 × 1019 eV, and the Xmax
is 756.4 g/cm2 .

Xmax is determined by maximizing the following likelihood, L, which is calculated byL=nexp,sum =

PMTs
∑
i
PMTs
∑

nobs
i log

nexp
i
nexp,sum,(5.2.6)

nexp
i ,

i

where nobs
and nexp
is the observed number of photo-electrons with i-th PMT.
ii
After the optimized Xmax is determined, the optimized value for Nmax is calculated asfollows,∑PMTs obs
nii.
(5.2.7)
Nmax = ∑PMTs
nexpii
Figure 5.7 shows an example of the result of the longitudinal profile reconstruction with
the IMC technique.
The primary energy of the EAS is calculated by integrating the Gaisser-Hillas function

5.3. Monte Carlo study

61

with optimized parameters as a result of the profile reconstruction.

5.2.4

Quality cuts

There are contaminations of inaccurately reconstructed events even if reconstruction process are successfully finished. For example, for the event of which true Xmax is outsideof the FOV of FDs the accurate estimation for Xmax is diﬃcult. Therefore, we need todiscard these inaccurately reconstructed events from analyzed event data set in order toavoid using these inaccurate events for further analyses. Requirements for selected events
to pass this quality cut are as follows:• The number of selected PMTs > 20
• The estimated Xmax is inside of the field of view of the FD• The estimated incident zenith angle < 55 degrees• The estimated core position is inside of the SD array• The minimum viewing angle (see text) > 20 degrees
The minimum viewing angle is the minimum angle between the shower axis and the viewing
angle from the FD to the each point of shower axis within the inside of the FOV. Thiscut reduces events which have photon component dominated by the Cherenkov photons,
because for such events, it is diﬃcult to reconstruct the longitudinal profile of fluorescencephotons with accurately.

5.3

Monte Carlo study

5.3.1

Monte Carlo data set

I prepared hybrid MC data sets for primary protons and primary photons with the procedures as described in Section 5.1. The parameters for making the MC data set are listed
in Table 5.2.

Table 5.2: Parameters for making the MC data set.
Zenith angle0◦ < θ < 60◦ , uniform
Azimuthal angle 0◦ < ϕ < 360◦ , randomCore position
within 30 km from CLF, random

It is diﬃcult to obtain enough MC event statistics in the whole energy range, 1018,0 ∼1020.5 eV, by throwing simulations with energy spectrum of E −3 . To overcome this situation, I made MC data with spectral index of −3.34 in the energy range from 1018.0 eV
to 1018.5 eV, and with induex of −1.0 in the range above 1018.5 eV. Then, each simulated

5.4. Data and Monte Carlo comparisons

62

event are weighted to fit the energy spectrum to the observed energy spectrum by the SD
measurement [6].
Statistics of the MC data set with each step of the quality cuts are listed in the Table
5.3.

Table 5.3: The statistics of the MC data for protons and photons. The parent
CORSIKA showers are repeatedly used for the simulations (see Section 5.1.2).
The number of survived events at each step
Protons
The number of parent CORSIKA showers
6500
The number of thrown MC events ∼ 3.5 × 106
The number of reconstructed events
119351
The number of selected PMTs > 20113907The estimated Xmax is inside of the FOV
103360The estimated incident zenith angle < 55 degrees
93873
The estimated core position is inside of the SD array
81752
The minimum viewing angle (see text) > 20 degrees
73765

5.3.2

Photons7280∼ 1.4 × 106
35634
31992
20770
17801
15999
14051

Resolutions

The resolutions of reconstructed geometry and shower profile for the hybrid reconstruction
technique are obtained from comparisons between simulated and reconstructed parameters. In order to evaluate the resolutions, the MC data set which passed the quality cuts
written in Section 5.2.4 is used.
The diﬀerence of core position between simulated and reconstructed values is shownin Figure 5.8, and the resolution of the core position is 111 m. That of the open angle isshown in Figure 5.9, the resolution is 0.7◦ .
The resolution and energy dependence of energy estimation is shown in Figure 5.10.
The energy resolution is about 9.5% with +2% systematic shift, and it has small energy
dependency. The Xmax resolution and energy dependency are shown in Figure 5.11. The
Xmax resolution is about 29 g/cm2 , and systematic shift is −14 g/cm2 . The energy dependence of the Xmax reconstruction is small, and that is negligible with in the resolution.
The systematic bias of Xmax may be come from reconstruction procedure, thus it also
aﬀect reconstruction results of real data.

5.4

Data and Monte Carlo comparisons

We use MC simulations not only for the event reconstruction procedures, but also getting
the physics results. Therefore, it is important to ensure that our MC simulations accu-

5.4. Data and Monte Carlo comparisons63

Figure 5.8: The histogram of diﬀerence between simulated and reconstructed core
position. The horizontal axis is the diﬀerence, and the vertical axis is normalized
entry. The resolution is calculated with region of 68%, and it is 111 m.

Figure 5.9: The histogram of diﬀerence between simulated and reconstructed open
angle of arrival direction. The horizontal axis is the diﬀerence, and the vertical axis
is normalized entry. The resolution is calculated with region of 68%, and it is 0.7◦ .

5.4. Data and Monte Carlo comparisons64

Figure 5.10: (Left) the histogram of energy ratio of simulated and reconstructedenergy. The systematic shift of reconstructed energy is about +2% and the resolution
is about 9.5%. (Right) the energy dependence of the energy ratio. Reconstructedenergies have small energy dependence.

Figure 5.11: (Left) the histogram of diﬀerence between simulated and reconstructed
Xmax . The systematic shift of reconstructed Xmax is about −14 g/cm2 and theresolution is about 29 g/cm2 . (Right) the energy dependence of the Xmax . Almost
no energy dependence can be found in Xmax reconstruction.

5.4. Data and Monte Carlo comparisons65

rately reproduce real observed data. For this purpose, in this section we compare various
observable values obtained from the analysis for MC simulated and real observed data.

5.4.1

Observed data set

The TA hybrid measurement has been continued since March 2008. The period of theobserved data set is from May 2008 to July 2013 to avoid unstable operation term at theearly term of the hybrid measurement. In this thesis, I used the data obtained by two
FD stations, BR and LR stations, and the SD array. The number of events which haveenergies grater than 1 EeV is 2970 after the quality cut.

5.4.2

Comparisons of observable parameters

In order to check the consistency between the observed data and MC data, comparisonplots for many observable parameters are made as shown in Figure 5.12-5.19. In thesefigures, the black crosses show measured values of the data, and the red and the blue
histograms indicate MC predicted values calculated with primary protons and primary
photons, respectively. The upper figures are histograms of each parameter, and lower
figures are the ratios of real data to MC data in each bin. For all the data in these figures,the quality cuts described above have been applied. The distributions of observables for
proton MC are in almost good agreement with that of the observed real data.

5.4. Data and Monte Carlo comparisons66

Figure 5.12: The Distributions of the number of event. Each MC events are weighted
to fit the energy spectrum to the obtained energy spectrum from the SD analysis [6].
The proton spectrum is in good agreement with that of real data without energy of1020 eV bin. The discrepancy is come from the choice of the bin size.

Figure 5.13: The Distributions of reconstructed Xmax . There is systematic shift
between proton MC and real data distribution. The amount of shift is about 20
g/cm2 .

5.4. Data and Monte Carlo comparisons67

Figure 5.14: The distributions of core positions along with the x-axis from west toeast direction centered at the CLF position. The MC data set distribution is ingood agreement with that of the real data set.

Figure 5.15: The distributions of core positions along with the y-axis from south tonorth direction centered at the CLF position. The MC data set distribution is ingood agreement with that of the real data set.

5.4. Data and Monte Carlo comparisons

68

Figure 5.16: The distributions of the zenith angle of EASs. The MC data set
distribution is in good agreement with that of the real data set with in statisticalerrors.

Figure 5.17: The distributions of the azimuthal angle of EASs. The MC data set
distribution is in good agreement with that of the real data set.

5.4. Data and Monte Carlo comparisons69

Figure 5.18: The distributions of the impact parameter of EASs. The MC data set
distribution is in good agreement with that of the real data set.

Figure 5.19: The distributions of the ψ angle, which is the angle between the shower
axis and the direction from the FD station to the shower core position on the SDP(see Figure 5.6). The MC data set distribution is in good agreement with that of
the real data set below 120◦ . Above 120◦ , there is systematic diﬀerence between thereal data set and MC data set.

Chapter 6
Search for ultra high energy
photons
In this chapter the results of the UHE photon search based on observed Xmax are presented.
First, I will describe the acceptance bias for detecting EASs induced by primary photons
and a data reduction method to reduce the acceptance bias. Second, a criterion for
selecting UHE photon-like events is discussed. Then, as a result I obtained photon fraction
upper limits, and I discuss the results based on comparison with theoretical predictions
and other experiments results.

6.1

Acceptance bias and data reduction

The averaged Xmax of photon induced showers are larger than nucleus induced showers
as discussed in Chapter 3. As discussed in Section 5.2.4 the quality cut based on thereconstructed parameters are applied in order to keep high reconstruction accuracies. One
of selection rules in the quality cut requires the reconstructed Xmax to be inside the FOVof FDs, thus deeply penetrated EASs are discarded by this criterion because Xmax of such
EASs are below the FOV of FDs. Since most of EASs of primary photons have significantly
deeper Xmax than that of primary protons at energies greater than 1 EeV, photon induced
EASs are more likely to be discarded with this criterion. This eﬀect is relatively low for
inclined showers, because inclined showers pass through more atmosphere before entering
the field of view of the FD. Thus, events with large Xmax are more likely to be remained
after the quality cut in inclined showers comparing with that of vertical showers.
Therefore, I applied an additional cut criterion discarding the events which have smaller
reconstructed zenith angle than 20◦ (photon enhance cut), to suppress the detection eﬃciency for protons, and comparatively to enhance the detection eﬃciency for photons.
The acceptance bias is calculated comparing the survival fraction, R, for photon primaries with that for proton primaries. The fraction R is defined as the ratio of thenumber of remained events after reconstructions and quality cut to the number of thrown70

6.2. Selection criterion for photon primary events

71

Figure 6.1: The survival fraction R. Red and blue points indicate R for protons
and that for photons, respectively. The lines show the fit results with an empiricalequation.

MC events, i.e.,
R=

Nreconstructed.
Nthrown(6.1.1)Figure 6.1 shows R for proton and photon primaries.
As a result, we get the acceptance bias as the ratio of R for photons to that for protons,
at several energy ranges as 0.512, 0.502, 0.482, and 0.445 at energy grater than 2, 3, 5,
and 10 EeV, respectively.

6.2

Selection criterion for photon primary events

As discussed above, EASs induced by primary photons have larger Xmax than that ofprotons on average. Thus, I use observed Xmax as a key to select photon primary events.
The selection criterion is determined from the MC simulations for photon primaries. I
calculated the average of reconstructed Xmax for photon MC, and adopted it for thephoton selection criterion. Figure 6.2 shows the averaged Xmax of the photon primary
EASs. The fitting result of the averaged Xmax for the photon primaries are the photon
selection criterion in this work. When the reconstructed Xmax of an event is larger than
the selection criterion at the energy of the event, the event is tagged as photon-like event.
The survival eﬃciencies for primary photons of this selection criterion are 0.498, 0.497,
0.532, and 0.515 at energies greater than 2, 3, 5, and 10 EeV.
The longitudinal development of simulated EASs of primary nuclei depend on highenergy hadronic interaction models. On the other hand, the interaction model are less

6.3. Data analysis

72

Figure 6.2: The blue dots are Xmax of each photon-induced EAS. The blue open
circles show the averaged Xmax at each energy bin, the size of each bin is 0.1 inlogarithm of energy. The lines show the fitting results of photon averaged Xmax .
There is a bend at 1019.6 eV for averaged Xmax of photons due to the preshower
eﬀect.

eﬀective on EASs , of primary photons, because EASs induced by photons are almost
purely electromagnetic cascades. In order to confirm this model dependence, averaged
Xmax for photons with QGSJETII–03, EPOS–LHC are simulated and compared eachother. Figure 6.3 shows the averaged Xmax with QGSJETII–03 and with EPOS–LHC.
The maximum discrepancy of the averaged Xmax between these two models are less than
10 g/cm2 . The maximum discrepancy is smaller than the Xmax resolution of the hybrid
analysis, thus, the model dependency is negligible for setting the photon selection criterion.
In Figure 6.4, the selection criterion and reconstructed Xmax of proton primary EASs
are plotted. In order to estimate the proton contamination with this photon selection
criterion, the photon selection is applied to the proton MC data set. Then, 3.3%, 2.7%,
1.8%, and 0.9% protons are tagged as photon-like events. Figure 6.5 shows distributionsof reconstructed Xmax of protons and photons at several energy ranges.

6.3

Data analysis

The observation period for data used in this analysis is the same as shown in Section 5.4.
The observed and the MC data are reconstructed and applied the quality cut describedin Section 5.2.4, and the photon enhance cut described in Section 6.1.
Figure 6.6 shows Xmax of observed events and the photon-like event selection criterion.
The events which have larger Xmax are selected as photon-like events. The numbers of the6.3. Data analysis

73

Figure 6.3: The interaction model dependence for Xmax of primary photon showers.
The filled circles indicates the averaged Xmax for each energy range with QGSJETII–
03, and the solid line is fitted to the filled circles. The open circles and dashed line
are calculated with EPOS–LHC. These values are calculated with CORSIKA.

Figure 6.4: The red dots are reconstructed Xmax for MC simulated events of primary
protons, and blue line is the selection criterion determined by Figure 6.2. Opencircles show the averaged Xmax of protons at each energy bin, the size of each bin
is 0.1 in logarithm of energy. We can find the contamination from protons to thephoton selection criterion.

6.3. Data analysis

74

Figure 6.5: Xmax histograms for several energy ranges in 1018.2 eV, 1018.5 eV, 1018.8
eV, 1019.1 eV, 1019.4 eV, and 1019.7 eV from upper to lower figures, respectively. The
red histograms and the blue histograms indicate protons and photons, respectively.
The blue arrow in each figure shows the photon selection criterion at each energy
range.

6.4. Upper limits of the photon fraction

75

Figure 6.6: Reconstructed Xmax of observed events and the selection criterion. The
circles show data points and the line is the selection criterion described in Section
6.2. The observed data are passed through the quality cut and the photon enhancecut.

photon-like events are 24, 11, 3, and 0 in the energy ranges with energies greater than 2,
3, 5, and 10 EeV, respectively. The numbers of all events are 1396, 773, 229, and 121 in
the same energy ranges.
Each of the photon-like event selected as above can not be distinguished whether real
photon or not, because of the proton contamination. It is expected that some protons
are able to have deeper Xmax than the selection criterion due to large fluctuation of alongitudinal development. Thus, we compared the number of photon-like events withexpected number of photon-like event from deeply penetrated proton EASs. In order toestimate the proton contamination, the same selection is applied to proton MC. As aresult, the number of photon-like events with proton MC are 46, 21, 6, and 1 with energy
ranges grater than 2, 3, 5, and 10 EeV, respectively. Therefore, the result from real data
analysis is not exceed the expectation from proton contamination, even though the shiftof Xmax between real data and proton MC (see Figure 5.13) is considered.

6.4

Upper limits of the photon fraction

We can not exclude all of the protons from the photon-like events, the upper limit of thephoton fraction which is the fraction of photon flux to all particle flux is calculated fromobserved data. In order to obtain a conservative upper limit, the photon-like events are
assumed to be truly induced by primary photons. This is because the Xmax distribution of
protons highly depend on which the high energy hadronic interaction model is used, thus
we can not estimate realistic contamination from protons without systematic uncertainty.

6.5. Systematic error from energy uncertainty

76

The upper limit of the number of photon-like events using the Poisson distribution is
calculated with assuming no background, i.e., all the photon-like events are true photons.
Poisson upper limit of the number of photon-like events with 95% confidence level isobtained as 33.8, 18.2, 7.8, and 3.0 in the range with energies above 2, 3, 5, and 10 EeV.
As discussed above the photon selection of this analysis discard not only nuclei but also
photons, thus the selection eﬃciency should be taken into account for calculating theupper limit. As a results we obtained 95% confidence level upper limit of the number ofphoton-like events, nγ |95% , as 67.8, 36.6, 14.6, and 5.8 in the range with energies above 2,
3, 5, and 10 EeV.
Finally, we obtained the upper limit of photon fraction with 95% confidence level,Fγ |95% , from following relation,
Fγ |95% =Iγ |95%nγ |95%=,
Iγ |95% + Ip
nγ |95% + (nobs − nγ |95% )Aγ /Ap(6.4.2)where Iγ and Ip are the fluxes of photons and protons, nobs indicates the number ofobserved events, nγ is the number of the photon-like events, Aγ and Ap are the aperturesof the experiment for primary photons and protons, respectively, and the notation of“95%” means 95% upper limit. The ratio of the apertures, Aγ /Ap , is substituted fromthe acceptance ratio derived in Section 6.1. Then, the upper limit of the photon fraction
with 95% confidence level are obtained as 9.1%, 9.0%, 8.5%, and 10.2% in energies above
2, 3, 5, and 10 EeV, respectively. The results of the calculations are shown on Table 6.1.
Figure 6.7 is the photon fraction upper limit comparing with other experimental results
and predictions from several models. In this figure, the result from the same processes ofanalysis as this work with monocular data set, which are reconstructed with monocularmethod [73], is also showed.

6.5

Systematic error from energy uncertainty

There are systematic uncertainties on energy determination with the hybrid analysis.
Table 6.2 shows a list of sources of the systematic uncertainties [71]. It consists of uncertainties of the detector calibration, the atmospheric attenuation, the fluorescence yield,
and the reconstruction error. The total systematic uncertainty on energy determination
is calculated with the quadratic sum of these uncertainties, and that is 21%.
In order to estimate the systematic uncertainty on the photon fraction upper limits,
the upper limits are recalculated with shifting the energy scales with −42%, −21%, +21%,
and +42%. Table 6.3 shows the results of the photon fraction upper limits calculated withshifting energy scales. Thus, I adopted the worst upper limits in each energy range 9.4%,
9.0%, 9.8%, and 26.6% in energies above 2, 3, 5, and 10 EeV.
Figure 6.8 shows the results of the photon fraction after the evaluation of the systematic
uncertainty.

6.5. Systematic error from energy uncertainty

77

Table 6.1: The results of the photon fraction upper limits with 95% confidence level.
See text for details.
> 2 EeVphoton-like
24

The number of
events, nγ
Poisson fluctuation with 95% c.l.
Photon selection eﬃciency
Upper limit of the number ofphoton-like events with 95% c.l.
nγ |95%
Acceptance ratio, Aγ /Ap
The number of all data, nobs
The upper limit of the photon
fraction with 95% c.l., Fγ |95% (%)> 3 EeV11> 5 EeV3> 10 EeV
0

33.8
0.498
67.8

18.2
0.497
36.6

7.8
0.532
14.6

3.0
0.515
5.8

0.512
1396
9.1

0.502
773
9.0

0.482
339
8.5

0.445
121
10.2

Figure 6.7: The photon fraction upper limits. Black arrows show the results of upper
limits from this work. The other arrows show other experimental results: gray is
the TA mono analysis result, red is the TA SD analysis result [74], light green is
the Auger hybrid result presented in 2009 [75], green is the Auger hybrid result
presented in 2011 [76], light blue is the Auger SD result [77], blue is the Haverah
Park [78, 79], yellow is the Yakutsk [80, 81], and pink is the AGASA result [82, 83].
The lines show the predictions from top-down models (see [84–86])

6.5. Systematic error from energy uncertainty

78

Table 6.2: The systematic uncertainties on energy determination [71].
Item
Detector sensitivity

Uncertainty (%)
10

Atmospheric collection
Fluorescence yield

11
11

Reconstruction
Total

10
21

Contributions
PMT(8%), mirror (4%), aging (3%), filter (1%)
aerosol (10%), Rayleigh (5%)
model (10%), humidity (4%), atmosphere (3%)
model (9%), missing energy (5%)

Table 6.3: The photon fraction upper limits calculated with shifting energy scales.
−42%# of photon-like events# of data
95% u.l. (%)−21% # of photon-like events# of data
95% u.l. (%)
Original # of photon-like events# of data
95% u.l. (%)+21% # of photon-like events# of data
95% u.l. (%)+42% # of photon-like events# of data
95% u.l. (%)> 2 EeV
8
620
8.7
17
1015
9.4
24
1396
9.1
23
1758
7.0
26
2073
6.6

> 3 EeV
2
315
7.7
5
535
7.6
11
773
9.0
12
1041
7.2
15
1300
6.9

> 5 EeV
0
157
7.2
2
237
9.8
3
339
8.5
4
464
7.4
5
600
6.6

> 10 EeV
0
42
26.6
0
79
15.2
0
121
10.2
0
165
7.6
0
200
6.3

6.5. Systematic error from energy uncertainty

79

Figure 6.8: The photon fraction upper limits derived from the worst case of upper
limits with shifting energy scale calculations. Notation in this figure is the same as
Figure 6.7.

Chapter 7
Discussion
The studies for the nature of extensive air showers induced by primary photons and theobservational results of searching for UHE photons using Xmax with hybrid reconstruction
technique are presented in this thesis. This is the first result which measured by the hybrid
detection technique in the northern hemisphere.
This analysis is insensitive to the hadronic interaction models, because the photon
primary EASs consist of almost purely electromagnetic cascade, and the eﬀect on Xmax
by the hadronic interaction model used in the analysis is suppressed.
As a result of the photon fraction upper limits without systematic uncertainty ofenergy, one of the SHDM models is constrained at energy range above 10 EeV. When
including the energy systematic uncertainty, the predicted photon fractions with topdown scenarios are not constrained. However, in consideration of combination with theTA surface detector result, super heavy dark matter models and topological defect modelof the UHECR generation are constrained with 95% confidence level, and Z-burst model
are survived with these upper limits in the norther hemisphere. It is important thatindependent search of UHE photons in the northern and the southern hemisphere, because
UHE photons can be expected to be anisotropically observed since sources or generation
points of UHE photons are limited to near from the Earth due to their mean free path,and also they are not deflected by the magnetic fields in the universe. Thus, the result
complements other results in the southern hemisphere, such as PAO results. In addition,
the uncertainties on other analyses caused by contamination of photon primary EASs,
such as composition analysis using average of the Xmax , are limited by this result.
If we get enough statistics of UHECR events (10 times) and there is no photon in thedata set, we would be able to achieve a few percent upper limit at energy range above
10 EeV (Figure 7.1), and almost all top-down scenarios, except Z-burst model, would be
excluded. In order to constrain the predicted photon fraction with GZK photons, however,
it is necessary that improving the analysis method.
In future analysis, we will use information of secondary particles at the ground, which
are derived as SD observable, and then it is expected that the discrimination power ofphotons from nuclei is improved, and a photon likeness of each event can be estimated.

80

Chapter 7. Discussion

81

Figure 7.1: The photon fraction upper limits with enough (10 times) statistics.
Notation in this figure is the same as Figure 6.7. Fraction upper limits at lower
energies are higher than that of 10 EeV, because separation of average Xmax between
protons and photons is worse at lower energy ranges.

We are now constructing the muon sensitive and countable detectors in the TA experiment, using segmented plastic scintillators, and lead as an absorber for electromagnetic
component. The information of muonic component in EASs are important to measure thechemical composition of UHECRs, especially in the photon search because primary photon EASs are include less muonic component. Then, the correlation analysis using photon
candidates will be available in future analysis, such as correlation between arrival timingsof photon candidates and timings of high energy phenomena in the universe. If UHECRs
are generated by transient phenomena, such as GRB, UHE photons and neutrinos can be
simultaneously detected with other optical experiment. The simultaneous detection canbe a strong clue to discover sources and to reveal the nature of UHECRs.
In order to detect UHE photons, composition sensitive measurement is very important,such as using the fluorescence detection technique, and also the statistics are importantto achieve more significant results.

Chapter 8
Summary and conclusion
The work in this thesis is focused on studying ultra high energy photons based on theobservations with the hybrid air shower detector, to reveal the nature of ultra high energy
cosmic ray origins.
An overview of recent progress in cosmic ray studies, especially focused on UHECR
physics, is given in Chapter 2. There are a lot of progress in the studies of UHECR
physics, however, the origins of UHECRs are still unknown. Since, the flux of cosmic rays
has steep power law index, it is diﬃcult to observe suﬃcient number of events at UHECR
regions to reveal the origins of UHECRs. Furthermore, the magnetic deflection also makes
it diﬃcult because almost all the UHECRs are charged particles. Therefore we need toattempt to other approaches for discovering UHECR origins.
Many of models of UHECR sources and propagation mechanisms predict the existenceof UHE photons. I reviewed these models in Chapter 3. In the top-down scenarios of
UHECR origins, relatively higher fraction of UHE photons to charged hadrons is expected
than that of the bottom-up scenarios. Thus, we can test these scenarios with measuring
photon fraction in observed UHECRs.
In order to distinguish UHE photons from UHECRs, the maximum point of the longitudinal development of EASs, Xmax , is a powerful discriminator and one of the composition
sensitive parameter. The averaged Xmax for primary photons is significantly larger thanthat of nuclei as it is shown in Figure 3.3.
TA is the largest hybrid detector in the northern hemisphere, which consists of 507
surface detectors with approximately 680 km2 area, and the three fluorescence detectors,
which viewing over the SD array. The TA experiment continues to the full operation since
May 2008, The TA hybrid detector observes the longitudinal development of EASs withfluorescence detectors and detects EAS particles at the ground with the surface detectors.
By using both of the data obtained from the SDs and from the FDs the geometry andthe primary energy of EASs are precisely estimated with the hybrid analysis, and also
the Xmax is precisely determined from the longitudinal development. Since Xmax is thecomposition sensitive parameter, the TA hybrid detector is suitable for searching UHE
photons from observed UHECRs. The hybrid data set simultaneously observed with the82

Chapter 8. Summary and conclusion

83

700 km2 surface detector array and with BR and LR fluorescence detector stations areused in this study searching UHE primary photons.
The data period is include the period from March 2008 to July 2013, and the quality
cut and photon enhance cut requirements (see Section 5.2.4 and 6.1) are applied to thedata.
The selection criterion to distinguish primary photons from primary nuclei used in thisstudy is determined through MC simulation study of air shower development and detector
response. We use de facto standard simulation package called CORSIKA to simulate
the longitudinal development of EASs, and GEANT4 package is used to simulate thesurface detector response. The response of the fluorescence detectors are simulated withthe software package developed by the TA collaboration. Finally, the average values ofreconstructed Xmax for MC simulated primary photon showers are adopted as the photon
selection criteria in this study.
The result of the photon selection we obtain 24, 11, 3, and 0 events as the photon-likeevents in the energy ranges with energies greater than 2, 3, 5, and 10 EeV, respectively.
Each of the photon-like event can not be distinguished whether real photon or not, becauseof the proton contamination. It is expected that Xmax of proton primaries can be larger
than the selection criterion due to large fluctuation of a longitudinal development. We
can not exclude all of the protons by using only Xmax information, thus the upper limitsof the photon fractions are derived from the numbers of photon-like events.
The upper limits are obtained from the data with taking into account the detection
biases for EASs by protons and by photons, and the selection eﬃciencies for the photon
selection. The systematic uncertainty on the photon fraction upper limits arising fromthe uncertainty of energy determination are also considered. Finally, the upper limits ofphoton fractions with 95% confidence level are obtained as 9.4%, 9.0%, 9.8%, and 26.6%
for the range of energies greater than 2, 3, 5, and 10 EeV, respectively. These upper limits
are the first result, which are measured with hybrid detectors technique in the northernhemisphere, and these upper limits also ensure that the uncertainties in other analyses
due to photon contamination, such as primary composition analysis, are reasonably small.
In conclusion, as a result of the upper limits the predicted photon fraction with topdown models which considered in this thesis are not constrained, but in consideration of
combination with the TA surface detector result, super heavy dark matter models andtopological defect model of the UHECR generation are constrained with 95% confidence
level, and Z-Burst model are survived with these upper limits in the norther hemisphere.
In near future, we will update the analysis method by adding the information of
secondary particles at the ground, which are derived as SD observable, then a photon
likeness of each event will be available. In addition, the directions and arrival timings ofphoton candidates are used in future analysis.

References[1] S. P. Swordy, Space Science Reviews, 99, pp. 85–94, 2001.
[2] J. Blümer, R. Engel, and J. R. Hörandel, Progress in Particle and Nuclear Physics,
63, pp. 293–338, 2009.
[3] T. Antoni et al., Astroparticle Physics, 24, pp. 1–25, 2005.
[4] K. Greisen, Physical Review Letters, 16, pp. 748–750, 1966.
[5] G. T. Zatsepin and V. A. Kuz’min, JETP Letters, 4, pp. 70–80, 1966.
[6] T. Abu-Zayyad et al., The Astrophysical Journal Letters, 768, L1, 2013.
[7] J. Abraham et al., Physical Review Letters, 101, 061101, 2008.
[8] J. A. Simpson, Annual Review of Nuclear and Particle Science, 33, pp. 323–381,
1983.
[9] E. Fermi, Physical Review, 75, 8, pp. 1169–1174, 1949.
[10] E. Fermi, The Astrophysical Journal, 119, 1, 1954.
[11] T. K. Gaisser, Cosmic Rays and Particle Physics, Cambridge University Press, 1990.
[12] L. Drury, Reports on Progress in Physics, 46, pp. 973–1027, 1983.
[13] P. O. Lagage and C. J. Cesarsky, Astronomy and Astrophysics, 118, pp. 223–228,
1983.
[14] T. Stanev, High Energy Cosmic Rays, Praxis Publishing Ltd., 2004.
[15] S. Yoshida and M. Teshima, Progress of Theoretical Physics, 89, 4, pp. 833–845,
1993.
[16] J. L. Puget et al., The Astrophysical Journal, 205, pp. 638–654, 1976.
[17] T. Yamamoto et al., Astroparticle Physics, 20, pp. 405–412, 2004.
[18] D. Kuempel, arXiv :1409.3129.

84

REFERENCES

85[19] J. W. Cronin et al., Nuclear Physics B (Proceedings Supplements), 138, pp. 465–
491, 2005.
[20] E. Waxman and J. Miralda-Escudé, The Astrophysical Journal Letters, 472, pp.
L89–L92, 1996.
[21] A. M. Hillas, Annual Review of Astronomy and Astrophysics, 22, pp. 425–444, 1984.
[22] M. B. Hindmarsh and T. W. B. Kibble, Reports on Progress in Physics, 58, p. 477,
1984.
[23] C. T. Hill, Nuclear Physics B, 224, pp. 469–490, 1983.
[24] T. J. Weiler, Astroparticle Physics, 11, pp. 303–316, 1999.
[25] M. Birkel and S. Sarkar, Astroparticle Physics, 9, pp. 297–309, 1998.
[26] W. Heitler, The quantum theory of radiation, 3rd edition, Clarendon Press, 1954.
[27] T. K. Gaisser and A. M. Hillas, Proceedings of 15th International Cosmic Ray Conference, 8, pp. 353–357, 1977.
[28] K. Kamata and J. Nishimura, Progress of Theoretical Physics Supplements, 6, pp.
93–151, 1958.
[29] S. Yoshida et al., Journal of Physics G: Nuclear and Particle Physics, 20, pp. 651–
664, 1994.
[30] F. Kakimoto et al., Nuclear Instruments and Methods in Physics Research A, 372,
pp. 527–533, 1996.
[31] R. Abbasi et al., Astroparticle Physics, 29, pp. 77–86, 2008.
[32] J. Linsley, Physical Review Letters, 10, pp. 146–148,1963.
[33] M. Lawrence et al., Journal of Physics G, 17, pp.733–757, 1991.
[34] V. Egorova et al., Nuclear Physics B, 136, pp. 3–11, 2004.
[35] D. Bird et al., The Astrophysical Journal, 240, pp. 410–428, 1985.
[36] M. Takeda et al., Physical Review Letters, 81, pp. 1163–1166, 1998.
[37] R. Abbasi et al., Astroparticle Physics, 32, pp. 53–60, 2009.
[38] R. Abbasi et al., Physical Review Letters, 100, 101101, 2008.
[39] The Telescope Array collaboration and the Pierre Auger collaboration et al.,arXiv :1310.0647[40] R. Abbasi et al., Astroparticle Physics, 64, pp. 49–62, 2015.

REFERENCES

86[41] J. Abreu et al., Proceedings of 32nd International Cosmic Ray Conference, 2011.
arXiv :1107.4804v1[42] The Pierre Auger Collaboration, Science, 318, pp. 938–943, 2007.
[43] C. Macolino et al., Journal of Physics: Conference Series, 375, 052002, 2012.
[44] T. Abu-Zayyad et al., The Astrophysical Journal, 777, 88, 2013.
[45] R. Abbasi et al., The Astrophysical Journal Letters, 790, L21, 2014.
[46] A. Aab et al., The Astrophysical Journal, 794, 172, 2014.
[47] M. Risse and P. Homola, Modern Physics Letters A, 22, pp. 749–766, 2007.
[48] G. B. Gelmini, O. E. Kalashev, and D. V. Semikoz, arXiv:astro-ph/0506128.
[49] G. B. Gelmini, O. E. Kalashev, and D. V. Semikoz, Journal of Experimental andTheoretical Physics, 106, pp. 1061–1082, 2008.
[50] T. Erber, Review of Modern Physics, 38, pp. 626–659, 1966.
[51] L. D. Landau and I. Y. Pomeranchuk, The Collected Papers of L. D. Landau, Pergamon, 1965.
[52] A. N. Cillis et al., Physical Review D, 59, 113012, 1999.
[53] S. R. Klein , Radiation Physics and Chemistory, 75, pp. 696–711, 2006.
[54] H. P. Vankov et al., Physical Review D, 67, 043002, 2003.
[55] T. Abu-Zayyad et al., Nuclear Instruments and Methods in Physics Research A,
676, pp. 54–65, 2012.
[56] T. Abu-Zayyad et al., Nuclear Instruments and Methods in Physics Research A,
689, pp. 87–97, 2012.
[57] S. Kawana et al., Nuclear Instruments and Methods in Physics Research A, 681, pp.
68–77, 2012.
[58] Masaaki Kobayashi, Takao Shinkawa, Takahiro Sato, Nuclear Instruments and Methods in Physics Research A, 337, pp. 355–361, 1994.
[59] H. Tokuno et al., Nuclear Instruments and Methods in Physics Research A, 601, pp.
364–371, 2009.
[60] Y. Tameda et al., Nuclear Instruments and Methods in Physics Research A, 609,
pp. 227–234, 2009.

REFERENCES

87[61] T. Tomida et al., Nuclear Instruments and Methods in Physics Research A, 654, pp.
653–660, 2011.
[62] M. Chikawa et al., “Cloud Monitoring with an Infra-Red Camera for the Telescope
Array Experiment”, Proceedings of 31st International Cosmic Ray Conference, 2009.
[63] D. Heck et al., CORSIKA: A Monte Carlo Code to Simulate Extensive Air Showers,
Forschungszentrum Karlsruhe Report FZKA 6019, 1998.
[64] B. T. Stokes et al., Astroparticle Physics, 35, pp. 759–766, 2012.
[65] B. T. Stokes et al., “Using CORSIKA to quantify Telescope Array surface detector”,
Proceedings of 31st International Cosmic Ray Conference, 2009.
[66] J. Alcaraz et al., Physics Letters B, 490, pp. 27–35, 2000.
[67] J. Alcaraz et al., Physics Letters B, 494, pp. 193–202, 2000.
[68] T. Nonaka et al., “Performance of TA Surface Array”, Proceedings of 31st International Cosmic Ray Conference, 2009.
[69] A. Taketa et al., “The trigger and DAQ system of the surface detector array of
the Telescope Array experiment”, Proceedings of 31st International Cosmic Ray
Conference, 2009.
[70] NOAA/ESRL Radiosonde Database. http://www-frd.fsl.noaa.gov/mab/raob/[71] T. Abu-Zayyad et al., Astroparticle Physics, 61, pp. 93–101, 2015.
[72] D. Ikeda, “Hybrid analysis of ultra-high energy cosmic rays observed with the Telescope Array”, PhD thesis, University of Tokyo, 2009.
[73] T. Fujii, “Measurements of the Energy Spectrum and the Mass Composition of
Ultra-High Energy Cosmic Rays with Telescope Array Fluorescence Detectors in
Monocular Mode”, PhD thesis, Osaka City University, 2012.
[74] T. Abu-Zayyad et al., Physical Review D, 88, 112005, 2013.
[75] The Pierre Auger Collaboration, Astroparticle Physics, 31, pp. 399–406, 2008.
[76] K. Kampert et al., “Highlights from the Pierre Auger Observatory”, Proceedings of
32st International Cosmic Ray Conference, 2011.
[77] J. Abraham et al., Astroparticle Physics, 29, pp. 243–256, 2008.
[78] M. Ave et al., Physical Review Letters, 88, pp. 2244–2247, 2000.
[79] M. Ave et al., Physical Review D, 65, 063007, 2002.
[80] A. V. Glushkov et al., JETP Letters, 85, pp. 131–135, 2007.

REFERENCES

88[81] A. V. Glushkov et al., Physical Review D, 82, 041101, 2010.
[82] K. Shinozaki et al., The Astrophysical Journal, 571, L117, 2002.
[83] M. Risse et al., Physical Review Letters, 95, 171102, 2005.
[84] G. Gelmini et al., Journal of Experimental and Theoretical Physics, 106, 11, 2008.
[85] J. Ellis et al., Physical Review D, 74, 115003, 2006.
[86] M. Risse et al., Physics Letters A, 22, 749, 2007.