Search for the end of the energy spectrum of primary cosmic rays

The experiments on the search for the end of the cosmic-ray energy spectrum with the particle detector array technique and the optical (Čerenkov and/or fluorescence) technique are summarized. In the highest observed energy region, the suppression of the flux above about 6×1019 eV, which corresponds to the so-called ‘Greisen–Zatsepin–Kuz'min (GZK) cutoff energy’, is confirmed by the new results from the High-Resolution Fly's Eye (HiRes) and the Pierre Auger Observatory (Auger). Summarizing the spectra observed from 1014 to 1020 eV with world wide experiments, there is still room for further investigations on the energy scale and the flux of primary cosmic rays. In order to establish the energy scale, it is necessary not only to understand the optical technique further, but also to determine the primary species in the highest energy region. Also, it is important to make more detailed measurements on the energy spectrum in the 1017–1018 eV region. To observe the super-GZK events and the anticipated recovery of the spectrum above the GZK cutoff, the search for the end of cosmic-ray energy spectrum will be continued even after a century since the discovery of cosmic rays in 1912.


Introduction
Cosmic rays are mainly protons and those in the highest observed energies may be the only samples of extra-galactic material directly detected. It has been of continuous interest to know from how far such samples are collected in relation to unknown physics at these energies as well as their exotic origin.
After our review on ultra-high-energy cosmic rays (UHECR) in 2000 [1], new results from the High-Resolution Fly's Eye (HiRes) experiment [2] and the Pierre Auger Observatory (Auger) [3,4] experiment were published. The most prominent feature of the results from both the experiments is the existence of the suppression of the flux above about 6 × 10 19 eV, which corresponds to the so-called 'Greisen-Zatsepin-Kuz'min (GZK) cutoff energy'. However, there is still plenty of room for further investigation about the energy scale and the flux of primary cosmic rays, especially in relation to the spectrum in the lower energy. Figure 1 shows a compilation of cosmic-ray spectra given by various groups over a very broad energy range. Only some results are plotted in the figure to highlight the main features of the spectrum. One of the mysteries of primary cosmic rays is that their energy spectrum extends over more than ten decades of energy almost continuously with only a few small changes in the slope in a power-law energy spectrum. It is hard to imagine any kind of astronomical accelerator that can cover such a broad energy range. Hence, we need to find changes, if any, in their composition, arrival directions and slope of energy spectrum with energy to solve the mystery. The energies where the slope changes are observed are around 4 × 10 15 eV (4 PeV), 6 × 10 17 eV (0.6 EeV) and 4 × 10 18 eV (4 EeV), which are popularly known as the knee, the second knee and the ankle, respectively.
In this paper, the experiments on the search for the end of the cosmic-ray energy spectrum, with the particle detector array technique and the optical (Čerenkov and/or fluorescence) technique are summarized. The energy spectrum in the highest energy region is discussed in relation to that in the lower energy region. A recent review of the composition and the anisotropy of UHECRs may be found in the Rapporteur paper of Teshima [5] presented at the 30th ICRC in Mérida, Mexico in 2007.

Historical background of UHECR experiments
There is a review on the history of cosmic-ray studies by Linsley [6], in which the early experiments to search for the composition and the end of the energy spectrum are described in detail, starting with the discovery of cosmic rays by Victor Hess in 1912. The early air fluorescence work at Cornell and in Japan is described by Tanahashi [7] in the proceedings of the same conference.
In figure 2 various experiments on UHECRs are listed along with their main observation periods. Here, we review first the history of the particle array technique in the early stages and

Particle detector array technique
During the period 1947-55, the MIT group [8] developed techniques of 'density sampling' and 'fast timing', in which the pattern of densities observed with an array of scintillation counters was used to locate the core of the extensive air shower (EAS) and its arrival direction was determined from the differences among the arrival times of shower particles at various counters. Between 1954 and 1957 an array of 15 counters, each 1.0 m 2 area, was operated at the Harvard Agassiz Station [9] and the energy spectrum up to 10 18 eV was derived. A liquid scintillator of about 25 gallons (about 100 l) of toluene was used as a counter. Unfortunately, one of the counters burst into flames due to lightning. Accordingly, a large-area plastic scintillator was 5 developed and the liquid scintillators were replaced by the plastic ones in 1956. The largest shower recorded in this experiment is about a few times 10 18 eV [8] as shown in figure 2 by an arrow. In this experiment, the feasibility of the fast timing method for determining the arrival direction of showers was established. Also, the use of not only the analog device (computer) but also digital computation for determining the core locations was demonstrated successfully. These methods are the basis of particle array experiments afterwards.
After the Agassiz experiment, Linsley et al [10] constructed a large array at Volcano Ranch in New Mexico and extended the primary energy spectrum to above 10 19 eV. The array had nineteen 3.3 m 2 plastic scintillation counters on the ground and those shielded by 10 cm of lead. The signals were displayed on an oscilloscope for the measurement of pulse amplitudes and relative arrival times of muons and all charged particles separately. With this array, Linsley [11] discovered the first candidate event with an estimated energy of 10 20 eV in 1962 during 3 years of full-scale operation. It should be noted that this event was observed before the discovery of the 2.7 K background radiation [12] in 1965 and the prediction of a cutoff (GZK cutoff) in the primary cosmic-ray energy spectrum by Greisen [13], and Zatsepin and Kuz'min [14] in 1966.
Soon after the prediction of this GZK cutoff, two large arrays, at Haverah Park in UK and at Narrabri in Australia Sydney University Giant Air-Shower Recorder (SUGAR) started their operations from 1967. In figure 2 are shown the years of steady operation in full scale for each experiment.
The main detectors at Haverah Park were waterČerenkov detectors laid out over an area of about 12 km 2 [15]. Each tank was 2.25 m 2 in area and 1.2 m in depth and 15 tanks (total area of 33.75 m 2 ) were deployed in each station. The communication of the coincidence trigger was made to the outer stations via a 7.0 GHz microwave link. This array was operated in stable mode for about 20 years without any deterioration in the water quality. This interesting observation may be one of the reasons for the selection of the waterČerenkov detector as the surface detector for the Auger project later. Four events whose estimated energy was above 10 20 eV were reported [15]. However, the energies of these events were reduced to below 10 20 eV after re-evaluation of shower energies with the newly developed simulation code [16].
SUGAR started its operation at Narrabri in 1968. This array had 56 stations, each of which consisted of a pair of 6 m 2 area scintillators buried 1.7 m below the ground [17]. Though the area covered was over 60 km 2 , the spacing between the stations, typically 1.61 km, was too large to determine accurately the core location, the energy and the arrival direction of showers with a large enough number of hit stations. Furthermore, after-pulsing in the 7-in diameter photomultipliers was a serious problem throughout the lifetime of the experiment, as the logarithmic pulse-width to pulse-height converters were used for estimating particle densities. Nine events whose energies were estimated to be nearly above 10 20 eV were reported [18]. Zenith angles of two events were about 70 • and their core locations were outside the boundary of the array. In three other events, the number of hit stations was only three and the detector response of one of the hit stations was saturated. These events may not be considered to be definitive candidates of events with energy above 10 20 eV. The primary energies of four other events were re-evaluated by Nagano using the muon lateral distribution determined with the Akeno Giant Air Shower Array (AGASA) [19]. Two events (SU14427 and SU6179) may be 10 20 eV candidates, if the extrapolation of AGASA muon lateral distribution to the atmospheric depth of individual SUGAR events is valid [20].
At Tokyo, Suga et al constructed the INS-LAS array in 1968, which was large compared with the previous INS array. However, it was relatively small compared with those at Haverah 6 Park and Narrabri. The primary purpose of this array was to use the decoherence method or the density spectra at different separations of detectors to estimate the primary energy spectrum. They reported the detection of a 4 × 10 21 eV cosmic ray in 1970 [21], whose energy was estimated from the total number of charged particles by fitting the lateral distribution given by Linsley [11]. If we use the AGASA lateral distribution of charged particles and muons [19] after taking into account the difference of the Mòliere unit between the INS and Akeno altitudes, and use the particle density per m 2 at 600 m from the core as an energy estimator, the energy of the INS-LAS event is reduced to about 2.5 × 10 20 eV [20]. However, it is still one of the highest energy events reported so far.
The Yakutsk array in Siberia is the most complex of the giant arrays. This array began taking data in 1970 and was expanded to cover an area of 18 km 2 in 1974 [22]. In addition to the surface detectors, many underground detectors are also deployed in this array. A particularly important feature of this installation has been the presence of 35 photomultiplier systems to measure the airČerenkov radiation associated with showers. This feature is useful for providing a calorimetric approach to the calibration of the primary energy experimentally, since the energy estimation of other particle array experiments has been based only on simulation studies. One event above 10 20 eV whose zenith angle is 58.9 • is reported [22]. In 1995 the Yakutsk array was re-arranged to cover an area of only 10 km 2 , so that detailed studies on shower structure could be made around 10 19 eV, where they have reported observing a change of shower characteristics [23]. The array is still in operation in 2008.
In the particle detector array technique, we should remark that the fluctuation of the density of shower particles far from the core is quite small irrespective of the hadronic interaction model or the chemical composition of the primary particle and hence the density far from the core is considered to be a good estimator of primary energy. This important insight was first pointed out by Hillas [24] in 1970 and it has been supported by many simulation studies afterwards [25]- [29]. Without this idea, the estimation of the energy of the primary particle has never been successful in the large ground array experiments, where detectors are deployed with large separation. The Haverah Park group first used this parameter and estimated primary energies of the observed showers [15].

