The energy spectrum and the chemical composition of primarycosmic rays with energies from 1014 to 1016 eV
Shoichi Ogio
Department of Physics
Tokyo Institute of Technology
February 2005

abstract
The energy spectrum of cosmic rays is well represented by a power law function of the primaryenergy. One of the significant features on the spectrum is a steepening around 1015 eV, whichis called the knee. The knee must be deeply related to the origin of the cosmic rays and/or
to the propagation mechanisms of primary cosmic rays in the Galaxy. Therefore, it is of crucial importance to obtain precise features of the primary energy spectrum as well as accurate
information on the primary chemical composition. In order to obtain a reliable conclusion forthe chemical composition of cosmic rays at the knee region, and to confirm our previous result,
we determined EAS longitudinal developments more directly with an equi-intensity method
analysis. We have measured extensive air showers with primary energies above 6 TeV at Mt.
Chacaltaya in Bolivia. The data were collected by an air shower array called the Minimum
Air Shower(MAS) array since March 2000. We applied the equi–intensity analysis method tothe extensive air showers extended over region of their maximum developments. We varied the
mixture of protons and iron nuclei in our simulations and compared the simulated equi–intensity
curves with the measured ones to determine the mixing ratio of protons as a function of the
primary energy. Moreover, on the basis of the obtained chemical composition, we derived theprimary energy spectrum from 1014 to 5 × 1016 eV. Consequently, we concluded that the powerlaw index of the spectrum is changing gradually around 1015.5 eV, and that the obtained proton
ratio decreases with increasing energy. We directly measured the longitudinal development of
air showers initiated by primaries with energies around the knee. We find that the average
mass number of primary cosmic rays increases with increasing energy above 1014.5 eV and thedominant component around the knee is not protons.
This result suggests that the steepening of the component spectra is caused at a single
rigidity of 1014.5 V. Moreover, it is consistent with our former result of the Čerenkov light observations. The measurement of Čerenkov light pulse shapes corresponded to an observationof the longitudinal developments of air showers at the earlier stages before the maximum developments. In contrast, with the equi–intensity method analysis we determined longitudinaldevelopments around their maximum and at the later stages. Therefore, with two different
observations, we measured the whole stage of EAS longitudinal developments and then we canobtain a trustworthy conclusion on the chemical composition.
On the basis of detailed simulation researches for the diffusive motions of charged particles
in the galactic turbulent magnetic field, I calculated the residence times of cosmic rays in thegalactic disk by solving one dimensional advective–diffusion equation. Consequently, this simple
model predicts the all–particle and the component spectra, the chemical composition, and theanisotropy, with assuming the natural values of the parameters for the equations such as themagnetic field strength and the galactic wind velocity.

Contents
1 Introduction

14

2 Cosmic rays: overview
18
2.1 Energy spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2 Chemical composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3 Propagation and acceleration of cosmic rays
3.1 Propagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Leaky box model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Acceleration of cosmic rays below 1014 eV . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Power . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 First order Fermi acceleration at supernova blast waves . . . . . . . . . .
3.3 Acceleration of cosmic rays above 1014 eV . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Acceleration of cosmic rays in a supernova shock which travels down astellar wind . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Reacceleration of cosmic rays in a termination shock of the galactic wind3.3.3 Reacceleration of cosmic rays with multiple shock encounters . . . . . . .
3.3.4 Production of cosmic rays in active galactic nuclei . . . . . . . . . . . . .
3.3.5 Acceleration by oblique shocks at SNRs . . . . . . . . . . . . . . . . . . .
3.3.6 Summary of acceleration models . . . . . . . . . . . . . . . . . . . . . . .
3.4 Advective diffusion propagation model for high energy cosmic rays above 1012 eV
3.4.1 Diffusive motion of charged particles in turbulent magnetic fields . . . . .
3.4.2 Advective–diffusion of cosmic ray particles . . . . . . . . . . . . . . . . . .
3.4.3 Residence time of cosmic rays in the galactic disk . . . . . . . . . . . . . .

22
22
24
25
26
26
30

4 Extensive Air Showers
4.1 Air showers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Interaction of electrons and gamma rays . . . . . . . . . . . . . . . . . . .
4.1.2 Longitudinal development . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.3 Lateral distribution of shower particles . . . . . . . . . . . . . . . . . . . .
4.2 Observations of extensive air showers . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Determination of arrival directions . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Determination of shower sizes . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.3 Muons and other components in air showers . . . . . . . . . . . . . . . . .
4.3 Equi–intensity method analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64
64
66
67
67
68
68
69
70
71

1

30
36
39
42
46
48
50
51
59
60

5 Review of instruments to measure chemical composition of cosmic rays
5.1 CRN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 RUNJOB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3 CASA–MIA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4 HEGRA–CRT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.5 CASA–BLANCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6 DICE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.7 BASJE–Čerenkov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.8 KASCADE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.9 Summary of recent measurements . . . . . . . . . . . . . . . . . . . . . . . . . . .

73
73
76
78
82
85
88
91
94
98

6 BASJE MAS array102
6.1 Site . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
6.2 Detectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.3 Electronics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
6.3.1 Local particle density measurements . . . . . . . . . . . . . . . . . . . . . 107
6.3.2 Fast timing measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
6.3.3 Trigger . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
6.3.4 Data acquisition system . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
7 Air shower analysis and characteristics of the array
114
7.1 Air shower analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
7.1.1 Determination of arrival direction . . . . . . . . . . . . . . . . . . . . . . . 114
7.1.2 Determination of shower size . . . . . . . . . . . . . . . . . . . . . . . . . 116
7.2 Performance of MAS array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
7.2.1 Effective area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
7.2.2 Resolutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
8 Experimental Results
127
8.1 Data selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
8.2 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
8.2.1 Parameterizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
8.2.2 Simulation scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
8.3 Size Spectra and Equi–Intensity Curves . . . . . . . . . . . . . . . . . . . . . . . 146
8.4 Chemical Composition and Energy Spectrum of Primary Cosmic Rays . . . . . . 151
8.4.1 Chemical Composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
8.4.2 Energy spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
9 Discussion
157
9.1 The systematic errors in the determination of chemical composition and primaryenergy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
9.2 Comparisons the results with the source models . . . . . . . . . . . . . . . . . . . 158
9.3 Comparisons the results with the advective diffusion propagation model . . . . . 163
9.3.1 All particle spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164
9.3.2 Spectra for the element groups . . . . . . . . . . . . . . . . . . . . . . . . 167
9.3.3 Average mass number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
9.3.4 Residence time and anisotropy . . . . . . . . . . . . . . . . . . . . . . . . 171

2

9.3.5

Summary of the advective diffusion model . . . . . . . . . . . . . . . . . . 173

10 Conclusions
177
10.1 Summary of observational results . . . . . . . . . . . . . . . . . . . . . . . . . . . 177
10.2 Cosmic ray source models and propagation models . . . . . . . . . . . . . . . . . 178
10.3 Future Prospects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178

3

List of Figures
1.1

1.2
2.1
2.2
2.3
3.1

3.2

3.3
3.4

3.5