Fluorescence technique
The possibility of using the earth's atmosphere as a vast scintillator was pointed out at the Norikura Symposium in 1958 [7]. The pioneering discussions were led by Suga and Chudakov at the 5th Inter-American Seminar at La Paz, Bolivia [30]. Based on the experimental measurement of the fluorescence efficiency in air produced by a deuteron beam [31], the first station was constructed on the hill of Mt Pleasant near Ithaca in 1964 by Greisen and his group [32]. With two additional stations, forming a triangle with distances of 11, 16 and 12 km between them, the experiment was continued till 1967. In 1965, Greisen proposed a new system to divide the celestial sphere into about 500 mosaic segments and to record the light intensity and its arrival time in each segment. This imaging system with 25 sided observatory was constructed in 1966 [33]. However, due to the severe climate of Ithaca in winter and the small statistical chance to observe very high-energy showers, the observations were unable to catch a clear signal of the fluorescent light generated by high-energy showers. However, this idea of a new detection technique may be the basis of the fluorescence detectors developed afterwards.

7
After returning from the Cornell experiment, Tanahashi started another experiment in 1968 at Mt Dodaira near Tokyo, using a Fresnel lens of 1.6 m diameter and 27 photomultiplier tubes (PMTs) at the focal plane [7]. Within 6000 events recorded during 90 h, his group luckily succeeded in detecting the first fluorescence light from an EAS of energy above 5 × 10 18 eV [34]. This event was clearly distinguishable from otherČerenkov light events, since the angular velocity of signals on each PMT was measured to be about 0.08 rad/0.45 µs and the light intensity on each PMT was almost constant. These features could be interpreted as due to a shower at a distance of 1.7 km if the shower arrived at right angles to the line of sight and the shower was nearly at its maximum development. Encouraged with this event, his group constructed another two types of telescopes. However, further development was not successful because the plastic lens newly manufactured was mixed one with UV stopper [7].
In 1974, the proposal of Fly's Eye was announced by the University of Utah group. In 1976, an unambiguous detection of fluorescent light emission from showers, in coincidence with the signal from a large particle detector array, was confirmed [35] at Volcano Ranch. This success led to the construction of the Fly's Eye detector, which started operation in 1981 [36] at Dugway near Salt Lake City. A few years later, two detectors, Fly's Eye I (FE I) and Fly's Eye II (FE II), separated by a distance of 3.4 km, became operational. The FE I consisted of 67 mirrors of 1.6 m diameter and a total of 880 photomultipliers was used to cover the whole sky. The observations continued for more than 10 years. However, a lot of pioneering work, such as the development of the optical calibration system [37] and the atmospheric monitoring methods [38], was required to obtain fruitful results on the energy spectrum of primary cosmic rays, measurement of the depth of the maximum shower development and so on. This may explain the fact that though the highest-energy event so far observed was detected in 1991, its details were published only in 1995 [39]. More significantly, the Fly's Eye experiment showed the importance of stereoscopic observations and laid the foundation of both the experimental and the analysis methods.

Highest energy region
The newly published energy spectra from the HiRes [2] and Auger [3,4] experiments are compared with those published earlier (after 1990) in figure 3, where the differential flux is multiplied by an energy-dependent power E 3 to reduce the steepness of the power-law spectrum. A prominent feature of the HiRes and Auger experiments is the suppression of the flux above the energy corresponding to the GZK cutoff with a statistical significance of 5.3 and 6 standard deviations, respectively. The HiRes energy spectrum is based on monocular data from HiRes-I and HiRes-II [2]. The Auger spectrum is the combined spectrum from 3 datasets, namely, vertical events observed by the surface detectors (SDs), inclined events observed by SDs and hybrid events detected by both the SDs and the fluorescence detectors (FDs). For every dataset the energy is calibrated using FDs.
The HiRes and the Auger experiments observe a flattening of the spectrum at an energy of 4.5 × 10 18 and 4 × 10 18 eV, respectively (ankle). It should be noted, however, that the ankle energy, E ankle , is not determined accurately, either by HiRes-I or HiRes-II independently. Similarly, in the case of the Auger experiment, the ankle energy cannot be derived for hybrid events or SD vertical events independently. They are determined by combining datasets from different measurements; in the case of HiRes, HiRes-I and HiRes-II, and in the case of Auger, hybrid events and SD vertical events.  Figure 3. Energy spectrum of primary cosmic rays from 10 17 to 10 20 eV. The differential flux in each bin is multiplied by an energy-dependent power E 3 .
In the same figure, Fly's Eye stereo spectrum [40] and HiRes stereo spectrum [41] are also plotted, which may have better energy resolution than the monocular data. These spectra cover the E ankle region with their own datasets. The ankle energy obtained from the Fly's Eye stereo data is 3.2 × 10 18 eV [40] and that from the HiRes stereo data is estimated from the figure to be 5.6 × 10 18 eV.
The AGASA spectrum [42] is plotted in the figure by reducing its energy by 10% based on the combined experiment with the Akeno array and the prototype AGASA array which is described later in section 4.1.3. This spectrum is in good agreement with the Akeno spectrum [43] in the overlapping energy region. The particle density at 600 m from the core, S(600), has been used as an energy estimator by AGASA and the conversion to primary energy is based on simulations [27]- [29]. The Akeno energy spectrum is based on the total number of shower particles N e and the conversion to primary energy is based on experimental data on the longitudinal development curves measured at Chacaltaya and Akeno [44].
The HiRes-MIA spectrum is determined with a hybrid detector consisting of the HiRes prototype detector and the Michigan muon array (MIA) [45]. By using the hybrid timing information, the geometrical reconstruction is improved and hence the energy determination from the longitudinal development observed by the HiRes prototype is more accurate.
The Haverah Park spectrum is the re-analyzed one using the QGSJET interaction model with the CORSIKA code [46]. In the figure, the case of mixed composition (34% protons and 66% iron nuclei) is plotted. A point at log E(eV) = 19.9 represents four re-calculated events whose energies were estimated to be larger than 10 20 eV in the original analysis [15]. The spectrum below 10 18.6 eV is in good agreement with the Fly's Eye and HiRes results. However, the points around 10 18.9 and 10 19.9 eV are quite high (nearly in agreement with Yakutsk spectrum) and the results from the re-analysis between those energies have not been reported. The Yakutsk spectrum reported in 2004 [47] is also plotted in figure 3. The energy parameter, S(600), the scintillator particle density at 600 m from the core, is calibrated experimentally with the calorimetric method by measuring the airČerenkov radiation. By adding the energy carried by the electromagnetic component and the muons below the ground and the unobserved portion (8%), the relation between S(600) and the primary energy is determined. The fluxes of the Yakutsk spectrum are quite high compared with the other results, however, they claimed the suppression of the spectrum above the E GZK cutoff in 1991 [48].
The flux values given by the Auger experiment are smaller compared with the other experiments. Though the spectrum is determined by the particle detector array data, the energy calibration is obtained from the fluorescence detector using the calorimetric method, by integrating the observed ionization energy deposition of shower particles in the atmosphere. Details of the comparison with other experiments are discussed later in section 4.2.2.
Spectral slopes below and above the E ankle of different experiments are listed in table 1. Slopes below E ankle are determined in the range between E th and E ankle . Though the differences in the spectral slopes among experiments are larger than the determination error of slopes in each experiment, the differences are relatively smaller than those in flux and E ankle . The average slopes of these experiments are 3.24 ± 0.07 and 2.75 ± 0.05 below and above E ankle , respectively, and may be accepted in the general discussions.
The E GZK value represents the suppression energy above which the flux decreases as expected from the GZK prediction.