The mean logarithmic mass hln Ai measured with balloon–borne detectors [7][8]and ground–based detectors [13][16][25][29][35][88][94], as a function of primaryenergy. Also the results of our former Čerenkov observations [82] are plotted.
A hatched region represents the result of other direct observations, which isaccumulated by Linsley [65]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Longitudinal development curves for primary proton and iron with the energies
1014 , 1015 and 1016 eV. Each curve is the average of 5000 simulated showers. . . 16
The cosmic ray energy spectrum over a very wide energy range. . . . . . . . . . 19
The cosmic ray elemental abundances(∼ 1010 eV) measured at the earth arecompared with the solar system abundances(all is relative to silicon). . . . . . . 20
Cosmic ray chemical composition below the knee measured with the balloon–borne instruments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Compiled data observed at the earth for (A) the B/C ratio and (B) the sub–Fe/Fe
ratio as a function of the energy(see Garcia–Munoz et al.,1984 [32] for references
to the experimental data.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Comparison between experimental measurements of the radioactive isotopes 26 Al
and 10 Be and calculated prediction [38]. (see Wefel,1988 [97] for references to theexperimental data.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Acceleration at the plane shock front. . . . . . . . . . . . . . . . . . . . . . . . .
The all particle spectrum and the various contributions from the different element
groups derived by Stanev et al. [86] are compared with the observational results.
Here the symbols denote the following experiments: Akeno(open circles), Haverah park(open squares),Yakutsk(open triangles, pointing down),Tien Shan(opentriangles, pointing up), Fly’s Eye(open hexagons),Proton 4(full squares), and
JACEE(full circles). The contributions of the six element groups are distinguished with the line type: H(solid), He(dots), CNO(dash), Ne–S(long dash),
Cl–Mn(dash dot), and Fe(long dash dot). . . . . . . . . . . . . . . . . . . . . . .
A comparison of the measured element abundances with the model calculationfor H(full squares), He(full circles), and Fe(full stars). . . . . . . . . . . . . . . .

4

23

25
27

34
35

3.6

(upper) Computed energy spectrum of accelerated iron nuclei at the equatorial
plane of the galactic wind, at radii of 10 kpc(lower curve) and at the shock(upper
curve). In this case the shock is at 300 kpc, and the outer absorbing boundary is
at 600 kpc. The wind velocity is 400 km s−1 , and the diffusion coefficient is taken
to be independent of energy, with a parallel diffusion mean free path of 14 times
the local gyro radius of a 1020 eV particle. The magnetic field is approximatedby an Archimedian spiral with a magnitude of 3 × 10−6 Gauss at a radial distanceof 10 kpc, and a rotational angular velocity of 2 × 10−15 s−1 . The particle drift
velocity is set equal zero. Injection is continuous, and uniform over the shock, atthe energy of 2 × 1017 eV. The spectrum is given at ∼ 7 × 109 years after the startof acceleration, when the distribution has reached a near steady state. (lower)
The computed spectrum as in the left figure, except that the mean free pathdepends on both the energy and position and is 25 times the local gyro radius of
the particles at any given energy. . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.7 The distinction between the multiple encounter acceleration model described in
the text(b), and the classical second order Fermi acceleration(a). . . . . . . . . .
3.8 The integral energy spectra is derived by using different values of the injection
rigidity(P0 ) and Nt . The background magnetic field is assumed to be 6µG. . . .
3.9 A synthetic spectrum of cosmic rays. The input parameters are P0 = 1014 eV and
B0 = 6µG. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.10 Spectra of produced neutron (full curves) and cosmic ray protons escaping froman AGN(dashed curve). Results are given for b = 10, x1 = 30, and LX =1042 (left–most curves), 1045 and 1048 ergs s−1 . The reduction at high energies isdue to interactions of neutrons with photons during the escape from the central
region, while the reduction at low energies is due to interactions of decay protons
with protons in the accreting plasma. . . . . . . . . . . . . . . . . . . . . . . . .
3.11 The possible contribution of cosmic rays accelerated in AGN to the observed
spectrum. Results are shown for b = 1(horizontal hatching), b = 10(thin obliquehatching) and b = 100(thick oblique hatching). The measured total cosmic rayintensity and the proton intensity are also shown in this figure. The dotted line
represents a single steady source contribution at 10 Mpc. . . . . . . . . . . . . .
3.12 The maximum energy Emax for a proton with three values of x versus magnetic
field inclination η. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.13 The predicted fluxes of total, proton, He, and other nuclear groups in the case of“Ours” parameter choice. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.14 Comparison between the all–particle spectra reported by various groups and thepredicted spectra with the oblique shock model. . . . . . . . . . . . . . . . . . . .
3.15 Energy dependence of chemical composition of primary cosmic rays in term of< ln A >. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.16 Left ordinate: the anisotropy amplitude; right ordinate: E 2.47 × dN/dE [45].
The markers and the solid line represent the observed anisotropy and the energy
spectrum, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.17 The assumed power spectrum of the magnetic field. . . . . . . . . . . . . . . . . .
3.18 Examples of the simulated field lines, (a) for Lirr = 10pc and (b) for 50pc. . . . .
3.19 An example of a simulated track of a teat particle (proton) with energies of
5 × 1014 eV for 5500 years. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

38
39
40
41

44

45
47
48
48
49

50
53
54
54

3.20 Snap shots of particle distributions after an explosion of a source(0,0) of particles
with energies 5×1014 eV. The magnetic field strength is assumed to be BL = 1µG
and δB/B = 1. The direction of the regular magnetic field is along with y–axis. .
3.21 Time dependences of the Fokker–Planck coefficients for the diffusive motions ofcharged particles with energies 5 − 9 × 1014 eV in the assumed galactic magnetic
field. Dxx and Dyy represent the coefficients for parallel and perpendicular to
the regular magnetic field, respectively. We assumed BL = 1.0µG, δB/B = 1
and Lirr = 10pc. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.22 Time dependences of the Fokker–Planck coefficients for the diffusive motions ofcharged particles with energies 5 − 9 × 1014 eV in the assumed galactic magnetic
field. Dxx and Dyy represent the coefficients for parallel and perpendicular to
the regular magnetic field, respectively. We assumed BL = 1.5µG, δB/B = 1
and Lirr = 10pc. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.23 The obtained diffusion coefficients depend on the Larmor radii. The diffusion
coefficients and the Larmor radii are normalized with cLirr and Lirr , respectively
and c is the light speed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.24 (upper) Illustration of the magnetic loop formation due to the instability and of
relativistic gas expulsion from the galactic disk [67]. (lower) A sketch of dark
filaments found in an optical image of spiral galaxy NGC253 [90]. . . . . . . . .
3.25 The residence time τR predicted with the advective diffusion model for cosmic ray
protons. This figure shows τR s for different advection velocities, i.e., the galactic
wind velocity, vg . For all the calculations in this figure Lirr is fixed to be 100 pc.
3.26 The residence time τR predicted with the advective diffusion model for cosmic
ray protons. This figure shows τR s for different maximum scale lengths of theirregularity of the galactic magnetic field, lirr . For all the calculations in thisfigure vg is fixed to be 5 × 10−4 pc/years = 490km/s. . . . . . . . . . . . . . . . .
3.27 The residence time τR predicted with the advective diffusion model for different
primary components of cosmic ray particles. For all the calculations in this figure
vg is fixed to be 5−4 pc/years = 490km/s and Lirr is fixed to be 50 pc. . . . . . .

55

56

57

58

60

62

62

63

4.1
4.2
4.3

A schematic view of an extensive air shower. . . . . . . . . . . . . . . . . . . . . 65
A cross–sectional view of a front plane of an air shower and detectors. . . . . . . 69
The determination of the development of EAS of different sizes through the atmosphere from observations at different zenith angles θ; n(N) is the rate ofoccurrence of shower size N. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

5.1
5.2

Schematic cross section of the CRN detector . . . . . . . . . . . . . . . . . . . .
(left): Differential energy spectra for the cosmic ray nuclei C, O, Ne, Mg, Si
and Fe measured with CRN(solid). Note that the fluxes multiplied with E 2.5 .
For other markers, see the reference [71]. (right): Abundances of the galactic
cosmic ray nuclei at the source(GCRS), relative to the local abundances(LG),are plotted against the values of first ionization potentials of the elements [71].
All abundances are relative to iron. . . . . . . . . . . . . . . . . . . . . . . . . .
Chamber structure in 1995 and 1996 experiments. . . . . . . . . . . . . . . . . .
The left figure is the measured spectra of the primary protons and primary He
nuclei. The right figure is the measured spectra of the heavy components. . . . .

5.3
5.4

6

74

75
76
77

5.5

The layout of CASA–MIA array. CASA is on the surface, and MIA is buried at
3 meter underground. The CASA stations are 15 m apart from each other. . . .
5.6 The energy parameter log(Ne + 60 × Nµ ) as a function of the energy for simulated
proton(open circle) and iron(filled triangle) vertical showers. . . . . . . . . . . .
5.7 The observed energy spectrum of the CASA–MIA experiment compared withthose of Tibet and of Akeno. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.8 The average composition hln Aiof CASA–MIA data(filled circle). In the simulations for this analysis the QGSJET interaction generator were used. The open
squares indicate the mean mass of the direct measurement of JACEE. . . . . . .
5.9 The energy spectra of each composition groups measured by the CASA–MIA
experiment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.10 Map of the HEGRA site showing only the scintillator stations (square dots) and
CRTs (circles). The dash-dotted line marks the area of shower core positions
accepted in the analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.11 Histograms of muon radial angles observed with CRT detectors in different intervals of core distances for 15000 < Ne < 50000. Note that particles coming fromthe shower axis have negative radial angles. . . . . . . . . . . . . . . . . . . . . .
5.12 Median radial angles of muons measured with CRT detectors as functions ofhorizontal distance of the CRT detectors from the shower cores. Simulations
are superimposed for pure protons and irons as well as the mixed compositions
corresponding to direct measurements of the JACEE collaboration both above
45 TeV(J45) and above 370 TeV(J370). . . . . . . . . . . . . . . . . . . . . . . .
5.13 The mean mass composition obtained by comparing HEGRA–CRT data withsimulations based on the VENUS/GHEISHA and DPM/Isobar interaction models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.14 The CASA–BLANCA array. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.15 The differential all particle cosmic ray flux measured by CASA–BLANCA. . . .
5.16 The mean logarithmic mass < ln A > measured by CASA–BLANCA as a functionof the energy. The four sets of symbols show the BLANCA data interpreted using
CORSIKA coupled with the indicated hadronic interaction model. . . . . . . . .
5.17 The spectra near the knee by the DICE group and by some other groups. . . . .
5.18 (upper) The result for the mean mass measurement by DICE and the model
predictions. (lower) The result of (p+α)/all. . . . . . . . . . . . . . . . . . . . .
5.19 The arrangement of the scintillation detectors of the SAS array and the Čerenkov
detectors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.20 Comparison of the cumulative fraction distribution of experimentally determined
T10−90 (dots) with the distributions expected for proton, carbon and iron primaries. These are result for (a) log Ne = 6.0 − 6.2, R = 140 − 150 m obtained in
the low energy mode observations, and (b) log Ne = 6.5 − 7.0, R = 150 − 160 m
obtained in the high energy mode observations. . . . . . . . . . . . . . . . . . .
5.21 Comparison of BASJE–Čerenkov results with the other direct and indirect experimental results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.22 A schematic layout of the KASCADE experiment. . . . . . . . . . . . . . . . . .
5.23 Electron size spectra (upper) and truncated muon size spectra(lower) in the intervals of zenith angle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.24 The total energy spectrum obtained with the KASCADE (Ne , Nµtr. ) analysis. .

7

78
79
80

80
81

82

83

84

84
85
86

87
89
90
91

92
93
94
95
96

5.25 The energy spectra of proton, helium and carbon primaries between 1 PeV to 10
PeV as a result of the KASCADE (Ne , Nµtr. ) analysis. . . . . . . . . . . . . . . . 97
5.26 The mean logarithmic mass by the KASCADE (Ne , Nµtr. ) analysis. . . . . . . . 98
5.27 Two examples of the hadronic observables measured with the KASCADE hadroncalorimeter for estimation of the primary cosmic ray mass composition. . . . . . 99
5.28 The mass parameters obtained from the six observables of the KASCADE experiment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.29 The mean logarithmic mass vs the primary energy reported by KASCADE(hadron
analysis) and by other experiments. . . . . . . . . . . . . . . . . . . . . . . . . . 101
6.1
6.2

BASJE MAS array. Only N–, F– and G–detectors are shown in this figure. . . . 103
The central area of MAS array. The 60m2 muon detector(the shielded detector)
is located at the center of the array (under the L–detectors). . . . . . . . . . . . 104
6.3 The schematic illustrations of the surface detectors of the MAS array. . . . . . . 105
6.4 The shielded muon detector has the fifteen 4 m2 scintillation detectors. . . . . . . 106
6.5 Block diagram of the MAS array electronics. . . . . . . . . . . . . . . . . . . . . 107
6.6 Log–amplification method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
6.7 The diagram of the high voltage adjustment. . . . . . . . . . . . . . . . . . . . . 108
6.8 The diagram of T0 measurement. . . . . . . . . . . . . . . . . . . . . . . . . . . 109
6.9 The diagram of τ measurement. . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
6.10 An example of τ measurement. The filter attenuation is logarithmically proportional to the pulse height. In this example, the τ value is determined 9.82 µs for
the F–1 detector. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.11 An example of the measurement for the linear relation between time differences
and TDC counts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
6.12 The block diagram of the data acquisition processes for the MAS array. . . . . . 113
7.1

Delay time of shower particles from the shower plane measured at Mt. Chacaltaya
with SAS array. The line in the figure is a calculated delay time with the formula(7.7). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
7.2 The fluctuation (σρ ) of shower particle density ρ. σρ is independent of air shower
size. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
7.3 The effective area when MAS array is under the trigger condition described in
the previous section. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
7.4 The expected energy distribution of the triggered events. The mode energy is
about 30 TeV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
7.5 Distributions of ∆θ of vertically incident primaries. All the events are selectedwith the condition r ≤ 20 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
7.6 Distributions of ∆θ of the simulated showers of θ = 44.4◦ . . . . . . . . . . . . . 120
7.7 Distributions of ∆θ for vertically incident primaries. All the events are selectedwith the condition r > 20 m. The angular determination errors of the events on
the outer region of the array is large by about 50 % compared to Figure 7.5. . . 121
7.8 An examples of the fitting ∆θ distribution with the formula (7.18). σθ is determined to be 1.94◦ for this distribution. . . . . . . . . . . . . . . . . . . . . . . . . 121
7.9 The plot of angular determination errors σθ as a function of primary energy. . . 122
7.10 Distributions of the core position difference for the vertically incident primaries.
All the events are selected with the condition r ≤ 20 m. . . . . . . . . . . . . . . 123
8

7.11 Distributions of the core position difference for the simulated showers of θ = 44.4◦ .
All the events are selected with the condition r ≤ 20 m. . . . . . . . . . . . . . . 124
7.12 Distributions of the size determination error for the vertically incident primaries.
All the events are selected with the condition r ≤ 20 m. . . . . . . . . . . . . . . 124
7.13 Distributions of the size determination error for the simulated showers of θ =44.4◦ . All the events are selected with the condition r ≤ 20 m. . . . . . . . . . . 125
7.14 The plot of the core distance resolution as a function of primary energy. The
markers indicate the systematic errors, and the error bar of each plot is thestatistical error. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
7.15 The plot of the size determination error as a function of primary energy. The
markers indicate the systematic errors, and the error bar of each plot is thestatistical error. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
8.1

Distribution of the reduced chi–square for the determinations of the arrival directions, (7.5). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
8.2 Distribution of the reduced chi–square for the determinations of the EAS sizes,(7.9). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
8.3 Scatter plots of χ2size vs EAS size. . . . . . . . . . . . . . . . . . . . . . . . . . . 129
8.4 Distributions: (first column) of the determination errors of core positions(|r−r0 |),(second column) of arrival directions(∆θ) and (third column) of EAS sizes(log(N/N0 ))for the full Monte Carlo simulated events with CORSIKA-QGSJET. . . . . . . . 130
8.5 Effective detection area as a function of the primary energy for primary protons
simulated with CORSIKA. The solid lines are drawn with the function (8.1). . . 132
8.6 Effective detection area as a function of the primary energy for primary irons
simulated with CORSIKA. The solid lines are drawn with the function (8.1). . . 133
8.7 The fitting parameters of the function (8.1) for the simulated primary protons.
The solid lines represent the polynomial functions of sec θ . . . . . . . . . . . . . 133
8.8 The fitting parameters of function (8.1) for the simulated primary irons. The
solid lines represent the polynomial functions of sec θ . . . . . . . . . . . . . . . 134
8.9 Effective detection area as a function of the primary energy for primary protons
simulated with CORSIKA. The solid lines are effective area calculated with thefunction (8.1) and with the polynomial functions of S1 , S2 and S3 . . . . . . . . . 134
8.10 Effective detection area as a function of the primary energy for primary irons
simulated with CORSIKA. The solid lines are effective area calculated with thefunction (8.1) and with the polynomial functions of S1 , S2 and S3 . . . . . . . . . 135
8.11 Angular resolutions as a function of the primary energy for primary protons
simulated with CORSIKA are fitted with a linear function. . . . . . . . . . . . . 136
8.12 Angular resolutions as a function of the primary energy for primary irons simulated with CORSIKA are fitted with a linear function. . . . . . . . . . . . . . . 137
8.13 The plots of the fitting parameters for the angular resolutions. . . . . . . . . . . 137
8.14 The plots of the reconstructed shower size as a function of the primary energy
for each sec θ for primary protons simulated with CORSIKA. . . . . . . . . . . . 138
8.15 The plots of the reconstructed shower size as a function of the primary energy
for each sec θ for primary irons simulated with CORSIKA. . . . . . . . . . . . . 139
8.16 The fitting parameters of log10 Nrec –log10 E relations for primary protons, forlow energy(line 1), middle energy(line 2) and high energy region(line 3). These
parameters are fitted with a polynomial function of sec θ. . . . . . . . . . . . . . 140
98.17 The fitting parameters of log10 Nrec –log10 E relations for primary protons, forlow energy(line 1), middle energy(line 2) and high energy region(line 3). These
parameters are fitted with a polynomial function of sec θ. . . . . . . . . . . . . . 141
8.18 Size resolutions as a function of the primary energy for primary protons simulatedwith CORSIKA are fitted with a linear function. . . . . . . . . . . . . . . . . . . 142
8.19 Size resolutions as a function of the primary energy for primary protons simulatedwith CORSIKA are fitted with a linear function. . . . . . . . . . . . . . . . . . . 143
8.20 The fitting parameters of σN –log10 E relations for primary protons and irons.
These values are fitted with a linear function of sec θ. . . . . . . . . . . . . . . . 143
8.21 The distributions of the simulated sizes for primary protons with the fixed several
primary energies and incident zenith angles. The solid lines are those simulatedwith CORSIKA–QGSJET no–thinning mode, and the dotted lines are those simulated with the parameterized functions derived in the text. . . . . . . . . . . . 145
8.22 The distributions of the simulated sizes for primary protons with the fixed several
primary energies and incident zenith angles. The solid lines are those simulatedwith CORSIKA–QGSJET no–thinning mode, and the dotted lines are those simulated with the parameterized functions derived in the text. . . . . . . . . . . . 145
8.23 The azimuth angle distributions of particular sec θ bins. . . . . . . . . . . . . . . 146
8.24 The observed integral EAS size spectra for each sec θ bin. . . . . . . . . . . . . . 147
8.25 The differential size spectrum for sec θ = 1.0 − 1.1 measured by BASJE–MAS
array. The intensity has been multiplied by N 1.5 to flatten the spectrum. . . . . 148
8.26 (a)The observed equi–intensity curves are compared with simulated ones for primary protons, iron nuclei and (50% protons + 50% irons), without any normalization. The attached value for each curve is the integral flux. In the analysis, we derived the curves for F (> N ) from 10−5.2 to 10−8.2 m−2 sr−1 s−1 with
10−0.25 step, but in this figure we show a part of these curves to avoid confusions.
(b)Simulated equi–intensity curves with the measured ones without any normalization. In the simulation, we took consideration of not only the obtained mixing
ratios of protons and irons, but also obtained all particle energy spectrum. . . . 149
8.27 The observed equi–intensity curves are compared with simulated ones. The attached value for each curve is the integral flux. The simulated curves are normalized at 578 g/cm2 with the measured curves. . . . . . . . . . . . . . . . . . . 150
8.28 The systematic and the statistical errors of np /(np + nFe ). The errors are estimated for simulated data with equivalent number of events as observed. . . . . . 153
8.29 The mean logarithmic mass hln Ai measured by BASJE–MAS array as a functionof the primary energy, are compared with the results of other experiments withballoon–borne detectors(JACEE [8], RUNJOB [7]) and ground–based detectors(CASA–
MIA [35], KASCADE(hadrons) [25], HEGRA–CRT [13], KASCADE(electrons)[94], CASA–BLANCA [29], DICE [88], Fly’s Eye [16]). Also the results of ourformer Čerenkov observations [82] are plotted. A hatched region represents theresult of other direct observations, which are accumulated by Linsley [65]. . . . 154
8.30 The relation between the primary energies and the reconstructed EAS sizes with
our array with 1.0 ≤ sec θ < 1.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

108.31 The differential all–particle cosmic ray flux measured by BASJE–MAS array.
Also plotted are the cosmic ray fluxes reported by JACEE [8], RUNJOB [7],SOKOL [52], proton satellite [36], KASCADE(hadrons) [25], KASCADE(electrons)[94], CASA–BLANCA [29], CASA–MIA [34], DICE [88], Tibet [4], EAS–TOP [1],and the dashed line represents the flux measured by Akeno group [72]. . . . . . 156
9.1

These plots are simulated results of the relation between the mean energies andthe atmospheric depths on the equi–intensity curve corresponding to a certain
rate of arrival. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
9.2 These plots are simulated results of the relation between the mixing ratio andthe atmospheric depths on the equi–intensity curve corresponding to a certain
rate of arrival. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
9.3 All–particle spectrum and the contributions of five components calculated withmodel (a)(upper panel) and with model (b)(lower panel) are compared with themeasured spectra. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
9.4 The predicted mean logarithmic mass hln Ai with the model (a)(solid line), and(b)(dashed line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
9.5 Calculated all particle spectra with the advective–diffusion models of J − 1(solidline), J − 2(thick dashed line) and J − 3(dotted line). The markers represent theobserved energy spectra, by SOKOL [52](¤), proton–satellite [36](°), RUNJOB[7](¥), JACEE [8](4), CASA–MIA [34](+), KASCADE–hadron [25](?), Tibet[4](∗), CASA–BLANCA [29](H), DICE [88](♦), EAS–TOP [1](F), KASCADE–
electron [94](N) and this work(BASJE-MAS)(•). Thin dashed line represents theflux measured by Akeno group [72]. . . . . . . . . . . . . . . . . . . . . . . . . . 165
9.6 Calculated all particle spectra with the models of J − 4(solid), J − 5(thick dashed)
and J − 6(dotted). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
9.7 Calculated all particle spectra with the models of RC − 1(solid), RC − 2(thickdashed) and RC − 3(dotted). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
9.8 Calculated all particle spectra with the models of RF − 1(solid), RF − 2(thickdashed) and RF − 3(dotted). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
9.9 The lines are the calculated energy spectra for the assumed five components withthe advective–diffusion models of J − 1(solid), J − 2(dashed), J − 3(dotted). The
marker plots represent the observed energy spectra of five groups of cosmic ray
species by the direct measurements of Ryan et al. [78](♦), SOKOL [52](°),
JACEE [8](•), CRN [71](+) and RUNJOB [7](H). . . . . . . . . . . . . . . . . . 167
9.10 The lines are the calculated energy spectrum for each component with the modelsof RC − 1(solid), RC − 2(dashed), RC − 3(dotted). . . . . . . . . . . . . . . . . . 168
9.11 The lines are the calculated energy spectrum for each component with the modelsof RF − 1(solid), RF − 2(dashed), RF − 3(dotted). . . . . . . . . . . . . . . . . . 168
9.12 The averaged mass number, < ln A >, calculated with the models of J − 1(solid),J − 2(dashed) and J − 3(dotted). The estimated values with various different
types of measurements [82]–[16] are also plotted in this figure. A hatched region represents the result of other direct observations, which are accumulated by
Linsley [65]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
9.13 The averaged mass number, < ln A >, calculated with the models of RC − 1(solid),
RC − 2(dashed) and RC − 3(dotted). . . . . . . . . . . . . . . . . . . . . . . . . 170

11

9.14 The averaged mass number, < ln A >, calculated with the models of RF − 1(solid),
RF − 2(dashed) and RF − 3(dotted). . . . . . . . . . . . . . . . . . . . . . . . . . 170
9.15 (a)The residence times calculated with the models of J − 1(solid), J − 2(dashedline) and J − 3(dotted). (b)The measured anisotropy amplitude compared withthe estimated values from τR on (a). The estimated anisotropies are calculatedwith a relation A ∝ 1/τR and a normalization factor obtained at 1012 eV. The
measured anisotropies are accumulated by Hillas [45](°), Linsley [65](•) andHayakawa [40](N). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
9.16 (a)The residence times calculated with the models of RC −(solid), RC − 2(dashedline) and RC − 3(dotted). (b)The measured anisotropy amplitude compared withthe estimated values from τR in (a). The estimated anisotropies are calculatedwith a relation A ∝ 1/τR and a normalization factor obtained at 1012 eV. . . . . 172
9.17 (a)The residence times calculated with the models of RF − 1(solid), RF − 2(dashedline) and RF − 3(dotted). (b)The measured anisotropy amplitude compared withthe estimated values from τR in (a). The estimated anisotropies are calculatedwith a relation A ∝ 1/τR and a normalization factor obtained at 1012 eV. . . . . 172
9.18 Calculated all particle spectrum and the contributions of five components withthe model J − 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174
9.19 Calculated all particle spectrum and the contributions of five components withthe model RC − 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174
9.20 Calculated all particle spectrum and the contributions of five components withthe model RF − 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
9.21 The schematic view of three cosmic propagation models. . . . . . . . . . . . . . 176

12

List of Tables
3.1

3.3
3.4

The maximum rigidity and the maximum energy of iron nuclei of cosmic rays
accelerated at the galactic wind termination shock. . . . . . . . . . . . . . . . .
Summary of the characteristic predictions of acceleration models for cosmic rays
> 1014 eV. †: these values depend on the assumption of the diffusive motions ofcosmic rays. ∗: it depends on the boundary condition, that is, the abundances
and the spectral indices for the components around 1 TeV. . . . . . . . . . . . .
The fitting parameters of equation (3.101) for the calculated diffusion coefficients.
Assumed five components of cosmic ray particles. . . . . . . . . . . . . . . . . . .

6.1

the characteristics of the each type of the detectors. . . . . . . . . . . . . . . . . 106

8.1

The integral fluxes, corresponding primary energies, the average EAS sizes forsec θ = 1.0 − 1.05 and estimated values of np /(np + nFe ) for each energy. The
unit of integral flux F (> N ) is m−2 sr−1 s−1 . . . . . . . . . . . . . . . . . . . . . 152
The systematic errors of np /(np + nFe ). . . . . . . . . . . . . . . . . . . . . . . . 152

3.2

8.2
9.1

9.2
9.3

37

49
53
60

The table of the integral fluxes, corresponding primary energies Erate , the average
EAS sizes and corresponding energy Esize The unit of integral flux F (> N ) is
m−2 sr−1 s−1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
Assumed relative abundance and three models of the spectral indices for fivecomponent of cosmic ray particles. . . . . . . . . . . . . . . . . . . . . . . . . . . 163
The model ID and the corresponding sets of the parameters pairs, Lirr and vg . . 164

13

Chapter 1

Introduction
The study of high energy cosmic ray particles and radiations is great interest in astrophysics.
Cosmic rays are believed to be concerned with the most energetic phenomena in the Galaxy
and the universe. The bulk of cosmic rays are nuclei, and the primary energies of observed
cosmic rays range up to ∼ 1020 eV. The origin of cosmic rays remains a mystery, however, there
is general agreement that cosmic ray sources are located in our Galaxy at least for cosmic rays
with energies below 1014 eV. Supernovae and neutron stars are eminently promising as potential
cosmic ray sources. The goal of cosmic ray physics is to clarify the cosmic ray origin and theproperties of the acceleration mechanisms, and those of the propagation in the Galaxy.
The origin of primary cosmic rays with energies above 1014 eV is still unknown. In orderto solve this mystery, we must obtain at least three accurate pieces of information on primarycosmic rays, that is, pieces are the energy spectrum, the anisotropy in the distribution of arrival
directions, and the primary chemical composition of cosmic rays. Many groups have reported on
the isotropy in arrival directions of cosmic rays. There is not any special regions nor astronomical
objects that emit cosmic rays intensely. Direct observations by balloon-borne experiments on
the energy spectrum and on the chemical composition have been performed because balloonborne experiments are efficient in the energy range below 1014 eV. Above 1014 eV, the steep
falling spectrum requires large detection area or long exposure times. A large detection area
can be realized with observations of extensive air showers(EASs) with ground-based detector
arrays.
The energy spectrum of cosmic rays is well represented by a power law function of the
primary energy, and one of the significant features on the spectrum is a steepening around 1015
eV, which is called the knee. The knee must be deeply related to the acceleration mechanismsof cosmic rays and/or to the propagation mechanisms of primary cosmic rays in the Galaxy.
In order to overcome the difficulty in the acceleration limit Emax ∼ 1014 eV for proton,
theorists have proposed many types of models, involving such a mechanism as post–accelerations
in supernova remnants [9] after the energy gain by the direct shock acceleration in supernova
remnants, or introducing a new source [74], and so on. Each model predicts a somewhat different
feature in the composition near the knee. For instance, Axford [9] proposed that the cosmic
ray components above the knee are primarily the same as those below the knee, and the energy
spectra of the latter with the cutoff at Z × Emax (Z is a charge of a particle) are boosted wellover the knee, due to the multiple collisions with large scale shock–waves in the interstellar
medium. Therefore the average mass of the primary elements does not change so drasticallyaround the knee. Alternatively, several authors propose that some new components might give14

< ln A >6

5

JACEE (direct)
RUNJOB (direct)
CASA BLANCADICECASA MIAKASCADE-electrons
KASCADE-hadrons
HEGRA (CRT)
Flys Eye
BASJE-Cherenkov (1996)

4

3

2

1

0 12
10

10

13

10

14

10

15

10

16

10

17

10

18

10

19

Energy (eV)Figure 1.1: The mean logarithmic mass hln Ai measured with balloon–borne detectors [7][8]and ground–based detectors [13][16][25][29][35][88][94], as a function of primary energy. Also
the results of our former Čerenkov observations [82] are plotted. A hatched region represents
the result of other direct observations, which is accumulated by Linsley [65].
rise to recovery in the energy spectrum above the knee. In this case, we expect the compositionchanges drastically either into heavier one or into lighter one above the knee.
Therefore, it is of crucial importance to obtain accurate information on the primary chemicalcomposition, because the chemical composition of cosmic rays and its variation with energy
reflect the acceleration mechanisms and the compositions at acceleration regions, and reflect
propagation mechanisms in the Galaxy. However, although many groups have reported theirresults as shown in Figure 1.1, these results are contradictory each other and the chemicalcomposition had not been determined definitely.
Unfortunately, the determinations of the energy spectrum and those of the composition
from EAS experiments are interdependent and affected by an adopted model of high energy
interactions. Hence, for EAS experiments, it is important to observe not only the EAS size atan array altitude, but also the EAS longitudinal development features, because longitudinal development curves reflect the composition as indicated in Figure 1.2. Since one of these methods
is a measurement of Čerenkov radiation induced by EAS particles, we observed Čerenkov lightpulse shapes for air showers with energies from 1015 to 1016.5 eV. Observations were carriedout from 1995 until 1997, and we determined the cosmic ray chemical composition [82]. Oneof the other approaches to determine the composition is to measure the number of muons andits correlation with air shower sizes. Thus most of the composition measurements are basedon the additional observables to the EAS size. Then, unfortunately, the determination of the15

shower size

7

10


6

10

10

5

10

4

10

3

10

2



10

1
0

200

400

600

2

800

1000

depth (g/cm )Figure 1.2: Longitudinal development curves for primary proton and iron with the energies
1014 , 1015 and 1016 eV. Each curve is the average of 5000 simulated showers.
chemical composition from EAS experiments is strongly affected by Monte Carlo simulation
code, particularly, models of nuclear interactions. This is one of the serious problems.
In the present work, we obtained the equi-intensity curves from the observed EAS size
spectra with various incident zenith angles, θ. The equi-intensity method analysis is based only
on air shower size spectra. According to Monte–Carlo simulations, using the CORSIKA code[41] with the QGSJET [55] hadronic interaction model, EASs initiated by protons of 1015.5 eVreach their maximum developments around 590 g/cm2 atmospheric depths and those initiatedby iron nuclei around 450 g/cm2 atmospheric depth [62]. According to Knapp et al. [62],there are systematic differences between interaction models on the order of 50g/cm2 for the
maximum development point(Xmax ) of protons, and the QGSJET model shows the most rapid
development among the major hadronic interaction models [42].
The cosmic ray observatory at Mt. Chacaltaya in Bolivia is located at an atmospheric depthof 550 g/cm2 so that we can observe EASs initiated by primary protons around the knee region
before their maximum developments. Thus, our site is most suitable to investigate the chemicalcomposition of cosmic rays around the knee. We installed an air shower array, which is called“Minimum Air Shower (MAS) array” in this observatory and started the observation of EASs
with energies above 6 TeV in 2000 [99].
The measurement of Čerenkov light pulse shape corresponds to observation of the longitudinal developments of air showers at the earlier stages before the maximum developments.
In contrast, with the equi-intensity method analysis we determine longitudinal developments

16around their maximum and at the later stages. Therefore, with two different observations,
we can measure the whole stage of EAS longitudinal developments and then we can obtain atrustworthy conclusion on the chemical composition.
The equi–intensity method analysis is based on a single observable value: air shower sizes,
and then, in comparison with the other methods for estimation of the chemical composition
with the amounts of muon component or Čerenkov photons. Therefore, this method is lessdependent on Monte–Carlo simulations.
In the following section I give a brief review of cosmic ray physics and unresolved problems
in it. Chapter 3 is a brief review of the accelerations and the propagation of cosmic rays
with energies below 1014 eV, and is a review of some proposed models of the accelerations ofcosmic rays with energies greater than 1014 eV. In Chapter 4, I describe general properties of air
shower phenomena produced by cosmic rays incident to the atmosphere and techniques for theirdetections. Before a detailed description of our air shower array, we give a brief review of theother cosmic ray instruments in Chapter 5. In Chapter 6 and 7 , I briefly describe details of theexperimental apparatus, the analysis procedure and the performance of the array. In Chapter
8, with the classifications of the selected EASs into corresponding sec θ bins, we obtain EAS
size spectra for different sec θ bins, and derive EAS longitudinal development curves with theequi–intensity method. Moreover we examine the primary chemical composition with comparing
the observed EAS longitudinal development curves with those calculated, and then derive theenergy spectrum of primary cosmic rays. Chapter 9 is discussions on our result about systematic
uncertainty, and comparisons between our results and the expectations with acceleration andpropagation models. Finally the present work is summarized in Chapter 10.

17

Chapter 2

Cosmic rays: overview
Cosmic ray particles hit the earth’s atmosphere at a rate of about 1000 per square meter per
second. They are primarily composed of ionized nuclei; about 90 % protons, 9 % alpha particles
and the rest heavier nuclei. The flux of cosmic ray electrons is about 1 % of the proton flux.
High energy gamma rays are also observed at the Earth, but their flux does not exceeds 10−3
compared to that of nuclear components. Most cosmic rays are relativistic, and a very few ofthem have ultra relativistic energies extending at least to 1020 eV.
A fundamental mystery of cosmic ray physics is their origins and acceleration mechanisms ofthem. The answer to the question of the origin of cosmic rays is not yet known well. However,
it is clear that nearly all of them come from the outside of the solar system, but from the insideof the Galaxy, because we know an anti–correlation between the solar activity and the flux ofcosmic rays at energies less than 100 GeV; they are more effectively excluded from the solar
neighborhood during the periods when expanding magnetized plasma from the sun – the solar
wind – is most intense. In contrast, ultra high energy cosmic rays with energies above ∼ 1018
eV may be of extra-galactic origin, because the gyro-radius of these cosmic rays in the galacticmagnetic field is equivalent to the galactic disk thickness, ' 300 pc.

2.1

Energy spectrum

Figure 2.1 is the cosmic ray energy spectrum over a very wide energy range [92]. This is fitted
well with an inverse power law function, and a formula of the differential flux is represented bydN∝ E −γdE(2.1)

The spectrum continues at least up to 1020 eV with a roughly constant index γ ' 2.8. However
we can see slight bends around 1015 eV and 1019 eV. These are called knee and ankle, respectively,
using analogy with a human leg.
In some theories, the knee is considered to reflect propagation nature of cosmic rays in the
Galaxy. The Larmor radius of a particle with an energy of 1016 eV and with the charge Z is
∼ 3(1/Z) pc in the galactic magnetic field, and this radius is equivalent to a typical size ofa magnetic field turbulence in the galactic disk. Therefore, above the energy of the knee theleakage of cosmic rays from the Galaxy may become faster with increasing primary energies.
In other theories, the knee is considered to reflect the size and the magnetic filed strength ofsources of cosmic rays. However, we expect the maximum energy of cosmic rays accelerated
18

Figure 2.1: The cosmic ray energy spectrum over a very wide energy range.
with a typical supernova remnant to be ∼ 1014 × Z eV. Although many detailed studies of theknee are carried out, there are many possible interpretations of the knee with including leakages
from the galactic disk, with energy limits of the acceleration mechanism or with new additional
components.
The ankle is interpreted as the point of the intersection of the galactic cosmic ray spectrumwith the extra-galactic one. The Larmor radius of cosmic rays with an energy E in the Galaxy
is calculated to be ∼ 3(E/1019 eV)(1/Z) kpc, thus cosmic rays propagate more freely, and thedeflection angle is less than one radian. Therefore, when we assume that these cosmic rays areaccelerated in galactic sources, we expect a very strong cosmic ray anisotropy(the anisotropy
amplitude is ∼ 100 %) with a single source. However, measured anisotropy amplitudes are small
above the ankle and this naturally leads that extra-galactic origins are dominant.

2.2

Chemical composition

The relative abundance of nuclear elements in cosmic rays provides important information tosolve the problem of origins and the acceleration sites of cosmic rays. A comparison of theelemental abundance of the low energy cosmic rays( E ∼ 1010 eV) and that of the solar system
is shown in Figure 2.2 [30]. There are two striking differences between two compositions. First,
nuclei with Z > 1 are much more abundant relative to protons in cosmic rays than those in
the solar system material. This is not explained clearly, but it could reflect genuine differencein composition at the source or the fact that hydrogens are relatively hard to be ionized for
injection into the acceleration process

19Figure 2.2: The cosmic ray elemental abundances(∼ 1010 eV) measured at the earth are compared with the solar system abundances(all is relative to silicon).

20

Figure 2.3: Cosmic ray chemical composition below the knee measured with the balloon–borne
instruments.
Second, the two groups of elements Li, Be, B and Sc, Ti, V, Cr, Mn are many orders ofmagnitude abundant in cosmic rays than in the solar system material. These elements are
essentially absent as the end products of stellar nucleo-synthesis. Li-, Be- and B-nuclei arenevertheless present in cosmic rays as the spallation products of the abundant nuclei: carbons
and oxygens. In the same way, Sc-, Ti-, V-, Cr- and Mn-nuclei are the spallation products ofirons; they are produced by collisions of cosmic rays in the interstellar medium. With the crosssections for the spallation, one can calculate total amount of matter traversed by cosmic rays
between a source and the Earth. (The details are described in Chapter 3.)
The cosmic ray mass composition around the knee has special importance. Although theorigin of the knee is a question in controversy, the escape probability from the Galaxy andfrom the acceleration region depends on the magnetic rigidity of a particle, p/Z, where p is
the momentum and Z is the charge of the particle. Accordingly the abundances of cosmic
ray nuclei as a function of energy reflect the properties of the cosmic ray propagation andtheir acceleration. Figure 2.3 [7] shows the cosmic ray composition below the knee measured
by the balloon–borne instruments. We can see the average mass of the cosmic rays increases
with increasing energy toward the knee. This supports the picture of the rigidity dependent
propagation and the acceleration.
For ultra high energy cosmic ray measurements, a depth of a shower maximum in theatmosphere(Xmax ) which is sensitive to mass of a primary particle is one of important observables to study the chemical composition. The studies of Xmax by Fly’s Eye suggested gradualchange of the chemical composition from the heavy components to light ones around 1018.5
eV. If this is true, it is a strong evidence for the cross over from galactic components to an
extra-galactic ones.

21

Chapter 3

Propagation and acceleration ofcosmic rays
Cosmic rays with energies from 1 to 105 GeV are considered to be accelerated at supernova
shocks, and to propagate diffusively in the turbulent magnetic fields in the Galaxy. Some of
the main aspects of diffusive propagations and shock accelerations are summarized in section
3.1 and 3.2. The descriptions in these sections are based on text books by Gaisser [30] and by
Longair [67], and on reviews by Cesarsky [20][21] and by Wefel [97]. Moreover, in the section
3.3 and 3.4, I will discuss some possible scenarios for accelerations and propagations of cosmic
rays with E > 1014 eV.

3.1

Propagation

Measurements of the chemical composition of cosmic rays below ∼ 1GeV/nucleon lead the factthat two groups of elements Li, Be, B and Sc, Ti, V, Cr, Mn are more abundant in the cosmic
radiation by many orders of magnitude than in solar system. These nuclei are almost absent as
the end products of stellar nucleo-synthesis. They are readily produced in cosmic radiation as
the spallation products of the abundant “primary” nuclei such as C, O and Fe. Measurementsof the “secondary” to “primary” ratios lead to following conclusions:
1. On average, cosmic rays in the GeV range traverse 5 ∼ 10g/cm2 equivalent of hydrogen
between the injection and the observation. This suggests that cosmic rays travel thousandsof times as much distances as the thickness of the galactic disk during their lifetimes,
because the amount of matter along a line of sight through the disk of the Galaxy is
about 10−3 g/cm2 .
2. The effective grammage decreases with the increasing energy, at least as far as observations
extend, as illustrated in Figure 3.1. This fact suggests that higher energy cosmic rays
spend less times in the Galaxy, and cosmic rays are accelerated before the propagations.
The propagation of cosmic rays is described by a transport equation written byZ∂vρ X dσi,k (E, E 0 )∂N= ∇ (Di ∇Ni )−[bi (E)Ni (E)]−∇uNi (E)+Qi (E, t)−pi Ni +Nk (E 0 )dE 0
∂t
∂E
m
dE 0
k≥i

(3.1)
22

Figure 3.1: Compiled data observed at the earth for (A) the B/C ratio and (B) the sub–
Fe/Fe ratio as a function of the energy(see Garcia–Munoz et al.,1984 [32] for references to theexperimental data.)
where Ni (E)dE is the density of particles of type i at position x with the energy between E
and E + dE. The first term on the right-hand side represents the diffusion, and the diffusion
coefficient can be interpreted as
1
D = λD v(3.2)
3
where v is a particle velocity and λD is a diffusion mean free path. With bi (E) ≡ dE/dt, thesecond term represents energy loss. The third term represents a convection with velocity u.
The source term is Qi (E, x, t) of the particle input of type i per cubic centimeter at position x
and time t with energies E to E + dE. The fifth term represents losses of nuclei of type i bycollisions and the decay, with
pi =

vρσi
1
vρ
1
+
=
+
m
γτi
λi
γτi

(3.3)

where γτi is a Lorentz dilated lifetime of the nucleus. Finally, the last term is a cascade term,
written here to include both feed–down from higher energy as in a nuclear cascade, and nuclear
fragmentation processes.

233.1.1

Leaky box model

The cosmic rays seem to propagate freely in a containment volume, with a constant probability
per unit time of escape, τesc ¿ c/h, here h is the half thickness of the galactic disk. In
this model, the diffusion term is replaced by −N/τesc . With the absence of collisions andother energy changing process and without convection, the solution for a delta function source,
Q(E, t) = N0 (E)δ(t), is
N (E, t) = N0 (E, t) exp(−t/τesc )(3.4)
Thus τesc is interpreted as the mean time spent by the cosmic rays in a confinement volume
and λesc ≡ ρβcτesc is the mean amount of the matter traversed by a particle of velocity βc. In
the case of neglecting energy gains and losses and convection, and in equilibrium, the function(3.1) is simplified to
¶
µ
βcρ
1
βcρ X
Ni (E)+Ni (E) += Qi (E) −σi,k Nk (E),(3.5)
τesc (E)
λi
γτi
m
k≥i

where σi,k is the spallation cross section.
A great deal of data on the chemical composition of low energy galactic cosmic rays wastreated within this framework, and the important result is that all the nuclei have the same
propagation history, i.e. by (3.5) with the single parameter λesc . Within this model, the energy
dependence of the secondary to primary ratios is attributed to the energy dependence of λesc .
Gupta and Webber [37] obtained λesc as(
10.8β
R < 4GV−2¡ 4 ¢δ
(3.6)
λesc (g cm ) = βcρτR =R > 4GV
10.8β R
with δ ' 0.6, where R is a rigidity in GV. This has an important implication for the source
spectrum, Qi (E). For a primary nucleus for which feed–down from fragmentation of heavier
nuclei can be neglected, the solution of (3.5) has the form
NP (E) =QP (E)τesc (R)
1 + λesc (R)/λP(3.7)

where the suffix P represents a primary nucleus.
For protons the interaction length λP is ∼ 55g/cm2 and λesc ¿ λP for all energies, so thatin this case the denominator of (3.7) is approximated with unity. Thus if the observed spectrumis
N (E) ∝ E −(γ+1) ,(3.8)
the source spectrum must be
Q(E) ∝ E −(γ+1−δ) .

(3.9)

For the observed spectral index γ + 1 ' 2.7, the spectral index at a source is given as ∼ 2.1.
This value is close to the predicted spectral index with first order Fermi acceleration at a strong
shock(See next the section).
The other major constraint on models of the propagation is based on the ratios of unstable
to stable isotopes of secondary nuclei. Unstable nuclei with lifetimes comparable to τesc , suchas 10 Be and 26 Al, can be used as “cosmic ray clocks”. For stable secondary nuclei, the solution

24

Figure 3.2: Comparison between experimental measurements of the radioactive isotopes 26 Al
and 10 Be and calculated prediction [38]. (see Wefel,1988 [97] for references to the experimental
data.)of (3.5) depends only on λesc and not on τesc and ρ. A measurement of the ratio of an unstable
to a stable isotope allows for as separate effects of the escape time and the density.
The most well studied example is Be. The isotope 10 Be is unstable with τS ∼ 3.9 × 106
years. Garcia–Munoz, Mason and Simpson [31] found τesc ∼ 2 × 107 years, with rather large
uncertainties. This result implies that cosmic rays propagate in a volume of mean density∼ 0.3 protons/cm3 (Figure 3.2), and therefore, this suggests that the containment volume inthe leaky box model is considerably larger than the disk of the Galaxy, and perhaps extends
into the galactic halo.

3.2

Acceleration of cosmic rays below 1014 eV

Supernova(SN) blast waves are considered as the major sites for accelerations of cosmic ray
particles. In this picture, particles are accelerated diffusively at outer supernova remnant(SNR)shocks which convert available hydrodynamic energy of a SN explosion into ultra-relativistic
particles with an overall efficiency of the order of 10%. This picture is generally accepted
because: (1) the efficiency is high(∼10%), (2) the distribution of cosmic rays in the Galaxy is25

well explained with the dispersion of shocks and the distributions of SNRs, and (3) this pictureleads the similarity between the observed and calculated energy spectra. In this mechanism the
maximum energy to which particles can be accelerated is the order of 1014 eV.

3.2.1

Power

The local energy density of cosmic rays is ρE ' 1eV/cm3 . If this values is typical throughout
the galactic disk, the power required to supply all the galactic cosmic rays is
LCR =

VD ρE
∼ 5 × 1040 ergs s−1
τR(3.10)

where the volume of the galactic disk, VD is
VD = πR2 d ∼ π(15kpc)2 (200pc) ∼ 4 × 1066 cm3(3.11)

and where τR is the residence time of cosmic rays in the galactic disk. The correct estimateof the residence time is τR ∼ 6 × 106 years by the leaky box model. Supernovae are the mostplausible sources to supply the required power for galactic cosmic rays, because the averaged
power released in the Galaxy with supernovae is
LSN ∼ 3 × 1042 ergs s−1 .

(3.12)

when we assumed a type II supernova occurs in every 30 years, and it ejects 10 M¯ with avelocity ∼ 5 × 108 cm/s. There are large uncertainties in these numbers, but it is plausible thatan efficiency of a few percent would be enough for SN blast waves to energize all the galactic
cosmic rays.

3.2.2

First order Fermi acceleration at supernova blast waves

The Fermi mechanism was first proposed by Fermi in 1949 as a stochastic means by which
particles colliding with clouds in the interstellar medium could be accelerated to high energy.
Here we consider the mechanism in rather simpler fashion, that is, the particle accelerations instrong shocks. If an energy gain per encounter is proportional to the primary energy, ∆E = ξE.
After n times encounters the energy En = E0 (1 + ξ)n , where E0 is an energy at injection into
the accelerator. Then, the number of encounters needed to reach energy E is,µ ¶E/ ln (1 + ξ) .
(3.13)
n = lnE0Thus, if the escape probability from the acceleration region per encounter is Pesc , the proportionof the particles accelerated to energies greater than E is
N (> E) ∝∞Xm(1 − Pesc )

m=n

with

µ
γ = ln1
1 − Pesc(1 − Pesc )n
1
=
=
Pesc
Pesc

¶
/ ln (1 + ξ) '

26

µ

E
E0

¶−γ

Pesc
1 Tcycle
= ×
ξ
ξ
Tesc

,

(3.14)(3.15)

where Tcycle and Tesc are the characteristic time for an acceleration cycle and that for an escapefrom the acceleration region, respectively. The ratio of these two times is Pesc . It is noticed
that the Fermi mechanism leads to the desired power law spectrum.
The diffusion of charged particles in the turbulent magnetic fields physically carried alongwith moving plasma is a mechanism for energy gains and losses. Here we consider a physicalsituation illustrated in Figure 3.3. In this situation, a large plane shock front moves with velocity−U1 , the shocked gas flows away from the shock with a velocity U2 relative to the shock front,
and |U2 | < |U1 |. Thus in the laboratory frame the gas behind the shock moves to left with

Figure 3.3: Acceleration at the plane shock front.
velocity V = −U1 + U2 . Here we consider a particle with energy E1 , it goes into the shockedgas where the particle begins to diffuse by scattering on irregularities of the magnetic fields. In
the rest frame of the moving gas the particle has a total energy,
E10 = γE1 (1 − β cos θ1 )(3.16)

where γ and β ≡ V /c are the Lorentz factor and the velocity of the cloud, respectively, and theprimes denote quantities measured in a frame moving with the shocked gas. All the scatterings
inside the cloud are the motions of the particles in the magnetic field and are therefore elastic.
Thus, the energy of the particle in the moving frame just before escapes is E20 = E10 . Thus, theenergy of the particle after its encounter is given as
E2 = γE20 (1 + β cos θ20 ).

(3.17)

Substituting (3.16) into (3.17), the energy change for the particular encounter is given by∆E
1 − β cos θ1 + β cos θ20 − β 2 cos θθ1 cos θ20=−1E1
1 − β2(3.18)

Averaging over cos θ20 and cos θ1 gives the average fractional energy gain per encounter, ξ. An“encounter” is one pair of back and forth across the shock. The distribution of cos θ20 isdn= 2 cos θ20 ,d cos θ20
27

0 ≤ cos θ20 ≤ 1(3.19)

and hcos θ20 i = 2/3. The distribution of cos θ1 for the plane shock is again the projection of anisotropic flux onto the plane, but this time with −1 ≤ cos θ20 ≤ 0, so that hcos θ20 i = −2/3. Thus
1 + 43 β + 49 β 2
44 U1 − U2ξ=−1∼ β =
1 − β233c(3.20)For the configuration of a large and plane shock, the rate of encounters is given by the
projection of an isotropic cosmic ray flux onto the plane shock front, that is
Z 1
Z 2π
cρCR
cρCR
d cos θdφ=
(3.21)
4π
400where ρCR is a number density of particles undergoing acceleration. The rate of convection
downstream away from the shock front is ρCR × U2 , thus the escape probability is given by
ρCR U2
4U2=
cρCR /4c(3.22)

Pesc3=
ξ
U1 /U2 − 1(3.23)

Pesc =
Thus the integral spectral index isγ=

It is noteworthy that the spectral index is independent of the absolute magnitude of velocitiesof plasma, and that it depends only on the ratio of an upstream velocity to downstream one.
The continuity of mass flow across the shock (ρ1 U1 = ρ2 U2 ) together with the kinetic theory ofgases gives(cp /cv + 1)M 2
U1
ρ2
=
=
.
(3.24)
U2
ρ1(cp /cv − 1)M 2 + 2
For a monoatomic gas the ratio of specific heats is cp /cv = 5/3, thus for a strong shock withM À 1,
4γ ' 1 + 2.
(3.25)
M
Not only the spectral index for the first order Fermi acceleration is universal, but it has anumerical value close to what is needed to describe the observed cosmic ray spectrum.
The finite lifetime of a supernova blast wave as a strong shock also limits the maximum
energy per particle that can be achieved with this mechanism. The acceleration rate is
ξE
dE
=
dt
Tcycle(3.26)

with the fractional energy gain per encounter, ξ. To integrate (3.26) and to estimate Emax , itis necessary to know the cycle time for one back and forth across the shock. This value can be
estimated with considering the number density distribution in the up and down stream regions.
The particle current with convection is given by
J = −D∇N + UN.

(3.27)In the upstream region, the fluid velocity U1 is negative relative to the shock front in equilibrium,
and there is no net current, anddN
κ1= −U1 N.
(3.28)dz
28

Then in the upstream region
N (z) = ρCR exp (−zU1 /κ1 ) ,(3.29)

where ρCR is the number density of cosmic rays at the shock. Therefore, the total number ofparticles per unit area in the upstream region is ρCR κ1 /U1 . The rate of encounter is cρCR /4,
thus the mean residence time of a particle in the upstream region is given by¶µ4κ1ρCR κ1 ³ cρCR ´−1=(3.30)U1
4
cU1
In the down stream region, this form is identical to that in the upstream region. Thus
µ
¶
4 κ1
κ2
Tcycle =+.
c U1 U2(3.31)The minimum diffusion coefficient gives a maximum possible acceleration rate and hence anupper limit to the energy accessible for an accelerator. The minimum diffusion coefficient isκmin =

rg c
1 Ec
∼
,
3
3 ZeB

(3.32)

where rg is the Larmor radius of a particle, rg = pc/ZeB, so that Tcycle ≥ 20E/(3U1 ZeB) fora strong shock with U2 = U1 /4. The resulting estimate of the maximum energy is
Emax ≤3 U1
ZeBU1 TSN
20 c

(3.33)

where TSN is the lifetime of an accelerator. The value of TSN is estimated as the time when asupernova has swept up its own mass, so that obtained from4
π(U1 TSN )3 ρISM = Mejecta .
3(3.34)For 10M¯ ejected at 5 × 108 cm/s into the normal ISM with 1 proton/cm3 , TSN ∼ 1000 years.
With an estimate of B ∼ 3 µG in ISM (3.33) gives
Emax ≤ Z × 3 × 1013 eV.

(3.35)

Consequently, the diffusive shock acceleration by supernova shocks appears to be the mechanism
for the acceleration of the bulk of cosmic rays, but seems to be limited to energies below∼ 1014 eV.

29

3.3

Acceleration of cosmic rays above 1014 eV

The simplest set of assumptions for the particle acceleration at supernova blast waves gives avery good explanation for origins of the bulk of cosmic rays. But it leaves open the origin of
those with energies greater than 1014 eV, although cosmic rays with energies up to 1020 eV are
observed with air shower experiments.
From (3.33) it is clear that the shock mechanism itself can achieve higher energy per particle
if the magnetic field is increased or if the time–scale for acceleration is increased. In this section,
we consider some of the possible scenarios for acceleration mechanisms of very high energy
cosmic rays.

3.3.1

Acceleration of cosmic rays in a supernova shock which travels downa stellar wind

Biermann [14] proposed a model that energetic cosmic ray particles are accelerated up to 3 ×
1018 eV in supernova shocks that travel down the density gradient of stellar winds. He assumed
a stellar wind which has the standard Parker spiral magnetic–field structure,
¶
µ 2
rs
rs2
2 1/2
(1 − µ ).
(3.36)(Br , Bφ ) = Bs,−
r2 rW rHere Bs is the surface radial magnetic field of a star, r is the distance from the star, rs is thesurface radius of the star, and rW = vW /Ωs where vW and Ωs are the wind velocity and theangular rotation rate of the star, respectively. Völk and Biermann [95] have argued that if theradial diffusion coefficient increases linearly with r, the adiabatic loss time and the acceleration
time in the first order Fermi theory have the same radial dependence leading to the preservationof the highest energy reached by particles until the shocks run into the stellar wind termination
shell. Biermann [14] estimated a maximum accelerated energy ofEmax =U2
ZerB2 = ZerB1 .
u1(3.37)

Here U1,2 and B1,2 are the velocity and the magnetic field strength on the both sides of theshock, and Ze is the particle charge. With the assumption that the product Br has the samevalue on the surface as that in the wind, and is of order 3 × 1014 cm Gauss for any OB stars andany Wolf–Reyet stars, a maximum energy of particles is inferred of
Emax (protons) = 9 × 1016 eV(3.38)

Emax (irons) = 3 × 1018 eV.

(3.39)

andIn this model, there are two critical postulates for radial diffusion based on the observational
evidence. The first critical assumption is that the convective random walk of energetic particles
perpendicular to the magnetic field can be described by a diffusive process with a downstream
diffusion coefficient κrr,2 which is given by the thickness of the shocked layer and the velocity
difference across the shock, and which is independent of energy:κrr,2 =

1 U2
r(U1 − U2 ).
3 U130(3.40)

The upstream diffusion coefficient can be derived in a similar way with the second assumptionthat the upstream length scale is just U1 /U2 times larger, i.e., r. Then the upstream diffusion
coefficient is1
κrr,1 = r(U1 − U2 ).
(3.41)
3
The energy gain due to the Lorentz transformations in one cycle of a particle remaining near
the shock and cycling back and forth from the upstream to the downstream can be written by
¶
µ
∆E
4 U1
U2
.
(3.42)
=
1−
E
3 c
U1
In addition, the energy gain associated with a drift is given by the product of the residencetime, the drift velocity and the electric field is obtained as
µ
¶
∆E
U2
U1 B
4r
1−
· Vd,θ · Ze=
(3.43)
E
3c
U1
c
µ
¶µ
¶
U2
U2
4 U1
fd 1 +
1−(3.44)=
3 c
U1
U1
where fd = 1/3(1 + U1 /2U2 ), fd = 1 for strong shocks, andµ¶1
U1 rg cVd,θ =
1+
.
3
2U2
r

(3.45)

The first term of Vd,θ is due to the gradient term and the second term is due to the curvature
which is increased by a lot of convective turbulences by a factor two. This additional energy
gain due to a drift changes the particle spectrum by
µ
¶
3U2
1
1−
(3.46)
U1 − U2
x
where

1
x=1+
3

µ
¶µ
¶
U1
U2
1+
1+
.
2U2
U1

(3.47)

Moreover, the acceleration time of a particle to reach a certain energy E is represented byµt(E) = t0
with
β=

E
E0

¶β

3U2 2 κrr,1.
U1 − U2 x rU1(3.48)(3.49)The particles injection rate is proportional to rb . This leads to a correction factor for theabundance of (E/E0 )−bβ . However, in a d–dimensional space, particles have rd more space
available to them than when they were injected. Thus this leads another correction factor
which is (E/E0 )−dβ . The combined effect of these two correction factors is a spectral change
by
κrr,1
3U1 2
−
(d + b).
(3.50)
U1 − U2 x
rU1
31

Hence the total spectral difference, as compared with the plane–parallel shock acceleration, isgiven with taking the minus sign by
µ
µ
¶
¶
κrr,1
3U1
U2 1
2
− 1 + (d + b).
(3.51)
U1 − U2 U1 x
x
rU1
For stellar winds we have b = −2 and d = +3, and so b + d = 1. Since the total spectral change
is given by 1/3 for strong shock(U1 /U2 = 4), the source spectrum is obtained as
Spectrum(source) ∝ E −7/3 .

(3.52)

After correction for the leakage from the galaxy, the spectrum at the earth isSpectrum(earth) ∝ E −8/3 ,(3.53), and it is very close to what is observed near the earth at particle energies below the knee.
In this model, the geometrical arguments lead to the knee energy and the spectrum of theenergetic particles beyond the knee. The critical assumption is about the diffusion coefficientin the lateral direction. The characteristic velocities of particles in θ are on average the drift
velocity Vd,θ , and the characteristic distance is the distance to the symmetry axis r sin θ. Thus
the diffusion coefficient is assumed to be1
κθθ,1 = |Vd,θ |r(1 − µ2 )1/2
3

(3.54)

where µ = cos θ. On the other hand, an upper limit to the diffusion coefficient in the lateral
directions is calculated from the residence time on any side of shock,µ¶4κrr,14κrr,24r
U2==
1−(3.55)
U1 c
U2 c
3c
U1
and the maximum velocity difference U1 − U2 , i.e.,4κθθ,max =
9

µ
¶ µ ¶
U2 3 U1 2
1−
rc.
U1
c

(3.56)

There is a critical energy when κθθ,1 reaches the maximum. It is another basic assumption inthis model that the latitude drift is dominated by the lateral convective motions induced by theradial convection above the critical energy. Thus the drift energy gain decreases by factor 2/3.
The critical energy is given byµ¶
3 U1 2
Ecrit =
Emax .
(3.57)
4 c
This critical energy corresponds to the knee energy Eknee . The reduced energy gain above theknee energy gives x = 11/6, and this leads to the overall spectrum at an injection,
Spectrum(source) = E −29/11(3.58)and the spectrum is deformed with the diffusive transport through the Galaxy,Spectrum(earth) = E −332

(3.59)

as observed. The bending points of the component spectra depend on their charge; with theincreasing energy the flux of protons drops off at the lowest energy, and that of heavy nuclei
does at high energy. Beyond the knee, irons and other heavy nuclei dominate in the high energy
galactic cosmic rays.
Stanev et al. [86] presented the all particle spectrum which is composed of following threecomponents:
1. The explosions of normal supernovae into an approximately homogeneous interstellar
medium drive blast waves which can accelerate protons to about 1014 eV. For these particles the spectral index is near –2.75 with taking account of the leakage from the Galaxy.
2. The explosions of Wolf–Rayet stars into their former stellar winds accelerate iron nuclei
up to energies of about 3 × 1018 eV. These particles have the slightly flatter spectra with
the index of −8/3 up to a rigidity dependent bending in the spectrum, and then beyond
the bending point the spectrum has rather steeper spectra with the index of –3 up to therigidity dependent cutoff energy.
3. The hot spots of radio galaxies produce particles with even higher energies, up to 1020 eV.
Their spectrum is approximately proportional to E −2 .
Using existing data on the chemical composition of cosmic rays near TeV energies as aconstraint, they showed the all particle spectrum and various contributions of the different
element groups, as shown in Figure 3.4,3.5. Thus, they conclude the followings: a) For the mostof energy range above 1014 eV the wind explosions can account for both the chemical composition
and the spectrum including the knee feature, b) the highest particle energies required from thestellar wind explosions imply a magnetic field in the preexisting stellar wind of at least 3 Gaussat a distance of 1014 cm, and c) the chemical abundance above 1014 eV is dominated by heavy
nuclei such as neon and higher ones.

33

Figure 3.4: The all particle spectrum and the various contributions from the different element groups derived by Stanev et al. [86] are compared with the observational results.
Here the symbols denote the following experiments: Akeno(open circles), Haverah park(opensquares),Yakutsk(open triangles, pointing down),Tien Shan(open triangles, pointing up), Fly’s
Eye(open hexagons),Proton 4(full squares), and JACEE(full circles). The contributions of thesix element groups are distinguished with the line type: H(solid), He(dots), CNO(dash), Ne–S(long dash), Cl–Mn(dash dot), and Fe(long dash dot).

34

Figure 3.5: A comparison of the measured element abundances with the model calculation for
H(full squares), He(full circles), and Fe(full stars).

35

3.3.2

Reacceleration of cosmic rays in a termination shock of the galactic
wind

The observed anomalous component in the cosmic ray flux in the inner solar system withenergies of the order of 1–10 MeV is explained with an acceleration model at the termination
shock of the solar wind. Jokippi and Morfill [53] suggested that cosmic rays up to 1014 eVenergy are accelerated by supernova shocks, but that particles with energies of 1015 eV andhigher are accelerated at the termination shock of the galactic wind in analogy with the solar
system. Although the existence of the galactic wind flowing outward from the Galaxy cannot
be regarded as established, they estimated the properties of the galactic wind with a simple
consideration. They took a supernova energy input rate into the interstellar medium ofĖSN = 1042 erg s−1(3.60)

and assumed that the energy input rate into the galactic winds is
2ĖW = Ṁ VW/2 = 5 × 1041 erg s−1 .

(3.61)

Moreover, they considered the available time for accelerations at the galactic termination shockis the age of our Galaxy Tgal = 1.5 × 1010 years, and assumed that the magnetic field strength
at the termination shock has the Parker spiral formµB = B0

R0
Rsh

!1/2
¶2 Ã
2 Ω2
Rsh
gal
1+
2
VW

(3.62)

where Rsh is the shock radius, B0 is the magnetic field at the reference radius R0 , and Ωgal =10−15 s−1 is the galactic rotation rate.
The mass loss rate is estimated as the diffusive galactic matter, 5 × 109 M¯ , divided by thelife time of the Galaxy, 1.5 × 1010 years, i.e., Ṁ = 2 × 1025 × F g s−1 , where F is a scaling
parameter. Then, this yields a wind velocity VW = 500F km s−1 .
They estimated the maximum rigidity of cosmic rays accelerated at the galactic wind using
the time scale for the particle acceleration which is given bytacc =

4κ
2
Vsh(3.63)

where κ is the diffusion coefficient and Vsh is the shock velocity. If the Bohm diffusion isassumed,1
κ = rg c(3.64)
3
where rg is a gyro radius. Thus the maximum rigidity is given byÃ!
Ω2gal R02 Ṁ F
450Tgal ĖW
2
Rmax (volts) =B0 η 1 +(3.65)c
Ṁ F
2ĖW η 2
where η = R0 /Rsh . In Table 3.3.2 are displayed Rmax for various η using R0 = 10kpc, B0 =3 × 106 Gauss, and F = 1.

36

η
0.2
0.1
0.05
0.02

Rsh (kpc)
50
100
200
500

Rmax (volts)
7 × 1018
3.2 × 1018
1.6 × 1018
6.4 × 1017

EF e (eV)
1.8 × 1020
0.8 × 1020
0.4 × 1020
0.2 × 1020

Table 3.1: The maximum rigidity and the maximum energy of iron nuclei of cosmic rays accelerated at the galactic wind termination shock.

They argued that particles accelerated at the terminus of the galactic wind will be subject
to modulation by the galactic wind in analogy with the solar modulation. The models of solar
problem suggest that the cosmic ray intensity behaves roughly as if the particles were decelerated
in a potential field equal to the electrostatic potential difference Epot between the heliosphericpole and the equator, which is
Epot = qV0 B0
qB0 Ṁ
R0 =cc 2πρ0 R0(3.66)

Using the standard values and setting ρ = 10−3 mp yields a value of 1015 eV. Here mp is themass of a proton. This suggests that cosmic ray particles with energies lower than this value
will be significantly modulated by the galactic wind.
They reported [54] a detailed analysis of numerical simulations of this model and represented
the expected energy spectra as shown in Figure 3.6.

37

Figure 3.6: (upper) Computed energy spectrum of accelerated iron nuclei at the equatorial
plane of the galactic wind, at radii of 10 kpc(lower curve) and at the shock(upper curve). In
this case the shock is at 300 kpc, and the outer absorbing boundary is at 600 kpc. The windvelocity is 400 km s−1 , and the diffusion coefficient is taken to be independent of energy, witha parallel diffusion mean free path of 14 times the local gyro radius of a 1020 eV particle. The
magnetic field is approximated by an Archimedian spiral with a magnitude of 3 × 10−6 Gauss ata radial distance of 10 kpc, and a rotational angular velocity of 2 × 10−15 s−1 . The particle drift
velocity is set equal zero. Injection is continuous, and uniform over the shock, at the energyof 2 × 1017 eV. The spectrum is given at ∼ 7 × 109 years after the start of acceleration, when
the distribution has reached a near steady state. (lower) The computed spectrum as in the left
figure, except that the mean free path depends on both the energy and position and is 25 times
the local gyro radius of the particles at any given energy.

38

3.3.3

Reacceleration of cosmic rays with multiple shock encounters

Axford [9] discussed about the merits and about the difficulties of different models, and drew
the conclusion that accelerations of cosmic ray particles beyond 1014 eV might be closely related
to the reacceleration of cosmic rays via multiple interactions with large scale structures in the
Galaxy, i.e., molecular clouds, SNRs, OB associations and rotating neutron stars, etc, andespecially with the associated shock waves. Ip and Axford [51] illustrated the the quantitative
aspects of this model taking into account the size frequency distribution of SNRs, the probability
distributions of energy gain/loss during the SNR encounters and the compositional abundancesof cosmic ray nuclei.
They considered three ways for cosmic rays to interact with SNR. First, in the case thatthe gyro–radius of a particle is considerably smaller than the size of a SNR and the particle
enters the SNR along a magnetic filed line, the energy will be reduced as a result of the adiabatic
cooling inside the expanding SNR. Second, if the initial pitch angle is large, particles are reflected
at encounters with a shock front, and the energies will be increased. Since these process may
be essentially considered to be of the nature of the first order Fermi acceleration, the absoluteenergy changes in these two types of encounters are on the order of 0.3% at the particle energy∼ 1016 eV. Finally, if particles encounter the shock with a guiding center impact parameterb within one gyro radius Rg , cosmic ray particles may be temporarily trapped and have thecycloidal motion across the shock front in the direction of the electric field E = −V × B, andgain relatively large energy (∆E/E up to 10%) at 1015 ∼ 1016 eV. This drift acceleration does
not operate when Rg reaches comparable to the radius of the SNR. They found the upper energy
limit Emax = 3 × 1017 eV for protons, and Emax = 3 × 1017 Z eV for heavy nuclei. In theircalculations they assumed the random encounters to be entirely uncorrelated, the life time of
SNR of the Sedov phase to be 2 × 105 years and the random explosion of supernovae in space
and time.

Figure 3.7: The distinction between the multiple encounter acceleration model described in thetext(b), and the classical second order Fermi acceleration(a).
The simulated cumulative spectra with an assumued galactic magnetic field of 6µG are shownin Figure 3.8. With their calculations, they concluded that the power law spectral index γ issensitive to the ratio of the dynamical escape time to the SNR encounter time, i.e., Nt = td /∆t,
and they found that a power law with γ = 3, and Nt is required to be 170 for B0 = 10µG, and39

300 for B0 = 3µG. Although the latter value is generally accepted to be the average value ofthe magnetic field strength in the galactic disk, the corresponding dynamical time estimatedtd = Nt ∆t ' 3 × 105 years for ∆t ' 103 years is too long. This is a matter of concern in theirmodel.

Figure 3.8: The integral energy spectra is derived by using different values of the injection
rigidity(P0 ) and Nt . The background magnetic field is assumed to be 6µG.
They also lead the resulting synthetic spectrum as shown in Figure 3.9 by assuming injection
rigidity P0 of 1014 V and Nt = 220. In their discussion, the discrepancy between the theoreticaland the observed spectra may be partly resolved by adjusting the relative abundances of theheavy nuclei at the injection energy of 1014 eV.

40

Figure 3.9: A synthetic spectrum of cosmic rays. The input parameters are P0 = 1014 eV and
B0 = 6µG.

413.3.4

Production of cosmic rays in active galactic nuclei

Protheroe and Szabo [74] proposed that a substantial fraction of cosmic rays with energies
between 1014 and 1019 eV may originate in active galactic nuclei(AGN). In this model, protons
are accelerated by the first order Fermi process at a shock in an accretion flow onto a supermassive black hole, and produce neutrons via the interactions with the photons at the central
region of AGN. These secondary neutrons escape from the central region and decay, and thus
the resulting protons are observed as cosmic rays with energies above 1014 eV.
In order to calculate the accelerations in an AGN central region, the authors assumed
three basic properties of AGN: a) Black hole masses to be proportional to the luminosity,
M ' x1 (LC /(1038 ergs s−1 ))M¯ , where x1 = R/rS is the ratio of the shock radius to the
Schwarzschild radius, and LC is the luminosity for the infrared to hard x–ray continuum, b) An
equi–partition between the energy density in the magnetic field and the radiation field at theshock, i.e., B 2 /8π = Urad ' LC /πR2 c to obtain the magnetic field and hence the gyro radiusrg at the shock, and c) The diffusion coefficient is larger by a factor b than the Bohm diffusion
coefficient, i.e.,
11(3.67)
D = b rg c ' 6.1 × 10−21 bx21 LC2 Ecm2 s−1
3
at the shock(r = R) where LC is measured in erg s−1 and E is measured in eV.
For the shock acceleration, the acceleration rate is dE/dt ' u21 E/20D where u1 is the
1/2
upstream flow velocity. Then the free–fall velocity onto the black hole u1 = x1 c, and thisleads an acceleration rate1dE
' 10−26 b−1 x1 LC2 Urad eV s−1(3.68)
dt
where the unit of Urad is eV cm−3 . Since proton energies reach the maximum when the acceleration rate equals to the total energy loss rate with the pion photo–production and the pair
production processes, they estimated the maximum energy Emax with a Monte Carlo simulation,
Emax ' E0 (x21 LC /b2 L0 )α , where E0 = 1.8 × 1016 eV, L0 = 2 × 1046 ergs s−1 , and α = 0.18(forx21 LC /b2 L0 < 1) or 0.52(for x21 LC /b2 L0 > 1).
As a result of the interactions during and after the accelerations, secondary particles which
include neutrons are produced. Although the neutrons themselves are subject to the pion-photoproduction interactions, the some fraction will escape from the central region. The authors foundthat for the energy range below ∼ 5 × 1016 eV neutrons traveling radially out from the shock
will escape from the intense radiation field of the central region and below 1015 eV all neutrons
will escape from the central region.
A relativistic neutron escaped from the intense radiation field of the central region will decay
on average after traveling a distance r0 ∼ 2.8 × 104 (E/eV)cm. The resulting proton will diffuse
in the magnetic field which is tied to the accreting plasma. The diffusion coefficient D in theaccreting plasma is assumed to be
5