Comparison with the spectrum in the lower energy region
Before discussing the differences among the experiments in the highest energy region, it is important to relate them to those in the lower energy region. In figure 4, the energy spectrum of primary cosmic rays from 10 13 to 10 20 eV is drawn. The differential flux in each bin is multiplied by an energy-dependent power E 3 as in figure 3.
In the energy region below the knee, the energy spectra from the satellite observations by Proton 3 (1970) [49,50] and SOKOL (1993) [51] and the balloon observations by JACEE (1995) [52] and RUNJOB (2005) [53] are plotted. In these experiments, the nuclear species of individual events are determined and the energies are measured in a relatively shallow calorimeter. However, the fraction of energy deposited in the calorimeter is subject to large fluctuation due to the fluctuation in the inelasticity in the nuclear interaction of the primary particle which determines the energy transferred to the electromagnetic component. Since the flux around the knee region is very small, the width of the energy bin in each experiment is quite large due to the fluctuations in the observed energies. The flux observed by RUNJOB is significantly lower than others. This is because the helium flux of RUNJOB is significantly lower compared with the other experiments.
The knee energy, E knee , and the spectral slopes below and above the knee are listed in table 2, together with the energy ranges over which the fits are made to determine the slopes. Though the fluxes and the slopes below E knee are nearly in agreement among the experiments and are smoothly connected to the satellite and balloon experiments, the differences in the values of E knee , the fluxes at E knee and the slopes above E knee are rather large, much beyond the experimental uncertainty. Considering the spectrum slope (∼ 3.1) above E knee and the slope (∼ 3.25) below E ankle , there must be a '2nd knee' between them. Spectral slopes below and above the 2nd knee from Akeno [43], Haverah Park [46], Fly's Eye and HiRes-MIA [45] are listed in table 3. Also listed in the table are the energy of the 2nd knee energy, E 2nd , and the energy range over which the fit is made to determine the slope. E 2nd may be between (4-6)×10 17 eV.
Though there are differences in E 2nd and the fluxes among the experiments, the spectra given by AGASA, Fly's Eye stereo and HiRes-II are smoothly connected to those in the lower energy region. On the other hand, the flux given by the Auger experiment seems to be lower and that of the Yakutsk is larger when trying to connect with the spectrum in the lower energy region. Energy spectra from 10 14 to 10 20 eV from world data from the fluorescence andČerenkov experiments. The Auger combined spectrum from the three datasets is drawn since the energy is calibrated using the FDs. Yakutsk spectrum is also shown here as its energy is calibrated by measuring the aiř Cerenkov radiation.
simulations codes, such as MOCCA [63], COSMOS [64], CORSIKA [65,66] and AIRES [67] have been developed and CORSIKA may be the most widely used code in recent experiments. Various hadronic interaction models, QCDJET [68], QGSJET [69]- [71], VENUS [72], DPMJET2 [73], DPMJET3 [74], SIBYLL [75]- [77], HDPM [65] and EPOS [78], are used for particle interactions in the high-energy region, and GHEISHA [79] and FLUKA [80,81] in the low-energy region in many experiments. Most of these codes have not remained fixed and have been frequently revised by the authors. The version of each code is not always mentioned in the discussion below and it may be necessary to refer to the citation for the experiment to know the code version and other details. The summary of energy spectra determined by particle detector arrays on the surface is shown in figure 5.
The MAS array of BASJE (BASJE-MAS) at Mt Chacaltaya is at 5200 m above sea level (a.s.l.) (550 g cm −2 ) and the showers with energies around the knee region are near their maximum development. The shower size at maximum, N max e , does not depend much on shower development fluctuations and/or different primary composition. By comparing the equi-intensity curves for various zenith angle bins with the simulated ones with five primary components and using CORSIKA code with the QGSJET model, the BASJE group determined both the energy spectrum and the energy dependence of the primary composition [60].  The Tibet array at 4300 m a.s.l. is also near the expected maximum for showers of energies around the knee. The Tibet spectra are obtained after converting the N e spectrum using the QGSJET01c or SIBYLL2.1 interaction models with heavy dominant (HD) or proton dominant (PD) composition model (the proportion of various nuclei changes with energy, refer to [62] for details). In figure 4, only the results obtained from a combination of QGSJET and HD are shown, but three combinations are shown in figure 5 and five combinations are listed in table 2.
There are differences of about 10% in energy between the BASJE and the Tibet spectra, however, they agree with the direct observations below the knee energy except the RUNJOB result.
In the case of HEGRA [58], the slope of theČerenkov light lateral distribution is used to estimate the height of the shower maximum above the observation level for each shower, then the conversion from N e to N max e is made. The conversion from N max e to primary energy is based on simulations using CORSIKA v5.20 with the QGSJET/GHEISHA option. The bestfit spectrum obtained from the spectra with the four energy reconstruction methods used are plotted in the figure.
EAS-TOP measures N e spectra for different zenith angle (θ ) bins and the conversion from N e (θ) to primary energy and mass has been made using the results from simulations with the CORSIKA-HDPM code [54].
CASA-MIA measures not only N e , but also the total number of muons N µ above 750 MeV (vertical). They found that the parameter N e + 60 × N µ is logarithmically linear with E and it is independent of the primary composition, based on simulations using the MOCCA code and the SIBYLL hadronic interaction model [55].
KASCADE measures N e and the truncated muon number N tr µ , which is the number of muons (>230 MeV) within core distances between 40 and 200 m [61]. Using the CORSIKA code with the high-energy interaction models QGSJET01 and SIBYLL2.1 and the low-energy interactions GHEISHA 2002, the primary energy spectra for five primary mass groups are determined.
At Akeno there are arrays of various detector spacings depending on the primary energy of interest, and the energy spectrum has been determined systematically over five decades in energy (10 15 -10 20 eV) using similar experimental procedures [43,44]. In the case of EAS hitting the compact arrays, N e is used to determine the primary energy. The relation for converting N e to N max e was obtained from the longitudinal development curves which were determined from equi-intensity cuts on integral N e spectra at various zenith angles at Chacaltaya and Akeno [44]. The conversion from N max e to primary energy was based on simulations by Hillas [82]. N e spectrum of the Moscow State University (MSU) is determined with the array of Geiger counters and hence the particle number observed is free from the definition of a single particle and the transition effect (electromagnetic cascade of the incident particles) in the detector. However, the conversion from N e spectrum to primary energy spectrum drawn in figure 5 [83] and the reasons for the change from the previous spectrum [84] have not been described in these references.
4.1.2. Conversion from particle density at a certain distance from the core. To estimate the primary energy of giant air showers with an array of detectors spread on the ground, one fits the observed particle densities to a lateral distribution function. Using this fit, the particle density expected at a certain distance from the core is determined and this density is used as an energy estimator, as described in section 2.1.
The detectors of the Haverah Park array were made of waterČerenkov tanks of 2.29 m 2 area and 1.2 m depth. Fifteen tanks were used at most of the stations. The signal density at a core distance of 600 m, ρ(600), was used as an energy estimator [15]. The energy spectrum was originally estimated from simulations using a shower model developed by Hillas [25,26]. However, the energy has been re-evaluated using the QGSJET98 interaction model in the CORSIKA code together with GEANT to simulate the detailed detector response to particles of different types and energies [46]. The energy spectra obtained for two different assumptions for the primary composition, namely, pure proton or pure iron, are shown in figure 5.
AGASA used the scintillation detectors of 2.2 m 2 area and the signal at 600 m from the core S(600) was used as an energy estimator. The energy conversion relation was initially derived using the COSMOS code and the QCDJET interaction model [27]. Later, the relation was re-calculated using the CORSIKA(v5.623) code with QGSJET98 and SIBYLL1.6 models [28] and also the AIRES(v2.2.1) code with QGSJET98 and SIBYLL1.6 models [29]. The summary of the results from these simulations are given in Takeda et al [42]. The systematic variation due to simulation codes, interaction models and mass composition is found to be within 10% among the combinations investigated. This small difference among the results for various interaction and composition models is clearly an advantage for the method of measuring S(600) using scintillators, where the observed particles are dominated by the electromagnetic component with only a small contribution due to muons.
The surface detectors of the Auger observatory are 10 m 2 area and 1.2 m depth wateř Cerenkov detectors, which are deployed on a 1500 m hexagonal grid. The covered area has been gradually increased and the results from an exposure about 4 times larger than AGASA are reported. The signal observed by the water tank at 1000 m from the shower core, which is due to the electromagnetic and muon components, is used as an energy parameter. Though the signal at 1000 m depends only slightly on the mass of the primary particle and the hadronic interaction model used for simulations, the fraction of signal due to muons is relatively large. Engel et al [85] studied the signal size at a zenith angle of 38 • and at 1000 m from the core (S 38 (1000)) through simulations with CORSIKA 6.5 and the hadronic interaction models QGSJET II.03 and FLUKA. By obtaining the best rate for the contribution of muons with the constant intensity cut method, they determined a measure of the energy scale of S 38 independent of the response from the fluorescence detector. This corresponds to assigning showers a ∼ 30% higher energy than the reported assignment based on the fluorescence detector at 10 19 eV (E SD = 1.3E FD ).
The spectrum shown in figure 5 is the one where the energy has been multiplied by 1.3 compared with the Auger combined energy spectrum in figure 3, assuming a constant factor 1.3 throughout the energy range measured. This modified spectrum is in good agreement with the re-evaluated Haverah Park spectrum below 10 19 eV and the AGASA result above 10 19 eV, except for energies beyond 10 20 eV. This may suggest that if we use the results from Monte Carlo simulations with the similar hadronic interaction model to estimate the primary energy from particle density such as S(600) or S 38 (1000), the primary energy is nearly in agreement among these experiments using surface detector arrays.
It should be noted, however, that using this factor 1.3 for the energy scale, the number of muons measured by Auger is about 1.15 times larger than the predictions for iron showers with hadronic interaction models used for simulations [85]. If a new hadronic interaction model EPOS is used for air shower simulations, many more muons are produced than QGSJET II.03 [86] and hence a factor 1.3 may be reduced. It is desirable to determine the energy scale of S 38 and its energy dependence from simulations using a hadronic interaction model which reproduce closely the measured muon number.
It is also remarked that the muon density at 1000 m from the core observed by AGASA could be interpreted by primaries between proton and iron with QGSJET01 or SIBYLL2.1 [87]. N e and S(600). A comparison has been made in the 10 18 eV energy region [43] between the energy determined from N e for the Akeno 1 km 2 array (E 1 ) and that using the S(600) for the 20 km 2 array (E 20 ), whose detectors were deployed with about 1 km separation and which was the prototype array of AGASA. The ratio E 1 /E 20 is found to be 1.10. Since the Akeno spectrum is in good agreement with the other experiments around the E knee , where the spectrum is smoothly connected to the spectra by the direct observations, the AGASA energy has been reduced by 10% by normalizing it to the Akeno energy.

Possible reasons for the difference in spectra among the experiments.
Though various methods for conversion from observables to primary energy are used at different altitudes of the experimental sites using different Monte Carlo simulation codes with different interaction/composition models, the differences in the fluxes and slopes above the E knee are not large among various experiments. This may be because the exposure (aperture × observation time) does not depend on the energy in most of the experiments and the simulation of the electromagnetic component near the shower maximum does not depend much on the hadronic interaction and/or primary composition models. However, there are some differences among experiments as follows: 1. E knee and the flux above E knee of CASA-MIA and BASJE-MAS are significantly lower than those of other experiments. 2. There are some differences in the energy and flux in the energy region, 10 18 -10 19 eV, between the AGASA and the re-evaluated Haverah Park spectra. 3. E 2nd is clearly observed around 6.3 × 10 17 eV by the AGASA, whereas E 2nd is not observed in the re-evaluated Haverah Park spectrum down to 3 × 10 17 eV. 4. There are some AGASA events above 10 20 eV.
In the following, possible reasons for these observed differences are discussed.

Single particle.
In most experiments the particle number incident on the detector is measured in units of a signal penetrating the scintillator or water tank vertically. The definition of 'single particle' depends on the experiment, an average value or a peak value of the pulse height distribution due to omni-directional particles or vertically incident particles selected by a particle telescope. Recently, the signal due to single particle defined in each experiment is evaluated using the well-known simulation code, GEANT4 [88,89], taking into account the transition effect in the detector and its surroundings as used in the experiment. Then the total number of shower particles N e or the lateral distribution can be determined in units of the single particle as defined in each experiment. However, it is not certain whether the effects due to the incidence of multiple particles or small showers, whose contribution to the pulse height distribution depends on the experiment, are taken into account in determining the 'single particle' signal. It should be remarked that there remains the effect of cascade showers from pair creation or bremsstrahlung of muons in the atmosphere even we select vertical muons with a telescope, whose proportion depends on the area of the detector, solid angle of the telescope and the altitude of the experimental site.
For example, the number of photoelectrons observed for a vertically passing particle through the water tank of Auger is about 60% larger than that obtained from simulations using the GEANT4 package [90]. Though this difference is explained by the authors as due to the PMT quantum and collection efficiencies and the optical coupling of the PMTs to the liner window and the contribution of electromagnetic component is relatively small in the WC detector, the effect due to multiple particle incidence should be measured experimentally since the detector area is large.

Shower front structure.
The arrival time distribution of signals on a detector is recorded, but the number of particles is integrated within some limited time. Therefore the overestimation or underestimation due to the shower front structure must be taken into account. The arrival time distribution of shower particles at distance r from the core is expressed by where the scaling parameter t 0 is a function of the distance at the observation level [91]. In the case of logarithmic amplifier with a decay constant of 10 µs as used in AGASA, S(600) may be overestimated by about 5±5% due to the broadening of the shower front structure [42]. Since this estimation is based on the extrapolation of a limited number of events near 10 20 eV, further study of the shower front structure is expected with the telescope array (TA) experiment [92] with accumulated events. On the other hand, S(600) may be underestimated in the case of Yakutsk, since the integration time of Yakutsk detectors is 1 or 2 µs on different detectors. This is too short to record the full signal at large axial distances [93]. Additionally, it seems probable that the short coincidence window of 1.2 µs between two 2 m 2 detectors spaced 90 cm apart may lead to a loss of triggers for thick shower front structure at larger core distance [93]. Though this criticism may explain the suppression of flux above 5 × 10 19 eV of Yakutsk, the systematically large flux of Yakutsk in the 10 18 -5 × 10 19 eV region cannot be explained.

Delayed particles.
One of the differences between a scintillation detector used by AGASA and the water tank used by Haverah Park and Auger is the response to the neutrons. The scintillator is sensitive to delayed neutrons, while water tank is not sensitive to scattered protons created by the neutrons. Since the AGASA group used the logarithmic amplifier, the signal amplitude is sometimes overestimated due to the arrival of delayed neutrons [42]. Considering the effect of the delayed neutrons on the energy estimation of 10 20 eV candidate events experimentally, the over-estimation is evaluated to be 0-10%, independent of primary energy [42]. If there remain effects of delayed neutrons in the signal, a possible flattening of the observed lateral distribution may be observed, however, such flattening is not observed up to 3 km from the core and up to 10 20 eV [42]. Since the scintillator is used as SD in the TA experiment [92], it is expected that the rate of delayed neutrons will be measured for events of several ×10 19 eV, whether the AGASA 10 20 eV events are due to the overestimation of S(600) caused by delayed neutrons.

4.1.4.4.
The dynamic range of the detector. The particle number incident on the detector increases with energy and relatively very high particle density occurs near the core where the shower front structure is very thin. Therefore, a large dynamic range of the detector within 10-100 ns including the PMT response, amplifier and its recording system is quite important for measuring high densities near the core as well as low densities far from the core.
In the case of AGASA, a logarithmic amplifier was used to cover a dynamic range from 0.3 to a few times 10 4 particles per detector, where the upper limit is due to the saturation of PMT response [94]. The number of saturated signals was negligibly small for the observed events up to the highest energy. On the other hand, the saturation of the FADC for the anode response in Auger SDs occurs below 10 3 particles per detector near the core [95]. Though the number of saturated detectors is only one or two near the core for events below 10 20 eV and the recovery methods from the saturated signal to real number are tried, the broadening of the dynamic range by more than a decade is highly expected to record an event above 10 20 eV unambiguously.

Optical technique
The energy spectra from 10 14 to 10 20 eV, in experiments where the energies are calibrated using the optical techniques are summarized in figure 6. In the energy region lower than a few 10 17 eV, the energy is determined by measuring the lateral distribution ofČerenkov light. The differences among the experiments above the E knee seems to be larger than that observed among the particle detector array experiments.
In the energy region above 10 17 eV, the fluorescence technique is used except for the Yakutsk experiment. Though the energy spectrum by the Auger is combined one from SD vertical, SD inclined and hybrid, the result is plotted in this figure, since the energy is calibrated by the fluorescence technique.
The energy spectrum of Yakutsk is derived from the S(600) spectrum. However, the relation of S(600) to primary energy is determined experimentally by the calorimetric method of measuring the airČerenkov radiation (70-80% of the total observed energy), the electromagnetic component and the muons on the surface. Therefore the Yakutsk spectrum is also plotted in figure 6.