−2
2 −1
2
D(E, r) ' 1.3 × 1061 bx−3
1 LC Er cm s

(3.69)

where E is in eV and r in cm, and protons will be trapped for a time tesc ' r2 /2D. Protons
traveling in the accreting plasma will be subject to the pp collisions and the collisions withphotons. The time scale for the pp collisions is given by tpp ' (nσpp c)−1 , where σpp ' 30mb is
the pp inelastic cross section, n is the number density of the nuclei in the accretion matter, andobtained by the authors as
1

1

3n(r) ' 1.3 × 108 Q(x1 )−1 x12 LC2 r− 2 cm−342(3.70)

where Q(x1 ) ' 1 − 0.1x10.31 is the conversion efficiency of the kinetic energy of in–falling matterto the relativistic particles. The time scale for pγ collisions, tpγ = E/(dE/dt)pγ , is estimatedwith Monte Carlo calculations. Thus, the probability of surviving for a proton is given by−1−1−1
Psurv ' t−1
esc /(tesc + tpp + tpγ ).

(3.71)

Finally, with an assumption that the 2–10 keV luminosity is proportional to LC , LX '0.005LC , they predicted the cosmic ray spectrum for the Einstein–de Sitter model with a function,
µ
¶
Z
Z
dICR
LX dLCR
1 c −1 Zmax g(z)
− 25
E
dz
=
(1 + z)dLX ρ0((1 + z)E, LX ), (3.72)
dE
4π H0
f (z)
f (z)
dE
0
where dLCR /dE is the differential luminosity of the escaping protons, ρ0 (LX )cm−3 (ergs s−1 )−1is the local x–ray luminosity function of AGN, and f and g describe the evolution of luminosity
and the number density in co-moving coordinate space, respectively. Using the AGN models to
describe the luminosity function and its evolution, and refining the estimation of the magnetic
field and photon spectrum in the central region and the treatment of the energy loss processes,
the authors presented [89] the spectra of produced neutrons, that of protons escaping from an
AGN(Figure 3.10), and the predicted cosmic ray spectrum due to acceleration in AGN(Figure3.11).
With these results, they concluded that:
1. The predicted contribution from AGN to the cosmic ray spectrum is the same orderof magnitude as the observed intensity in the region of the knee and the higher energies.
Minor adjustment, for example, by increasing the ratio of LC to LX by ∼ 2 in the assumed
AGN continuum, could give better agreement with the observations at 1016 eV.
2. If this model is correct, any extra galactic component in the region of the knee will be 100%
protons, and then one would expect to observe an enhancement in the relative abundanceof protons in the cosmic rays at ∼ 1016 eV.

43

Figure 3.10: Spectra of produced neutron (full curves) and cosmic ray protons escaping from anAGN(dashed curve). Results are given for b = 10, x1 = 30, and LX = 1042 (left–most curves),
1045 and 1048 ergs s−1 . The reduction at high energies is due to interactions of neutrons withphotons during the escape from the central region, while the reduction at low energies is due tointeractions of decay protons with protons in the accreting plasma.

44

Figure 3.11: The possible contribution of cosmic rays accelerated in AGN to the observed spectrum. Results are shown for b = 1(horizontal hatching), b = 10(thin oblique hatching) andb = 100(thick oblique hatching). The measured total cosmic ray intensity and the proton intensity are also shown in this figure. The dotted line represents a single steady source contribution
at 10 Mpc.

45

3.3.5

Acceleration by oblique shocks at SNRs

Kobayakawa, Sato and Samura [63] explained the knee behavior of the energy spectrum withthe acceleration mechanism caused by oblique shocks. They examined that oblique shocks canaccelerate particles more efficiently, because the reflections at the shock front are more importantand rapid accelerations occur comparing to parallel shocks. Thus the maximum energy Emax
is raised by a couple of orders compared to those in parallel shocks. Emax is mainly limitedby the finite lifetime of the supernova blast wave. The effective shock acceleration can occur
until the swept up mass reaches to the ejected mass ∼ 10M¯ . The lifetime of the shock tsh is
estimated to be some hundred years for ejecta expanding with the velocity of ∼ 108 cm/s into
the medium of the average density, 1 proton/cm3 . Balancing the acceleration time scale withtsh , the authors lead,
³
´ −1

2 α + r2 sin2 α
2
r
cos
1
1
Rsh (r − 1)sin α1x2
+ ¡(3.73)
Emax =
U1 eB1 Z cos2 α1 +
3
¢
2
rcx
x
2
2
2
2
cos α + r sin α11where α1 is the angle between the shock normal and the field directions, r is the compression
ratio of the shock, i.e., r = U1 /U2 , and x is the square root of the ratio of κk to κ⊥ , that isx2 = κk /κ⊥ . The definitions of Rsh , U1 , e, B and Z were mentioned before. Substituting thevarious constants for plausible values, Emax is given as,h¡¢i −1¶µ¶µ¶µ2 + r2 1 − η 2 
2
r
η
Rsh
U1
1−η
B1
x2
η2 +
Emax = 2.5 × 1016
+
eV
3
7
2
30µG
3pc
10 m/s x[η 2 + r2 (1 − η 2 )] 2 (3.74)
where η = cos α1 . As shown in Figure 3.12, Emax is strongly dependent on the field inclination η.
For example, Emax in quasi–perpendicular shocks is larger by two or three orders of magnitude
for each value of x than that in parallel shocks. The maximum energy for parallel shocks(η = 1),
Ecrit can be estimated in the case of strong shocks(r = 4) with x = 30 and with typical valuesof B1 , Rsh and U1 ,Ecrit = 1.25Z × 1014 eV,(3.75)and they defined Emax (η = ηmin ) ≡ Ecut .
In their model, the field inclination η is assumed to distribute uniformly, and then theprobability in the width dη can be written asf (η)dη =dη
1 − ηmin(3.76)

where ηmin = U1 /c ∼ 1/30. Moreover, they assumed the injection efficiency ²(η) depending on
η, and took simply,²(η) = η.
(3.77)
The first order Fermi acceleration mechanism gives a power–law energy spectrum, so that thedifferential spectrum dJ/dE ∝ E −q , with q = (r + 2)/(r − 1). The authors considered thedependence of the spectral index on the obliquity, according to a test particle simulation workby Naito and Takahara [73], assuming that q is proportional to η such as,q(η) = aη + b
46(3.78)

Figure 3.12: The maximum energy Emax for a proton with three values of x versus magnetic
field inclination η.
where a = 0.68, b = 1.41 and ηmin = 1/30. This dependence of power indices on the obliquity
enhances intensities of accelerated particles by the oblique shocks. For Ecrit < E < Ecut , taking
into account the injection efficiency and the chance probability of η, a correction factor is givenas,µ¶−q(η0 )+q(ηa )Z ηE2η − ηmin0ηdη 0(3.79)
g(η) =
2
1 − ηmin 1 − ηmin
E
crit
ηmin
Finally, the differential energy spectra can be expressed as,(E < Ecrit ) CE −γdJ
CE −γ g(η) (Ecrit ≤ E ≤ Ecut )=dE 0(Ecut < E)

(3.80)

here C and γ are fixed parameters for each chemical component and are determined with source
conditions. The authors assumed two different sets of parameters, C and γ. According to theexperimental results at E = 1TeV they predicted the energy spectra of primary cosmic raycomponents as shown in Figure 3.13 ,3.14 and the energy dependence of chemical composition
as shown in Figure 3.15, and the parameters with subscript “HEGRA” have been taken fromthe analysis by the HEGRA air shower experiment, and the values with suffix “Ours” have been
determined by the authors. They fixed the value x = 30 in their numerical calculations. The
curves of the total flux in Figure 3.14 have a smooth bending around 1016 eV. The curve labeled“HEGRA” is consistent with DICE and CASA–MIA data, while the curve “Ours” reproduces
well Tibet and Akeno data up to several times 1017 eV. In Figure 3.15, they showed the average

47

value of the logarithm of mass number A as an indicator of the composition,P
f ln AiPi
hln Ai ≡,fi(3.81)

and both curves of the model predictions monotonically increase with the energy.

Figure 3.13: The predicted fluxes of total, proton, He, and other nuclear groups in the case of“Ours” parameter choice.

Figure 3.14: Comparison between the all–particle spectra reported by various groups and thepredicted spectra with the oblique shock model.

3.3.6

Summary of acceleration models

Here I summarize the acceleration models for cosmic rays with energies > 1014 eV. The characteristics of the model predictions are listed in Table 3.2. In the galactic wind termination48

Figure 3.15: Energy dependence of chemical composition of primary cosmic rays in term of< ln A >.
shock model, the chemical abundance depends on the properties of the diffusive motions ofcosmic rays in the galactic disk and the galactic halo. For example, when we assume an energy
independent diffusion coefficient and a leakage from the galactic disk, the chemical composition
is independent of primary energies.
The predicted composition by the oblique shock model depends on the boundary condition,that is, the abundances and the spectral indices for the components around 1 TeV. However, theysuggested that the predicted averaged mass number is greater than 2 when we take boundary
conditions being consistent with the measured component spectra. For example, in this model,
the proton dominant composition determined by KASCADE is inconsistent with the observed
component spectra in TeV energy region.

model/mechanism

Emax (proton)
9 × 1016 eV

Dominant component
at the kneeincreasing

Energy dependence ofhln AiNe–Si

SNR shock+stellar windtermination shockof the galactic windmultiple shock encounter
AGN origin
SNR oblique shock∼ 1020 eV
3 × 1017 eV∼ 1018 eV∼ 1018 eV–†∼ constant∼ constantincreasing–†
protonprotonmiddle∗Table 3.2: Summary of the characteristic predictions of acceleration models for cosmic rays
> 1014 eV. †: these values depend on the assumption of the diffusive motions of cosmic rays. ∗:it depends on the boundary condition, that is, the abundances and the spectral indices for thecomponents around 1 TeV.
49

3.4

Advective diffusion propagation model for high energy cosmic rays above 1012 eV

Here I consider a cosmic ray propagation model based on the simple galactic magnetic field modeland on detailed Monte Carlo simulations of charged particles in the Galaxy with assumptionsof the galactic winds and the magnetic field structures.
In some models, it is explained that the knee of the energy spectrum of cosmic rays isdue to the energy dependence of the acceleration processes in a source or due to differences inacceleration sites, as described in the previous section. On the other hand, some models attempt
to explain the cosmic ray spectrum around the knee energy with the propagation process of
cosmic ray particles in the Galaxy [75]. Many models of acceleration/propagation lead a criticalenergy depending on the charge of particles, so that these models predict that the mean mass
number increases with the increasing energy around the knee. Moreover, the evidences to specify
models and to be informative for constructions of a propagation theory, such as the primary
to secondary ratio around the knee region, have not been measured because of experimental
difficulties. Therefore, the detailed acceleration and propagation models are not establisheduntil now.
A propagation model must explain important facts, i.e., the anisotropy of the arrival directions of cosmic rays and its energy dependence. It is known that the observed first harmonic
anisotropy amplitude increases with the increasing energy above 1015 eV [65]. The plot reported
by Hillas [45] shows that the amplitude has the same energy dependence as E 2.47 ×dN/dE(Figure3.16).

Figure 3.16: Left ordinate: the anisotropy amplitude; right ordinate: E 2.47 × dN/dE [45].
The markers and the solid line represent the observed anisotropy and the energy spectrum,
respectively.
In this section, I will briefly describe about previous researches for the diffusive motionsof charged particles in the turbulent magnetic fields by Senda, Ogio and Kakimoto [79], andpresent a new model of propagations of cosmic ray particles energies more than 1012 eV. Thismodel predicts the observed features such as the energy spectrum, the anisotropy, the mean
mass number and their energy dependence.
50

3.4.1

Diffusive motion of charged particles in turbulent magnetic fields

The magnetic field lines in the Galaxy are nearly parallel to the spiral arms and the strength(|B|) is estimated to be 1 − 3µG with the rotation measures for pulsars and galactic radiosources [22][39][50][85]. In addition to this regular component, the irregularities of the magnetic
field exist, and have roughly the same strength as the regular component. The maximum scale
length of the irregularities(Lirr ) is estimated to be 10 − 100pc with a large uncertainty.
We had derived the Fokker–Planck coefficients of the diffusive motion of cosmic ray particles
along and perpendicular to the regular magnetic field lines with the test-particle simulations[79]. Here I briefly explain this calculation, the details are reported in reference [80]. Our
calculations are based on the results by Giacalone and Jokipii [33].
The equation of the charged particle motion in a magnetic field is
dp= ev × B(3.82)
dt
where the galactic magnetic field B is expressed by the sum of a regular component BL and an
irregular component δB,
B = BL + δB.
(3.83)The equation (3.82) is expected as the following differential equations at the position x,
dx
dt
dv
dt

= v
=

e
(v × B)γm(3.84)(3.85)

These equations are simplified with defining a vector y = (x, v) as followsdy= F(t, y).
dt(3.86)The equation (3.86) is numerically solved with the fourth order Runge–Kutta method. A
position and a velocity of a test particle at a time t + ∆t is given as,
k(1) = ∆ · F(t, y)(3.87)∆tk(2) = ∆ · F(t +, y + k(1) )(3.88)
2∆tk(3) = ∆ · F(t +, y + k(2) )(3.89)
2
k(4) = ∆ · F(t + ∆t, y + k(3) )(3.90)´
1 ³ (1)y(t + ∆t) = y(t) +k + 2k(2) + 2k(3) + k(4)(3.91)
6
The time step is taken as
1 2πrg(3.92)∆t =50 chere rg is the Larmor radius of a particle and c is the light speed.
The irregular component of the magnetic field is express as a sum of the transverse waves,each of which has a wave number kn as follows
δB =

NX
max

£
¤
£ ¡
¢¤
A(kn ) cos αk ex0 ± i sin αk ey0 exp i kn z 0 + β

n=1

51

(3.93)

here A(kn ) is the amplitude, αk is a random angle, and β is a random phase. ex0 , ey0 and ez 0 are
the unit vectors of the co-ordinate axis of the space of r0 = (x0 , y 0 , z 0 ). The expression (3.93))
satisfies divδB = 0. The space of r0 , in which the wave number vectors point in the directionof (0,0,1), is related with the real space of r, as follows
r0 = Rz (φ)Ry (θ)r
cos φ − sin φ 0
Rz (φ) =  sin φ cos φ 0 
0
0
1


cos θ 0 sin θ010 
Ry (θ) = − sin θ 0 cos θ(3.94)(3.95)(3.96)

Thus, in our calculations, the five values (θk , φk , αk , βk and the sign in (3.93)) are randomly
defined for each space.
Many observations of the rotation measure suggest that the regular and the irregular components of the galactic magnetic field have the same order in the magnitude of their energy
density, and thus we assumed
N

max
B2L
δB2
1 X
=
=
A2 (kn )8π8π8π(3.97)

n=1

moreover we assumed that the energy density spectrum P (kn ) has the Kolmogorov spectrum(Figure 3.17)4πkn2 ∆kn.
(3.98)
P (kn ) ∝1 + (kn Lirr )11/3
The interval of the wave number kn is calculated with a relation ∆kn /kn = 0.05 in our
simulations. In Figure 3.18 it is shown that the examples of magnetic field lines are shownfor Lirr = 10pc and 50pc, and Figure 3.19 is the track of a test particle proton with energiesof 5 × 1014 eV during 5500 years in the assumed galactic magnetic field of BL = 3µG, andLirr = 10pc.
We calculated the diffusive motions of 1000 test particles with energies of 5 − 9 × 1014 eV,assuming a single source of particles and Lirr = 10 pc for two cases of the magnetic field strengthof 1.0 and 1.5 µG(Figure 3.20). By the calculated position ∆x and the time after a explosion∆t , we obtained the Fokker–Planck coefficients,Dij=< ∆xi ∆xj > =< ∆xi ∆xj >
2∆t
N
1 X
∆xi ∆xjN(3.99)(3.100)k=1

The obtained time dependences of the Fokker–Planck coefficients are shown in Figure 3.21,
3.22. As shown in these plots the Fokker-Planck coefficients become constant during several
thousands years after a SN explosion. Thus, the motion of high energy particles in the galactic
disk is diffusive, and the Fokker–Planck coefficients are considered to be the diffusion coefficientsfor the diffusive motion in the galactic magnetic field.
52

Figure 3.17: The assumed power spectrum of the magnetic field.
We plotted the diffusion coefficients versus the normalized Larmor radii rg /Lirr as shown in
Figure 3.23. These linear relations on logarithmic graph is expressed by [48]
µ
¶ µ
¶
rg α
Dii
Dii
=
(3.101)cLirr
cLirr 0 Lirr
The fitting parameters of equation (3.101) for the calculated diffusion coefficients are listed inTable 3.3.
parameterαlog(Dii /cLii )0

parallel
0.6018 ± 0.0218
0.1893 ± 0.0281

perpendicular
0.6627 ± 0.0142−1.481 ± 0.018

Table 3.3: The fitting parameters of equation (3.101) for the calculated diffusion coefficients.

53

Figure 3.18: Examples of the simulated field lines, (a) for Lirr = 10pc and (b) for 50pc.

Figure 3.19: An example of a simulated track of a teat particle (proton) with energies of
5 × 1014 eV for 5500 years.

54

Figure 3.20: Snap shots of particle distributions after an explosion of a source(0,0) of particles
with energies 5×1014 eV. The magnetic field strength is assumed to be BL = 1µG and δB/B = 1.
The direction of the regular magnetic field is along with y–axis.

55

Figure 3.21: Time dependences of the Fokker–Planck coefficients for the diffusive motions ofcharged particles with energies 5 − 9 × 1014 eV in the assumed galactic magnetic field. Dxx
and Dyy represent the coefficients for parallel and perpendicular to the regular magnetic field,
respectively. We assumed BL = 1.0µG, δB/B = 1 and Lirr = 10pc.

56

Figure 3.22: Time dependences of the Fokker–Planck coefficients for the diffusive motions ofcharged particles with energies 5 − 9 × 1014 eV in the assumed galactic magnetic field. Dxx
and Dyy represent the coefficients for parallel and perpendicular to the regular magnetic field,
respectively. We assumed BL = 1.5µG, δB/B = 1 and Lirr = 10pc.

57

Dxx/cLirr

1

-1

10
10

-2

10

-1

10

Dyy/cLirr1

rg/Lirr

-1

-2

10

-3

10
10

-2

10

-1

1

rg/Lirr

Figure 3.23: The obtained diffusion coefficients depend on the Larmor radii. The diffusion
coefficients and the Larmor radii are normalized with cLirr and Lirr , respectively and c is thelight speed.

58

3.4.2

Advective–diffusion of cosmic ray particles

When we assume that the regular component of the galactic magnetic field is completely parallelto the spiral arms of the Galaxy, the leakage of cosmic ray particles is dominated by D⊥ . For
example, from equation (3.101), D⊥ for protons with energies 1014 eV is 0.005 pc2 /year for
BL = 3µG and Lirr = 100pc. With this value, we can estimate the residence time (τR ) in thegalactic disk with the half thickness l ∼ 150pc,
τR =l2∼ 5 × 106 years.
D⊥(3.102)On the other hand, an extrapolation of the leaky box model to 1012 eV leads a residence
time of 2 × 105 years. In a different way, a residence time is estimated from the anisotropy. The
observed anisotropy amplitude around 1012 eV is 10−3 . This amplitude is the order of the ratioof the diffusion velocity vD to the particle velocity v, so that vD ' 300 km/s. Then the averaged
residence time is estimated as l/vD ' 5 × 105 years. Thus, the estimated residence time withthe leaky box model and that from the observed anisotropy amplitude are small comparing withthe expectation (3.102), and therefore, an idea of the diffusive leakage of cosmic ray particles
perpendicular to the galactic magnetic field apparently contradicts to the observations.
With many observations for the magnetic field of the spiral galaxies, the loops and filament
structures are found [69], which are considered to be due to expansions of the magnetic field via
supernova explosions or the magnetic buoyancy [90](Figure 3.24). Thus in those structures wecan expect the existence of the magnetic field perpendicular to the galactic disk with open fieldlines, in which cosmic ray particles are possible to escape more rapidly with Dk . Moreover, thegalactic wind which has analogy with solar wind supposed to be exist. The galactic wind is theoutflow of the interstellar medium and that of the magnetic field irregularities. The existence ofthe galactic wind is not established, but the theoretically predicted velocity of the galactic winddepends on the assumptions of the boundary conditions between the Galaxy and the galactic
halo. In a model of the galactic disk surrounded with the dark matter halo, vg is predicted tobe ∼ 300 km/s [96].
Cosmic ray particles flow out along with the galactic wind because the Larmor radii of cosmic
rays with energies lower than 1016 eV are much smaller than the scale length of the irregularities.
Since, the cosmic ray particles leak out through the loops and the filament structures of thegalactic magnetic field, this process is considered to be not only diffusive motions with Dk , butalso the outflows with the galactic wind of velocity vg . Thus the cosmic ray leakage is considered
to be an advective–diffusion process with Dk and vg .
Here we assume that the direction of the leakages of cosmic rays is perpendicular to thegalactic disk surface, and that there is leakage along neither with spiral arms nor with radial
direction. Moreover, when we assume that the strength of magnetic field and the cosmic raydensity n are uniform in the galactic disk, the leakage of cosmic rays is expressed with thefollowing one dimensional advective–diffusion equation,
µ
¶
∂n
∂n
∂
∂n
+ vg
=
Dk
(3.103)∂t∂x∂x∂xwhere vg is a galactic wind velocity and Dk is a diffusion coefficient calculated with (3.101) fora particular charge and energy. From the advective–diffusion equation(3.103), I estimated theresidence time τR which is defined as the time, during which the number of cosmic ray particles
in a confinement volume reaches 1/e of a initial value for cosmic ray particles with energies
1012 − 1017 eV and for various pairs of the parameters Lirr and vg .
59