4.2.1.Čerenkov lateral distribution.
The BLANCA array consisted of 144Čerenkov light detectors arranged with an average separation of 35-40 m. It measured the lateral distribution of Cerenkov light in coincidence with the CASA scintillation detector array at Dugway, USA [57]. The energy of each shower was derived usingČerenkov density at 120 m from the shower core and the slope parameter of theČerenkov lateral distribution between 120 and 350 m. The results from the simulations with the CORSIKA code and four hadronic interaction models, QGSJET, VENUS, SIBYLL and HDPM were used to estimate the energy spectrum. In figure 6, the results based on the QGSJET model are plotted. DICE data are from two telescopes located at the CASA-MIA site, also at Dugway, Utah [56]. Each consisted of a 2 m diameter spherical mirror with a focal plane detector of 256 photomultipliers (PMTs). The telescopes were separated by 100 m. Cosmic-ray events within the field of view produced a focal plane image at the PMT cluster corresponding to the direction and intensity ofČerenkov light coming from the air shower. A combination of the amount of Cerenkov light and the depth of the shower maximum was used to estimate the primary energy. The CORSIKA code with VENUS interaction model was used to derive theČerenkov size function [56].
Inside the HEGRA particle detector array, 49 open PMTs were arranged on a grid with 30 m spacing called AIROBIC array [58]. TheČerenkov density at core distance of 90 m (L 90 ) was used as an energy parameter. The energy spectrum was reconstructed by L 90 alone without using a particle detector array. Since L 90 spectrum is within the systematic and statistical uncertainty of the best-fit spectrum shown earlier in figure 5, that spectrum is also drawn in figure 6.
At Yakutsk the smallČerenkov array of ∼3 × 10 5 m 2 area was added in the central part of the large Yakutsk array [59]. EAS with energies above 10 15 eV were recorded for about 2716 h. Cerenkov light density at 100 m from the core Q(100) was used as an energy parameter. The Yakutsk spectrum above 10 17.5 eV is described in section 3.1. The spectrum obtained from the smallČerenkov array is in agreement with the HiRes-MIA, the Fly's Eye and HiRes spectra in the overlapping energy region. However, the Yakutsk spectrum is quite high at energies above 10 17.5 eV compared with other spectra. The reasons for the difference in the Yakutsk spectra between the low-and high-energy regions have not been described in any paper by the Yakutsk It is true, in principle, that the energy estimation by the fluorescence technique is most reliable since the primary energy is obtained experimentally by integrating the energy deposited in the atmosphere. Therefore the energy estimation does not depend on the hadronic interaction model or the primary species. However, there are several other factors in the determination of the energy experimentally in optical methods. This may be why there are some differences in the energy spectrum between Fly's Eye/HiRes and Auger.

Important items relevant to the optical technique.
To compare the Fly's Eye, HiRes and Auger results, the important items related to the estimation of the deposited energy by the fluorescence technique are compared in table 4.

Single photon.
Relative number of photons observed by the system are carefully calibrated routinely in every experiment.
To compare the energy spectra determined by the optical methods, the determination of the signal size due to a single photon may be the most important. In the laboratory experiment, a photomultiplier whose thermal noise is small enough may be selected for the necessary discrimination and the single photon counting method, irrespective of the signal amplitude, 20 can be applied in many experiments, for example [97]. Even in such an arranged experiment, the absolute photon flux is subject to a systematic uncertainty of more than 10% at present.
In fluorescence experiments in field, the number of photons incident on the photomultiplier may be counted in units of the average of the single photon amplitude distribution, where the noise level may be larger than the threshold level of single photon amplitude. The absolute number of photons observed by the photomultiplier is calibrated with a silicon photodetector, whose efficiency has been calibrated at the National Institute for Standard in USA, for both Fly's Eye/HiRes [96] and Auger [98]. However, the sensitivity and the size of the photodetector calibrated are quite different from the photomultipliers used in the experiment. The average depends on the photomultiplier used and the definition of a single photon is supposed to be different in each experiment by some amount. This may be a part of possible difference among the fluorescence experiments. It is highly desirable to calibrate the photomultipliers used in the field by measuring the amplitude distribution from nitrogen fluorescence photons in the laboratory.
InČerenkov measurement of Yakutsk and theČerenkov experiments in the knee region, the definition or the determination of the absolute number of photons is not described in detail. Since the observedČerenkov signal is quite large compared to the noise level, the absolute value seems to be determined less carefully than in the fluorescence experiments. This may be one of the reasons for rather large discrepancies amongČerenkov experiments.

Fluorescence yield.
The fluorescence yield from the EAS particles was quantitatively studied and summarized by Bunner in his PhD thesis [99] and his results have been used for the Fly's Eye and HiRes analysis. The altitude dependence of the fluorescence yield and its energy loss dependence were examined by Kakimoto et al [100] and their results have been used by the HiRes. Since the wavebands measured were limited only to 337, 381 and 391 nm in the studies of Kakimoto et al and Nagano et al [97] re-measured the fluorescence yield in 14 wavebands. The attenuation for each wavelength must be taken into account for the distant EAS, since the attenuation by Rayleigh scattering depends on the wavelength λ −4 .
Afterwards, the importance of more accurate measurements on the fluorescence yield for energy determination of EAS by fluorescence technique has been recognized and its dependence on pressure, temperature and humidity have been measured by several groups. Recent developments were reported at the 5th Fluorescence Workshop held at El Escorial in 2007 and their summary by Arqueros et al [101] is found in the Proceedings of the Workshop.
In determining the primary energy with the fluorescence detector at Auger, the absolute fluorescence yield of the 337 nm band is taken from Nagano et al [97] and the relative yields at other bands and their pressure dependence are taken from Ave et al [102].
The fluorescence efficiency, defined as the radiated energy divided by the energy loss in the observed medium, summarized by Bunner [99] and measured by Kakimoto et al [100] and Nagano et al [97], is compared for three main wavebands in table 5. It should be noted that Bunner estimated the values as a set of weighted average of three experiments, Bunner [31], Davidson and O'Neil [103] and Hartman [104]. These experiments used different methods with different operating conditions and an accuracy of each of these experiments is estimated to be not better than ±30% [99]. On the other hand, Kakimoto et al and Nagano et al studied β particles from the 90 Sr source penetrating the air and used fixed filters which accepted some contributions from the small side bands.  In order to compare the reported yields given in different units for different experimental conditions, Arqueros et al [105] normalized them to photons per deposited energy at 293 K and 1013 hPa as listed in the last column in table 6. Only some values relevant for the following discussion are shown here. Though the new AIRFLY result [106] is preliminary, its systematic error is estimated to be less than 10%, which is smaller than those of Kakimoto et al (>10%) and Nagano et al (13%).
Though there are significant differences in photon yield among the experiments, their exact differences on the energy determination are not simple to estimate. The contributions of all the lines passing through the filters used and the wavelength dependence of the attenuation of the number of photons must be taken into account in each event. Nevertheless, the following remarks should be made when comparing the experiments.
1. According to Arqueros et al [105], the yield at 293 K and 1013 hPa for 337 nm band of Nagano et al is estimated to be 5.5 ph MeV −1 for deposited energy taking into account the escape of delta rays from the field of view of the experiment. The Auger group uses the reported value of Nagano et al whose fluorescence efficiency is determined under the assumption that the energy lost by the electron is fully deposited in the field of view. This value is about 10-15% smaller than the value listed in the table and hence about 5.0 ph MeV −1 [107]. Therefore the yield for 337 nm by Nagano et al used by Auger is about 25% larger than that of preliminary AIRFLY. If the preliminary AIRFLY yield is confirmed, the estimate of primary energy for the Auger FD data based on Nagano et al is likely to increase by about 25%.

Reconstruction method.
Shower size profile is used by the Fly's Eye and HiRes experiments, while the energy deposit profile is used by Auger to reconstruct the shower development. As discussed by Unger et al [108], it seems to be rather straightforward to use the energy deposit profile for the determination of the calorimetric energy. This is mainly because a different rate of fluorescence photons per electron or positron has to be assumed for the early and late stages of shower development, since the energy spectrum of particles in EAS changes in the course of its developmment and the ionizaiton energy deposit depends on the particle energy.
It is desirable to reconstruct the shower development of the same event by both methods to estimate the differences in order to understand the difference between the HiRes and Auger energy.

4.2.3.4.Čerenkov light subtraction.
Fluorescent light is contaminated by directČerenkov light, and the scatteredČerenkov light by air molecules and aerosols. TheČerenkov light is iteratively subtracted from the observed light to fit the Gaisser-Hillas function [109] in the case of the Fly's Eye [36]. In HiRes experiments, using the fit to the Gaisser-Hillas function, the number ofČerenkov photons scattered into the detector acceptance is calculated. The relative numbers of photoelectrons from fluorescence andČerenkov sources are determined after applying the filter transmission efficiency, quantum efficiency and the atmospheric attenuation. Then the charged particle shower curve to fit the Gaisser-Hillas function is recalculated after subtracting theČerenkov contribution. This process is iteratively continued until a satisfactory agreement is achieved with the observed light at the detector [96].
In case of the Auger experiment, the total energy deposit for all particles of all energies is determined as a function of the depth by solving a series of equations, in which both the fluorescence and theČerenkov light are treated as signal. Then the energy deposit profile is fitted to the Gaisser-Hillas function [109] and the total energy is determined [108].
Since the quantity of direct and scatteredČerenkov light depends much on the viewing angle to the shower axis from the detector, a comparison of the energy spectra for different azimuthal directions relative to the direction of the fluorescence detector may be helpful for understanding the validity of the subtraction.

4.2.3.5.
Field of view and the estimation of the aperture. Except for FE I, the field of view of the detector is limited to a narrow window. Since only a limited number of showers, whose shower maxima are within the field of view and their observed track length is larger than some criterion, are selected, the aperture depends on the fluctuations in X max , the distance to the shower axis and the arrival direction of the EAS. The estimation of the fraction of real showers, which are not selected, requires the real X max distribution.
The apertures of Fly's Eye and HiRes depend on the primary energy and composition and increase with energy. Detailed Monte Carlo simulations have been done for estimating them, including the experimental details. However, it should be remarked that the slope of the energy spectrum is quite sensitive to the slope of the energy dependence of the aperture, which varies with the composition.
On the other hand, a constant aperture over the whole energy range concerned is used for the analysis of the combined datasets of Auger. The aperture for the surface detector array events is much more reliable than that for the fluorescence events only.

4.2.3.6.
Atmosphere. The '1976 US standard' atmospheric model with the desert aerosol profile for altitudes less than 4 km is used for the Fly's Eye analysis [37]. For the HiRes stereo analysis, the profile obtained from a daily 00:00 UTC balloon flight at the Salt Lake City International Airport is used [110]. At the Auger site, the atmospheric profiles are measured with balloons and the average atmosphere in each month is determined [111]. These average atmospheres are used in the analysis, including the temperature dependence of the fluorescence yield in the form described by Kakimoto et al [100].
According to Keilhauer et al [112], the expected light from an EAS with a given energy deposit is reduced by 7-11% depending on the season, if the temperature-dependent collision cross sections and water vapor quenching are used. This reduction value is obtained by applying to the measured atmosphere at Pierre Auger Observatory. This means that the energy of the Auger energy spectrum must be increased by 7-11%.
The effect of these temperature-dependent collision cross sections must be taken into account for the Dugway atmosphere also, even though the effect due to water vapor quenching may be less at the HiRes and Fly's Eye site than at the Auger site.

4.2.3.7.
Relation to the depth of shower maximum. In figure 7, the energy dependence of the average depth of shower maximum, X max , determined by Auger hybrid and HiRes stereo are compared (from Kampert [113]). Both experiments are in agreement below 3 × 10 18 eV, where the aperture depends on the primary composition for the HiRes [113], whereas there are differences of about 20 g cm −2 above 3 × 10 18 eV, which is nearly the same as the accuracy of X max determination in each experiment. Since the relative values of the observed photons are enough for X max determination, items related to the absolute number of photons described earlier are not important in explaining the difference. Therefore, there must be some differences in selection criteria, which are related to the limited field of view, the estimation of aperture, the subtraction ofČerenkov photons and so on.
Not only the errors related to the determination of photon numbers, but also other errors relevant to X max determination may lead to an accumulated error in the determination of energy. Therefore, more detailed studies are desirable for each group.

Normalized spectrum of world data
Berezinsky et al [114] have proposed to normalize the observed energy spectrum at the theoretical pair-production dip, expected around 1 × 10 18 and 1 × 10 19 eV. These positions and the shape of the dip are robustly fixed by the interaction of protons with CMB and can be used for energy calibration of the detectors. By fitting each spectrum to the theoretical dip with minimum χ 2 method, they obtain the values of the factor λ, required to shift the observed energy spectrum, to be 0.9, 0.75 and 1.2 for AGASA, Yakutsk and HiRes, respectively [115]. This may be a good method as far as primary cosmic rays are mainly protons. After multiplying Average X max as a function of the primary energy observed by Auger Hybrid (red closed circles) and HiRes stereo (blue closed squares). From Kampert [113]. the primary energy by λ, the AGASA, Yakutsk and HiRes spectra are in good agreement with the Akeno spectrum in the 10 18 eV region and the Akeno spectrum may be used in the lower energy region as a good reference. Since other spectra cannot be fitted to the theoretical dip, the approximate shift by factors λ, listed as E * λ in the figure, are made, if necessary. The modified results are shown in figure 8 over the broad energy range, 10 14 -10 20 eV. This composite primary energy spectrum may be characterized with the knee energy as 4 × 10 15 eV, the 2nd knee as 6 × 10 17 eV and the ankle energy as 5 × 10 18 eV.
The energy dependence of the predicted X max based on four different hadronic interaction models are also shown for pure proton and iron primaries in figure 7 [113]. The red lines are for protons and the blue ones are for iron nuclei. The curves labeled A and B are predictions of X max based on the QGSJET01 model for the galactic-extragalactic transition model by Allard et al [116] and the dip-model by Berezinsky et al [117], respectively. If we use the X max from the HiRes stereo result, the conclusion about the primaries being mainly protons at the highest energy along with the results shown in figure 8 may be acceptable. The correlation of the highest energy cosmic rays with nearby extragalactic objects, suggested by the Auger observations, is also compatible with proton primaries [118]. The autocorrelation function calculated from the 27 highest Auger events seems to be quite significant, with angular scale of 5 • over the random distribution [119], which again suggests protons to be the primary particles at the highest energy.
On the other hand, if we use the Auger hybrid X max result, the primary composition seems to be mixed in the energy region above E ankle and the normalization to the pair production dip may not be appropriate. In this case, the flux of the HiRes energy spectrum is likely to be changed by a factor whose value corresponds to a shift of the average X max by −20 g cm −2 . Nevertheless, the primary energy spectrum normalized to the HiRes spectrum is shown in figure 9, since it can also be connected smoothly with the spectra in 10 18 eV region except with the Akeno spectrum. The factors used for multiplying the energy bins for each experiment are listed in the figure. In this figure, the E 2nd may be near 2 × 10 17 eV, which is significantly lower than the case shown in figure 8. It is quite important to determine not only the composition in the 10 19 eV region, but also the E 2nd and the spectral slopes below and above E 2nd , in order to determine the most reliable energy spectrum over a broad energy range.