Figure 3.24: (upper) Illustration of the magnetic loop formation due to the instability and of
relativistic gas expulsion from the galactic disk [67]. (lower) A sketch of dark filaments foundin an optical image of spiral galaxy NGC253 [90].

3.4.3

Residence time of cosmic rays in the galactic disk

The equation (3.103) is calculated numerically with the Crank–Nicholson scheme. The equationis translated to the finite difference equations with an operator splitting method. In thesecalculations we took a one-dimensional confinement volume (with the length of 300 pc) where a
cosmic ray source Q = Q0 δ(t) is located at the middle point. The step sizes for the calculations
are 12.5 pc for distance, and 100 years for time, respectively. Here we assumed that the strengthof the galactic magnetic field is 3µG for the regular and the irregular components, and decreases
exponentially with a scale length of 1 kpc. The cosmic ray particles are assumed to consist of
five nucleus groups as listed in Table 3.4.
IDcomponents
mass number
charge1
proton
1
1

2
He
4
2

3
CNO
14
7

4
Ne–Si
24
12

5
Fe
56
26

Table 3.4: Assumed five components of cosmic ray particles.
In these calculations, all cosmic ray particles escape from the galactic disk with advective–
diffusive motions with Dk and with vg , i.e., all the galactic magnetic field lines are assumed
to be perpendicular to the galactic plane. This is a limiting case of the galactic magnetic field

60

configuration. It is natural to consider that a proportion of magnetic field lines are perpendicular
to the disk and the rest are along with the spiral arms. This proportion affects the estimatedresidence time of cosmic rays to extend with the same factor. However, the absolute value of τR
shown in Figure 3.25 is not much different from the residence time estimated from the observed
results, so that the proportion is not much smaller than unity.
The calculated residence times of cosmic ray protons for different advection velocities areshown in Figure 3.25. Also, the residence times for different Lirr and for different componentsof cosmic rays are shown in Figure 3.26 and 3.27, respectively.
Apparently, the energy dependence of τR is power law (∝ E −γ ) and the curves bend at∼ 1014 eV, that is, the flat curves (index: −γ1 ) are bending to the steep ones (index: −γ2 )with the increasing energy. Moreover, while γ1 clearly depends on vg , γ2 is approximately
independent of vg . This result shows that the leakages of lower energy cosmic rays are dominatedby the advection because the diffusion coefficient is too small to be efficient. This conclusion is
supported by the result shown in Figure 3.27 in which τR for low energy cosmic rays and γ1 is
not strongly depend on the charge of the particles.
In contrast, it is concluded that, for high energy cosmic rays, the leakages of cosmic rays are
dominated by the diffusion, and τR has the same power law energy dependence as the diffusion
coefficient in equation (3.101). Moreover, the energies of bending points on the curves of τR in
Figure 3.27 obviously depend on the charge of particles, and τR are smaller for lighter nuclei,
and thus for high energy cosmic rays, the averaged mass number of cosmic ray particles is
expected to increase with the increasing energy.
The power law indexes of τR are about −0.1 at 1012 eV and −0.6 at 1017 eV. The steepening
is gradual and the bending energy is E ∼ 1014 − 1015 eV. Since the energy spectrum N (E) =Q(E)τR (E), if we assume the power law index of a source spectrum Q(E) is constant of −2.6,
the indexes of N (E) are expected to be about −2.7 at 1012 eV and −3.2 at 1017 eV, and thesteepening of the spectrum at ∼ 1015 eV. These expected features are remarkably consistent
with those of the observed energy spectrum of cosmic rays.

61

τR
10

6

10

5

10

4

12

13

14

15

16

17

log10E

τR

Figure 3.25: The residence time τR predicted with the advective diffusion model for cosmic ray
protons. This figure shows τR s for different advection velocities, i.e., the galactic wind velocity,
vg . For all the calculations in this figure Lirr is fixed to be 100 pc.

10

6

10

5

10

4

12

13

14

15

16

17

log10E

Figure 3.26: The residence time τR predicted with the advective diffusion model for cosmic ray
protons. This figure shows τR s for different maximum scale lengths of the irregularity of thegalactic magnetic field, lirr . For all the calculations in this figure vg is fixed to be 5 × 10−4
pc/years = 490km/s.

62

τR
10

6

10

5

10

4

12

13

14

15

16

17

log10E

Figure 3.27: The residence time τR predicted with the advective diffusion model for different
primary components of cosmic ray particles. For all the calculations in this figure vg is fixed tobe 5−4 pc/years = 490km/s and Lirr is fixed to be 50 pc.

63Chapter 4

Extensive Air Showers
The cosmic ray spectrum is well described by a inverse power law in the energy, ∝ E −γ . Thus
for observations of high energy cosmic rays one requires long observation times and a large
effective detection area. The rate of arrivals of cosmic rays with energies above 1012 eV is very
few, thus one cannot observe directly with a detector provided on a balloon or mounted ona satellite. Accordingly high energy cosmic rays are observed through detections of extensive
air showers(EAS) developed in the atmosphere. This chapter is a description of air shower
phenomena and of their properties. Detailed discussions are found in a book by Rao andSreekantan [76].

4.1

Air showers

An extensive air shower is a swarm of particles produced by iterative interactions between
the atmosphere and primary and secondary cosmic ray particles. In 1930’s, Auger and hiscolleagues found many events of coincident signals on separated detectors. They concluded
that the coincidences are due to showers of secondary particles induced by high energy primarycosmic rays at a top of the atmosphere.
The inelastic interaction cross section of cosmic ray protons with energies of a few hundred
GeV on an air nucleus is about ∼ 250 mb, which corresponds to the interaction length ∼ 90
g/cm2 , where the mean mass number of the air is assumed ∼ 14. While the total thickness ofthe earth’s atmosphere is ∼ 1000 g/cm2 , so that, a high energy primary cosmic ray interacts
with the air and induces a cascade shower of secondary particles.
At the first interaction of a cosmic ray nucleus with the air, many pions and kaons areproduced. These secondary particles interact further with the air, whereas some of them decay
because of their finite lifetimes. A charged pion gives rise to a muon component in an air showerascosmic ray + air nucleus → π ± + π 0 + K ± + · · ·
π ± → µ± + ν(ν̄).

64

A kaon also contributes to a muon component via various decay modes to pions, such as,K±→ µ± + ν(ν̄)

K±→ π± + π0K±→ π± + π+ + π−··· .

A neutral pion decays into two gamma rays practically before any other interactions because of
their short lifetime, ∼ 10−16 s, in spite of the relativistic time dilation. A gamma ray initiates
the iterations of the electron–positron pair production and the bremsstrahlung, which induces
electromagnetic cascades,
π0

→ γ+γ

γ

→ e+ + e−e±→ e± + γ··· .

primary
particle

Most of air showers are mixture of the nuclear and the electromagnetic cascades. Figure 4.1
shows a schematic view of an extensive air shower.

π

γ



π

γ





nuclearcascadeπelectromagnetic
cascade

γ

νµ
µ

µ


νµ

Figure 4.1: A schematic view of an extensive air shower.
As a shower travels down in the atmosphere, a number of particles, which is called a shower
size, increases and reaches the maximum development. For example, a shower size at the
maximum development amounts to ∼ 105 for a shower of primary energy E = 1014 eV. Afterthe maximum development, a number of particles decreases with absorptions of particles byatmospheric atoms and molecules. These processes are called the longitudinal developments ofextensive air showers.
65

4.1.1

Interaction of electrons and gamma rays

An air shower is predominantly composed of gamma-rays, electrons and positrons. Various
approximate formulae to describe the electromagnetic cascade were given by many authors andhere I present a brief summary. This result is also helpful to discuss the nuclear cascade.
A relativistic electron (and positron) loses its energy mainly by the bremsstrahlung. Its
cross section for a relativistic electron is given by the Bethe–Heitler formulae, which is written
as,
·µ
¶ ³
¸
´ 1
4Z 2 re2 dv
22− 13
σbrems. (E, v)dv =1 + (1 − v) − (1 − v) ln 184Z
+ (1 − v) , (4.1)
137 v
3
9
where v = hν/E, Z is the atomic number of a target nucleus, re is the classical electron radius.
Here we define the radiation length X0 of electrons in a medium as,´
1
4Z 2 re2 N ³
1=
ln 184Z − 3 ,
(4.2)
X0
137 A
where N is the Avogadro’s number and A is the mass number of a target nucleus. Using X0
and integrating (4.1) by v = 0 to 1, the energy loss rate due to bremsstrahlung is
¶
µ
E
dE
'− .
(4.3)
dX brems.
X0
The value of X0 for the air is ∼ 38 g/cm2 .
Gamma rays with energy greater than 1 MeV lose their energies mainly via electron–position
pair productions. The cross section for pair production, σpp (hν, u) at energy hν is given by
·µ
¶ ³
¸
´ 1
4Z 2 re2222− 31
σpp (hν, u)du =du u + (1 − u) − (1 − u) ln 184Z
− (1 − u) ,
(4.4)
137
3
9
where u = E/hν and E is the energy of a produced electron. The interaction length of gamma
rays for pair productions is obtained by integrating (4.4) by u = 0 to 1, and we can write theenergy loss rate for gamma rays as
¶
µ
7 E
dE
=−
.
(4.5)
dX pair
9 X0Thus the interaction length for pair production is comparable to that for bremsstrahlung.
In an extensive air shower, when the average energy of electrons is less than the critical energy
Ecrit , the ionization loss process becomes dominant and the shower size begins to decrease. The
ionization energy loss rate is given by
!
#
" Ã
µ
¶
dE
πN Zre2π 2 m2e
−a ,
(4.6)
= −2
me ln
3
dX ion.
A(1 − β 2 ) 2 I (Z)2
where I(Z) is the mean ionization potential: I(Z) ∼ 80 eV, and a constant a ∼ 3. The
critical energy Ecrit is the energy at which an electron energy loss due to ionization processes
comes to comparable to that due to bremsstrahlung. If we assume the minimum ionization loss−(dE/dX) = 2 MeV/gcm−2 , we obtained the critical energy of electrons in the air; Ecrit ∼ 80
MeV.
For muons in an EAS, the practical energy loss process is only ionization loss, therefore,
they do not suffer significant energy losses even after a maximum development. Thus, most ofmuons survive and reach the ground before they lose their energy entirely or decay.
664.1.2Longitudinal development

A total number of electrons in an electromagnetic cascade Ne at an atmospheric depth t isapproximately given by
· µ
¶¸
0.31
3
Ne ∼ √ exp t 1 − ln s ,
(4.7)
y
2
where t is in the unit of the radiation length of electrons in the atmosphere, X0 ' 38 g/cm2 .
The variable y is defined asµ¶E0
y = ln
,
(4.8)
Ecrit
where E0 is the primary energy. The age parameter s is defined ass=

3t.
t + 2y(4.9)

From this definition, s = 0 at the first interaction point, and s = 1 at the maximum developmentof an EAS. The number of particles at the maximum development is obtained by substituting
1 for s in (4.7),
0.31
E0
.
(4.10)
Nmax (E0 ) '
1/2 E
crit
[ln (E0 /Ecrit )]Fenyves et al. carried out Monte Carlo studies for air shower developments induced byprotons and iron nuclei [27]. They showed that the air shower developments induced by primary
nuclei are expressed by the same formula as electromagnetic cascades, (4.7), with replacing t
with t0 = t + β(tmax − t)tγ , where tmax is the depth at the maximum development point, andβ(= 0.2 ∼ 0.4) and γ(∼ 0.5) are constants.

4.1.3

Lateral distribution of shower particles

As shower particles gain transverse momentum by the multiple scattering, the swarm of shower
particles expands laterally. The root mean square of a scattering angle of electrons after traversing the amount of matter, x g/cm2 , is
µ ¶
Es
2
<θ> =
x,
(4.11)
E
where

r

4π
me c2 = 21 MeV.
α
For electrons with energy E = Ecrit at x = X0 , we have
Es =

rm =

Es
X0 ,
Ecrit

(4.12)(4.13)

and rm ∼ 9.5 g/cm2 for the air. This is called Molière length, which gives the characteristic length of lateral spread of electromagnetic cascades. A structure function of the lateral
distribution is represented with the variable r/rm . For an electromagnetic cascade of age s,
the structure function f (r/rm , s) is well expressed with the Nishimura–Kamata–Greisen(NKG)
function:µ ¶s−2 µ¶rr s−4.5f (r/rm , s) = C(s)1+.
(4.14)
rm
rm
67

The normalization factor C(s) is taken to satisfy the following relation,Z ∞2πxf (x, s)dx = 1(x ≡ r/rm ).

(4.15)

0

Thus,
· Z ∞
¸−1
s−4.5
s−1
C(s) = 2πx(1 + x)dx= [2πB (s, 4.5 − 2s)]−10=

Γ (4.5 − s)2πΓ (s) Γ (4.5 − 2s)(4.16)

where B is Beta function, and Γ is Gamma function.
Using f (r, s), a local particle density ρ(r, s) at a radial distance r from the shower axis isgiven as,Ne
(4.17)
ρ(r, s) = 2 f (r, s)
rm
where Ne is a shower size.
Fenyves et al. shows that for air showers induced by primary nuclei, the lateral distributionsof shower particles can be written as the formula (4.14) with replacements of rm → rm /2 andt → t0 described in the previous section.
A shower disk has a finite thickness and a curvature, not a flat structure. Agnetta et al.
reported the time structure of the air shower front, and they showed the delay time τ from aflat shower plane [2],³ r ´2.1
ns,(4.18)
τ = 7.7 + 9.5
79 m
where r is the radial distance from a shower axis. Linsley gave an empirical formula of theshower disk thickness as a dispersion of arrival times of shower particles [66],στ = 2.6

4.2

³ r
´1.5
ns.
+1
30 m

(4.19)

Observations of extensive air showers

One technique for detections of air showers is the coincidence method with an air shower
array, which is a set of ground–based detectors(Figure 4.2). The total number of particles
in an air shower, so called the shower size, is determined from the local particle densities
sampled by detectors which are placed on the ground. The arrival direction of an air shower
is determined from relative arrival times of shower particles incident on detectors. This fast–
timing technique can be used because shower disks are thin enough comparing with their lateral
spread as mentioned above.

4.2.1

Determination of arrival directions

Here we assume that the surface of shower disks is a plane. The equation of a shower front isexpressed bylx + my + nz − ct = 0(4.20)

68

xis
ra
sh
ow
e

sh

ow
e

rd

detector 1

detector 2

detector 3

isk

detector 4

shower coreshielded muon
detector

Figure 4.2: A cross–sectional view of a front plane of an air shower and detectors.
where c is the speed of light, t is an arrival time of a shower front, (x, y, z) is a position of adetector, (l, m, n) is a set of direction cosines of a shower axis with a zenith angle θ and with aazimuth angle φ,l = sin θ cos φ(4.21)

m = sin θ sin φ(4.22)

n = cos θ(4.23)The vector (l, m, n) is determined by minimizing χ2 which is defined asXχ2 =wi [lxi + myi + nzi − c (ti − t0 )]2(4.24)iwhere (xi , yi , zi ) is the position of the detector i, and ti is the arrival time of a shower front at
the detector i. The weight wi is determined from errors of timing measurements. The equation
to be solved is
∂χ2
∂χ2
∂χ2
=
=
=0
(4.25)∂l∂m∂nunder the condition l2 + m2 + n2 = 1. The plane approximation for a shower front is valid
within distances less than 30 m from a shower axis.

4.2.2

Determination of shower sizes

A shower size, Ne , is determined by fitting the observed particle densities to an assumed lateral
distribution function, e.g., Eq.(3.17) in the shower plane using the least square procedure. In
this case χ2 is defined as
´2
X ³
exp
wi ρobs,(4.26)
χ2 =i − ρiiwhere ρobs
and ρexp
are observed and expected local particle densities at the detector i; wi is theiiweight factor for detector i. One practical method to find the core position (x0 , y0 ) is obtained69

by calculating density–weighted means of detector coordinates, such as,
X
X
x0 =
ρobs
x
/
ρobs
i
i
i
i

y0 =

X

(4.27)iρobs
i yi /

i

X

ρobs
i ,

(4.28)iwhere (xi , yi ) is the position of the detector i. The grid–search method for searching theparameters (x0 , y0 ) to minimize χ2 is also deployed. The actual method which we used is
described in Chapter 7.
The shower size is a good estimator for the primary energy. An approximated relation
between shower sizes Ne and primary energies E0 is
E0 ∼ 2 × 109 Ne eV.

(4.29)

The actual energy–size relation depends on the mass composition of primary cosmic rays
and on the hadronic interactions. The detailed measurement of primary composition and theassumptions for the hadronic interactions in the air are very important for the determinationof the averaged relation between primary energies and shower sizes.

4.2.3

Muons and other components in air showers

Muons
Electrons are the most dominant component in air showers. However, an air shower contains
also muons as decay products of charged pions and kaons. Since muons do not strongly interact,they loose their energies due to only ionizaion or decay. A total number of muons in an air
shower strongly reflects the nature of particle interactions and is sensitive to primary species.
The number of muons, Nµ , with energies greater than Eµ is expressed as [76],Nµ (> Eµ ) ∝ Neα ,(4.30)

where Ne is the electron size; the exponent α is ∼ 0.9 for Eµ ∼ 1 GeV, and α ∼ 0.7 for Eµ ∼ 200
GeV. The lateral distribution of muons in an air shower is expressed as,³r ´−2.5
ρµ (r) = CNµ r−0.75 1 +,R(4.31)

where C is a normalization constant.
Shielded detectors are used to detect muons in order to avoid background electrons. Somegroups discussed cosmic ray composition through measurements of air shower muons [25]. Furthermore, the muon component in air showers is important in high energy gamma ray astronomy.
Some groups measure muous in order to distinguish air showers induced by primary gamma raysfrom those induced by primary nuclei.
Čerenkov lightČerenkov photons are emitted by a charged particle which has the velocity greater than the
speed of light in a medium. In a medium with a refraction index n, particles with β(= v/c)

70

above 1/n only can emit Čerenkov light, thus n is setting a threshold energy. Čerenkov light isemitted the half angle θ around the direction of the particle. The angle θ is given asθ = cos−11,nβ(4.32)

where β = v/c. The reflectivity of the atmosphere is approximately expressed as a following
function, at depth t and temperature T [44],¶µµ¶273.2 Ktn = 1 + 0.000292.
(4.33)
1030 gcm−2
T
The minimum energy for electrons which emit Čerenkov photons in the atmosphere is ∼ 21
MeV at sea level, and the emission angle is about 1.4◦ .
Fluorescence light
Charged particles in an air shower ionize and excite air molecules, and then, the molecules inparticular excited states release energies as emissions of fluorescence light. Since the photon
yield is very small [57], the measurements of the fluorescence light are possible for a large air
shower induced by a very high energy primary. The Fly’s Eye experiment [10] is the pioneer ofsuch observations of extremely high energy cosmic rays.

4.3

Equi–intensity method analysis

The equi–intensity method is one of the techniques to measure the longitudinal development
of EAS. This method is based on the measurements of air shower sizes and rate of arrivals. In
the equi–intensity method it is assumed that extensive air showers induced by primary particles
with a same energy develop in the same way and that rate of arrivals of cosmic rays does notvary with time and direction. We can then plot the frequency distribution of the sizes of EASs
at the different slant depth, D1 · · · Dn . Since showers at large zenith angles have penetrated
large slant depths to reach the detector, the size that corresponds to a given intensity(showers
per m2 · sr · s) will be smaller for larger zenith angle. When it is assumed that showers which
are observed with the same frequency of occurrence at different zenith angles are initiated byparticles of the same energy, the developments of showers with depth through the atmospherecan be derived directly from the distributions observed at different zenith angles(Figure 4.3).
This procedure provides an average picture of the development of EASs for a particular primaryenergy.
This method has a great advantage that the obtained equi-intensity curves are less independent of Monte–Carlo simulations than any other method, because the simulations depend only
on the detection efficiencies of an air shower array, and because this method is not affected by
the elementary processes for the detailed features of muons, hadrons, Čerenkov photons andatmospheric fluorescence photons. Moreover, observations are not restricted with weather andthe phase of the moon, and there is no systematic uncertainty due to atmospheric monitoring.
The depth of shower maximum(Xmax ), which is one of the shower parameters used for
inferring composition of cosmic rays, is obtained from the equi–intensity method. It is necessary
to do a very high altitude experiment in order to measure the EAS development curves forprimary energies of interest. An EAS initiated by a proton with an energy of 1015 eV reaches
71zenith angle=θ1

zenith angle=θ2

depth=D 1

depth=D 2

F(N)size spectra

N1

N2

N(D)
equi−intensity curve

N1

F(N)

energyN2

D1

D2

Figure 4.3: The determination of the development of EAS of different sizes through the atmosphere from observations at different zenith angles θ; n(N) is the rate of occurrence of shower
size N.
the maximum development at about 600 g/cm2 and that initiated by an iron nucleus at about
470 g/cm2 . Thus the only site at such high altitude is Mt. Chacaltaya in Bolivia located at theatmospheric depth of 550 g/cm2 .
Note that one cannot compare equi-intensity curves with longitudinal development curves
immediately. A measured size with air shower array does not translate directly into a numberof electrons and positrons in the air shower, in part because the detected number of particles
is contaminated by muons and in part because the determined size is affected by the thresholdenergies of the detectors. Moreover these effects have primary energy dependences, and the equiintensity curves are affected by the energy-size relations, and energy spectra of components.

72

Chapter 5

Review of instruments to measurechemical composition of cosmic rays
In order to discuss the origins, the accelerations and the propagations of cosmic rays, it is themost important to measure the energy spectrum, the mass composition as a function of energy
and the degree of isotropy of their arrival directions. While the arrival direction distribution
can be directly measured, the other have to be deduced from measured EAS parameters. The
primary energy spectrum is derived through measurements of shower size spectra, Čerenkov
photon size spectra, muon size spectra, or shower development curves obtained from fluorescencelight measurements, etc. The composition of the primaries is deduced through measurementsof shower parameters such as the mean number as well as the fluctuations of the number ofmuons, the fraction of delayed hadrons, the high energy end of the hadron energy spectrum,or the depth of the shower maximum, etc., which are sensitive to mass number of a primaryparticle. Before we discuss the results of EAS studies with the BASJE MAS array, we shallbriefly summarize the other past and present experiments on the primary energy spectrum andon the chemical composition.

5.1

CRN

The detector referred to as Cosmic Ray Nuclei(CRN) Detector on a Space Shuttle measured
the nuclear composition of cosmic rays with energies well beyond TeV per amu. Figure 5.1
shows a schematic cross section of the instrument [47][87]. The instrument employed plastic
scitillators(T1, T2), gas Čerenkov counters(C1, C2) and a transition radiation detector(TRD).
This instrument measured nuclear charge Z, energy E over the range 40 GeV/amu to several
TeV/amu and the trajectories of particles. Owing to the limitation in the dynamic range of the
instrument and owing to the inherent magnitude of fluctuations in the transition radiation and
the Čerenkov signals, the CRN instrument was not designed to observe protons and He nuclei.
Nuclear charges Z are measured with the scintillators utilizing the Z 2 dependence of theionization loss. Their signals are nearly independent of the energy over the range. The scintillators also provide the coincidence triggers and the signals for the TOF measurements. The
gas Čerenkov counters are filled with a N2 /CO2 mixture at 1 atm, and provide energy measurements 40–150 GeV/amu. The TRD consists of six pairs of plastic fiber radiators and multiwire proportional chambers(MWPCs); the detection threshold for the transition radiation isreached at about 500 GeV/amu. For the energy below 500 GeV/amu, the MWPCs measure73

Figure 5.1: Schematic cross section of the CRN detector
the ionization losses of the particles, and is used for the coarse energy assignments. Above
500 GeV/amu, the yield of transition radiation X–rays rises rapidly with the energy up to 10
TeV/amu and provides the energy measurements. In addition to the energy measurements, theMWPCs also provide the position information to reconstruct the trajectory of a traversing particle. Consequently, the performance of the CRN instrument is characterized as; (1) the charge
resolution δZ = 0.2 for oxygen(Z = 8), δZ = 0.35 for iron(Z = 26), (2) the energy resolutions
δE/E at 100 GeV/amu and 1TeV/amu are 35% and 11% for oxygen, respectively, and thoseare 13% and 8% for iron, respectively, and (3) the trajectory resolution is about 1◦ .
The CRN instrument was flown in 1985 July/August in the Spacelab–2 mission of the Space
Shuttle Challenger. While the total flight duration was about 8 days, the net observation timeat the full aperture amounted 78 hours. Müller et al [71] reported the particle fluxes up to
2 TeV/amu for major primary nuclei as shown in Figure 5.2(left). In the same report, theauthors estimated the source abundances of the element groups from their measurements withassuming a simple leaky box model of the particle propagation in the Galaxy, and then they
showed the relative abundances of cosmic ray nuclei against the first ionization potentials(FIP)of the elements as shown in Figure 5.2(right). They pointed out the data which indicate thesudden suppression of particles when the FIP exceeds 10 eV. Thus, they concluded that the firststage injections of cosmic ray particles take place in stellar photo spheres at temperatures around104 K [84]. Moreover they concluded that the bulk of the cosmic rays must be accelerated frommaterial that originated in stellar eruptions or stellar winds and is enriched with ions formed
in stellar photo spheres, and also they suggested that some of this material might be ejected by
objects such as pre–supernova red giant stars or Wolf–Rayet stars.

74

Figure 5.2: (left): Differential energy spectra for the cosmic ray nuclei C, O, Ne, Mg, Si and Fe
measured with CRN(solid). Note that the fluxes multiplied with E 2.5 . For other markers, see
the reference [71]. (right): Abundances of the galactic cosmic ray nuclei at the source(GCRS),
relative to the local abundances(LG), are plotted against the values of first ionization potentialsof the elements [71]. All abundances are relative to iron.

75

5.2

RUNJOB

RUNJOB(RUssia–Nippon JOint Balloon–program) is a joint collaboration on the observationsof primary cosmic rays with the balloon–borne emulsion chambers. They launched six long
duration balloon on the trans–Siberian trajectories from 1995 until 1999. They reported experimental results of four flights in 1995 and 1996 [7]. Total exposure was 231.5m2 hr at the
average altitude of 32 km, and the cut–off rigidity was about 3 GV along the trajectories. The
detector on board had two chambers, each of which had an area of 40×50 cm2 and the schematic
illustrations are shown in Figure 5.3. The detector covered an energy range of 10–500 TeV for
protons, 3–70 TeV/amu for He, and 1-5 TeV/amu for Fe–group.

Figure 5.3: Chamber structure in 1995 and 1996 experiments.
The charge of a primary particle was determined from a darkness of a primary track on the
nuclear emulsion plate of the top layer [6]. The resolution of the charge determination, which
depends on the path length, was determined 0.167 for He and 1.089 for Fe in the case of 200µm
path length.
The primary energy E0 was determined from the shower energy ΣEγ with the linear relation
E0 = Cγ ΣEγ . The calorimetric determinations of ΣEγ was difficult, because the thickness ofthe calorimeter layer was only four radiation length. Thus, for relatively low E0 they estimatedΣEγ from maximum darkness in the transition curves of the spot darkness on X-ray films, andfor high E0 , ΣEγ were estimated from the emission angles of γ–rays. The energy resolution,∆Eγ /Eγ , was about 15%.
Within four flights in two years, they obtained 483 tracks of primary cosmic rays, andreported the primary particle flux of five element groups as shown in Figure 5.4. The remarkablefeatures in their results are:
• The proton spectrum is consistent with those reported by other groups but He spectrum
is not consistent with JACEE and SOKOL results. The spectra of proton and He spectra
have the same index ∼ 2.8.
• The CNO spectrum is slightly steeper than that of other groups.
• The slopes of the energy spectra of heavy components seem to become gradually harderwith the increasing mass number, i.e., ∼ 2.70 for CNO-group, and ∼ 2.55 for Fe-group.
76• Average mass is nearly constant over the wide energy range 20–100 TeV.

Figure 5.4: The left figure is the measured spectra of the primary protons and primary He
nuclei. The right figure is the measured spectra of the heavy components.

77

5.3

CASA–MIAThe CASA–MIA detector [18] was the ground–based array of 1089 surface particle detector
stations (Chicago Air Shower Array: CASA) and 1024 underground muon detectors (MIchigan
Array: MIA) located at Dugway, Utah. The layout is shown in Figure 5.5. The mean atmo-

Figure 5.5: The layout of CASA–MIA array. CASA is on the surface, and MIA is buried at 3
meter underground. The CASA stations are 15 m apart from each other.
spheric overburden is 870 g/cm2 . Each station contained four counters, each of which consistsof a two inches diameter PMT glued to the square sheet of an acrylic scintillator with 61 cm
on a side and 1.5cm thick. Each station had 1 radiation length lead sheet on its top. The array
was triggered when any three stations report at least three hit counters each.
The MIA counters were arranged in 16 patches, each of which contained 64 individual
counters. Each counter was a 1.6m by 1.9m acrylic scintillator viewed by a five inches PMT.
The MIA patches were buried at 3 meters underground and registered signals induced by muons
with energies exceeding 750 MeV at the surface.
In order to estimate the primary energy independent of the type of primary particles theyused a combination of the muon size Nµ and the electron size Ne . They reported that theempirically determined optimal value of the energy parameterization (Ne + 60 × Nµ ) (Figure5.6) has the most compositional insensitivity, and the best value of the factor of Nµ increases
slightly with zenith angle reaching about 65 for the zenith angle of 45 degrees. This methodof the energy determination yielded systematic differences between iron and proton energy
assignments with less than 5%. The average absolute values of the energy reconstruction errors
78

Figure 5.6: The energy parameter log(Ne + 60 × Nµ ) as a function of the energy for simulated
proton(open circle) and iron(filled triangle) vertical showers.
decreased from about 25% near 1014 eV to about 16% at 1015 eV.
The CASA–MIA group reported the all particle energy spectrum as shown in Figure 5.7[34], and they concluded that the knee of the energy spectrum is rather a smooth transition
over energies from 1015 eV to 3.0 × 1015 eV.
Moreover, this group studied the mass composition of cosmic rays [35]. They classified theobserved events into “iron-like” group(called Heavy group) and “proton-like” group(called Light
group) by parameters: ρe , ρµ , and α with “K”Nearest Neighbor (KNN) test, where ρe is thedensity of the surface particle and α is the slope of the lateral distribution near the core, and ρµis the muon density at large core distance. They showed < ln A >(Figure 5.8) and the energy
spectra of the data groups classified according to the KNN–identified mass(Figure 5.9). In
their results, the composition of cosmic rays agrees with the direct measurements near 1014 eV,and evenly consists of both the light elements and the heavier elements. Moreover the average
mass increases with increasing energy, in the energy range from 1015 eV to 1016 eV, whichcoincides with the bending point of the energy spectrum.

79

Figure 5.7: The observed energy spectrum of the CASA–MIA experiment compared with thoseof Tibet and of Akeno.

Figure 5.8: The average composition hln Aiof CASA–MIA data(filled circle). In the simulations
for this analysis the QGSJET interaction generator were used. The open squares indicate themean mass of the direct measurement of JACEE.

80

Figure 5.9: The energy spectra of each composition groups measured by the CASA–MIA experiment.

81

5.4

HEGRA–CRTThe HEGRA experiment is a multi component air shower installation with an air shower array;a Čerenkov light detector array and a independent Čerenkov telescope system [3] on La Palma,
Canary Islands. The particle detector array consisted of 221 scintillator stations of 1m2 areaeach(Figure 5.10). The Cosmic Ray Tracking(CRT) detectors were installed and operated in the
HEGRA air shower array, and had the lower threshold energy and the good angular resolution
even for small air showers. One CRT detector consisted of two 2.5 m2 circular TPC type drift
chambers, and a 10 cm thick iron plate was placed as a muon filter between both chambers[12]. In total, ten CRT detectors were installed and operated in HEGRA site for about three
years(1993–1996).

Figure 5.10: Map of the HEGRA site showing only the scintillator stations (square dots) and
CRTs (circles). The dash-dotted line marks the area of shower core positions accepted in theanalysis.
For studying the mass composition of cosmic rays, the analysis of CRT data were based ona statistical method using angular distributions of particles in air showers [11]. The authors
reported that a measured average radial angle was closely related to an average longitudinal
shower development [11]. Most muons are produced near the shower axis and rarely deflectedby more than a few tenth of a degree due to multiple scattering and the geomagnetic field.
Thus, muon radial angles can, in principle, be transformed into muon production heights for in82

a given core distance region.
Figure 5.11 shows histograms of muon radial angles for different core distance intervals [13].
Resulting measurements of median muon radial angles for 10 m core distance bins and for six
different intervals of air shower sizes are shown in Figure 5.12. They estimated the cosmic
ray mass composition by comparing the measured data with proton and iron nuclei simulations.
The resulting average mass composition is shown in Figure 5.13. They concluded that the resultshows an increase of hln Ai in 1014 − 1015 eV and it agrees very well with the highest energy
points of the direct measurements made by the JACEE collaboration, and an extrapolation of
the directly measured spectrum with an additional constraint with a kink at a fixed rigidity of
about 200 TV is consistent with the measured flux and composition.

Figure 5.11: Histograms of muon radial angles observed with CRT detectors in different intervalsof core distances for 15000 < Ne < 50000. Note that particles coming from the shower axis
have negative radial angles.

83

Figure 5.12: Median radial angles of muons measured with CRT detectors as functions of horizontal distance of the CRT detectors from the shower cores. Simulations are superimposed forpure protons and irons as well as the mixed compositions corresponding to direct measurementsof the JACEE collaboration both above 45 TeV(J45) and above 370 TeV(J370).

Figure 5.13: The mean mass composition obtained by comparing HEGRA–CRT data withsimulations based on the VENUS/GHEISHA and DPM/Isobar interaction models.

84

5.5

CASA–BLANCA

Observations of atmospheric Čerenkov photons emitted by charged particles in air showers
are considered as one of the integrated measurements of the longitudinal developments of air
showers. For this purpose one approach is to measure lateral distributions of the Čerenkov
photon densities on the ground. Slopes of the lateral distributions are related to depths of the
maximum shower developments, and hence, they are related to the composition of the primarycosmic rays. However, a total number of the Čerenkov photons is a good estimator for a primaryenergy. Therefore, the detections of a large number of Čerenkov lateral distributions can provide
information on how the composition changes with the energy.
To obtaine high quality Čerenkov lateral distribution data, the BLANCA(Broad Lateral
Non–imaging Cherenkov Array) detectors were built in CASA installation. Using CASA as
the cosmic ray trigger, the BLANCA detectors were operated at clear and moon less nights.
In the analysis, they used the CASA data to determine the shower core positions and thearrival directions, and used the BLANCA data for the precise measurements of Čerenkov lateral
distributions [28][29].
BLANCA [19] consisted of 144 angle–integrating detectors which recorded the lateral distributions of atmospheric Čerenkov photons. The BALNCA detectors were not uniformly spaced
but had an averaged separation of 35–40m as shown in Figure 5.14. Each BLANCA detector
contained a large Winston cone which concentrated the light striking an 880cm2 entrance aperture onto a PMT. The effective concentration ratio is 15, the effective half–angle is about ∼ 10◦ ,and a typical BLANCA unit had a detection threshold of approximately one blue photon per
cm2 .

Figure 5.14: The CASA–BLANCA array.
CASA–BLANCA were operated in 90 moon less nights between January 1997 and May
85

1998. After removing the periods of hazy or cloudy weather, approximately 460 hours of theČerenkov observations remained. Events were selected with the arrival directions, with the core
locations and with conditions of measurements. The geometrical and temporal cuts resulted in
an exposure to cosmic rays of 1.83 × 1010 m2 sr.
The obtained lateral distributions of photon densities were fitted to an empirical function
with three fitting parameters, C120 , s and β:½
C120 exp [s(120m − r)] 30m < r ≤ 120m
C(r) =(5.1)
C120 (r/120m)−β120m < r ≤ 350m
The energy of each air shower was determined from C120 using the relation between C120 andprimary energy obtained with the Monte Carlo simulation studies. As a result, they reported
the differential all particle cosmic ray flux as shown in Figure 5.15, and concluded that the datashows a knee with one–half decade wide, and best–fit knee energy is 2.0 PeV with power law
indices of −2.72 ± 0.02 for lower energies, and −2.95 ± 0.02 for higher energies above the knee.

Figure 5.15: The differential all particle cosmic ray flux measured by CASA–BLANCA.
The mean primary mass was directly determined from the Čerenkov lateral distribution
slope s. With their Monte Carlo studies, they insisted that, at fixed energy, s depends linearlyon hln Ai, and the determined < ln A > is shown in Figure 5.16. Their result shown in the
figure indicates that the cosmic ray composition is lighter near 3 PeV than that around either
300 TeV or 30 PeV, and becomes heavier with the increasing energy above 3 PeV.

86

Figure 5.16: The mean logarithmic mass < ln A > measured by CASA–BLANCA as a functionof the energy. The four sets of symbols show the BLANCA data interpreted using CORSIKA
coupled with the indicated hadronic interaction model.

87

5.6

DICE

The DICE(Dual Imaging Cherenkov Experiment) project was another experiment of atmospheric Čerenkov light measurements at the CASA–MIA site. The DICE project had twotelescopes and each telescope consisted of a 2m diameter f/1.16 spherical mirror with a focal
plane detector of 256 close packed 40mm hexagonal PMTs which provide ∼ 1◦ pixels with an
overall field of view 16◦ × 13.5◦ centered about the zeith. The telescopes were on fixed mounts
separated by 100 m.
A cosmic ray event in the field of view produced a focal plane image on the PMT cluster,and the shape of the image reflected the emission point distribution and the intensities ofČerenkov light. When the direction of an air shower axis and the distance of an air shower corefrom the telescopes were known, one can reconstruct the longitudinal development of the air
shower with a simple geometry. Essentially, this procedure is geometrical and is independentof numerical simulations except for calculations which determine the angular distribution ofČerenkov photons around shower axes. This feature is an aim of this experiment, and then
DICE was designed to be independent of the details of the air shower simulations.
The DICE group determined the location of shower maximum in the atmosphere(Xmax ) byfitting the shape of a shower image in each of the DICE PMT clusters. An accurate determination of energy was derived from a combination of the amount of Čerenkov photons and Xmax .
A total shower energy E and a primary particle mass A were obtained from the geometry,Čerenkov size (Ch) and Xmax . The exact forms of these fitting functions were obtained byMonte Carlo simulations, and given as follows:
Ch = C0 E(PeV)γ A−² e−βr(m)(5.2)

Xmax = X0 + Xm (log E − log A),

(5.3)

where
C0 = 1.89 × 106 ,(5.4)

γ = 1.144 + 0.0905 log A,(5.5)

β = 0.0161(1 − 0.128 log A) + 0.124 log E(PeV)(1 + 0.322 log A),(5.6)² = 0.186,

(5.7)
2

(5.8)

Xm = 80g/cm .

(5.9)

X0 = 560g/cm ,2

The events were collected by DICE in coincidence with CASA–MIA over a period from mid
1994 to early 1996 [88]. The all particle flux by the DICE experiment is shown in Figure 5.17.
They calculated the mass composition with two sets of measure values: (1) the location of Xmax
and the fitted shower energy, (2) the muon and the electron sizes in combination with the fitted
energy. The results are shown in Figure 5.18. These results do not support a simple “rigidity
steepening” which would lead to steady increase in the mass composition across the knee region.
The dotted curve is a prediction with a model where a source spectrum is proportional to E −3
above the knee region and the energy dependent escape from the Galaxy was assumed to reacha plateau near the knee region. The dashed line, to which the data seem more consistent, was amodel which introduces a proton source to compensate the lost flux above the SNR acceleration
cutoff rigidity. This model was similar to that of Szabo and Protheroe described in Chapter
88Figure 5.17: The spectra near the knee by the DICE group and by some other groups.
3 [89] where an extra-galactic component of particles produced by AGNs provides the cosmic
rays at the knee region.
The DICE group made reanalysis with considering the systematic errors including the mirror
spot sizes, the dead spaces between the PMTs and the electronics saturations [64]. The masscomposition as the results of their reanalysis have not be published yet.

89

Figure 5.18: (upper) The result for the mean mass measurement by DICE and the model
predictions. (lower) The result of (p+α)/all.

90

5.7

BASJE–Čerenkov

The BASJE experiment observed air showers at Mt. Chacaltaya (5200 m a.s.l.) in Bolivia.
Shirasaki et al. [82] reported the result of their observation of atmospheric Čerenkov photons
associated with air showers detected by Čerenkov detectors and the Small Air Shower(SAS)array [56](Figure 5.19). For estimation of the cosmic ray composition, they observed Čerenkov
light pulse shapes which strongly reflect longitudinal developments of air showers.

Figure 5.19: The arrangement of the scintillation detectors of the SAS array and the Čerenkov
detectors.
One Čerenokov light detector consists of seven 5–inch PMTs. In 1995 and 1996, they used
two Čerenkov detectors separately at 150 m from the center of the array(point C1 and C2 on
Figure 5.19) to achieve wide effective detection area and enough statistics at the energy range
1015.8 − 1016.5 eV. In the next year, four Čerenkov detectors were install at same position(C3)and the signals of the detectors were added with fast buffer amplifiers for enhancement of signals
to noise ratio and to decrease the threshold energy to 1015 eV.
They used the parameter T10−90 , which is defined as the rise time of a time–integratedČerenkov pulse, for estimating the mass composition. Figure 5.20 shows the composition dependence of this parameter. With comparing the observed T10−90 distributions with the simulated
ones for three different mixtures of components(ex. Figure 5.20), they estimated the mean mass
number as shown in Figure 5.21.
They discussed four different source models of cosmic ray particles, and concluded thatconsidering additional heavy matter ejected from pre–supernova, a stellar wind model proposed91

Figure 5.20: Comparison of the cumulative fraction distribution of experimentally determined
T10−90 (dots) with the distributions expected for proton, carbon and iron primaries. These areresult for (a) log Ne = 6.0 − 6.2, R = 140 − 150 m obtained in the low energy mode observations,and (b) log Ne = 6.5 − 7.0, R = 150 − 160 m obtained in the high energy mode observations.
by Biermann [15] shows a good agreement with both the all particle flux and the mean mass
number.

92

Figure 5.21: Comparison of BASJE–Čerenkov results with the other direct and indirect experimental results.

93

5.8

KASCADEThe KASCADE site [61] is located at the site of the Forschungszentrum Karlsruhe, Germany
at an altitude of 110 m a.s.l., and they observe EASs especially in the knee region. Their air
shower array consists of three major detector components(Figure 5.22): the scintillator detector
surface array, the muon tracking detectors and the central detector complex.

Figure 5.22: A schematic layout of the KASCADE experiment.
The surface array covers an area of 200 × 200m2 and consists of 252 detector stations for thedetection of electron and muon components. Each station of the inner part of the array contains
four liquid scintillation counters to detect electromagnetic components of an EAS. Each outerstation contains only two of these detectors, but it has plastic scintillators to detect muons below
an iron–lead absorber. Using the information of these detectors they reconstruct the shower
core positions, the arrival directions of EASs, the total numbers of electrons, Ne , and the socalled the truncated muon numbers, Nµtr. . Nµtr. is the number of muons inside the concentricring around a shower core with the inner radius of 40 m and the outer radius of 200 m. The
hadronic components of an EAS is studied with a 320m2 large iron calorimeter. The impact
point, the energy, and the direction for hadrons are reconstructed above a threshold energy ( =50 GeV ).
Ulrich et al. [94] determined the primary energy spectra for different mass groups by analyzing the shower size spectra of Ne and Nµtr. . The data used in the analysis consist of Ne and Nµtr.
spectra in three zenith angle bins; 0◦ < θ < 18◦ , 18◦ < θ < 25.9◦ and 25.9◦ < θ < 32.3◦ (Figure945.23). Their analysis was based on the relation between the measured shower size spectrum

Figure 5.23: Electron size spectra (upper) and truncated muon size spectra(lower) in the intervals of zenith angle.
dJ/d log N and the primary energy flux dJA (log E)/d log E for primary particles with mass
number of A,
X Z +∞ dJA (log E)
dJ=
pA (d log N |d log E)d log E(5.10)
d log Nd log E−∞
A=1

Here pA (d log N |d log E) is the probability for a primary particle of type A with an energy
E to be reconstructed as an air shower with shower size N . The measured size spectrum
transfered to a vector of numbers of events, yi (i = 1, · · · , n), binned with EAS sizes (bin width,∆ log N = 0.05) The corresponding vector for the energy spectrum, xAj (j = 1, · · · , m), binned
with the energies (bin width, ∆ log E = 0.1). They obtained the vector xAj with solving therelation
m
X
A A
Rij
xj .
(5.11)
yi =
j=1

95

A is an element of the response matrix RA , which defined as,Here Rij

R
A
Rij
=

i d log N

R

dJA (log E)
d log E pA (d log N |d log E).
R
dJA (log E)
j d log E d log E

j d log E(5.12)In the analysis, the data vector y consists of 6 spectra. The mass composition of primary cosmic
rays was assumed to have four components; protons, He, C and Fe. The response matrix RAof the four components are arranged in one response matrix R consisting of 24 sub matrices,and the vector x represents the four unknown energy spectra.
They obtained the elements of RA with the simulation calculations by CORSIKA–QGSJET
model. The consequence of the analysis is shown in Figure 5.24 as the all particle spectrum andthe energy spectra for four components. They summarized about their results as follows:

Figure 5.24: The total energy spectrum obtained with the KASCADE (Ne , Nµtr. ) analysis.
1. The energy of the knee on the all particle spectrum was about 5 PeV.
2. Comparing with the results of other experiments, the all particle spectrum showed a good
agreement with the data of Tibet and HEGRA.
3. From the differential energy spectra multiplied by E 3 as shown in Figure 5.25, the indexesof the power law describing the individual energy spectra below the knee was steeper for
proton than for He, and was flatter for C than for He.
4. From the differential spectra, the factors of the bending point energies of the individual
elements relative to that of proton were 2–3 for He, 6–9 for C and ∼ 20 for Fe. Thissuggest a rigidity dependence of the individual knees.

96Figure 5.25: The energy spectra of proton, helium and carbon primaries between 1 PeV to 10
PeV as a result of the KASCADE (Ne , Nµtr. ) analysis.
5. They obtained the mean logarithmic mass hln Ai as shown in Figure 5.26. It slowly started
to increase at about 2.5 PeV with increasing energy, and hln Ai becomes a constant above
45 PeV.
The KASCADE collaboration also reported the estimated mass composition using thehadronic component of EASs measured by the hadron calorimeter [25]. In this analysis theshowers were classified according to their truncated muon size Nµtr. measured by the surface
detector array. They estimated a mean mass parameter λ as a function of Nµtr. from the analysis of the six observable values, namely, the lateral hadron density distributions, the distance
distributions in minimum–spanning trees, the lateral energy density distributions, the energy
spectra of hadrons, the energies of the most energetic hadron, and the distributions of the fractional energy of hadrons with respect to the most energetic hadron in the showers. Figure 5.27
shows the distributions of two observables comparing with the simulations. λ is linearly relatedwith hln Ai as λ = hln Ai / hln 56i. It indicates the distance to the two extreme compositions:
λ = 0 for a pure proton and λ = 1 for a pure iron composition. The result of their estimationsof λ is shown in Figure 5.28. Consequently, they combined the mean values < λ > of all theobservable values and for two interaction models, and determined the corresponding < ln A >as shown in Figure 5.29.
The hadron analysis result follows the trend indicated by the RUNJOB data and by themeasurement by Shirasaki et al [82]. However this result shows a too heavy composition comparing their (Ne , Nµtr. ) analysis result. The authors pointed out that the conflict between two
results is due to the simulations, and that some investigations have revealed that the simulations generate too many hadrons at observation level [5], [77]. They mentioned that the missing
hadrons are interpreted by the hadron analysis as a heavy composition, and that generating
too few electrons accounts for the light composition found in the analysis using the electron
97

Figure 5.26: The mean logarithmic mass by the KASCADE (Ne , Nµtr. ) analysis.
component only.

5.9

Summary of recent measurements