Super GZK events
After the normalization of the energy described above, there remain some AGASA events which seem to exceed the GZK energy. There are a few events which clearly exceed the GZK energy, Fly's Eye event [39] and AGASA event [120]. One strong candidate above 10 20 eV has been reported by the HiRes, albeit on a night of variable weather conditions [121]. Excepting the reported 1.6 × 10 20 eV event, an event of 2.0 × 10 20 eV, assuming typical atmospheric clarity at the Auger site, is also reported [122] though this event is outside the acceptance criteria used for the Auger spectrum. The depths of the shower size maximum for the HiRes and Auger events are around 850 g cm −2 . These are deeper than the average of the shower maximum obtained from an extrapolation from the lower energy region.
Are they the result of an overestimation of the energy due to any experimental problem or due to statistical effect? If we look at the AGASA spectrum carefully, there seems to be a suppression of the flux between 10 19.6 and 10 20 eV and the super 10 20 eV events seems to be a different component.
The longitudinal development of showers initiated by gamma rays above 10 19 eV is suppressed by the Landau-Pomeranchuk-Migdal (LPM) effect [123,124] in the atmosphere. For primary gamma rays above ∼ 10 19.5 eV, pair creation may occur in the geomagnetic field and the shower development in the atmosphere is similar to that without the LPM effect, resulting from a superposition of lower energy gamma-ray subshowers [125]. As a result, the shower energies are underestimated by a muon sensitive detector such as a water tank, but not by an electromagnetic component sensitive detector such as a scintillator. The maximum position for showers due to gamma-ray primaries is deeper than that for the hadronic showers and their detection may be inefficient in a shallow field of view such as in the HiRes-I. The gamma-ray possibility for the highest energy AGASA events is examined by Shinozaki et al [126] and there is no indication that most of the observed events are gamma-ray showers. Similar results are obtained by Risse et al [127] and Rubtsov et al [128]. The possibility of a gamma-ray origin for the Fly's Eye highest energy event is examined by Halzen et al [129] and found to be unacceptable.
Though no indication of gamma-ray primaries has been found above the GZK energy as yet, the endeavor to find a definitive signature should be continued with gamma-ray sensitive detection methods.

Conclusion and outlook
New results on the primary cosmic-ray energy spectrum in the highest energy region have been reported. The spectrum has been found not to be extended with a constant slope beyond about 5 × 10 19 eV. Since the composition of the primary flux in this energy region has not yet been confirmed, an interpretation of the ankle and the suppression of the spectrum around 5 × 10 19 eV as due to the pair creation dip and the GZK cutoff energy, respectively, remains to be confirmed with more studies.
Summarizing the spectrum observed from 10 14 to 10 20 eV with world wide experiments, the E 2nd is around 6 × 10 17 eV, if the energy calibration with the pair production dip and the Akeno spectrum are accepted. The E 2nd is around 2 × 10 17 eV, if the HiRes spectrum is accepted. The energy scale and/or the flux of the Auger spectrum seems to be too low when compared with the spectrum in the lower energy region.
The following are expected to minimize or eliminate some of the existing problems in UHECR experiment: 1. Though the fluorescence technique measures, in principle, the total energy deposit in the atmosphere which is a good measure of energy, there are possibilities of accumulating errors at several stages. Some possibilities to increase the primary energy determined by the present fluorescence detectors are mentioned below.
• The energy of the Auger FD events must be increased by 7-11% depending on the season, if the temperature-dependent collision cross sections and the actual humidity profiles at the site are used [112]. It is necessary to examine how much correction is necessary for the Fly's Eye and HiRes experiments for the Dugway atmosphere. • The preliminary result of AIRFLY experiment for the 337 nm band [106] is smaller by about 25% than that of Nagano et al [101] presently used by Auger. If the AIRFLY result is confirmed, the Auger energy is likely to be increased by a similar amount. Combining both of these amounts, there is a possibility of an increase of about 35% (10 + 25%) for the Auger FD energy and hence the Auger spectrum is similar to that shown in figure 5. In the case of Fly's Eye and HiRes, the increases may be approximately 10 and 35%, respectively, if absolute values are normalized at 337 nm.
Further study on the determination of the absolute yield by the AIRFLY group is going on at GeV and MeV energies [106]. A confirmation of the absolute value for the 337 nm band is highly expected. Comparing to the novel but complicated procedure in absolute yield determination of AIRFLY, the photon counting method within a gated window by electrons is quite simple. Its systematic error is only limited by the quantum and collection efficiencies of the PMT used. Their improvement is also highly expected. 2. A new energy calibration, using an electron linear accelerator, which is under preparation at the TA site [130], may help in reducing systematic errors of most items described in section 4.2.3. Observing an air shower generated by a beam of electrons of known energy, the calibration for the fluorescence yield (in units of a single photon signal defined in the experiment), the reflectivity of the mirror, the transmission through the filter and the quantum and collection efficiencies of the photomultiplier will be made altogether under controlled conditions. That is, most important calibration problems such as absolute fluorescence yield and a definition of single particle are possible to solve at the same time within the experiment. 3. If we use the results from Monte Carlo simulations with available hadronic interaction models to estimate the primary energy, the conversion from surface array data to primary energy seem to be in agreement within possible experimental errors. However, there remain discrepancies in energy scale from the fluorescence measurement, and a deficit in the number of observed muons in EAS by Auger even with iron primary assumption [85].
It is desirable to determine the energy scale of S 38 or S(600) and its energy dependence from simulations using a hadronic interaction model which reproduce the measured muon numbers. 4. The experiment which covers the energy range 10 17 -10 18 eV is so far limited to the Akeno [43]. It is quite important for the confirmation of the energy spectrum in the highest energy region to examine whether E 2nd is near 10 17 or 10 18 eV and also to determine the spectrum slope below and above E 2nd . The results from experiments whose energy range may cover 10 17 -10 18 eV, such as KASCADE-Grande [131], HEAT (high elevation Auger telescopes) [132], AMIGA (Auger muons and the infill for the ground array) at the Auger site [133] and TALE (telescope array low energy extension) at the TA site [134] shall be eagerly expected. 5. Some candidate events beyond 10 20 eV are not still excluded as due to overestimation of primary energy or simply statistical effect. The experimental results with the scintillation detector array by the Telescope Array project [92], which started its operation in the northern sky, may give a clue whether the AGASA super-GZK events are due to density overestimation caused by delayed neutrons or shower front structure at highest energy region observed. 6. The accumulated events by the Auger South within a few years by hybrid methods and further studies on the absolute calibration of FD detectors are expected to resolve the differences in the results available presently on the detailed shape of the primary energy spectrum at the highest energy and the existence of super-GZK events.
The spectrum with the E GZK cutoff is anticipated to recover beyond a few 10 20 eV if the primary cosmic-ray spectrum extends to energies far beyond 10 21 eV [135]. To open a new window for charged particle astronomy covering the whole sky, the Auger North [136] and the JEM-EUSO projects [137] are under preparation. The JEM-EUSO is expected to explore the anticipated recovery of the spectrum. Therefore our search for the end of the cosmic-ray energy spectrum will be continued even after a century since the discovery of cosmic rays in 1912.