Here I summarize the recent ground–based measurements of the chemical composition. All thegroups detected not only the electron components of EAS particles, but also other components,
namely Čerekonkov photons, muons and/or hadrons, because the angular spreads of Čerenkov
photons, muons and hadrons relative to a shower axis and the pulse shapes of Čerenkov photons
reflect the longitudinal developments of EASs. In addition, the multiplicity of muons andhadrons strongly depend on the mass of the primary particles. The results of the measurements
can be classified into the three groups; (1) hln Ai increases with the increasing primary energy,
and the averaged mass is heavy (hln Ai ' 2.8) at the knee energy (∼ 1015.5 eV), (2) hln Ai
increases with the increasing primary energy, and the averaged mass is light (hln Ai ' 1.6) at
the knee energy, and (3) hln Ai is almost constant up to 1016 eV.
From the classification of the results and of the observable values, it is concluded that:
1. In spite of using the similar observable parameters for estimations of chemical composition,
the results often seem to conflict between different experiments.
2. At the same sites, the results are inconsistent when the different observable are used.
3. The results with the detections of Čerenkov lights except the BASJE group show the sametendency that hln Ai is energy independent.
The first two conclusions reflect that simulations are imperfect as descriptions of the nature of
EASs. The KASCADE group pointed out that the simulations used for their analysis seem togenerate too many hadrons at the observational level. The last conclusion shows the difficultiesof the corrections both of the efficiencies and the atmospheric attenuation for photons. Moreover,98

Figure 5.27: Two examples of the hadronic observables measured with the KASCADE hadroncalorimeter for estimation of the primary cosmic ray mass composition.
we consider as these two effects lead to the observational biases through the parameterizationsof the lateral distributions of Čerenkov photons.
Finally, we conclude that the measurement of chemical composition must be based on theobservable independent on simulations and the observational conditions. Thus, I chose theequi–intensity method analysis for the observed EASs at the high altitude equivalent to the
maximum development of EASs.

99

Figure 5.28: The mass parameters obtained from the six observables of the KASCADE experiment.

100

Figure 5.29: The mean logarithmic mass vs the primary energy reported by KASCADE(hadron
analysis) and by other experiments.

101

Chapter 6

BASJE MAS array
This chapter is detailed description about BASJE1 MAS array, including the experimental
details about the site, the detectors and the electronics for the experiment.

6.1

Site

The BASJE MAS(Minimum Air Shower) array is located at Chacaltaya Cosmic Ray Observatory (16◦ 200 5200 S, 68◦ 070 5700 W ) in the Republic of Bolivia. This observatory was established in
1962. The altitude of the site is 5200 m above sea level, which is the highest altitude among
the constantly operated observatories all over the world. The site is at about 35 km north–eastof the capital city of Bolivia: La Paz.
The altitude of 5200 m a.s.l. corresponds to the atmospheric depth of 550 g/cm2 . According
to many Monte Carlo simulations, EASs with primary protons of 1015 eV reach their maximum
developments at about 600 g/cm2 and those with primary iron nuclei at about 470 g/cm2 .
Thus, we can observe EASs induced by primary protons before their maximum developments
at this altitude, and this site is most suitable to investigate cosmic rays around the knee region.
The MAS array constructed as an improvement of the SAS(Small Air Shower) array [57]combined with the detectors of the LMC array [98] which were specially designed to observe
the supernova SN1987A. The aim of the construction of the MAS array were, with arranging
the detectors higher density than the SAS array, to achieve the lower threshold energy and1

The BASJE(Bolivia Air Shower Joint Experiment) collaborators:
O.Burgoa, A.Furuhata, D.Harada, F.Kakimoto, Y.Kurashina, S.Ogio, H.Tokuno and K.Yamada: Department of
Physics, Tokyo Institute of Technology, Meguro, Tokyo, 152-8551, Japan
H.Furukawa, M.Kubo, A.Morizawa and H.Yoshii: Department of Physics, Ehime University, Matsuyama, Ehime,
790-8577, Japan
E.Gotoh, H.Nakatani, K.Nishi, S.Shimoda, N.Tajima and Y.Yamada: The Institute of Physical and Chemical
Research, Wako, Saitama, 351-0198, Japan
T.Kaneko: Department of Physics, Okayama University, Okayama, 700-8530, Japan
Y.Matsubara: Solar-Terrestrial Environment Laboratory, Nagoya University, Nagoya, A ichi, 464-8601, Japan
P.Miranda and A.Velarde: Instituto de Investigaciones Fisicas, Universidad Mayor de San Andres, La Paz, 8635,
Bolivia
Y.Mizumoto, Y.Shirasaki, and Y.Tsunesada: National Astronomical Observatory, Mitaka, Tokyo, 181-8588,
Japan
K.Murakami: Nagoya University of Foreign Studies, Nissin, Aichi, 470-0197, Japan
Y.Toyoda: Faculty of General Education, Fukui University of Technology, Fukui, 910-8505, Japan102

the higher accuracy of size determinations. At present, the MAS array consists of sixty-eight
unshielded scintillation detectors and a shielded scintillation detector as shown in Figure 6.1
and 6.2. These detectors were deployed over an almost flat field of approximately 104 m2 . The
average spacing between the detectors is about 5 m at the central area, and 30 m at the outerarea. The maximum distance in the elevations among the G–detectors is 20.9 m, and among
the N– and F–detectors, which are widely spread around the array, that is 59.8 m.

(m)

150

100

50

0

-50

-100

-150
-150

-100

-50

0

50

100

150

(m)

Figure 6.1: BASJE MAS array. Only N–, F– and G–detectors are shown in this figure.
The signals from all the detectors are sent to the central electronics station through coaxial
cables for being processed and for being recorded. The details of the detector equipment andthe electronics will be described in the following sections.

103



G13

G9

G1

G5

+30m

Detector:
G


NT


L


S

NT7
NT4
S16

S11

S5

S13

G10

G14

NT1

S9

G6

G2

S3

L7
L11

L3

S17
−30m

NT8

S20

L4

L12

L1

L9
NT5

S1

S6

L8

NT2

S15

S2

S21
L10

G15
G11

L6

S19

NT6

−30m

S4

G7

G3

S8

S12

NT9

G12

S7

L2

S18
S14

G16

+30m

L5

S10

G8

NT3

G4

Figure 6.2: The central area of MAS array. The 60m2 muon detector(the shielded detector) is
located at the center of the array (under the L–detectors).

104

6.2

Detectors

As mentioned in the previous section, the MAS array consists of the sixty–eight surface detectors
and one shielded muon detector, and each detector is comprised by a set of plastic scintillatorsand photo–multiplier tubes(PMTs). The schematic cross sectional illustrations of the surface
detectors are shown in Figure 6.3, and that of the shielded muon detector is shown in Figure
6.4.
R1828
(timing)

scintillator(5.0cm)

RCA8055(density)L

R877(density)G
scintillator(5.0cm)
R1250(timing)S

R1512(density)

scintillator(5.0cm)
R329(timing & density)

NT
XP2040(timing & density)

F
N

scintillator(7.5cm)

RCA8055(density)

scintillator(7.5cm)

scintillator(7.5cm)

Figure 6.3: The schematic illustrations of the surface detectors of the MAS array.
The detectors are classified into G, NT, L, S, N, F and Mu detectors. The specifications ofeach type of the detectors are listed in Table 6.1. Each of G and L detectors has two PMTs; one
PMT measures the number of particles and the other PMT is a fast-timing channel to measurearrival times of particles. Each of S and NT detectors has only one PMT, and an output signalof it is divided to the density channel and the fast-timing channel. Each of N and F detectors
has only one PMT for the density channel.
We also installed a shielded muon (Mu) detector of 60 m2 area at the center of the array todetect muons in air showers. This detector is the matrix of 3 × 5 scintillation detectors of 4 m2
area each as shown in Figure 6.4. The threshold energy of this detector is 600 MeV for verticallyincident muons. The characteristics of the Mu detector are reported in reference [58]. Signalsof all the detectors are sent to the central station through coaxial cables. The low voltages forthe pre–amplifiers and the high voltages for the detectors are supplied from the central station.

105















































































































































































































































































































































































































































































sheilding material
plastic scintillator(5cm thickness, 4m 2)

PMT(K1328 or R1512)

Mu3

Mu6

Mu9

Mu10 Mu13

Mu2

Mu5

Mu8

Mu11 Mu14

Mu1

Mu4

Mu7

Mu12 Mu15

Figure 6.4: The shielded muon detector has the fifteen 4 m2 scintillation detectors.

detector
G
NT
L
S
N
F
Mu

number of
detectors
16
9
12
21
4
6
15

area(m2 )
1.0
0.87
4.0
1.0
0.83
0.83
4.0

thickness ofa scintillator(cm)
5
7.5
5
5
7.5
7.5
5

PMT
densitytiming
R1512
R1828
XP2040
R1512
R1250
R329
RCA8055
R877
R1512

Table 6.1: the characteristics of the each type of the detectors.

106

6.3

Electronics

The signals of the detectors are sent to the central station and processed to be recored. A block
diagram of the data acquisition system of the MAS array is shown in Figure 6.5.
Particle Detectors

L4,L5,L8,L9
Timing(58ch)

Discriminators

Density(74ch)

GATE:4µ s

exponentialpulse

Log output

Main Amplifiers

NT density(9ch)

trigger

Coincidence

ADC
GPS Clock

Conputer with DVD−RAM

TDC

ADC

IGOR

LAM
Linux

Start/Stop/Clear

Crate Controller
CAMAC

Figure 6.5: Block diagram of the MAS array electronics.

6.3.1

Local particle density measurements

Output signals of each density channel are processed by pre–amplifiers. The pre–amplifiers
are installed at the vicinity of each PMT in the detector box to transmit signals with a low
impedance to reduce noises and distortions of pulse profiles while signals are fed through cables.
For the measurements of air showers with very high energy cosmic rays around the energiesof 1015 eV, the required dynamic range for local particle density measurements is four orders ofmagnitude. So that, a linear amplifier is not suitable to use because of saturations of signals.
Therefore, in our experiment, we use following log–amplification methods to achieve the wide
dynamic range. The signals of a PMT are shaped into exponential pulses by the RC circuit in
the pre–amplifier, and fed into a main amplifier at the central station. After that, an amplified
signal with the height VN is discriminated at the particular level Vth and converted to a rectangle
output pulse with the width of TN (Figure 6.6). The pulse width TN of the output pulse islogarithmically proportional to the pulse height VN which corresponds to the number of particles
N . The relation of TN and VN is written by
Vth = VN exp (−TN /τ )(6.1)

= V0 exp (−T0 /τ )(6.2)

where τ is the decay constant of the exponential pulse, and depends on the constants of the RC
circuit. T0 is the pulse width of the output which corresponds to single particle signal V0 . Thus107the number of shower particles N at a detector is obtained byµ¶VNTN − T0
N=
= exp
V0
τ

VN

V0

N particles

single particle

Vth

Vth

T0

TN
Figure 6.6: Log–amplification method.

The constants, V0 , T0 and τ , are calibrated for all the local particle density channels. In our
experiment, the pulse height of the single particle, V0 , is 2.7 Volts at the linear–output of themain amplifiers for all the density channels. V0 s for all the density channels are monitored daily
with a pulse height analyzer(PHA), and adjusted with changing the applied high voltages for
the PMTs(Figure 6.7). The applied high voltages range from 700 to 1200 V, and the counting
scintillator
PMT
PHA
Pre. Amp.

Inverter

Pulse Shaping Amp.

single peak
G,N,F 1400ch
L,S 1900ch

(L,S)
G NF
S
S
L

main amplifier linear out

detector selecter(G,N,F)

Figure 6.7: The diagram of the high voltage adjustment.
rates of the signals of greater than the single particle level are about 550 Hz/m2 .
The main amplifiers are designed for the T0 values to be of 1 µs for L and S detectors, andto be 8 µs for G, N, F and Mu detectors. However, T0 depends on the manufacturing tolerancesand variations of the constants of the circuit, and T0 s are measured by the calibration system
shown in Figure 6.8. The single particle input is generated with a set of a LED light source,108the PMT and the pre–amplifier. The brightness of the LED is adjusted to a PHA peak being
2.7 V at the linear–output of the main amplifier which is calibrated. The corresponding pulse
width of the log–output pulse is T0 and it is measured with a pulse width analyzer(PWA). The
error for T0 is about 0.05 µs for L, S detectors, and is 0.1 µs for other detectors.
Pulse Gen.
Tektronics 2101

PHA
Inverter

Pulse Shaping Amp.

single peak
G,N,F 1400ch
L,S 1900ch

(L,S)

LED

G NF
S
S
L

detector selecter

PMT(G,N,F)

PWA

linear out

Pre. Amp.

resolution G,N,F 100ns
L,S 50ns

Pulse WIdth Distribution
main amplifier

log out
48inputs

Figure 6.8: The diagram of T0 measurement.
The time constant τ is calibrated using the relation between the input pulse height and thelog–output pulse width. The block diagram of the τ measurement is shown in Figure 6.9. The
input pulse is generated as same as T0 measurements, and is variable with transmittance filters.
Figure 6.10 shows an example of the measured relations between the input pulse height and thelog–output width. The slope of this figure gives a τ value. The averaged τ value for L and S
detectors is 1 µs, for G, N, F and Mu detectors is 10 µs, and the linearities of the relations are
confirmed up to more than 104 particles for every density channels. The error for τ is about
0.05 µs. The log–outputs of the main amplifiers are processed to digital values with the ADC
modules which were developed by the Advanced Engineering Center, RIKEN. The resolutionof the ADC corresponds to 10 % for the measured number of particles. Including the errors of
T0 and that of τ and ADC resolution, the determination error for particle densities is about 3%
for L detectors, 11% for S detector and 10% for other detectors.
Since we do not use the log amplification system for the measurements of the local densitywith NT detectors, the dynamic range of the number of particle is very narrow, and thus thedata of NT density channels were not used for air shower size determinations.

6.3.2

Fast timing measurements

The accuracy in the determination of arrival directions of air showers depends on the time
response of each detector. Then, for fast–timing channels, we use PMTs with fast rise times(less than 10 ns), and the fast–timing signals are fed into the central station without pre–amplifications in order to avoid the distortions of the leading edges on signals. However, weused the pre–amplified signals for S detectors’ fast–timing channels because the pre–amplifiersfor S detectors [93] have Burr Brown BUF600 buffer amplifiers and have a fast time response,109height: 8V(MAX)width: 20ns(MIN)

Pulse Gen.
Tektronics 2101

filter selection dialLED+
filters
PMT
PWA
Pre. Amp.
Pulse WIdth Distribution
main amplifier

log out
48inputs

Figure 6.9: The diagram of τ measurement.
in which the typical rise time of outputs is 12 ns for the PMT(HAMAMATSU R329, typicalrise time = 3.2 ns) output signals.
The fast–timing signals are discriminated and are fed through the delay cables. Finally, therelative arrival times of these signals are measured with CAMAC TDC(LeCroy 2228A). The
discrimination level is set at the pulse height of 0.8 particle equivalent and the typical pulse
frequency is 700 Hz/m2 . The transit times of the PMTs and the cable delays for the fast–
timing channels are calibrated with a reference detector. The measured fluctuations on total
delay times are about 1 ns for G, 2 ns for L and 3 ns for S and NT detectors.
The time resolutions of the TDCs are 50 ps for the L detectors, 250 ps for the S detectors
and 500 ps for the other detectors. The linearities between the time difference and the TDC
counts are calibrated with a CAMAC TDC tester(REPIC RPC–070). Figure 6.11 shows anexample of the TDC calibration. The typical error in the linearity of the TDC is 0.001 %, andis sufficiently small for our requirement.

6.3.3

Trigger

A trigger signal of MAS array is generated with a 4–fold coincidence of the undelayed fast–
timing signals of the central four detectors, i.e., L4, L5, L8 and L9 within the gate time windowof 4 µs. The discrimination levels for the trigger channels are set to the pulse height of 0.8
particle equivalent. The triggering rate is about 8.5 Hz.
When one of the undelayed hit signals of the trigger detectors is occurred, the coincidence
module generates the common start signal for starting digitization in the CAMAC TDCs. Consecutively, if the four trigger detector signals coincide within 4 µs, a triggering signal is to ADC
modules for density measurements, and to sent the GPS clock module to register the event
date/time. This trigger signal is also sent to a CAMAC LAM(Look At Me) module in order togenerate an interrupt signal for starting the data acquisition procedure with the CAMAC crate
controller and the personal computer. If a coincidence is not occurred, the coincidence module
sends a clear signal to all the modules(Figure 6.5).
110Figure 6.10: An example of τ measurement. The filter attenuation is logarithmically proportional to the pulse height. In this example, the τ value is determined 9.82 µs for the F–1
detector.

6.3.4

Data acquisition system

A local time when an air shower hits the array is measured with a digital clock calibrated
with GPS. The stability of this clock is better than 10−7 . This clock is also developed by the
Advanced Engineering Center, RIKEN. This GPS clock and the ADC modules are accessible
through a CAMAC IGOR(Input Gate and Output Register, Kinetic 3063). When the dataacquisition system is triggered, a set of data, i.e., the output of TDC and ADC and the local
time of the event, is sent to the CAMAC crate controller(TOYO CC7700) and is transfered to
the DAQ computer through the PCI board(TOYO CC/PCI).
The DAQ computer is an IBM compatible PC with the cpu of MMX Pentium(166MHz)
and the Linux OS(kernel version is 2.0.35). In our experiment the size of data is 348 bytes per
events, and thus the total amount of data is reached 0.3 GB for a continuous operations in a day.
In order to store this huge size of data, we adopted DVD–RAM drive(Panasonic LF–D102) anddisks(2.6 GB) with considering its capacity, robustness, reliability and cost. The data storage
has a sufficiently fast transfer speed, which is measured about 600 kB/s on the Linux OS.
A set of application programs to control the data acquisition procedures is developed. Thisis comprised by three programs and these processes are communicate each other through thefunctions on the Liunx called the shared memories and the signals(Figure 6.12).
The program named “eco” is a data acquisition program with the CAMAC functions. The
CAMAC functions included in the program were developed by the KEK online group. This
program reads out data from CAMAC modules and stores on a local hard disk drive(HDD).
One data file contains 150,000 events(about 52.2 MB). Each data file is transfered to DVD–RAM with the program named “datast”. The forty data files are written in one DVD–RAM.
These two processes are controlled by the program named “xmas”. This is X–window systemapplication written in Tcl/Tk, and is graphical user interface to control all the process and to111

Figure 6.11: An example of the measurement for the linear relation between time differences
and TDC counts.
provide the monitoring functions for the acquired data and for the acquisition processes.

112

user

controlxmas

monitorstatuscommand/agreementcommand/agreement

Shared memory

statusstatuseco

datast

data

CAMAC

control

HDD

DVD
Figure 6.12: The block diagram of the data acquisition processes for the MAS array.

113

Chapter 7

Air shower analysis andcharacteristics of the array
7.1

Air shower analysis

7.1.1

Determination of arrival direction

The arrival directions of air showers are determined by the relative arrival times measured
through the TDCs and with the positions of the fast-timing detectors. Each TDC value isconverted to a relative time, and corrected with the calibrated delay time. When we assume
that the surface of a shower disk is a plane, each set of the detector position(x, y, z) and thearrival time of a shower front(t) should satisfy the following relation,lx + my + nz − ct = 0(7.1)

where c is the speed of light. (l, m, n) are the direction cosines of the shower axis, and have thefollowing relations with the zenith angle θ and the azimuth angle φ,l = sin θ cos φ(7.2)

m = sin θ sin φ(7.3)

n = cos θ(7.4)

The vector (l, m, n) is determined by minimizing the χ2 which is defined asXχ2 =wi [lxi + myi + nzi − c (ti − t0 )]2

(7.5)

i

where (xi , yi , zi ) is the position of the detector i and ti is an arrival time of a shower front at
the detector i. The weight wi is determined from the error for the timing measurements. The
equation to be solved is
∂χ2
∂χ2
∂χ2
=
=
=0
(7.6)
∂l
∂m
∂n
under the condition l2 + m2 + n2 = 1. In our analysis, this procedure is iterated by three times
with excluding the data sets which have the relative times separated more than 15 ns from theestimated shower plane in order to avoid the determination errors based on accidental hits on
the detectors by local muons.
114

The actual shower front has a corn–like structure, not a plane perpendicular to the shower
axis. When we require an improvement of the accuracy of arrival direction determinations, thecorrections of the corn structure should be applied for relative arrival times.
The corn structure of shower front is described with delay times of shower particles fromthe approximated shower plane as shown in Figure 7.1. The mean delay time is expressed by
the empirical formula obtained through the measurements of the SAS array [57] as follows,
Td = (−0.120 ± 0.164) − (1.334 ± 0.337) log r + (2.125 ± 0.156)(log r)2 ns
with the fluctuation³
σt = 0.777 1 +

r ´
,
30m

(7.7)
(7.8)

where r is the distance from a shower axis.

Figure 7.1: Delay time of shower particles from the shower plane measured at Mt. Chacaltaya
with SAS array. The line in the figure is a calculated delay time with the formula (7.7).
The accuracy in the determination with the corn structure is 20 % better than that withthe plane approximation. However, the calculation time of the cone structure fitting becomes
longer that of the plane fitting, for example, by about 10 times with one present calculation
code. Then we did not apply the procedure with the corn structure corrections for the following
analysis.

115

7.1.2

Determination of shower size

The shower core position and the shower size of an air shower are determined with minimizing
the following χ2 ,
´2
X ³
exp
χ2 =
wi ρobs
−
ρ
.
(7.9)
i
i
i

is expected and observed particle densities at the detector i, respectively. When weand ρobs
ρexp
ii
measured local densities at the same core distance with more than two detectors, the densities
for them. When theare averaged over these detectors and the averaged density is used as ρobs
i
particles in a locally high particle density region near a shower axis hit a detector i, the weight
wi for (7.9) becomes too large, and the determined core position is drawn to the detector. Since
the determined core positions are localized near the detector positions, the density averaging
procedure for the equidistant detectors are performed in order to avoid this localizations. For
this averaging calculation, the resolution of core distances is 1 m, because this is a typical sizeof the detectors. This procedure is efficient for the events with large air shower sizes.
An expected density, ρexp
is calculated with the following functions,iNef (ri , s),2rm¶µ ¶s−2 µri s−4.5ri
1+,
f (ri , s) = C(s)rmrm
Γ (4.5 − s).
C(s) =
2πΓ (s) Γ (4.5 − 2s)ρ(ri , s) =(7.10)(7.11)(7.12)

here r is the core distance and rm is the half of the Molière length. In our site at 5200 m a.s.l.,
the length is 155m [57], and thus rm = 77.5 m. Ne is a shower size, and it is estimated fromthe observed particle densities by
P 0 obs 2
wi · ρi · rm /f (ri )
P 0
Ne =
.
(7.13)wi
The weights wi of (7.9) and wi0 are defined wi = 1/σi2 and wi0 = f 2 (ri )/σi2 with the standard
deviation σi of ρexp
i . From the detailed analysis of the observed air shower particles with the
SAS array, we found that the standard deviation σ is expressed well with the following empiricalformula [57],σ = (1.25 ± 0.68)ρexp + (7.08 ± 0.24)(ρexp )2 .
(7.14)The observed standard deviation as a function of expected density is plotted in Figure 7.2 withthe calculated curves by (7.14).
We used a simple grid–search method to find a minimum chi–square solution. This calculations are performed in orthogonal co–ordinate systems with shower axes as the z co–ordinates.
1. The density-weighted average of detector positions is chosen as a starting values of agrid–search, (x0 , y0 ).
2. Place a square with sides of L meters, and the starting point(x0 , y0 ) is centered in thesquare. The initial length of L is 120 m.
3. Search the point(xm , ym ) minimizing χ2 among 7 × 7 = 49 grid points in the square.
4. Replace L and (x0 , y0 ) with L/6 and (xm , ym ), respectively. Then repeat the step 2–4
until L < 1 m.
116

Figure 7.2: The fluctuation (σρ ) of shower particle density ρ. σρ is independent of air shower
size.

7.2

Performance of MAS array

7.2.1

Effective area

The effective area of the array is estimated with the triggering efficiency, which is defined as,η(r, E) =number of triggered events,
number of simulated events(7.15)

for air showers of a primary energy E and the core distance of r ∼ r + dr. So that, the effective
area Seff (E) is given by,Z∞

Seff (E) =2πrη(r, E)dr.

(7.16)

0

Seff (E) for the array shown in Figure 7.3 [99] was calculated for the events of the zenith angle
θ < 60◦ with the simulation code developed by Shirasaki [83] and Tsunesada [93].
The expected energy distribution of observed cosmic rays with the array is calculated with
Seff (E) and primary energy spectrum. When we assume a prevalent spectrum with the index
of 2.7, the expected energy distribution of cosmic rays is obtained with the following formula,µ¶µ¶
Z E+dEdN
dN
dE =Seff (E)dE(7.17)
dE triggereddE primaryE
117

Figure 7.3: The effective area when MAS array is under the trigger condition described in theprevious section.
as shown in Figure 7.4. As the result, the mode energy is found about 30 TeV. The integrationof the figure over the whole energy gives the expected trigger frequency. This is about 10 Hz
and is consistent with the measurements.

7.2.2

Resolutions

In order to evaluate the resolutions and the systematic errors of the shower sizes and of thearrival directions, we used an air shower simulation program, CORSIKA [41]. It is a detailedMonte Carlo simulation program to study the evolution and the properties of extensive air
showers in the atmosphere. It was developed for the KASCADE experiment, and is one of thestandard simulation code that is used by many experiments all over the world. We developed
the program called “cor2mas “ to convert the outputs of CORSIKA to simulated data whichhave the identical format with the measured data of the MAS array. On the translations by
this program, we took account of the responses of the detectors and the electronics, the cabledelays, and the fluctuations (of number of incident particles, arrival times and delays). The
simulated data were processed with the completely identical analysis programs which are used
for the actual observed data.
The total number of simulated events with CORSIKA no–thinning mode is 3 × 106 forprimary protons and same number for primary irons. The simulation calculations covered
primary energies from 1012 to 6 × 1015 eV, and core locations from 0 to 40 m from the center of
the array. For the simulated data we applied the identical data selections with those for actual
measured data, and accordingly the selected events have the reduced chi–square values for thedeterminations of EAS arrival directions less than 3 and the determined core position within
118

Figure 7.4: The expected energy distribution of the triggered events. The mode energy is about
30 TeV.
20 m from the center of the array.
Angular resolutions
We examine the accuracy in determination of shower arrival directions by comparing both
directions of the incident particles and of the corresponding reconstructed showers. Figure 7.5
and 7.6 are the distributions of the opening angle ∆θ for the several primary energies for the
protons whose incident zenith angles are 0◦ and 44.4◦ . These figures are the histograms for thesimulated events which have the core positions(r) less than 20 m from the center of the array.
For the events outside of this circle, the angular resolution are about 50 % large as shown in
Figure 7.7. Thus in the following descriptions of the resolutions we apply the data selections
with r < 20 m. The distribution is well expressed with the following formula;
µ
¶
dN
1
∆θ2
= ∆θ exp − 2(7.18)d∆θσθ
2σθ
Figure 7.8 shows an example of the fitting of a opening angle distribution with this formula.
We take the σθ in the formula (7.18) as the angular resolution of the array, and plot σθ as afunction of energies in Figure 7.9. For the vertically incident primary particles, the angularresolution is 1.4◦ for EASs with primary energy of 1014 eV, and is 0.8◦ for those of 1015 eV.

119

60
200
40
100

20
0
0

2
4
6
8
opening angle ∆θ(degree)

0
0

10

2
4
6
8
opening angle ∆θ(degree)

10

2
4
6
8
opening angle ∆θ(degree)

10

2
4
6
8
opening angle ∆θ(degree)

10

80
100

60
40

50
20
0
0

2
4
6
8
opening angle ∆θ(degree)

0
0

10

60

40
30

40

20
20

10
0
0

2
4
6
8
opening angle ∆θ(degree)

0
0

10

Figure 7.5: Distributions of ∆θ of vertically incident primaries. All the events are selected withthe condition r ≤ 20 m.

60
6
40

4

20

2
0
0

2
4
6
8
opening angle ∆θ(degree)

0
0

10

2
4
6
8
opening angle ∆θ(degree)

10

2
4
6
8
opening angle ∆θ(degree)

10

2
4
6
8
opening angle ∆θ(degree)

10

80
40
60

30

40

20

20

10

0
0

2
4
6
8
opening angle ∆θ(degree)

0
0

10

30
15
20
10
10

5
0
0

2
4
6
8
opening angle ∆θ(degree)

10

0
0

Figure 7.6: Distributions of ∆θ of the simulated showers of θ = 44.4◦ .

120

250
40

200

30

150

20

100

10

50

0
0

2
4
6
8
opening angle ∆θ(degree)

0
0

10

2
4
6
8
opening angle ∆θ(degree)

10

2
4
6
8
opening angle ∆θ(degree)

10

2
4
6
8
opening angle ∆θ(degree)

10

300
200
200
100

100

0
0

2
4
6
8
opening angle ∆θ(degree)

0
0

10

200
100

150

75
100

50

50

25
0
0

2
4
6
8
opening angle ∆θ(degree)

10

0
0

Figure 7.7: Distributions of ∆θ for vertically incident primaries. All the events are selected withthe condition r > 20 m. The angular determination errors of the events on the outer region of
the array is large by about 50 % compared to Figure 7.5.

90

80

70

60

50

40

30

20

10

0
0

1

2

3

4
5
6
opening angle ∆θ(degree)

7

8

9

10

Figure 7.8: An examples of the fitting ∆θ distribution with the formula (7.18). σθ is determined
to be 1.94◦ for this distribution.
121

angular resolution(degree)

12

10

8

6

4

2

0
10

10

2

energy(TeV)

Figure 7.9: The plot of angular determination errors σθ as a function of primary energy.

122Size determination errors
The core locations and the sizes of air showers are calculated with the grid–search proceduredescribed in the previous section. The core position differences, the size determination error∆ log(Nrec /Nsim ), the core distance resolutions are shown in Figure 7.10–7.15. Figure 7.15
shows that the EAS size resolution is 35% for vertically incident EASs with primary energy of1014 eV, and is 20% for 1015 eV. Also, it is shown there are about 30% systematic differences
between the input and the determined sizes for the primaries of E > 30 TeV.

200

60

150

40

100
20
0
0

50
20
40
core position diff. ∆R(m)

60

0
0

60

20
40
core position diff. ∆R(m)

60

20
40
core position diff. ∆R(m)

60

20
40
core position diff. ∆R(m)

60

100

40
50
20
0
0

20
40
core position diff. ∆R(m)

0
0

60

40

80

30

60

20

40

10

20

0
0

20
40
core position diff. ∆R(m)

60

0
0

Figure 7.10: Distributions of the core position difference for the vertically incident primaries.
All the events are selected with the condition r ≤ 20 m.

123

60
10
40
5

20

0
0

20
40
core position diff. ∆R(m)

0
0

60

20
40
core position diff. ∆R(m)

60

20
40
core position diff. ∆R(m)

60

20
40
core position diff. ∆R(m)

60

40

60

30

40

20
20

10

0
0

20
40
core position diff. ∆R(m)

0
0

60

15

20

10

15
10

5

5

0
0

20
40
core position diff. ∆R(m)

0
0

60

Figure 7.11: Distributions of the core position difference for the simulated showers of θ = 44.4◦ .
All the events are selected with the condition r ≤ 20 m.

40

200

20

100

0

0

-2
-1
0
1
2
size determination error ∆log(Nrec/Nsim)

-2
-1
0
1
2
size determination error ∆log(Nrec/Nsim)

100
150

75
50

100

25

50

0

0

-2
-1
0
1
2
size determination error ∆log(Nrec/Nsim)

-2
-1
0
1
2
size determination error ∆log(Nrec/Nsim)

60
40
40
20

0

20

-2
-1
0
1
2
size determination error ∆log(Nrec/Nsim)

0

-2
-1
0
1
2
size determination error ∆log(Nrec/Nsim)

Figure 7.12: Distributions of the size determination error for the vertically incident primaries.
All the events are selected with the condition r ≤ 20 m.
124

10


8

60

6
40
4
20

2
0

0

-2
-1
0
1
2
size determination error ∆log(Nrec/Nsim)

-2
-1
0
1
2
size determination error ∆log(Nrec/Nsim)



80

60

60
40
40
20

20
0

0

-2
-1
0
1
2
size determination error ∆log(Nrec/Nsim)

-2
-1
0
1
2
size determination error ∆log(Nrec/Nsim)

30
20
20

15
10

10

5
0

-2
-1
0
1
2
size determination error ∆log(Nrec/Nsim)

0

-2
-1
0
1
2
size determination error ∆log(Nrec/Nsim)

core distance resolution(m)Figure 7.13: Distributions of the size determination error for the simulated showers of θ = 44.4◦ .
All the events are selected with the condition r ≤ 20 m.

20

15

10

5

0

-5

-10

-15

-20
10

10

2

10

3

energy(TeV)

Figure 7.14: The plot of the core distance resolution as a function of primary energy. The
markers indicate the systematic errors, and the error bar of each plot is the statistical error.

125

size determination error ∆log(Nrec/Nsim)

2

1.5

1

0.5

0

-0.5

-1

-1.5

-2
10

10

2

10

3

energy(TeV)

Figure 7.15: The plot of the size determination error as a function of primary energy. The
markers indicate the systematic errors, and the error bar of each plot is the statistical error.

126

Chapter 8

Experimental Results
In order to estimate the chemical composition of cosmic rays, we measured EAS longitudinaldevelopment with the equi–intensity method analysis. In the analysis, we used the angulardependences of EAS sizes, the rate of arrivals of EASs and the resultant size distributions. Thus
the equi–intensity method analysis is based only on air shower size spectra. The simulations andthe results of the chemical composition and the energy spectrum are presented in this chapter.

8.1

Data selection

The details of the analysis is described in the previous chapter. The observed data were analyzed with the following procedure: (1)The arrival direction of each EAS is determined with aleast square fit, assuming that the EAS front is flat, (2)The positions of the unshielded detectors are projected onto the EAS front plane, (3)The core location, the EAS size, and the age
parameter are determined on this plane by minimizing the chi–square. Through this analysis,
we obtained a set of the parameters characterizing each EAS event, i.e., an EAS size(Nobs ), an
age parameter(s), a zenith angle(θ), an azimuthal angle(φ), and a core location.
The data used for the present analysis were collected between March and November of
2000. The total observation time was 8.9 × 106 s. During this period, the array was in nearlycontinuous operation. There were, however, short interruptions for the detector calibrations
and the occasional power failures. Since these interruptions can give rise to non–uniformities
on the arrival direction distributions in the observed data, we checked the continuity of theobservation in each sidereal day. The data in the sidereal days which had a time gap of twoconsecutive events larger than 8 minutes (2◦ in hour angle) were excluded from the analysis.
In order to choose EAS events whose arrival directions and air shower sizes are determinedin good accuracy, we apply the selections of EAS events with the reduced chi–square values
for the determinations of EAS arrival directions(χ2dir ) and with that for determinations of EAS
sizes(χ2size ) to the data. Figure 8.1 and 8.2 are distributions of the reduced chi–square values,
χ2dir and χ2size for a 20000 event sample. We selected the events with χ2dir less than 3. However,χ2size have a size dependence shown in Figure 8.3, and then we did not select with χ2size in order
to avoid the systematic selection bias.
We selected the events of which the core positions were within distances less than 20 m fromthe center of the array. As shown in Figure 8.4, the events fallen outside of the region havelarge determination errors of arrival directions and those of core positions, and consequently,
have large size determination errors.
127

6000

5000

4000

3000

2000

1000

0

0

1

2

3

4

5

χdir2

6

7

8

9

10

Figure 8.1: Distribution of the reduced chi–square for the determinations of the arrival directions, (7.5).
Moreover, we selected the events of which zenith angles are less than 60◦ . As the result, thetotal observation time is 8.9 × 106 s and 7.5 × 107 events are selected.

128

3000

2500

2000

1500

1000

500

0

0

1

2

3

4

5

χsize2

6

7

8

9

10

χsize2

Figure 8.2: Distribution of the reduced chi–square for the determinations of the EAS sizes,(7.9).

10

9

8

7

6

5

4

3

2

1

0

3

4

5

6

7

8

log(N)

Figure 8.3: Scatter plots of χ2size vs EAS size.

129

100
200

75

100 TeV
1 PeV

50

100

25
0



0

20

40

60

80

100

0



0

20

40

|r-r0|

60

80

100

|r-r0|



80
200
60

150

40

100

20
0

50
0

5

10

15

0

0

5

∆θ
150

300

100

200

50

100

0

-2

-1

10

15

∆θ

0

1

2

0



-2

-1

0

1

2



log10(N/N0)

log10(N/N0)

Figure 8.4: Distributions: (first column) of the determination errors of core positions(|r − r0 |),(second column) of arrival directions(∆θ) and (third column) of EAS sizes(log(N/N0 )) for thefull Monte Carlo simulated events with CORSIKA-QGSJET.

130

8.2

Simulations

In the analysis of the experimental data, in order to estimate the chemical composition ofcosmic rays with comparing between the observed and the simulated size distributions, thenumber of simulated events must be at least the same as that of the collected data. However,making full Monte-Carlo simulations for such a large data set is a huge task. Therefore, we
simulated events with the size–energy relations, the energy dependent size fluctuations, theangular resolutions, and the detection efficiencies. These relations were obtained with the fullMonte–Carlo simulations of CORSIKA [41] of no–thinning mode with the QGSJET hadroninteraction model. The Monte Carlo program CORSIKA is the most used EAS simulation
program, and simulates the 4–dimensional evolution of EASs in the atmosphere initiated byphotons, hadrons and nuclei. It has the full Monte Carlo simulation mode called ’no–thinningmode’. In the full Monte-Carlo simulations, we took account of the characteristics of the MAS
array detectors including the pulse height distributions, the time delay and the fluctuations forthe minimum ionizing particles. The simulated data were analyzed with the identical procedure
with that for the real observed data.
We generated EAS events with CORSIKA for the EAS size from 103 to 107 , and the events
had the uniform core position distribution around the center of the array within 40 m radius.
The number of simulated EAS with CORSIKA is about 6 × 106 .

8.2.1

Parameterizations

The simulated data were selected with the same criteria as that for the measured data, andused to obtain (1)effective area, (2)angular resolution, (3)shower size and (4)the fluctuation of
shower size for fixed primary energy as parameterized functions of the primary energy and sec θ.
effective area
The plots of the effective detection areas as a function of the primary energy for each sec θ areshown in Figure 8.5 for primary protons and Figure 8.6 for primary irons. These plots werefitted with the following empirical formula,
µ
¶
10S2
(8.1)
log S = S1 − S3 log 1.0 + 2 .
EHere, S1 , S2 and S3 are the fitting parameters and E is a primary energy in TeV. The parameters
are plotted in Figure 8.7 and 8.8 as functions of sec θ. These plots were fitted with polynomial
functions of sec θ, and then the results for primary protons are
S1 = 2.6897 + 0.64697 sec θ − 0.25758(sec θ)2

(8.2)
2

S2 = −0.30813 + 2.3568 sec θ + 0.10936(sec θ)(8.3)

S3 = 2.6021 − 0.83939 sec θ,(8.4)

and for primary irons are
S1 = 2.8983 + 0.44159 sec θ − 0.20076(sec θ)2

(8.5)

2

S2 = 1.0792 + 1.5149 sec θ + 0.18284(sec θ)(8.6)

S3 = 7.0.

(8.7)
131

2

S[m ]

2

S[m ]

2
0
-2
-1

0

1

2

3

2
0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

S[m ]

2
0
-2
-1

0

1

2

3

2
0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

S[m ]

3

log10(E[TeV])



2
0
-2
-1

0

1

2

3

2
0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

3

log10(E[TeV])



S[m ]

3

log10(E[TeV])



2

0
-2
-1

0

1

2

3

2

0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

S[m ]

2
0
-2
-1

3

log10(E[TeV])



0

1

2

3

4

2
0
-2
-1

log10(E[TeV])


0

1

2

3

4

log10(E[TeV])

Figure 8.5: Effective detection area as a function of the primary energy for primary protons
simulated with CORSIKA. The solid lines are drawn with the function (8.1).
The parameters S1 , S2 and S3 well represent the effective area for the possible energy and sec θrange as shown in Figure 8.9 and 8.10.

132

2

S[m ]

2

S[m ]

2
0
-2
-1

0

1

2

3

2
0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

S[m ]

2
0
-2
-1

0

1

2

3

2
0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

S[m ]

3

log10(E[TeV])



2
0
-2
-1

0

1

2

3

2
0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

3

log10(E[TeV])



S[m ]

3

log10(E[TeV])



2

0
-2
-1

0

1

2

3

2

0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

S[m ]

3

log10(E[TeV])



2
0
-2
-1

0

1

2

3

4

2
0
-2
-1

0

1

2

log10(E[TeV])

34log10(E[TeV])





S1

Figure 8.6: Effective detection area as a function of the primary energy for primary irons
simulated with CORSIKA. The solid lines are drawn with the function (8.1).

3.2
3.1
3
2.9

1

1.2

1.4

1.6

1.8

2

S2

secθ
5
4
3
2
1

1.2

1.4

1.6

1.8

2

S3

secθ
2
1.5
1
0.5
0

1

1.2

1.4

1.6

1.8

2

secθ

Figure 8.7: The fitting parameters of the function (8.1) for the simulated primary protons. The
solid lines represent the polynomial functions of sec θ .
133

S1

3.2
3.1
3
2.9

1

1.2

1.4

1.6

1.8

2

S2

secθ
5
4
3
2
1

1.2

1.4

1.6

1.8

2

secθ
S3

8
7
6
5

1

1.2

1.4

1.6

1.8

2

secθ

2

S[m ]

2

S[m ]

Figure 8.8: The fitting parameters of function (8.1) for the simulated primary irons. The solidlines represent the polynomial functions of sec θ .

2
0
-2
-1

0

1

2

3

2
0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

S[m ]

2
0
-2
-1

0

1

2

3

2
0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

S[m ]

3

log10(E[TeV])



2
0
-2
-1

0

1

2

3

2
0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

3

log10(E[TeV])



S[m ]

3

log10(E[TeV])



2

0
-2
-1

0

1

2

3

2

0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

S[m ]

2
0
-2
-1

3

log10(E[TeV])



0

1

2

3

4

2
0
-2
-1

log10(E[TeV])


0

1

2

3

4

log10(E[TeV])

Figure 8.9: Effective detection area as a function of the primary energy for primary protons
simulated with CORSIKA. The solid lines are effective area calculated with the function (8.1)
and with the polynomial functions of S1 , S2 and S3 .
134

2

S[m ]

2

S[m ]

2
0
-2
-1

0

1

2

3

2
0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

S[m ]

2
0
-2
-1

0

1

2

3

2
0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

S[m ]

3

log10(E[TeV])



2
0
-2
-1

0

1

2

3

2
0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

3

log10(E[TeV])



S[m ]

3

log10(E[TeV])



2

0
-2
-1

0

1

2

3

2

0
-2
-1

4

0

1

2

log10(E[TeV])

4


2

S[m ]

2

S[m ]

2
0
-2
-1

3

log10(E[TeV])



0

1

2

3

4

2
0
-2
-1

log10(E[TeV])


0

1

2

3

4

log10(E[TeV])

Figure 8.10: Effective detection area as a function of the primary energy for primary irons
simulated with CORSIKA. The solid lines are effective area calculated with the function (8.1)
and with the polynomial functions of S1 , S2 and S3 .

135

angular resolutionσθσθThe plots of the angular resolutions as a function of primary energy for each sec θ are shown in
Figure 8.11 for primary protons and Figure 8.12 for primary irons. These plots are fitted witha linear function of sec θ, that is A × sec θ + B. The fitting parameters are shown in Figure
8.13. In the simulations, the angular resolutions are calculated with the natural cubic splineinterpolation of these plots.

4
2
0

4
2
0

3
4
log10(E[TeV])

1

2

3
4
log10(E[TeV])

1

2

3
4
log10(E[TeV])

1

2

3
4
log10(E[TeV])

1

2

3
4
log10(E[TeV])

1

2

3
4
log10(E[TeV])

σθ

2

σθ

1

4
2
0

4
2
0

3
4
log10(E[TeV])
σθ

2

σθ

1

4
2

4
2

0

0
3
4
log10(E[TeV])
σθ

2

σθ

1

4
2
0

4
2
0

3
4
log10(E[TeV])

4

σθ

2

σθ

1

2
0

2
0

4

1

2

3
4
log10(E[TeV])

Figure 8.11: Angular resolutions as a function of the primary energy for primary protons simulated with CORSIKA are fitted with a linear function.

136

σθ

σθ

4
2
0

4
2
0

3
4
log10(E[TeV])

1

2

3
4
log10(E[TeV])

1

2

3
4
log10(E[TeV])

1

2

3
4
log10(E[TeV])

1

2

3
4
log10(E[TeV])

1

2

3
4
log10(E[TeV])

σθ

2

σθ

1

4
2
0

4
2
0

3
4
log10(E[TeV])
σθ

2

σθ

1

4
2

4
2

0

0
3
4
log10(E[TeV])
σθ

2

σθ

1

4
2
0

4
2
0

3
4
log10(E[TeV])

4

σθ

2

σθ

1

2
0

2
0

4

1

2

3
4
log10(E[TeV])

Figure 8.12: Angular resolutions as a function of the primary energy for primary irons simulatedwith CORSIKA are fitted with a linear function.

7
-0.6

6.5

-0.8

6

-1

5.5
5

-1.2

4.5

-1.4

4
-1.6

3.5

-1.8

3

-2

2.5

1

1.25

1.5

1.75

2

1

1.25

1.5

1.75

2

1

1.25

1.5

1.75

2

7
-0.6

6.5

-0.8

6

-1

5.5
5

-1.2

4.5

-1.4

4
-1.6

3.5

-1.8

3

-2

2.5

1

1.25

1.5

1.75

2

Figure 8.13: The plots of the fitting parameters for the angular resolutions.

137

shower size

log10Nrec

log10Nrec

The plots of the reconstructed shower sizes(Nrec ) as a function of the primary energy for eachsec θ are shown in Figure 8.14 for primary protons and Figure 8.15 for primary irons. These

6

6

4

4

2
10

10

10
E[TeV]

4

2
1

10

10

1

10

10

1

10

10

1

10

10

1

10

10

2

10

3

10
E[TeV]

4

log10Nrec10

3log10Nrec

1

2

6

6

4

4

2
10

10

10
E[TeV]

4

2

2

10

3

10
E[TeV]

4

log10Nrec10

3log10Nrec

1

2

6

6

4

4

2
10

10

10
E[TeV]

4

2

2

10

3

10
E[TeV]

4

log10Nrec10

3log10Nrec

1

2

6

6

4

4

2
10

10

10
E[TeV]

4

2

2

10

3

10
E[TeV]

4

log10Nrec10

3log10Nrec

1

2

6

6

4

4

2
1

10

10

2

10

3

10
E[TeV]

4

2

2

10

3

10
E[TeV]

4

Figure 8.14: The plots of the reconstructed shower size as a function of the primary energy foreach sec θ for primary protons simulated with CORSIKA.
plots are fitted well with lines,
log Nrec = A + B log E(8.8)

here E is primary energy in TeV. The parameters of these lines are shown in Figure 8.16 andFigure 8.17, as a function of sec θ. These plots are fitted with polynomial functions of sec θ,and thus the results for primary protons are
A = 1.6989 + 0.90041 sec θ= 2.2141 − 0.10137 sec θ2= 5.9711 − 3.9656 sec θ + 0.74167(sec θ)
B = 1.5139 − 0.65675 sec θ= 1.5899 − 0.35788 sec θ2= −0.31205 + 1.8017 sec θ − 0.52547(sec θ)138(line1)(8.9)(line2)(8.10)(line3)(8.11)(line1)(8.12)(line2)(8.13)(line3)(8.14)

log10Nrec

log10Nrec

6

6

4

4

2
10

10

10
E[TeV]

4

2
1

10

10

1

10

10

1

10

10

1

10

10

1

10

10

2

10

3

10
E[TeV]

4

log10Nrec10

3log10Nrec

1

2

6

6

4

4

2
10

10

10
E[TeV]

4

2

2

10

3

10
E[TeV]

4

log10Nrec10

3log10Nrec

1

2

6

6

4

4

2
10

10

10
E[TeV]

4

2

2

10

3

10
E[TeV]

4

log10Nrec10

3log10Nrec

1

2

6

6

4

4

2
10

10

10
E[TeV]

4

2

2

10

3

10
E[TeV]

4

log10Nrec10

3log10Nrec

1

2

6

6

4

4

2
1

10

10

2

10

3

10
E[TeV]

4

2

2

10

3

10
E[TeV]

4

Figure 8.15: The plots of the reconstructed shower size as a function of the primary energy foreach sec θ for primary irons simulated with CORSIKA.
and for primary irons are
A = −0.56604 + 1.7134 sec θ= 2.3537 − 1.2348 sec θ2= 5.8099 − 4.9532 sec θ + 1.1001(sec θ)
B = 2.6532 − 1.2265 sec θ= 1.5409 − 0.058930 sec θ2= 0.19496 + 1.4247 sec θ − 0.44074(sec θ)139(line1)(8.15)(line2)(8.16)(line3)(8.17)(line1)(8.18)(line2)(8.19)(line3)(8.20)

B

A

4
3

2

1.5

2
1
1
0.5

0
1

1.25

1.5

1.75

0

2
secθ

4

B

A

-1

3

1

1.25

1.5

1.75

2
secθ

1

1.25

1.5

1.75

2
secθ

1

1.25

1.5

1.75

2
secθ

2

1.5

2
1
1
0.5

0
1

1.25

1.5

1.75

0

2
secθ

4

B

A

-1

3

2

1.5

2
1
1
0.5

0
-1

1

1.25

1.5

1.75

2
secθ

0

Figure 8.16: The fitting parameters of log10 Nrec –log10 E relations for primary protons, for low
energy(line 1), middle energy(line 2) and high energy region(line 3). These parameters are fittedwith a polynomial function of sec θ.

140

B

A

4
3

2

1.5

2
1
1
0.5

0
1

1.25

1.5

1.75

0

2
secθ

4

B

A

-1

3

1

1.25

1.5

1.75

2
secθ

1

1.25

1.5

1.75

2
secθ

1

1.25

1.5

1.75

2
secθ

2

1.5

2
1
1
0.5

0
1

1.25

1.5

1.75

0

2
secθ

4

B

A

-1

3

2

1.5

2
1
1
0.5

0
-1

1

1.25

1.5

1.75

2
secθ

0

Figure 8.17: The fitting parameters of log10 Nrec –log10 E relations for primary protons, for low
energy(line 1), middle energy(line 2) and high energy region(line 3). These parameters are fittedwith a polynomial function of sec θ.

141

The size fluctuations
The size fluctuations(σN ) for fixed primary energies are obtained as the standard deviationsof the size distributions. σN are plotted as a function of the primary energy for each sec θ in
Figure 8.18 and 8.19. These plots are fitted well with the following linear relation as,

σN

1

σN

1
0

0

-1
0

-1
0

1

2

2

3
4
log10[E(TeV)]

1

2

3
4
log10[E(TeV)]

1

2

3
4
log10[E(TeV)]

1

2

3
4
log10[E(TeV)]

1

2

3
4
log10[E(TeV)]

σN

0
-1
0

1

1

σN

1

3
4
log10[E(TeV)]

0
1

2

-1
0

σN

1

σN

1

3
4
log10[E(TeV)]

0

0

-1
0

-1
0

1

2

σN

1

σN

1

3
4
log10[E(TeV)]

0

0

-1
0

-1
0

1

2

σN

1

σN

1

3
4
log10[E(TeV)]

0

0

-1
0

-1
0

1

2

3
4
log10[E(TeV)]

Figure 8.18: Size resolutions as a function of the primary energy for primary protons simulatedwith CORSIKA are fitted with a linear function.
σN = A + B log E(8.21)

here E is the primary energy in TeV. The determined values for A and B are shown in Figure
8.20, and fitted with the linear function of sec θ. The results for primary protons are
A = −0.023634 + 0.32377 sec θ(8.22)

B = −0.045893 − 0.029957 sec θ,(8.23)

A = −0.17884 + 0.42024 sec θ(8.24)

B = 0.029940 − 0.089345 sec θ.

(8.25)

and for primary irons,

142

σN

1

σN

1
0

0

-1
0

-1
0

1

2

3
4
log10[E(TeV)]234
log10[E(TeV)]

1

2

3
4
log10[E(TeV)]

1

2

3
4
log10[E(TeV)]

1

2

3
4
log10[E(TeV)]

1

2

3
4
log10[E(TeV)]

σN

1

σN

1

1

0

0

-1
0

1

2

3
4
log10[E(TeV)]

σN

1

σN

1

-1
0

0

0

-1
0

-1
0

1

2

3
4
log10[E(TeV)]

σN

1

σN

1
0

0

-1
0

-1
0

1

2

3
4
log10[E(TeV)]

σN

1

σN

1
0

0

-1
0

-1
0

1

2

3
4
log10[E(TeV)]

0.7

-0.05

0.6

-0.1

0.5

-0.15

A

0.8

B

Figure 8.19: Size resolutions as a function of the primary energy for primary protons simulatedwith CORSIKA are fitted with a linear function.

0

0.4

-0.2

0.3

-0.25

0.2

-0.3
1

1.25

1.5

1.75

2

1

1.25

1.5

1.75

0.7

-0.05

0.6

-0.1

0.5

-0.15

0.4

-0.2

0.3

-0.25

A

0.8

0.2

2

secθ
B

secθ
0

-0.3
1

1.25

1.5

1.75

2

secθ

1

1.25

1.5

1.75

2

secθ

Figure 8.20: The fitting parameters of σN –log10 E relations for primary protons and irons.
These values are fitted with a linear function of sec θ.
143

8.2.2

Simulation scheme

The simulation procedure for the EAS events observed with MAS array has a following scheme.
1. Determine the primary energy randomly. However, the frequency of the occurrence is
proportional to the cosmic ray energy spectrum.
2. The primary particle type is determined. The program chooses randomly between a proton
and an iron nucleus. The probability is proportional to the assumption of the averaged
mass number of primary particles.
3. The incident direction of the air shower is determined randomly.
4. Calculate a detection efficiency with the primary energy and with the incident direction.
The event which satisfies the event selection conditions is randomly chosen according to
the calculated detection efficiency.
5. Determine the air shower size and its fluctuation with (8.8) and (8.21).
6. Output the shower parameters and go back to the top of the loop.
The data format of the output of the simulation program is identical to the output of
the actual analysis program for MAS array, and then the both of observed and simulated
shower parameters,i.e., the air shower size, arrival directions, etc, are processed with the sameprocedures of the equi–intensity method analysis described in the next section.
In Figure 8.21 and 8.22, I show the simulated EAS sizes for primary protons compared withthose simulated with CORSIKA no–thinning mode for several primary energies and incident
zenith angles. The calculations with the fitted parameters well reproduce the size distributions
with the full Monte Carlo simulations.

144

90
100
80

80

12

70

10

60
8

50

60

6

40
40

30

4

20

20

2

10
0

4

5

0

4

5

log10Nrec

40

0

18

14

16

30

14

25

12

10

10

8

12

8

15

5
log10Nrec

35

20

4

log10Nrec

6

6
10
5
0

4

4

2

2
5

5.5

0

6

5

5.5

log10Nrec

6

0

5

5.5

log10Nrec

6
log10Nrec

Figure 8.21: The distributions of the simulated sizes for primary protons with the fixed several
primary energies and incident zenith angles. The solid lines are those simulated with CORSIKA–QGSJET no–thinning mode, and the dotted lines are those simulated with the parameterizedfunctions derived in the text.

40

180

45

160

35

140

30

120

25

30

20

25

40

100

35

80

20

15
60

15
10

40

10

5

20
0

5

6

0

5
5

5.5

log10Nrec

6

5.5

6

6.5
log10Nrec

25

18
50

0

log10Nrec

16
20

14
40

12

15

10

30

8
20

10

6
4

10

5

2
0

4

0

5

5

5.5

log10Nrec

6
log10Nrec

0

5.5

6

6.5
log10Nrec

Figure 8.22: The distributions of the simulated sizes for primary protons with the fixed several
primary energies and incident zenith angles. The solid lines are those simulated with CORSIKA–QGSJET no–thinning mode, and the dotted lines are those simulated with the parameterizedfunctions derived in the text.
145

8.3

Size Spectra and Equi–Intensity Curves

eventseventsThe observed EASs were classified into two dimensional bins of EAS size(N obs ) and sec θ (θ:zenith angle of arrival direction). The size of each bin is 0.2 in logarithm of N obs , and 0.1 insec θ. We then obtained the size distributions for each sec θ interval. It it possible that there
is some deformation in these distributions from the true distributions due to the uncertainties
in size and in θ. Therefore, the correction factors for the number of events in each bin wereobtained with the EAS simulations.
Moreover the non–uniformity of azimuth angle distributions must be taken into consideration. Since the array is located on the side of the mountain, the array slopes slightly down tothe valley. Thus the detection efficiency of the array is not uniform in azimuth directions. The
azimuth angle distributions of each sec θ bin is shown in Figure 8.23. These non–uniformities are
due to the geometry of the array and are confirmed to be independent of EAS sizes. Thus, thecorrection factors were calculated from the observed azimuth angle distributions partly shownin Figure 8.23 which bin sizes are 0.1 in sec θ and 3◦ in azimuth angle, and the number of eventswere normalized to that in a maximum efficient direction.

5

10

4
3

10

2

10

5

10

3

10



10

4

10

2

10

10

10

1

1
0

100

200

300

0

100

200



10

5

10

4

10

3
2

10



10

5

10

4

10

3
2

10

10

10

1

1
0

100

200

300

0

100

200

azimuth[degree]


10

6

10

eventsevents



5



10

4

10

5
4

10

3

10

3

10

2

10

2

10

300

azimuth[degree]


6

10
10

300

azimuth[degree]
eventseventsazimuth[degree]

10

1

1
0

100

200

300

0

azimuth[degree]

100

200

300

azimuth[degree]

Figure 8.23: The azimuth angle distributions of particular sec θ bins.
After the corrections for measured size distributions, finally, we obtained the integral N obs
spectra for each sec θ bin, as shown in Figure 8.24. We selected EAS with N obs > 104.5 forfurther analysis, because our simulation shows that the detection efficiencies for N obs > 104.5
are approximately 100% and independent of an adopted primary component.
A differential size spectra dF/d log N obs for sec θ = 1.0−1.1 is shown in Figure 8.25, in which
the intensity have been multiplied by N 1.5 to flatten the spectrum and enhance its steepening.
The size spectrum is well represented by the power–law, and steepens at around N obs = 106.0 ,and there is not any other bending points. Thus the measured size spectrum is fitted with the
146

-3

-2

-1 -1

F(>N) (m sr s )

-2

-4

-5

-6

-7

-8

-9

-10

-11

-12

3.5

4

4.5

5

5.5

6

6.5

7

7.5

8

shower size N
Figure 8.24: The observed integral EAS size spectra for each sec θ bin.
following relation,dF
log([m−2 s−1 sr−1 ]) =d log N½

(3.85 ± 0.03) − (1.57 ± 0.01) × log N(6.80 ± 0.87) − (2.07 ± 0.14) × log N(log N < 5.9)(log N ≥ 5.9)

(8.26)

where N is equal to N obs .
If it is assumed that air showers initiated by primary cosmic rays with a certain energy arriveat the earth isotropically with a constant rate, the air shower size plots with the same flux as afunction of θ, i.e., the atmospheric depth represents the corresponding longitudinal development
curve. For this equi–intensity method analysis, we classified again the observed data into thebins of sec θ of which bin size = 0.05 in order to draw the details of longitudinal development.
We used the events only with sec θ < 1.5 because the angular determination errors becomelarger than the bin size∆ sec θ = 0.5 for sec θ ≥ 1.5. After correcting the size distributions again
with same procedures as described before, we derived the equi–intensity curves for air showers
with the integral of rate F (> N obs ) from 10−5.2 to 10−8.2 m−2 sr−1 s−1 for 10−0.25 steps. The
results are shown in Figure 8.26.
Please note that one cannot compare the curves in Figure 8.26(a) with the longitudinaldevelopment curves immediately. The measured EAS size with the air shower array is nottranslated directly into the number of electrons and positrons in the air shower, in part because
the detected number of particles is contaminated by muons and in part because the determined
size is affected by the threshold energy of the detectors. Moreover the equi–intensity curves areaffected by the energy-size relations, and by the energy spectra of components, as described147

N1.5dF/dlogN (m-2sr-1s-1)

10

4

10

3

4

4.5

5

5.5

6

6.5

7

7.5

shower size N

Figure 8.25: The differential size spectrum for sec θ = 1.0 − 1.1 measured by BASJE–MAS
array. The intensity has been multiplied by N 1.5 to flatten the spectrum.
above. The maximum development points show disagreement with the Xmax of the longitudinal development curves simulated with CORSIKA/QGSJET [62]. For example, for 4 PeV one
obtains Xmax ∼ 590g/cm2 for primary protons with CORSIKA/QGSJET, but at corresponding flux, 10−7.2 m−2 sr−1 s−1 , the maximum point of the equi–intensity curve is ∼ 660g/cm2 forprimary protons as shown in Figure 8.26(a).
In Figure 8.26(a), we show equi–intensity curves obtained with the same analysis proceduresfor the simulated data as were applied to the real observed events. The simulation was performed
with four different models in which the composition of the primaries was varied: pure protons,
pure iron nuclei, and two different mixtures of protons and iron nuclei. The air shower sizes
at certain fluxes depend on the adopted all particle spectrum and on the adopted mixing ratio
of protons and iron nuclei. Only when we simulate air showers with the actual all particle
energy spectrum and with the real mixing ratio, the equi–intensity curves will coincide with themeasured ones. Otherwise, a measured curve will not agree with a simulated one for a samerate of arrival.
The knowledge of the cosmic ray chemical composition is essential to the determination of
the energy spectrum from a measured air shower size spectrum as in the case of our analysis. On
the other hand, unless an adopted energy spectrum for the simulations is totally different fromthe actual one, the correlations between a mixing ratio and the shape of an equi–intensity curve
are less dependent on an adopted energy spectrum. Therefore, as the first step of determinationof the mixing ratios, we compared shapes of the observed and the corresponding simulated equi–148

shower size N

7

7

6.5

6.5

6

6

5.5

5.5

5

5

400

500

600

700


800
2


depth(g/cm )

400

500

600

700


800
2


depth(g/cm )Figure 8.26: (a)The observed equi–intensity curves are compared with simulated ones for primary protons, iron nuclei and (50% protons + 50% irons), without any normalization. The
attached value for each curve is the integral flux. In the analysis, we derived the curves for
F (> N ) from 10−5.2 to 10−8.2 m−2 sr−1 s−1 with 10−0.25 step, but in this figure we show a partof these curves to avoid confusions. (b)Simulated equi–intensity curves with the measured ones
without any normalization. In the simulation, we took consideration of not only the obtained
mixing ratios of protons and irons, but also obtained all particle energy spectrum.
intensity curves. Next, we calculated the relation between primary energies and air shower sizes
with the Monte Carlo simulations under the requirement of the obtained mixing ratios. The
observed size spectrum is then converted to the all particle energy spectrum. Using the resultant
energy spectrum, we examined the mixing ratios and the energy spectrum again, iteratively. In
total, we repeated this process four times to converge with the final results shown below.
The slope of equi–intensity curves becomes flatter with increasing proton ratio. However,as shown in Figure 8.27, the flattening is not linearly proportional to proton ratio. So that, inorder to determine a mixing ratio, the measured equi–intensity curves are compared with mixed
composition simulated ones. This is described in more detail in the next section.
The simulated curves of primary protons with an integral rate less than 10−6.2 m−2 sr−1 s−1
definitely show that the depths of shower maximum points should be deeper than the atmospheric depth of our observation site. However, in our results, these maximum points are notseen. Therefore, the major component of primary cosmic rays is considered to be heavier thanproton even in this energy range.

149

shower size N

6.6

6.4

6.2

6

5.8

5.6

400

450

500

550

600

650

700

750

800


2


depth(g/cm )Figure 8.27: The observed equi–intensity curves are compared with simulated ones. The attached value for each curve is the integral flux. The simulated curves are normalized at 578
g/cm2 with the measured curves.

150

8.4

Chemical Composition and Energy Spectrum of Primary
Cosmic Rays

8.4.1

Chemical Composition

We determined the primary chemical composition by comparing the observed equi–intesity
curves with the calculated ones. In our calculation, we assumed that the primary chemicalcomposition is a mixture of protons and iron nuclei, and that both the energy spectra of primary
protons and those of iron nuclei are proportional to E −γ with γ = 2.67 up to E = 3.9 × 1015 eV
and with γ = 3.41 above this energy [59][98]. We examined not only the primary composition
but also all particle energy spectrum. We simulated EASs with varying the proportion of protons
to the total (sum of protons and iron nuclei, (np /(np + nFe )). The proportion was varied from
0.0 to 1.0 with steps of 0.1.
The longitudinal development curves with the simulated data sets are obtained for each
proton mixing ratio. Comparing the observed curves with the simulated ones, we estimated themixing ratios by the least squares method with the following observation equation,vi = l(zi , x0 ) − lobs (zi ) +∂l(zi , x0 )∆x∂x

(8.27)

Here lobs (zi ) is an observed equi–intensity curve and l(zi , x0 ) is a calculated equi–intensity curvefor the mixture x0 , zi is the atmospheric depth for i-th sec θ bin. The best estimated x0 is determined with a minimum chi–squares method and the correction value of ∆x is determined withthe least squares method. Notice that each primary energy listed in Table 8.4.1 is determinedfrom the integral flux, and not from the EAS size, because the size–energy relations are strongly
depend on the chemical composition of primary cosmic rays.
Here we used the following energy spectrum measured with SAS array [81], 6.44 · 1018 × E −2.62 (14.1 < log E < 14.6)
6.91 · 1020 × E −2.76 (14.6 < log E < 15.7)(dF/dE)SAS =(8.28)
29−3.30
2.31 · 10 × E(15.7 < log E)
Moreover, for primary protons and iron nuclei, the primary energies and the energy bin widths
corresponding to an identical EAS size and a identical size bin width are different. So that,the mixing ratios as a function of primary energy were corrected by considering the energy
spectrum, and the energy–size relations of primary protons and iron nuclei. Finally, the resultsof the estimation of the mixing ratios x0 + ∆x are listed in Table 8.1.
The errors in the table are the fitting errors of 68% confidence level and the statisticalerrors. The systematic and statistical errors for the the values of np /(np + nFe ) are estimatedas the result of the identical analysis of simulated data. For the estimation, we simulated 50
independent data sets, each of which has the identical number of events to the observed data.
The differential flux of the simulated data has the identical energy dependences to the measured
all–particle spectrum described below, and the proton mixing ratios of the simulated data wereequal to those of the measured ones. The systematic and the statistical errors for the equi–intensity curves with the several different fluxes and their mixing ratio dependences are shownin Figure 8.28.
The systematic errors were estimated from the primary energies of protons and of iron nuclei
to contribute to each flux. A difference between the estimated energy of protons and that ofiron nuclei gives a systematic error of an estimated energy. We obtained the systematic and151log F (> N )
−5.20
−5.45
−5.70
−5.95
−6.20
−6.45
−6.70

log E[eV]14.55 ± 0.10
14.70 ± 0.09
14.85 ± 0.07
15.00 ± 0.06
15.15 ± 0.05
15.29 ± 0.04
15.42 ± 0.02

log N
5.41
5.56
5.71
5.86
6.00
6.13
6.26

np /(np + nFe )0.444 ± 0.029
0.369 ± 0.021
0.334 ± 0.023
0.260 ± 0.026
0.291 ± 0.024
0.286 ± 0.085
0.259 ± 0.027

log F (> N )
−6.95
−7.20
−7.45
−7.70
−7.95
−8.20

log E[eV]15.54 ± 0.01
15.65 ± 0.00
15.76 ± 0.01
15.87 ± 0.02
15.98 ± 0.05
16.08 ± 0.04

log N
6.38
6.49
6.59
6.70
6.87
6.98

np /(np + nFe )0.284 ± 0.039
0.312 ± 0.034
0.220 ± 0.080
0.148 ± 0.080
0.077 ± 0.110
0.095 ± 0.101

Table 8.1: The integral fluxes, corresponding primary energies, the average EAS sizes for sec θ =1.0−1.05 and estimated values of np /(np + nFe ) for each energy. The unit of integral flux F (> N )is m−2 sr−1 s−1 .
log F (> N )
−5.20
−5.45
−5.70
−5.95
−6.20
−6.45
−6.70

systematic bias
+0.0233+0.0011-0.0019-0.0121-0.0100+0.0066+0.0231

statistical error±0.0102±0.0128±0.0172±0.0094±0.0115±0.0158±0.0200log F (> N )
−6.95
−7.20
−7.45
−7.70
−7.95
−8.20

systematic error+0.0028-0.0048-0.0186-0.0502+0.0473+0.1451

statistical error±0.0339±0.0206±0.0702±0.0652±0.1052±0.0971

Table 8.2: The systematic errors of np /(np + nFe ).

statistical errors with the simulations for the mixing ratios as listed in Table 8.2. We calculated
the systematic errors in primary energy and the values appear in Table 8.1.
By using these result, the mean logarithmic mass number hln Ai (where A is atomic number)is calculated and is compared with those of other measurements as shown in Figure 8.29. Our
result shows good agreement with the results of the direct measurements by JACEE [8] andby RUNJOB [7] at primary energies from 1014 to 1015 eV. Moreover, hln Ai obtained with our
result is increasing with the primary energy above the knee up to 1016 eV.
Our result also agrees well with the results of our Čerenkov radiation observations [82],of CASA–MIA [35], of KASCADE with the analysis using hadrons and muons [25], and of
HEGRA–CRT [13]. Clearly, the observed hln Ai increases with the primary energy around theknee. This feature is apparently inconsistent with the preliminary results of KASCADE with
unfolding analysis using electrons and muons [94], and of CASA–BLANCA [29].
Our result combined with the direct measurement of hln Ai shown in Figure 8.29 indicatesthat hln Ai is constant up to about 1014.5 eV. Above this energy, hln Ai increases with the energy
up to 1016 eV. The factor between these two characteristic energies is about 30, and this is equal
to the charge of iron, i.e., Z = 26. Thus one possible explanation of this feature of the measured
hln Ai is that the energy spectrum of each cosmic ray component is steepen at a fixed rigidity.
As mentioned in the previous section, the present result shows that the dominant component
in the energy region around the knee becomes heavy such as iron nuclei. In our previous observations of Čerenkov radiation induced by EASs, we observed EAS longitudinal development at
the stages before shower maxima. With the present analysis, we determined the longitudinal152

output proton ratio

1

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

input proton ratio
Figure 8.28: The systematic and the statistical errors of np /(np + nFe ). The errors are estimatedfor simulated data with equivalent number of events as observed.
developments at the later stages. Nonetheless, both measurements of the chemical composition with two different and the independent observations are consistent each other. Thus, we
successfully measured whole longitudinal development of EASs with the two observations andthereby reach an estimation of the chemical composition.
The present result is consistent with the results both CASA–MIA and KASCADE(hadrons),
but inconsistent with those of KASCADE(electrons) and CASA–BLANCA. The validity of
our result is shown in the observed longitudinal development curves by comparisons with thesimulated curves of the primary protons. While the calculated EAS longitudinal development
curves are dependent on the hadron interaction model, our adopted QGSJET model shows themost rapid developments among the major models. Therefore, it is not possible to explain
our observed development curves with any hadronic interaction models with proton dominantcompositions.

153

< ln A >6

5

JACEE (direct)
RUNJOB (direct)
CASA BLANCADICECASA MIAKASCADE-electrons
KASCADE-hadrons
HEGRA (CRT)
Flys Eye
BASJE-Cherenkov (1996)
BASJE-MAS

4

3

2

1

0 12
10

10

13

10

14

10

15

10

16

10

17

10

18

10

19

Energy (eV)Figure 8.29: The mean logarithmic mass hln Ai measured by BASJE–MAS array as a function of the primary energy, are compared with the results of other experiments with balloon–borne detectors(JACEE [8], RUNJOB [7]) and ground–based detectors(CASA–MIA [35], KASCADE(hadrons) [25], HEGRA–CRT [13], KASCADE(electrons) [94], CASA–BLANCA [29],
DICE [88], Fly’s Eye [16]). Also the results of our former Čerenkov observations [82] are plotted. A hatched region represents the result of other direct observations, which are accumulatedby Linsley [65].

1548.4.2

Energy spectrum

Taking into account the determined proton mixing ratio, the average relation between the
primary energies and the observed EAS sizes for EASs with 1.0 ≤ sec θ < 1.1 was calculatedwith the simulations. The results are plotted in Figure 8.30.

primary energy E(eV)

MAS array : Size Energy Relation
18

17

16

15

14

13

12

3

4

5

6

7

8

shower size N
Figure 8.30: The relation between the primary energies and the reconstructed EAS sizes with
our array with 1.0 ≤ sec θ < 1.1
The size–energy relation was obtained with fitting of the plots, and is derived as follows,½
0.95 × log Nrec + 9.65 (log E < 15.7)
log E =
(8.29)
0.85 × log Nrec + 10.36 (log E ≥ 15.7)
Finally, using this relation, we derived the primary energy spectrum.
The resultant differential energy spectrum, dF/dE is shown in Figure 8.31 and is well
expressed by the following,½
dF −2 −1 −1 −1(19.26 ± 0.07) − (2.66 ± 0.01) × log E[eV ] (log E < 15.4)
log([m s sr eV ]) =(27.44 ± 1.77) − (3.19 ± 0.11) × log E[eV ] (log E ≥ 15.4)
dE(8.30)In Figure 8.31, we also plot the results reported by other groups. At energies around 1014
eV the present result is consistent within uncertainty with the direct measurements and with
CASA-MIA [34] and DICE [88] measurements. Comparing our result with those by Tibet andby KASCADE, both of the absolute intensity and the energy of the knee in our spectrum arelow. These differences could be due to the systematic difference of energy estimation procedures
155

10

5

10

4

10

3

10

2

2.5

-2

-1 -1

1.5

E dF/dE (m sr s GeV )



10

12

10

13

10

14

10

15

10

16

10

17

E (eV)
Figure 8.31: The differential all–particle cosmic ray flux measured by BASJE–MAS array.
Also plotted are the cosmic ray fluxes reported by JACEE [8], RUNJOB [7], SOKOL [52],proton satellite [36], KASCADE(hadrons) [25], KASCADE(electrons) [94], CASA–BLANCA[29], CASA–MIA [34], DICE [88], Tibet [4], EAS–TOP [1], and the dashed line represents theflux measured by Akeno group [72].
as mentioned by the CASA–MIA group [35]. According to their discussion, the required energy
shift is small. Lowering the Tibet energy scale by 20% would reduce the discrepancy between
the experiments.
We show simulated equi–intensity curves with the measured ones without any normalizationin Figure 8.26(b). In the simulation, we took into consideration not only the obtained mixing
ratios of protons and iron nuclei, but also the all particle energy spectrum(8.30). In this plot weshow only the statistical uncertainties. When we consider the uncertainties in the energy spectrum determination, which strongly affect the uncertainties in the rate of arrival, we concluded
that the simulation gives good agreement with the experimental results.

156

Chapter 9

Discussion
9.1

The systematic errors in the determination of chemical composition and primary energy

In the present analysis, we assumed the primary cosmic ray component is a mixture of two
nuclei, i.e., protons and iron nuclei. To see how the longitudinal developments reflect thecomposition, here we consider (4.8) and (4.9) in the superposition approximation for incident
nuclei. The superposition model leads the mean depth of maximum development, Xmax , for a
primary energy E0 as follows,
¶
µ
X
E0
(9.1)
Xmax = X0
wi ln
Ai Ecritiµ¶
E0= X0 ln− hln Ai ,
(9.2)
Ecrit
where X0 is the radiation length for the air, Ai is the mass number, wi is the relative abundanceof the component i, and Ecrit is the critical energy for the air. In the relation (9.2), it is shownthat the longitudinal development curves depend only on hln Ai. That is to say, in the superposition approximation, through the measurement of the averaged properties of the longitudinaldevelopment curves, we can determined a correct hln Ai for any mixture of components in the
superposition approximation.
With the detailed Monte–Carlo studies, we can estimate the tolerance of the superposition
approximation. For example, the mixture of 50% proton + 50 % iron nuclei and the mixture of
10% protons + 90% carbon nuclei show the equivalent Xmax , 655 g/cm2 and 653 g/cm2 by thesimulations with CORSIKA and QGSJET model(log E = 16.8 − 17.0). On the other hand, eachhln Ai for the different mixtures considered above is 2.01 and 2.24, respectively. The differenceof hln Ai for these mixture is about 0.2, and this value is a systematic error based on the choiceof the components for the analysis.
One can estimate primary energy from the two different observable values, i.e., the rates ofarrivals and the EAS sizes. The primary energies estimated from the rates of arrivals (8.30),
Erate and the energies determined from the EAS sizes with the size–energy relation (8.29), arelisted in Table 9.1. From this table, the averaged difference log Esize − log Erate is about 0.2.
This is the systematic energy determination errors for the present result totally including theeffects of the equipments, the composition uncertainty, the analysis and the simulations.
157

log F (> N )
−5.20
−5.45
−5.70
−5.95
−6.20
−6.45
−6.70

log Erate [eV]
14.61
14.76
14.91
15.05
15.19
15.32
15.44

log N
5.41
5.56
5.71
5.86
6.00
6.13
6.26

log Esize [eV]
14.79
14.93
15.07
15.22
15.35
15.47
15.60

log F (> N )
−6.95
−7.20
−7.45
−7.70
−7.95
−8.20

log Erate [eV]
15.55
15.66
15.78
15.89
16.00
16.12

log N
6.38
6.49
6.59
6.70
6.87
6.98

log Esize [eV]
15.78
15.88
15.96
16.06
16.20
16.30

Table 9.1: The table of the integral fluxes, corresponding primary energies Erate , the average
EAS sizes and corresponding energy Esize The unit of integral flux F (> N ) is m−2 sr−1 s−1 .

In the equi–intensity method analysis it is postulated that showers which are observed withthe same frequency of occurrence at different zenith angles are initiated by particles of thesame energy. Here the validity of this assumption is checked with Monte Carlo simulations.
Figure 9.1 shows the relation between the mean energies and the atmospheric depths on theequi–intensity curve corresponding to a certain rate of arrival. In the Monte Carlo simulation
I used the mixing ratios shown in Table 8.1. The averaged energies over all the components
are practically constant as shown in Figure 9.1. However, the mean energies of each component
vary with the atmospheric depths because of the difference of the longitudinal developments
between air showers induced by protons and those induced by iron nuclei. Therefore, the mixing
ratios on the equi–intensity curve also vary with the atmospheric depths as shown in Figure
9.2. Practically, the fact that the mean energies and the mixing ratios are not constant on anequi–intensity curve does not causes the systematic errors on the determinations of the mixing
ratios. Because the measured equi–intensity curves are compared with the calculated ones withthe air shower simulations including the assumptions of the all–particle energy spectrum andthe chemical composition.
The difference of the indexes of component spectra between protons and iron nuclei is apossible cause of a systematic error, because in equi–intensity method analysis we supposedthat all the spectral indexes are identical in a range of the energies of cosmic rays contributing
to a equi–intensity curve. However, the indexes are not identical as shown in Figure 9.18. In
order to estimate the systematic error, I simulated the equi–intensity curves for two cases. In
case 1, the indexes of the component spectra are not identical in a range of the energies of thecontributing cosmic rays, i.e., hln Ai is increase with the energy in the range. In case 2, theindexes are identical, i.e., hln Ai is constant in the energy range. The comparisons between
case 1 and case 2 are shown in Figure 9.1 and 9.2. The differences between the two cases arevery small, and less than the measurement errors, and then the systematic errors based on the
different indexes are negligible.

9.2

Comparisons the results with the source models

We compare the results of our all–particle flux and hln Ai to those of a composition model withfive cosmic ray components (protons, He, CNO, Ne–Si, Fe). In this model, it is assumed thateach component spectrum has the index that measured by the RUNJOB collaboration and thecomponent spectra are steepened at the fixed rigidity 1014.5 V. The calculated flux of each
158

15.5

15.4

15.3

15.2

15.1

15

14.9

14.8

14.7

14.6

14.5

1

1.05

1.1

1.15

1.2

1.25

1.3

1.35

1.4

1.45

1.5

Figure 9.1: These plots are simulated results of the relation between the mean energies and theatmospheric depths on the equi–intensity curve corresponding to a certain rate of arrival.
component is summed each other according to the relative abundances measured by SOKOL[52] at 1012 eV and the total flux is normalized to the all–particle spectrum obtained also by
SOKOL at the same energy. Moreover, we examined two different cases for the model. In
the first case (a), each spectral index is steepened by 0.6, corresponding to the rigidity. Thisis expected in the case that the energy dependence of the diffusion coefficient dominates thecosmic ray propagation processes changes at a fixed rigidity. In the second case (b), the spectral
index changes, irrespectively of A, to −3.2 in the energy corresponding to the same rigidity.
This is expected in the case that the dominant acceleration process of cosmic rays is changed
above the rigidity. The values of 0.6 in model (a) and -3.2 in model (b) are assumed with ourmeasured all–particle spectrum. The calculated spectra and the resultant hln Ai are shown in
Figures 9.3 and 9.4, respectively. Although both the calculated fluxes in Figures 9.3(a) and9.3(b) are slightly less than the measured one at 1014.7 − 1015.7 eV, the all–particle fluxes at theother energy range and the predicted hln Ai of models (a) and (b) are consistent with present
results. This suggests that iron nuclei are the dominant component at the primary energies
greater than 1015 eV. The model predictions do not fit the measured spectrum between 1014.7
and 1015.7 eV and result in two knees at 1014.7 eV and at 1015.7 eV. Therefore, the simple models
described here are not sufficient to produce the measured spectrum and composition.
In the report of HEGRA–CRT group [13], they suggested that the each spectrum of the
primary components is steepened at a fixed rigidity and that the dominant component at theknee energy is CNO. They also showed an increase in hln Ai with the energy. Their simple
model is consistent with our present result up to 1015 eV. However, the hln Ai in their model
saturates around this energy and does not fit to our result at the higher energies. Hörandel [49]proposed the model with introducing the charge dependent cut–off energies and the ultra–heavy
nuclei(Z=30–92) components, this is inconsistent with our result because the model predicts
protons are dominant at the knee.
159

1

0.95

0.9

0.85

0.8

0.75

0.7

0.65

0.6

0.55

0.5

1

1.05

1.1

1.15

1.2

1.25

1.3

1.35

1.4

1.45

1.5

Figure 9.2: These plots are simulated results of the relation between the mixing ratio and theatmospheric depths on the equi–intensity curve corresponding to a certain rate of arrival.
The model of the particle acceleration by oblique shocks proposed by Kobayakawa, et al.
[63] predicts the knee and the gradual increase of hln Ai with the increasing energy from 1014 to1016 eV without any assumption of a rigidity dependent cut–off. Their prediction on hln Ai is
consistent with our result, but the predicted absolute value of hln Ai is smaller than our result.
In the model of Völk and Biermann [95], cosmic rays from 1013 eV to the knee are mainly
accelerated during explosions of massive stars. Biermann [14] developed this model further andexamines explosions of Wolf–Rayet stars. He concluded that at the knee, the particles segregatewith particle energy according to their charge, that protons drop off first, then the C–N–O
elements, next Mg, Si, etc., and finally iron nuclei. At the surfaces of Wolf–Rayet stars helium
and heavier elements are enhanced rather than protons. This can be reflected effectively to thechemical composition of primary cosmic rays. As discussed in our previous paper by Shirasaki etal. [82], the measured hln Ai suggests that the accelerated particle abundance must be heavier
than the stellar winds of Wolf–Rayet stars. Since the accelerated particles are a mixture of
the stellar wind particles and ejected matters, Biermann’s model seems to be very promisingfavourable to our former and present results.

160

4

10

3

2.5

E dN/dE

10

2

10
12
10

10

13

10

14

10

15

10

16

10

17

E (eV)

4

10

3

2.5

E dN/dE

10

2

10
12
10

10

13

10

14

10

15

10

16

10

17

E (eV)

Figure 9.3: All–particle spectrum and the contributions of five components calculated withmodel (a)(upper panel) and with model (b)(lower panel) are compared with the measured
spectra.

161

< ln A >6

5

JACEE (direct)
RUNJOB (direct)
CASA BLANCADICECASA MIAKASCADE-electrons
KASCADE-hadrons
HEGRA (CRT)
Flys Eye
BASJE-Cherenkov (1996)
BASJE-MAS

4

3

2

1

0 12
10

10

13

10

14

10

15

10

16

10

17

10

18

10

19

Energy (eV)Figure 9.4: The predicted mean logarithmic mass hln Ai with the model (a)(solid line), and(b)(dashed line).

162

9.3

Comparisons the results with the advective diffusion propagation model

In the previous section, I presented a very simple models to explain the chemical composition
and the energy spectrum. In that model the fixed rigidity at the bending points of the spectra
and the index jump of the steepening are determined empirically. The advective diffusion
model described in section 3.4 for the calculations of the residence time in the galactic disk withthe advective–diffusion function naturally predicted the observed features such as the energy
spectrum, the anisotropy, the mean mass number and its energy dependence.
Using the calculated residence times τR in section 3.4.3, the component spectra, the allparticle spectrum, the chemical composition and the anisotropy amplitude are estimated witha simple assumed relation N (E) = Q(E)τR (E). The expected flux of i–th component Ni andall particle flux, Nall , were given by
Ni = τRi × ci AE −γi
XNi
Nall =

(9.3)
(9.4)

i

where ci is relative abundance of i–th component at a source, A is a normalization factor and
γi is the spectral index of i–th component. For this calculation, ci were determined from theconsequences of the direct observations for primary cosmic rays of less than 1012 eV, because
the observed spectra for low energies are expected to be equal to source spectra in the advective
diffusion model. Here ci are calculated from the measured relative abundance by SOKOL on
COSMOS satellite above 2.5 TeV [52], and these are listed in Table 9.2, and consistent with thereported values by CRN for more lower energy cosmic rays [71].
component
proton
He
CNO
Ne–Si
Fe

abundance(E > 1012 eV)
0.39
0.27
0.13
0.10
0.11

model J
2.7
2.5
2.4
2.6
2.5

spectral index
model RC model RF
2.6
2.7
2.7
2.7
2.4
2.6
2.6
2.6
2.5
2.4

Table 9.2: Assumed relative abundance and three models of the spectral indices for five component of cosmic ray particles.
The spectral index γi for i-th component is estimated from the direct measurements of
SOKOL, Ryan et al. [78], Sanriku [60], JACEE [8] and RUNJOB [7]. Since the measured
spectra by SOKOL, JACEE and RUNJOB are different each other, I took the three different
sets of γi listed in Table 9.2. Model J is chosen to fit the results of JACEE, and this has rather
flatter spectra for He and CNO compared with the proton spectrum. Model RC is chosen to
fit the He spectrum by RUNJOB and the spectra of other component are consistent with theresult of JACEE and with that of SOKOL. This model predict a flat CNO spectrum. Model
RF is chosen to fit the spectra by RUNJOB, except the Fe spectrum. This model has a slightly
flat Fe spectrum.
The calculated all particle spectrum was normalized to that of observed at 1012 eV withthe factor A. The other parameters in this model are vg and Lirr . These parameters are not163

determined with the current observations as mentioned in section 3.4. Thus, I used the six pairof these parameters as listed in Table 9.3 to calculate the energy spectrum, the composition
and its energy dependences.
Model ID
Lirr (pc)vg (km/s)(pc/yrs)

110

2
3
50 100490
0.0005410

5
650 100196
0.0002

Table 9.3: The model ID and the corresponding sets of the parameters pairs, Lirr and vg .

9.3.1

All particle spectrum

The all particle spectra calculated with the models J − 1(spectral index model J + model ID=1),
J − 2 and J − 3 are shown in Figure 9.5, and those with the models J − 4, J − 5 and J − 6 areshown in Figure 9.6. The spectrum of J − 1 is too flat and too large flux at energies less than1014 eV, but this curve fits well with the Tibet and KASCADE–electron results. The curvewith the model J − 2 is gradually steepening around 1015 eV and the spectral index is changed
from −2.7 to −3.1. These features of the spectral shape and the flux which are predicted withthe model J − 2 are consistent with the direct measurements below 1014 eV and our results of
BASJE–MAS, CASA–MIA and the other measurements up to 1017 eV. The predicted flux withthe model J − 3 is slightly less than the observed results in the energies greater than 1015 eV.
The curves of J − 5 and J − 6 are too steep to fit the measured spectrum, because theadvection velocity vg = 0.0002pc/years is too slow, so that the diffusion process is dominatedfor the whole energies regions 1012 − 1017 eV. Moreover, the curve of model J − 4 resembles thatof J − 2, thus only the models of 1, 2 and 3(vg = 0.0005pc/years) are considered for further
discussions.
The spectra predicted by the model RC and RF which are shown in Figure 9.7 and 9.8. These
curves are not much different from that of the model J. From Figure 9.8, the spectra of the model
RF have the bends with index of −2.7 to −3.0 at the knee, because the dominant component
above the knee is iron nuclei, which has rather flatter spectrum of the index −2.4 than the othermodels. The curves of RF − 1 and RF − 2 are not inconsistent with the measurements, but theflux predicted by RF − 2 at 1015 − 1016.5 eV is less than the observed one.

164

4

10

3

10

2

2.5

-2

-1 -1

1.5

E dN/dE (m sr s GeV )

10

10

12

10

13

10

14

10

15

10

16

10

17

E (eV)

Figure 9.5: Calculated all particle spectra with the advective–diffusion models of J − 1(solidline), J − 2(thick dashed line) and J − 3(dotted line). The markers represent the observed
energy spectra, by SOKOL [52](¤), proton–satellite [36](°), RUNJOB [7](¥), JACEE [8](4),CASA–MIA [34](+), KASCADE–hadron [25](?), Tibet [4](∗), CASA–BLANCA [29](H), DICE[88](♦), EAS–TOP [1](F), KASCADE–electron [94](N) and this work(BASJE-MAS)(•). Thindashed line represents the flux measured by Akeno group [72].

4

10

3

10

2

2.5

-2

-1 -1

1.5

E dN/dE (m sr s GeV )

10

10

12

10

13

10

14

10

15

10

16

10

17

E (eV)

Figure 9.6: Calculated all particle spectra with the models of J − 4(solid), J − 5(thick dashed)
and J − 6(dotted).
165

4

10

3

10

2

2.5

-2

-1 -1

1.5

E dN/dE (m sr s GeV )

10

10

12

10

13

10

14

10

15

10

16

10

17

E (eV)

Figure 9.7: Calculated all particle spectra with the models of RC − 1(solid), RC −2(thick dashed)
and RC − 3(dotted).

4

10

3

10

2

2.5

-2

-1 -1

1.5

E dN/dE (m sr s GeV )

10

10

12

10

13

10

14

10

15

10

16

10

17

E (eV)

Figure 9.8: Calculated all particle spectra with the models of RF − 1(solid), RF − 2(thickdashed) and RF − 3(dotted).

1669.3.2

Spectra for the element groups

The most significant differences of the observed spectra by JACEE and RUNJOB are appearedin the He and CNO flux. The spectra for the element groups are shown in Figure 9.9 which arecalculated with the model J, 9.10 with RC and9.11 with RF . The spectra of model J fits well with the JACEE data within the error as
shown in Figure 9.9. The He spectrum of the model RC , and the He and CNO spectra of
the model RF are based on the RUNJOB data. For every models, the calculated spectra are
well consistent with the observed spectra assumed in whole energy range in which the direct
measurement data exist.
4

10

3

10

2

2.5

E dN/dE

10

10
12
10

10

13

10

14

10

15

10

16

10

17

2.5

E dN/dE

E (eV)
10

3

10

2

10
1
-1

10
10

12

10

13

10

14

10

15

10

16

10

17

E (eV)

Figure 9.9: The lines are the calculated energy spectra for the assumed five components withthe advective–diffusion models of J − 1(solid), J − 2(dashed), J − 3(dotted). The marker plots
represent the observed energy spectra of five groups of cosmic ray species by the direct measurements of Ryan et al. [78](♦), SOKOL [52](°), JACEE [8](•), CRN [71](+) and RUNJOB[7](H).

167

4

10

3

10

2

2.5

E dN/dE

10

10
12
10

10

13

10

14

10

15

10

16

10

17

2.5

E dN/dE

E (eV)
10

3

10

2

10
1
-1

10
10

12

10

13

10

14

10

15

10

16

10

17

E (eV)

Figure 9.10: The lines are the calculated energy spectrum for each component with the modelsof RC − 1(solid), RC − 2(dashed), RC − 3(dotted).

4

10

3

10

2

2.5

E dN/dE

10

10
12
10

10

13

10

14

10

15

10

16

10

17

2.5

E dN/dE

E (eV)
10

3

10

2

10
1
-1

10
10

12

10

13

10

14

10

15

10

16

10

17

E (eV)

Figure 9.11: The lines are the calculated energy spectrum for each component with the modelsof RF − 1(solid), RF − 2(dashed), RF − 3(dotted).

1689.3.3

Average mass number

As an indicator of the composition, here we calculate the average value of the logarithm of mass
number A given byP(Ni × ln Ai )P
< ln A >≡(9.5)
Ni< ln A >where Ni is obtained with equation (9.3). The results of the model calculations are shown in
Figure 9.12, 9.13 and 9.14 with the observed results. The difference among the models of 1, 2
and 3 is not significant, and every model has the same tendency, that is, hln Ai increases with
the increasing energy up to 1017 eV. This feature is consistent with many experimental results
including the BASJE–MAS result, but not consistent with some measurements. In the model J
and the model RC , which show the same tendency, hln Ai are saturated below 3, so that thesemodel are inconsistent with our results. On the other hand, the curve of the model RF shows agood agreement with the direct measurements by JACEE, BASJE–MAS and some observationsin the energy range 1012 eV − 1017 eV.
65

JACEE (direct)
RUNJOB (direct)
CASA BLANCADICECASA MIAKASCADE-electrons
KASCADE-hadrons
HEGRA (CRT)
Flys Eye
BASJE-Cherenkov (1996)
BASJE-MAS

4

3

2

1

0 12
10

10

13

10

14

10

15

10

16

10

17

10

18

10

19

Energy (eV)Figure 9.12: The averaged mass number, < ln A >, calculated with the models of J − 1(solid),J − 2(dashed) and J − 3(dotted). The estimated values with various different types of measurements [82]–[16] are also plotted in this figure. A hatched region represents the result of otherdirect observations, which are accumulated by Linsley [65].

169

< ln A >6

5

JACEE (direct)
RUNJOB (direct)
CASA BLANCADICE

KASCADE-electrons
KASCADE-hadrons
HEGRA (CRT)
Flys Eye

CASA MIABASJE-Cherenkov (1996)
BASJE-MAS

4

3

2

1

0 12
10

10

13

10

14

10

15

10

16

10

17

10

18

10

19

Energy (eV)< ln A >Figure 9.13: The averaged mass number, < ln A >, calculated with the models of RC − 1(solid),
RC − 2(dashed) and RC − 3(dotted).

6

5

JACEE (direct)
RUNJOB (direct)
CASA BLANCADICE

KASCADE-electrons
KASCADE-hadrons
HEGRA (CRT)
Flys Eye

CASA MIABASJE-Cherenkov (1996)
BASJE-MAS

4

3

2

1

0 12
10

10

13

10

14

10

15

10

16

10

17

10

18

10

19

Energy (eV)Figure 9.14: The averaged mass number, < ln A >, calculated with the models of RF − 1(solid),
RF − 2(dashed) and RF − 3(dotted).

170

9.3.4

Residence time and anisotropy

The calculated residence times with the model J, RC and RF are shown in Figure 9.15(a),9.16(a), 9.17(a). The differences of the residence times between the models are not significant,τR ∼ 105 years at 1012 eV and ∼ several × 104 years at 1017 eV. The calculated τR at 1012 eV
is of the same order of magnitude as the estimated value, 5 × 105 years, from the anisotropy
amplitude described in Section 3.4.
The predicted anisotropy amplitudes estimated with the relation of A ∝ 1/τR are shownin Figure 9.15(b), 9.16(b), 9.17(b). The normalization factors were chosen to fit with themeasurements at 1012 eV. The model estimations have the tendency for the anisotropy amplitude
to be almost constant up to ∼ 1014 eV , and to increase with the increasing energy above∼ 1014 eV. This feature and the estimated values are consistent with the measured amplitudes.

6

10

5

10

4

τR

10

12 12.5 13 13.5 14 14.5 15 15.5 16 16.5 17

log10E

10

11

12

13

14

15

16

17

log10E

Figure 9.15: (a)The residence times calculated with the models of J − 1(solid), J − 2(dashedline) and J − 3(dotted). (b)The measured anisotropy amplitude compared with the estimatedvalues from τR on (a). The estimated anisotropies are calculated with a relation A ∝ 1/τR
and a normalization factor obtained at 1012 eV. The measured anisotropies are accumulated by
Hillas [45](°), Linsley [65](•) and Hayakawa [40](N).

171

6

10

5

10

4

τR

10

12 12.5 13 13.5 14 14.5 15 15.5 16 16.5 17

log10E

10

11

12

13

14

15

16

17

log10E

Figure 9.16: (a)The residence times calculated with the models of RC −(solid), RC − 2(dashedline) and RC − 3(dotted). (b)The measured anisotropy amplitude compared with the estimatedvalues from τR in (a). The estimated anisotropies are calculated with a relation A ∝ 1/τR anda normalization factor obtained at 1012 eV.

6

10

5

10

4

τR

10

12 12.5 13 13.5 14 14.5 15 15.5 16 16.5 17

log10E

10

11

12

13

14

15

16

17

log10E

Figure 9.17: (a)The residence times calculated with the models of RF − 1(solid), RF − 2(dashedline) and RF − 3(dotted). (b)The measured anisotropy amplitude compared with the estimatedvalues from τR in (a). The estimated anisotropies are calculated with a relation A ∝ 1/τR anda normalization factor obtained at 1012 eV.
172

9.3.5

Summary of the advective diffusion model

Here we summarize the advective–diffusion model. The comparisons of the measured all particle
spectrum with the model calculated ones showed that the most optimum pair of the parameters
is vg = 0.0005 pc/years and Lirr = 50 pc. These values agree with the observations of thegalactic magnetic fields and a model of the galactic wind.
In Figure 9.18, 9.19 and 9.20, we show the all particle spectrum and the contributions of
the element groups with the model 2 with the spectral index model J, RC and RF . The model
RF − 2 is only model of which hln Ai fit the satellite–, the balloon borne– and our measurements.
In RF model the He and CNO spectral indices are assumed to be same as the RUNJOB data.
When the flatter spectra for these species as measured by JACEE will be established, it isnecessary to exist an extra contribution of heavy nuclei in cosmic rays with energies greaterthan 1015 eV to fit with the observed results.
The advective–diffusion model predicts the energy independent escape time for cosmic rays
with energies below about 1013 eV. This is inconsistent with the prediction of the leaky box
model, τesc ∝ R−0.6 , where R is a rigidity. But this contradiction can be solved by considering
the nested leaky box model [23]. In this model, it is assumed that small confinement regions
are surrounding the sources with relatively high density in which particles diffuse for a shorttime. The energy dependence of the secondary to the primary nuclei ratio is attributed to anenergy dependent leakage from the source regions, and is characterized by the energy dependent
escape length, λ1 (E). The Galaxy is considered as an outer volume in which the nuclei fromshrouded sources may traverse a further small amount of matter with the energy independent
escape length, λ2 (Figure 9.21). In this model, the earth is not inside the inner volume, so that,
we observe the source spectrum. Thus cosmic ray accelerators will need to have a differential
spectrum with γ ' 2.7. This prediction agrees with the assumption for the advective diffusion
model. Therefore, we can draw a following scenario: the low energy cosmic rays propagate in
the galactic halo with a process of the nested leaky box model, and for the higher energy cosmic
rays, the confinement volume is the disk and such cosmic rays are leaking with the advective
diffusion process, as shown in Figure 9.21.
The escape times predicted with the nested leaky box model and with the advective diffusion
model are not depend on the particle energies, so that the spectral index is not changed up to
the knee. However, the outer confinement volume of the nested leaky box will have to be larger
than that in the simple leaky box model to allow the 10 Be/Be ratio to fall to its observed
value with keeping small λ2 . Since the outer confinement volume is not the galactic disk as theadvective diffusion model, the size of confinement volume should change at the intermediate
energy, which may be less than 1012 eV, and this evidence would appear in an observed primaryenergy spectrum. Thus, we need more detailed investigations of the properties of the galactic
halo and the energy spectrum above 1012 eV.

173

4

10

3

10

2

2.5

-2

-1 -1

1.5

E dN/dE (m sr s GeV )

10

10

12

10

13

10

14

10

15

10

16

10

17

E (eV)

Figure 9.18: Calculated all particle spectrum and the contributions of five components with themodel J − 2.

4

10

3

10

2

2.5

-2

-1 -1

1.5

E dN/dE (m sr s GeV )

10

10

12

10

13

10

14

10

15

10

16

10

17

E (eV)

Figure 9.19: Calculated all particle spectrum and the contributions of five components with themodel RC − 2.
174

4

10

3

10

2

2.5

-2

-1 -1

1.5

E dN/dE (m sr s GeV )

10

10

12

10

13

10

14

10

15

10

16

10

17

E (eV)

Figure 9.20: Calculated all particle spectrum and the contributions of five components with themodel RF − 2.

175

leaky box model

λ esc (Ε)

Galaxy
Earth

Halo?

nested leaky box modelλ2
λ 1 (Ε)sourceinner volume
Earthouter volume

confinement volume

advective diffusion modelλL
λH
inner volume

inner volume
Earth

Earth

outer volume

outer volume
higher than ~1012eV

12

lower than ~10 eV

Figure 9.21: The schematic view of three cosmic propagation models.

176

Chapter 10

Conclusions
10.1

Summary of observational results

Using the equi–intensity method we obtained the mean equi-intensity curves of EASs with theprimary energies from 1014 to 1016 eV observed with the MAS array at Mt. Chacaltaya, 5200
m a.s.l.(550 g/cm2 ). The data used in the analysis consists of 7.5 × 107 selected EAS events.
The apparent maximum development points which are expected to be found with a proton
dominant composition model are not found in the measured equi-intensity curves around theknee energy. By comparing the measured curves with those calculated with the Monte–Carlosimulation, the mean logarithmic mass hln Ai was obtained as a function of the primary energy.
The measured hln Ai increases with the increasing energy in the range of 1014.5 − 1016 eV, andthis is consistent with our former Čerenkov light observations and with the measurements ofsome other groups. The observed hln Ai suggests that the steepening of the component spectra
is caused at a single rigidity of 1014.5 V.
The measurement of Čerenkov light pulse shapes corresponds to an observation of the longitudinal developments of air showers at the earlier stages before the maximum developments. In
contrast, with the equi-intensity method analysis we determined the longitudinal developments
around their maximum and at the later stages. Therefore, with two different observations,
we measured the whole stages of EAS longitudinal developments and then we could make thetrustworthy conclusion on the chemical composition based.
From the observational results, we obtained a power–law form EAS size spectrum for themeasured EAS events with the zenith angle range of sec θ = 1.0 − 1.1,½
dF
(3.85 ± 0.03) − (1.57 ± 0.01) × log N (log N < 5.9)
−2 −1 −1[m s sr ]) =(10.1)
log((6.80 ± 0.87) − (2.07 ± 0.14) × log N (log N ≥ 5.9)
d log N
where N is the measured size.
With the investigated chemical composition we obtained the energy spectrum of primarycosmic rays in the energy range of 1014.5 − 1016 eV,½
dF −2 −1 −1 −1(19.26 ± 0.07) − (2.66 ± 0.01) × log E[eV ] (log E < 15.4)
log([m s sr eV ]) =(27.44 ± 1.77) − (3.19 ± 0.11) × log E[eV ] (log E ≥ 15.4)
dE(10.2)The present result of cosmic ray flux is consistent with other experiments, and the obtained
all–particle spectrum has a gradual steepening around 1015.5 eV with the spectral index jumpfrom −2.66 to −3.19. Around the energies of 1014 eV the present result is consistent with those177

of the direct measurements within uncertainties and consistent with those reported by CASAMIA and by DICE. Comparing our result with those by Tibet and by KASCADE, both of theabsolute intensity and the knee energy in our spectrum are inconsistent. These differences could
be due to the systematic difference of energy estimation procedures.

10.2

Cosmic ray source models and propagation models

The supernova acceleration model with the stellar winds and with ejected matters of Wolf–
Rayet stars is one of the plausible models to explain our results. We conclude that the actualmodel should predict that the dominant component above 1015 eV is iron nuclei, and that thehln Ai increases with the increasing energy and it becomes about 3.5 at 1016 eV.
With the detailed research for the diffusive motions of charged particles in the galactic
turbulent magnetic fields, I calculated the residence times of cosmic rays in the galactic disk
with solving an one dimensional advective–diffusion equation. Consequently, this simple model
predicts all the observed features, i.e., the energy dependences of the flux for the all–particle
and the components, the chemical composition and the anisotropy with assuming the naturalvalues of the parameters for the equation, such as the magnetic field strength and the galactic
wind velocity. However, in this model we assumed that the regular magnetic fields are open
toward the galactic halo, that the maximum energy of cosmic ray accelerations is at least 1016
eV and that the source spectrum index γ ' 2.7. Thus, we need more detailed investigations forthe galactic magnetic fields and the galactic winds and new discoveries and investigations forcosmic ray accelerators with observations of cosmic ray anisotropies [93] and very high energy
gamma rays [26].

10.3

Future Prospects

The results of our measurement and other observations suggest that cosmic rays with energies
less than 1016 eV originate in this galaxy. Up to the present date cosmic rays with energies
greater than 1020 eV are observed with air shower experiments. Ultra high energy cosmic rays(> 1019 eV) may be of extra galactic origin because the arrival direction distribution of ultrahigh energy cosmic rays does not show any enhancement of cosmic rays from the galactic plain.
Therefore, the maximum accelerated energy for galactic cosmic ray sources is expected to be lessthan 1020 eV, and the slight bend in the energy spectrum at 1018.5 eV (called ankle) is believed
to show the transition from galactic sources to extra galactic ones with the increasing energy.
The difference between origins of cosmic rays will be observed as the variation of the chemicalcomposition with th primary energy, so that, the study of the chemical composition of cosmic
rays with energies greater than the knee gives some crucial information on possible galacticand extra galactic sources and acceleration sites of cosmic rays. The BASJE collaboration ispreparing the observations of air showers with primary energies from 1015 to 1017 eV. The
particle detectors of the MAS array will be re–deployed in an array with the area of 400 m ×
700 m. The rearrangement will be finished in 2005, and the observations will be continued for
3 years.
For air showers with primary energies greater than 1017 eV, the measurement of Xmax by the
Fly’s Eye group showed that the averaged mass number is decreasing with the increasing primaryenergy. On the other hand, the Akeno group reported that they cannot observe any evidenceof the variation of the chemical composition, so that, both of the results are contradictory to178each other. The chemical composition of this energy region will be stringently verified by the
proposed TALE (Telescope Array Low energy Extension) project using the hybrid detection
method with fluoresces light telescopes and surface particle detectors [46].
When an air shower is detected with a surface detector array, the primary energy and thespecies are not determined independently each other. This is the main drawback to measuringcosmic rays with air shower arrays. If we determine a primary energy of a measured air shower,
we need information of the stage of the longitudinal development of the shower. However, thelongitudinal developments depend on primary species. Therefore, the size–energy relation id
determined with Monte Carlo simulations taking a assumed chemical composition into consideration. For estimations of the chemical composition, muons or Čerenkov photons are measured
simultaneously with air shower sizes, or we use the equi–intensity method. However, primary
species are not determined on event by event in any methods, only the averaged mixing ratioof species can be estimated. In order to obtain the energy spectra of each component groupwith air shower observations, we need a new method for the determinations of primary species
independently with primary energies. The most promising method is the hybrid observationsof air shower, i.e., the simultaneous detections with air fluorescence telescopes and with surface
detector arrays. The arrival directions and the core positions of air showers are determined by
the surface array. On the other hand, the fluorescence telescopes are analogous with hadroncalorimeters. The primary energies are determined from atmospheric fluorescence photon yields,
and we can separately obtain longitudinal development curves. Therefore, in the hybrid observations, the estimations of species are independent with the determinations of primary energies,
and then it is possible to estimate the primary species on event by event. With simple Monte
Carlo simulations of air shower developments, it is demonstrated that we can classified the primary particles according to the longitudinal developments into three different types of species,
i.e., light, middle and heavy nuclei, and can obtain the energy spectra of each component withhybrid measurements [70].
Since accelerated high energy particles interact with magnetic fields and the ambient mattersof sources, sources or their environments radiate high energy gamma rays. So that, we believethat the astronomical objects which accelerate cosmic rays are observed as high energy gamma
ray sources. Many sources of gamma rays with energies less than 1013 eV have been discovered
with satellites and ground–based Čerenkov telescopes. However, for all the very high energy
gamma ray sources, the process are explained with the synchrotron radiations and subsequent
inverse Compton scattering, i.e., the emissions of high energy electrons. SN1006 is only instanceof the source which radiate π 0 decay gamma rays induced by high energy nuclei [91]. However
the emissions from this object also can be explained with the radiations from electrons. Up
to now, acceleration sites of cosmic rays with energies greater than 1014 eV have not been
discovered yet. Even for lower energy cosmic rays which are accelerated in supernova remnants,
source objects are not identified. We hope discoveries of the sources by observations with thenext generation Čerenkov telescopes, such as HESS [43], MAGIC [68], etc.
Since cosmic charged particles move on complex trajectories in galactic magnetic field, wecannot trace back to their sources directly from their arrival directions. However, if there are
nearby cosmic ray sources or a gradient in the source distribution, we expect a small anisotropy
oriented to their directions. We found a significant anisotropy in the galactic longitude regionl = 220 ∼ 260◦ from the observations with the BASJE MAS array, and we concluded that thisanisotropy is due to the nuclear component of cosmic rays [93]. The calculations of cosmic ray
propagations in the galaxy clarify the contributions of the local sources to the cosmic ray flux atthe earth as well as the resultant anisotropies. It is shown that the Vela SNR and Monogem can179give rise to the anisotropy found in our observations. In order to identify the contributing source
to the anisotropy we need detailed measurements of the anisotropy and its energy dependence
and need to reveal the structures of the galactic magnetic field.
Cosmic rays with energies less than 1014 eV are accelerated in supernova remnants andpropagate in the galactic halo. This model had been confirmed by observations of the flux ratio
between the primary nuclei and the secondary nuclei which are fragmentation of the primary
nuclei, and by measurements of the abundance ratios of radio isotopes. We need to measure
the amount of matter traversed and the escape time of cosmic rays with energies greater than1014 eV in order to definitively identify the acceleration and propagation mechanisms. For thesemeasurements we have to directly detect primary cosmic rays with satellites or balloon–borne
equipment. The CREAM experiment [24] using the transition radiation detectors challenges todirect measurements of knee region cosmic rays, but we need improved detectors having morelarger effective area and the higher charge resolution.

180

Acknowledgments
The authors would like to thank to Prof. Fumio Kakimoto for his support to over all mycosmic ray researches. I also acknowledge Prof. Hisashi Yoshii(Ehime University), Prof. Yutaka Matsubara(Solar–Terrestrial Environment Laboratory, Nagoya University), Prof. Kazuaki
Murakami(Nagoya University), Prof. Tatsunosuke Kaneko(Okayama University), Dr. Yuji Shirasaki(National Astronomical Observatory of Japan), Dr. Yoshiki Tsunesada(National Astronomical Observatory of Japan), Prof. Yoshihiko Mizumoto(National Astronomical Observatoryof Japan) and Prof. Yoshio Toyoda(Fukui Institute of Technology) for their useful comments
and discussions. I acknowledge to Mr. Katsuo Nishi, Mr. Susumu Shimoda, Mr. Norio Tajima
and Mr. Yutaka Yamada for their support for the observations. I am grateful to all othercolleagues in the BASJE Collaboration. I would like to appreciate helpful support by all thestaff of Instituto de Investigaciones Fisicas, Universidad Mayor de San Andres, La Paz, Bolivia.
Especially I wish to acknowledge Ing. Pedro Miranda and Prof. Alfonso Veralde.
I express thanks to Dr. Kenichi Kadota(Musashi Institute of Technology), Mr. Daijiro
Senda(Tokyo Institute of Technology), Mr. Yutaka Kurashina(Tokyo Institute of Technology),
Dr. Hisao Tokuno(University of Tokyo) and Mr. Dai Harada(Tokyo Institute of Technology)
for useful discussion.
My thanks are due to all the colleagues in the Kakimoto Laboratory for their encouragement.
The author also would like to appreciate for the continued support of our experiment at Mt.
Chacaltaya by the staff of Institute for Cosmic Ray Research, University of Tokyo.
This work is supported in part by a Grant–in–aid for Scientific Research from the Ministryof Education, Culture, Sports, Science and Technology of Japan.

181

Bibliography[1] Aglietta, M., et al., Astropart. Phys.,10, 1(1999)[2] Agnetta, G., et al., Astropart. Phys., 6, 301(1997)[3] Aharonian, F.A., and Heizelmann, G., Nucl. Phys.(Proc. Suppl.), 60B, 193(1998)[4] Amenomori, M., et al., Ap. J.,461, 408(1996)[5] Antoni, T., et al., Proc. 26th ICRC(Salt Lake City), 1, 131(1999)[6] Apanasenko, A.V., et al., Proc. 25th ICRC, Durban, 5, 13(1997)[7] Apanasenko, A.V., et al., Astorpart. Phys., 16, 13(2001)[8] Asakimori, K., et al. Proc. 24th ICRC(Roma),2, 707(1995); Astrophys. J.,502,278(1998)[9] Axford, W.I., Proc. ICRR int’l symp., Astrophysical Aspects of the Most Energetic Cosmic
Rays, Eds., M.Nagano and F.Takahara(World Scientific), 134 (1991)[10] Baltrusaitis, R.M., et al., Nucl. Inst. Meth., A240, 246(1985)[11] Bernlöhr, K., et al., Astropart. Phys., 5, 139(1996)[12] Bernlöhr, K., et al., Nucl. Inst. Meth., A369, 284(1996)[13] Bernlöhr, K., et al., Astropart. Phys., 8, 253(1998)[14] Biermann, P.L., Astron. & Astrophys.,271,649 (1993)[15] Bierman, P.L., Astron. Astrophys., 271, 649(1993)[16] Bird, J.D. Phys. Rev. Lett., 71, 3401(1993)[17] Boothby, K., et al., Proc. 24th ICRC(Roma), 2, 697(1995)[18] Borione, A., et al., Nucl. Inst. Meth., A346, 329(1994)[19] Cassidy, M., et al., Proc. 25th ICRC(Durban), 5, 189(1997)[20] Cesarsky, C.J., Ann. Rev. Astron. Astrophys.,18,289 (1980)[21] Cesarsky, C.J., Proc. 20th ICRC, Moscow,8,87 (1987)[22] Chi, X. and Wolfendare, A.W., J. Phys. G, 16, 1409 (1990)
182

[23] Cowsik, R. and Wilson, L., Proc. 13th ICRC(Denver), 1, 500 (1973); Meneguzzi, M., Nature
Physical Science, 241, 100(1973)[24] CREAM collaboration, http://cosmicray.umd.edu/cream/cream.html[25] Engler, J, et al., Proc. 26th ICRC(Salt Lake City), 1, 349(1999); Kampert, K.–H., et al.,Proc. 26th ICRC(Salt Lake City), 3, 159(1999)[26] for exmaple, Enomoto, R., et al., Nature, 416, 823(2002)[27] Fenyves, E.J., et al., Phys. Rev. D, 22, 1839(1996)[28] Fortson, L.F., et al., Proc. 26th ICRC(Salt Lake City), 3, 125(1999)[29] Fowler, J.W., et al., Astropart. Phys., 15, 49(2001)[30] Gaisser, T.K., Csomic Rays and Particle Physics (Cambridge University Press) (1990)[31] Garcia–Munoz, M., Mason, G.M. and Simpson, J.A., Ap.J.,217, 859 (1977)[32] Garcia–Munoz, M., et al., Ap.J.,280, L13 (1984)[33] Giacalone, J. and Jokipii, J.R., Ap. J., 430, L137 (1994), J. Geo. Res., 101, 11095 (1996)[34] Glasmacher, M.A.K., et al., Astropart. Phys., 10, 291(1999)[35] Glasmacher, M.A.K., et al., Proc. 26th ICRC(Salt Lake City), 3, 129(1999)[36] Grigorov, N.L., et al., Proc. 12th ICRC(Hobart),5,1760(1971)[37] Gupta, M. and Webber, W.R., Ap.J.,340, 1124 (1989)[38] Guzik,, T.G., et al., Proc. 19th ICRC(La Jolla),2,76 (1985)[39] Han, J.L. and Qiao, G.J., Astron. & Astrophys., 288, 759 (1994)[40] Hayakawa, S., Cosmic Ray Physics(John Wiley & Sons)(1969)[41] Heck, D., et al., Report FZKA(Forschungszentrum Karlsruhe), 6019(1998)[42] Heck, D., et al.,Proc. 27th ICRC(Hamburg),1,233(2001)[43] HESS collaboration, http://www.mpi-hd.mpg.de/hfm/HESS/HESS.html[44] Hillas, A.M., J. Phys. G., 8, 1475(1982)[45] Hillas, A.M., Ann. Rev. Astron. Astrophys., 22, 425 (1984)[46] HiRes collaboration, http://www.hires.physics.utah.edu[47] L’Heureux, J., et al., Nucl. Inst. Meth., A295, 246(1990)[48] Honda, M., Ap. J., 319, 836 (1987)[49] Hörandel, Proc. 27th ICRC(Hamburg),1,71(2001)

183

[50] Inoue, M. and Tabara, H., Pub. Astron. Soc. Japan, 33, 603 (1981)[51] Ip, W.H., and Axford, W.I., Particle Acceleration in Cosmic Plasmas, Eds., G.P.Zank and
T.K.Gaisser(A.I.P. NewYork), 400 (1992)[52] Ivanenko, I.P., et al., Proc. 23th ICRC(Calgary),2,17(1993)[53] Jokipii, J.R., and Morfill, G.E., Ap.J.,290, L1 (1985)[54] Jokipii, J.R., and Morfill, G.E., Ap.J.,312, 170 (1987)[55] Kalmykov, N.N. and Ostapchenko, S.S., Yad. Fiz.,56, 105 (1993)[56] Kakimoto, F., et al., Nucl. Inst. Meth., A373, 282(1996)[57] Kakimoto, F., et al., Nucl. Inst. Meth., A372, 527(1996)[58] Kakimoto, F., et al., J. Phys. G., 9, 339(1983)[59] Kakimoto, F., et al. Proc. 23th ICRC(Calgary),2,103(1993)[60] Kamioka, E., Astropart. Phys., 6, 155(1997)[61] Klages, H.O., et al., Nucl. Phys. B(Proc. Suppl.), 52B, 92(1997)[62] Knapp, J., et al., Astropart. Phys., 19, 77 (2003)[63] Kobayakawa, K., Sato, Y. and Samura, T., Proc. 26th ICRC(Salt Lake City),4,467(1999);
Astropart. Phys., submitted; private comunication[64] Larsen, C.G., Kieda, D.B., and Swordy, S.P.,Proc. 27th ICRC(Hamburg), 1, 134(2001)[65] Linsley, J., Proc. 18th ICRC, Bangalore, 12, 135 (1983)[66] Linsley, J., Proc. 19th ICRC(La Jolla), 3, 461(1985)[67] Longair, M.S., High Energy Astrophysics (Cambridge University Press) (1994)[68] MAGIC collaboration, http://magic.mppmu.mpg.de[69] Mathewson, D.S. and Ford, V.L., Mon. Not. R. Astron. Soc., 74, 139 (1970)[70] Matsuzaki, E., Master thesis(Tokyo Tech.), (2004)[71] Müller, D., et al., Ap. J., 374, 356(1991)[72] Nagano, M., et al., J. Phys. G,10, 1295(1984)[73] Naito, T., and Takahara, F., Prog. Theor. Phys.,93, 287 (1995), Mon. Not. R. Astron.
Soc.,275, 1077 (1995)[74] Protheroe, R.J., and Szabo, A.P., Phys. Rev. Lett.,69, 2885 (1992)[75] Ptuskin, V.S., et al., Astron. & Astrophys., 268, 726 (1990)[76] Rao, M.V.S., and Sreekantan, B.V, Extensive Air Showers (World Scientific), (1998)
184

[77] Risse, M., et al., Proc. 26th ICRC(Salt Lake City), 1, 135(1999)[78] Ryan, M.J., et al., Phys. Rev. D, 40, 729(1989)[79] Senda, D., Ogio, S. and Kakimoto, F., Proc. 25th ICRC, Durban, 2, 269 (1997)[80] Senda, D., Master thesis, Tokyo Tech. (1997)[81] Shirasaki, Y., Mater thesis,Tokyo Tech.(1994)[82] Shirasaki, Y., et al., Astropart. Phys., 15, 357(2001)[83] Shirasaki, Y., and Kakimoto, F., Astropart. Phys., 15, 241(2001)[84] Silberberg, R., and Tsao, C.H., Ap. J., 352, L49(1990)[85] Sofue, Y. and Fujimoto, M., Ann. Rev. Astron, Astrophys., 24, 459 (1986)[86] Stanev, T., Astron. & Astrophys.,274,902 (1993)[87] Swordy, S.P., et al., Ap. J., 349, 625(1990)[88] Swordy, S.P., and Kieda, D.B., Astropart. Phys., 13, 137(2000)[89] Szabo, A.P., and Protheroe, R.J., Astropart. Phys.,2, 375 (1994)[90] Tajima, T. and Shibata, K., Plasma Astrophysics(Addison–Wesley), (1997)[91] Tanimori, T., et al., Astrophys. J. Lett., 497, L25-L28(1998)[92] The Telescope Array Project Design Report (2000)[93] Tsunesada, Y., Doctor thesis(Tokyo Tech.), (2002)[94] Ulrich, H., et al., Proc. 27th ICRC(Hamburg), 1, 97(2001)[95] Völk, H.J., and Biermann, P.L., Ap.J.,333, L65 (1988)[96] Wang, B., Ap. J., 444, 590 (1995)[97] Wefel, J.P., Genesis and Propagation of Cosmic Rays Eds., M.M.Shapiro and J.P.Wefel(NATO ASI series, Reidel), C220, 1 (1988)[98] Yoshii, H., et al., Proc. 24th ICRC(Roma),2,703(1995)[99] Yoshii, H., et al., Il Nuovo Cimento,24C,507(2001)

185