Measuring extensive air showers with Cherenkov light detectors of the Yakutsk array: the energy spectrum of cosmic rays

The energy spectrum of cosmic rays in the range E∼1015 eV to 6×1019 eV is studied in this paper using air Cherenkov light detectors of the Yakutsk array. The total flux of photons produced by the relativistic electrons (including positrons as well, hereafter) of extensive air showers in the atmosphere is used as an energy estimator of the primary particle initiating a shower. The resultant differential flux of cosmic rays exhibits, in agreement with previous measurements, a knee and ankle feature at energies of 3×1015 and ∼1019 eV, respectively. A comparison of observational data with simulations is made in the knee and ankle regions in order to choose the models of galactic and extragalactic components of cosmic rays that describe well the energy spectrum measured.


Introduction
The Cherenkov light emitted in the atmosphere by relativistic electrons of extensive air showers (EASs) of cosmic rays (CRs) carries important information about the shower development and the primary CR particles.Well-known application of the air Cherenkov light technique is the γ-ray astronomy.But in this article we aim at another application of the technique -namely to measure Cherenkov photons from EAS initiated by CR of energy above 10 15 eV.These giant showers supply plenty of light so that one can detect it with unarmed photomultiplier tubes (PMT) triggered within pulse duration.
Since the first observation by Galbraith and Jelley [1] and a systematic measurement of air Cherenkov light properties in the Pamir experiment [2], a number of EAS arrays have been equipped with Cherenkov light detectors.It seems that the most durable and plentiful in Cherenkov light data is the Yakutsk array experiment [3,4].The total flux of light is used to estimate the primary energy in a model independent manner and the radial distribution of the light intensity at ground level is used to infer the position of shower maximum, X max , in the atmosphere [3].
In this article we focus on the experimental data obtained in Yakutsk with Cherenkov light detectors aiming at the CR energy spectrum.The paper is structured as follows.In Section 2 we outline the general characteristics of the Yakutsk array experiment, while in Sections 3 and 4 particular properties of the Cherenkov light detectors are given: detector design and calibration of the signal.In Section 5 the measurement and monitoring of the atmospheric extinction of light is described.The measured lateral distribution of the Cherenkov light intensity at the observation level is given in Section 6.The method used to estimate the energy of the particle initiating the EAS is described in Section 7. Resultant energy spectrum of CRs is discussed in Section 8. Our conclusions are set out in Section 9.In two Appendices additional material is given essential for the subjects considered.

The Yakutsk array
The Yakutsk array is located in Oktyomtsy near Yakutsk, Russia (61.7 0 N, 129.4 0 E), 100 m above sea level (1020 g/cm 2 ).At present it consists of 58 ground-based and 6 underground scintillation detector stations to measure charged particles (electrons and muons) and 48 detectors -PMTs in shuttered housing to observe the atmospheric Cherenkov light.During more than 30 years of lifetime the Yakutsk array has been re-configured several times, the total area covered by detectors was maximal about 1990 (S ef f ∼ 17 km 2 ), now it is S ef f ∼ 10 km 2 .In the central part of the array there is a denser domain with 100-250 m detector spacing.During the whole observation period approximately 10 6 showers of the primary energy above 10 15 eV are detected; the three highest energy events selected with axes within the array area and zenith angle θ ≤ 60 0 have an energy E > 10 20 eV.
The actual detector arrangement of the array is shown in Figure 1.Charged particle detectors of 2 m 2 area are built in stations in couples; the Cherenkov light detectors -PMTs of 176 cm 2 and 3 × 176 cm 2 acceptance area, forms the medium, C 1 (∼ 500 m spacing), and the autonomous, C 2 (50 to 200 m spacing), subsets.The latter was added in 1995 with the aim to study air showers in the energy range 10 15 − 10 17 eV via the Cherenkov light measurements [5].All detectors/controllers and data processing units of the array are connected into the data handling network shown in figure 2.

Air Cherenkov light detector design
Charged particle detectors of the array were described in [3].Here we will detail a Cherenkov light detector unit.It consists of a vertically mounted PMT (FEU-49B, 15 cm diameter) with amplifier in a metal container blackened inside [6].An upper hole provides θ ≤ 55 0 aperture (figure 3).To protect the photocathode from sunlight the motorized light-proof lid is set.At night, all lids of the array can be commanded remotely to open.PMTs and amplifiers are powered around-the-clock to guarantee a stability of performance.When the lid is open, a fan blows with warm air to keep snow and dust out of the photocathode surface.
There is a variant of detector with three PMTs in a housing which can operate independently or in summation of signal in order to increase acceptance area at the shower periphery.In addition, dedicated detectors were used to measure the shape and width of the Cherenkov signal from the shower.As an example, the pulse shape of the Cherenkov signal at the shower periphery is shown in figure 4, while the halfwidth of the signal as a function of radial distance is given in figure 5.
The spectral sensitivity of the PMT used is shown by the dashed line in figure 6 together with air Cherenkov light spectrum (solid line) and the atmospheric

transmission curve.
Two kind of triggers were used to select the showers from the background at the Yakutsk array: produced by the scintillators and by Cherenkov light detectors.In the first case a coincidence signal (if the particle density, ρ, is greater than 0.5 m −2 in two scintillators of each station within 2 µs) passes on to the central controller.Trigger-500 is then produced in the case of a coincident signal (in 40 µs) from three or more stations with ∼ 500 m spacing (≥ sixfold coincidence).Similarly, trigger-1000 is produced by ∼ 1 km spacing stations.After 1992 when 18 new stations were added, the array area is increased from 2.5 km 2 to 7.2 km 2 where trigger-500 operates.That is why we can deal with EAS in the energy range from 3 × 10 16 eV to 3 × 10 19 eV using the same trigger.
The C 1 subset of Cherenkov light detectors has no hardware trigger of its own.Instead, the scintillator trigger controls the detection of Cherenkov signal in this subarray.The amplitude of PMT signal is recorded within 5 µs after the coincident signal in two scintillators of the host station.In order to exclude detection of an accidental Cherenkov signal when there is no trigger-500/1000, a double trigger condition is in use in data analysis: • three or more scintillator stations have detected the charged particles (ρ > 0.5 m −2 ); • the Cherenkov light intensities in three or more PMTs are greater than 2.4 × 10 5 /4.5 × 10 5 m −2 depending on the acceptance of PMT.
The thresholds have been chosen to keep the signal/noise ratio greater than 3.
The second case is C 2 subset which has an independent Cherenkov light trigger formed by three or more PMTs having detected light intensities above a given threshold within 10 µs.The signal integration time of the individual PMT is 0.5 µs.

Detector calibration
In order to find the absolute response of the array to Cherenkov light we have considered the sensitivity of all the PMTs in the detectors and have measured the conversion factors from light intensity into ADC value.The Cherenkov light emitted by relativistic muons in a transparent medium with known refraction index (optical radiator) was used in the detector calibration, as was suggested in [8].In this paper we are keeping along the description given in [6,9].Distilled water and plexiglas are used to measure the output signal in a geometry shown in figure 7. A telescope composed of two PMTs above and below the radiator, spaced 1 m from the detector, bottom one under the lead shield (10 cm), selects high-energy muons from background CRs.On the photocathode of the detector PMT is placed (above a thin layer of glycerin) a plexiglas disk of 5.5 cm thickness and 15 cm diameter with polished side surface.As another version of radiator, distilled water was used without any pan layer.The Cherenkov radiation angle in water is 35 0 , while it is greater in plexiglas (48 0 ) and in glycerin (47 0 ).
Integrating it above the threshold energy (E µ > 0.25 GeV) with the muon spectrum measured at sea level we have N γ = 2850.Then the ratio of the air Cherenkov light intensity Q(R) at distance R from the shower core to N γ is where A s , A r are the detector signal amplitudes from the shower and radiator; S is photocathode area; R s ≤ 0.04 is the reflectivity factor of air-glass junction.
where K r (λ) is a radiator transparency of effective thickness 4 cm (taking into account a radiation cone).The spectral characteristics of PMT, plexiglas and glycerin are measured with the spectrophotometer.Similar values were derived in the case of water radiator [9].Optical radiators used in the laboratory for absolute calibration of the detector produce photon number insufficient for the routine ratio calibration of a multitude of array detectors in the field.Instead, the scintillator disk (D=15 cm, h=5 cm) is in use as the light source in this case.Intercalibration of different radiators has been performed using the experimental setup (figure 7).

Atmospheric extinction of light
Cherenkov light undergo extinction in the atmosphere because of absorption on molecules by Rayleigh scattering and Mie scattering by aerosols and absorption.
During the observation periods in winter (average temperature −40 0 C) Cherenkov light absorption in the atmosphere is negligible and only molecular and aerosol scattering of photons are taken into account.Molecular scattering is almost constant, while aerosol concentration in the boundary layer above the surrounding terrain is of diurnal and seasonal variability.We are using the event rate of the showers in the energy range 10 15 − 10 16 eV detected with a PMT subset in hour and 15 min intervals in order to monitor the atmospheric transparency for the light generated in EAS [3].The method is based on the Cherenkov light flux proportional to the primary energy: where N (> Q thr ) is integral number of events with the light intensity above the threshold detected in i-th and basic periods; τ is atmospheric extinction coefficient; κ is the energy spectrum index.
Resultant extinction coefficient is a product of variable τ a due to current aerosol concentration and basic coefficient τ 0 caused by the molecular scattering and minimal aerosol extinction.Rayleigh scattering parameters evaluated for the Yakutsk array conditions give τ R = 0.9 [6] and the corresponding attenuation length λ R = 300 km, averaged over wavelength range under consideration.
The integral spectrum index of the Cherenkov light flux around the knee is [6]: Evaluating the coefficient τ 0 as a function of X max in the shining point approximation [11] by assuming the attenuation length due to aerosol scattering constant above 1.5 km height and decreasing below [12], one can derive the average extinction coefficient shown in figure 8.The approximation is where X max is in g/cm 2 .Finally, we have found τ = 0.765 ± 0.015 at θ = 20 0 , E = 7 × 10 16 eV averaged over all observation periods, and inserting spectrum index in equation (2).

Lateral distribution of air Cherenkov light
During ∼ 15000 hours of observation ∼ 60000 showers of energy above 6 × 10 16 eV were detected by the medium C 1 subset.The autonomous C 2 array data consist of ∼ 200000 showers with E > 1.2×10 15 eV detected during ∼ 3200 hours of observation.The Cherenkov trigger condition only is used to select these showers.
Relative dispersion of the light intensity observed with different detectors at R ∼ 500 m from the shower core has been evaluated selecting showers with 3 detectors hit in the interval 400 < R < 500 m not far from each other (≤ 75) m.Intensity dispersion contains, in addition to inherent fluctuations instrumental errors, core location error, intercalibration error, etc. Supposing chance variation of the signal around an average Q(R) ∝ R −2.35 [13], we have found δQ/Q = 0.25 ± 0.13.In order to minimize the Q(R) uncertainty due to EAS axis location error, one has to select showers with axes within R < R opt m, (80 eV) where stations are spaced closely.
In the highest energy domain measurements of Cherenkov light exist, carried out by the Haverah Park group [14].In figure 9 our results are given in comparison with [14].The data are normalized at 200 m core distance and exhibit consistency in the lateral distribution function (LDF) shape.The solid lines are an approximation to our data given in [9].
Model simulations give a variety of lateral distribution functions, some of them are close to our data, as illustrated in figure 10.The lateral distribution depends on the attenuation lengths, multiplicity of secondaries in the interactions, primary mass composition etc., which can be parameterized with X max .Another considerable influence is the angular distribution of electrons in the shower.The combination of Our cumulative results on LDF measurements (zenith angle θ < 30 0 ) are given in figure 11.The data of both subsets C 1 , C 2 are parameterized by the intensity at 150 m from the shower core, Q 150 , the only core distance really present in the shifting range of measurements when the primary energy is rising from E 0 ∼ 10 15 to 10 19 eV.The data are consistent with the previous results of the Yakutsk array concerning the Cherenkov light and can be described by the suitable EAS model simulation [15,16].
No abrupt change of LDF parameters is seen, so we choose a rather smooth approximation curve to fit the experimental data in the whole energy range [4]: where

Energy estimation
The total flux of Cherenkov light emitted, Q tot , is our main estimator of the primary particle energy.In order to derive the relation between Q tot and ionization loss of the shower electrons in the atmosphere we have used formula (1) for the number of photons induced, N γ , in ratio to that from optical radiator, N r : where v is electron velocity; n, n r are the refraction indexes in air and radiator, respectively; c = 1.Emission is possible when vn > 1 or The number of photons induced in the depth interval dx is where ; ρ r is radiator density; ρ 0 , n 0 are air density and refraction at the observation level.
The total light flux is given using the electron differential energy spectrum in the shower at the depth x: where τ (x 0 − x) is the extinction coefficient of light along the path x 0 − x.
The number of electrons in a shower is approximately N (x, E, E 0 ) N (x, E 0 )χ(E), where N (x, E 0 ) is the number of electrons with E > 0; χ(E) is the universal electron spectrum at the shower maximum.Then leads to the relation with the ionization loss of electrons in the atmosphere if we assume the extinction τ = 1 [19].
In the real case of the Yakutsk array conditions T = −30 0 C, P = 754 Torr and the extinction of light described in Section 5, it was found that the relation is [19]: for X max ∈ (500, 1000) g/cm 2 ; the accuracy of the relation is ∼ 5%.It demonstrates the advantage of the air Cherenkov light measurement technique: the relation is determined by X max and extinction of light only; interaction model dependence is parameterized by means of X max .
In Appendix A we analyze the shower parameters governing the energy fractions transferred to EAS components.

Experimental evaluation of the energy transferred to EAS components
The energy fractions of the main EAS components can be estimated using the Yakutsk data.The ionization loss of electrons is measured when detecting the total flux of the Cherenkov light at ground level.The detector disposition of the Yakutsk array is appropriate to measure Q tot in the range above about Q 150 = 10 7 m −2 as is shown in table 1.In each Q 150 bin the LDF extrapolation formula (4) is used to calculate the total flux.
Conversion of the measured Q tot to E i is carried out along equation ( 9) with parameters relevant to the particular observation period, taking into account detector calibration and atmospheric extinction of light.
Other portions of the energy carried out by electromagnetic and muonic components beyond sea level is evaluated via the total number of electrons: where 0 , t 0 are the critical energy and radiation length of electrons in air; attenuation length λ e is derived from zenith angle dependence of N e [20]; and muons measured at the ground level: where the average energy of muons E µ is taken from the MSU array data [21].
Residuary energy fractions transferred to neutrinos E ν , nucleons E h etc., unmeasurable with this array, are estimated using computational modeling [19].The resulting apportioning of the primary energy E 0 is given in table 2.
The energy fraction carried by electromagnetic component (E em = E i + E g ) appears to be the basic contribution to the total energy of the shower, and its energy dependence (measured with Cherenkov light detectors + scintillators of the Yakutsk array [22]) is illustrated in figure 12 in comparison with CORSIKA/QGSJET estimation [23].
Due to the air Cherenkov total light flux and the electron and muon number which are experimental values measured at ground level, only about 5% of the primary energy in the interval E 0 ∈ (10 18 , 10 19 ) eV is calculated by using model assumptions.So we consider the energy estimation used to be model-independent within these bounds.
Moonless nights, when air Cherenkov light measurements are possible, constitute ∼ 10% of the observation period.In order to evaluate the primary energy of the bulk of showers, the correlation S 600 = 1.56 × 10 −8 Q 1.01 150 is used between the charged particle density at 600 m from the shower core, S 600 , and the light intensity at 150 m from the core (figure 13) which, in turn, is related Table 3. Parameters of the relation between Q 150 and the total flux of air Cherenkov light and the energy of the primary particle initiating EAS.
2.39 × 10 6 0.90 (6.87 ± 1.44) × 10 10 0.87 ± 0.02 10 6 − 10 to the total flux of the Cherenkov light in the atmosphere Q tot = a 1 Q b1 150 , parameters are given in first columns of the table 3. Using derived relation between the primary energy and Q tot (equation 9) the final formula was found to link Q 150 with E 0 in the interval θ < 15 0 [19]: where the numerical values of a 2 , b 2 are summarized in the last two columns of table 3 for different Q 150 intervals.The observed densities S 300 /S 600 at various zenith angles are connected to the 'vertical' one (θ = 0 0 ) along attenuation curve [24].In order to measure the attenuation length of these densities for fixed energy, we have used two different methods -well-known equi-intensity cut method, and fixing the Cherenkov light intensity at 400 m from the core as the equivalent of the primary energy, taking into account the light absorption in the atmosphere.In figure 14 the results are given.Experimental points are consistent with each other for the two methods used and can be described by the sum of two components -a soft component (electrons, attenuation length λ e = 200 g/cm 2 ) and a hard component (muons, λ µ = 1000 g/cm 2 ) [20]: where β 300 is the hard-component fraction.Attenuation curve for S 600 is the same but β is different: Experimental uncertainties when estimating the EAS component energies are summarized in the last column of the table 2. The main contribution arise from δE i which is governed by uncertainties in the atmospheric transparency (15%), detector calibration (21%) and the total light flux measurement (15%).Errors in estimation of N e , λ e , N µ determine the next two items (for ionization loss in the ground and δE µ ).Resultant uncertainty in energy estimation is the sum of all errors weighed with the second/third columns of the table 2: δE 0 ∼ 32% [19].Extra 20% are added due to a Q 150 -S 600 conversion uncertainty.
To illustrate the energy estimation method used, four showers detected at the Yakutsk array in the range E 0 ≥ 10 20 eV, θ < 60 0 , and axes within the array area, are given in table 4 (one event is added slightly below the threshold because the energy estimation error is larger than the tiny difference).

The energy spectrum of cosmic rays derived from air Cherenkov light measurements
The Cherenkov light detector subsets of the Yakutsk array give us the opportunity to reconstruct the energy spectrum of cosmic rays in the energy range from E 0 ∼ 10 15 to 6 × 10 19 eV.The total number of electrons and muons measured at the ground level are used to estimate the additional energy fractions carried by EAS components but the final relation (10) comprises the light intensity alone.The intensity of CRs has been evaluated using the number of EAS events derived from the showers detected with PMTs, and aperture S ef f T Ω of the sub-arrays in a particular energy interval, where the array area bounded by the perimeter, S ef f , depends on the primary energy and zenith angle; Ω is the acceptance solid angle; and T is the sum of observation periods.Under the condition of the array configuration changing due to rearrangements and non-active detectors, S ef f has to be calculated for a given period using the Monte Carlo technique.
The acceptance area has been simulated as a function of Q 150 averaged in the zenith angle interval (0 0 , 30 0 ) for the two Cherenkov light detector subsets shown in figure 1.A lateral distribution fit accord to Eq. ( 4) is used together with instrumental and statistical errors (Gaussian with 25% relative deviation) to model the trigger of each detector subset with 100000 fake showers.In the case of the autonomous sub-array, the Cherenkov trigger is simulated, while in the C 1 case the double trigger for scintillator and PMT signals has been modelled.The average (S ef f = S 0 n triggered /10 4 , where S 0 is the array area inside the perimeter; n triggered is triggered number of events) is shown in Figure 15 versus the parameter Q 150 .Corrections for inoperative detectors are not shown here.
The shower data gathered after the latest array re-configuration were used to work out the spectrum.Namely, 1993-2007 for the medium C 1 subset, and 1995-2007 for autonomous subset.In order to evaluate the intensity of the primary flux we have collected data during observation periods with a light extinction better than 0.65 and shower axes within area of the corresponding sub-array.
The resulting differential all-particle spectrum of cosmic rays is shown in figure 16 in comparison with the data from other Cherenkov detector arrays, namely, BLANKA [25] and Tunka [26].The present data (given in tabular form in Appendix C: tables C1, C2) exhibit the spectrum irregularity near E ∼ 10 19 eV, the 'ankle', seen by all arrays in the area [27]; at lower energies the results of all three arrays are compatible with a 'knee' at E ∼ 3 × 10 15 eV revealed in the pioneering works of the MSU array [28].Below E = 10 18 eV there is a transition region between the two subsets, C 1 , C 2 , so the 'second knee' visible here may be due to the data sewed together.
Hereinafter, two energy regions below and above E = 10 18 eV are analyzed separately in order to consider the two irregularities in the spectrum.The region below 10 18 eV is shown in detail in figure 17 where CR intensity is multiplied by E 2.75 in order to emphasize the spectrum irregularity.In addition to spectra around the knee measured by Cherenkov detector arrays (BLANKA, Tunka and Yakutsk taken from the previous Figure ), the data of scintillation detector arrays (Akeno [29], KASCADE [30] and Tibet [31]) are given for comparison.Due to air Cherenkov light signal proportional to the number of electrons in the shower, no difference is expected in the spectrum shape measured with scintillator or PMT arrays.Measured CR intensities are corrected in the case of Cherenkov detector arrays along the algorithm in Appendix B. Integral energy spectrum indices below and above the knee are assumed to be 1.67 and 2.1, respectively, as measured by the Tibet-III array, while energy evaluation errors are estimated as 0.12 (Blanka), 0.2 (Tunka) and 0.25 (Yakutsk).Original intensities from Akeno, KASCADE and Tibet arrays are not changed because the energy/intensity reconstruction procedures from N e , N µ measurements include the conversion factor (B.1) needed.
For the reconstructed energy and intensity there is an interaction model/primary composition dependence in all the data from arrays.For instance, the Tibet-III results indicate a ∼ 20% systematic error due to chemical composition and ∼ 10% discrepancy between QGSJET01c and SIBYLL2.1 interaction models below E = 10 16 eV.The uncertainty may by much worse if to use interaction models not carefully tuned to the measured EAS observables.
We have used energy correction factors to compare the measured differential As examples of model calculations, four predicted spectra are shown together with the Cherenkov light data of the Yakutsk array in figure 18.Only the shape of the all-particle spectra can be compared with the data because of the free parameters in models -the intensity of CRs and, to a lesser degree, a knee position due to the magnetic field uncertainty in the sources and in the interstellar medium.
All the model spectra are compatible with experiment, especially, in view of the dispersion in the results of arrays (e.g.figure 17) which is greater then the shape variation due to models.At energy E > 0.1 EeV the anomalous diffusion model spectrum is harder then that observed in Yakutsk, but the points are within experimental errors.The lack of CRs in SNR acceleration and diffusion models above energy 0.1 EeV can be filled up by the extragalactic component.
There is seen to be a hint of the fine structure in the energy spectrum measured with the Yakutsk array.Although it is possible that the undulations around the knee are caused by the instrumental errors ‡, but on the other hand, there are 40 size spectra and 5 Cherenkov light spectra measured before 2005, which demonstrate the second excess ('peak') at lg E = lg E knee + 0.6 besides the knee itself [36].The second peak in our data is approximately at this energy.However, more Cherenkov light data are needed to scrutinize the subject.Our autonomous sub-array is able to supply with a sufficient sample of EAS events within the next few years.
Recently, Lagutin et al. [37] examined the contribution of a nearby SNR-type source to the energy spectrum of CRs produced in anomalous diffusion model.They found several sequential peaks caused by H, He and CNO nuclei around the knee in all-particle spectrum, confirming results of the single source model [34] in the case of the 'background' anomalous diffusion model.‡ vertical bars show statistical errors only ) Figure 18.Our energy spectrum and the results of galactic cosmic ray simulations.Curves: SNR acceleration model [32] (BV); diffusion model [33] (KP); single source model [34] (EW) and anomalous diffusion model [35] (LNU).Our conclusion concerning the part of the energy spectrum below 10 18 eV is that while all four models considered are compatible with our measurements in the knee region, e.g. the intensities below and above the knee, the only model able to describe the fine structure in the spectrum is the single model.So this model combining a recent nearby SNR with the background (anomalous) diffusion of CRs in Galaxy is the best fit for the Yakutsk array data.
When comparing the upper half of the spectrum measured for E > 10 18 eV, with the energy spectra observed by giant EAS arrays, including that of the Yakutsk array, the somewhat different energy estimation has to be taken into account [3,20,38].Due to the comparatively small acceptance area of the Cherenkov detector subsets (figure 15), our data are reliable up to ∼ 10 19 eV; this includes the ankle region.Another purpose is to compare UHECR intensities measured with different techniques.
In the table 5 are given the energy estimation errors of the AGASA [39], HiRes [40], PAO [41] and Yakutsk [19] experiments and the intensity conversion factors calculated as described in Appendix B for energy bins where the index is constant and the variation of instrumental errors is negligible.The integral energy spectrum index of CRs is assumed to be 2.3 and 1.9 below and above the ankle in the spectrum, and κ = 4 at lg E > 19.8, along the results of HiRes fitted by the broken power law [40].
Observed UHECR spectra are given in the left panel of the figure 19 with the intensity correction factors applied: observed intensities are decreased by R J and  spectra are displaced over 0.434σ along lg E. § This is preliminary, crude procedure to reconstruct the spectrum breaks, but for the purpose of the intensity comparison it may be sufficient.Two spectra of the Yakutsk array are compatible within errors above 10 19 eV but diverge at lower energies.The possible reason of a discrepancy can be systematic errors in primary energy and CR intensity estimation near the scintillator threshold of the trigger-500.The work is in progress to surmount the divergence.
Comparison of UHECR spectra measured using different detectors and energy reconstruction methods infer the existence of systematic differences between resultant energies obtained at the arrays.The average values, Êi , should be corrected, too.While the true primary energy is unknown, cross calibration of energy estimation methods can be carried out adjusting correction factors, R E , for the pairs of observed spectra, converging it together.Then the resulting spread of factors elucidates the confidence interval for the CR energy estimated.
Table 6 demonstrates the variety of correction factors for the energy estimation methods of the different groups.
To illustrate the result of corrections R j × R E applied to estimated energies and § This may be unwarranted in the PAO case due to [42].for the given primary energy E intensities, the measured spectra are re-plotted in the right panel of figure 19.Energy scale factors are used here (arbitrarily, but preferably close to the median) from the last column of table 6, although any other column may be used as well.
The shape of the UHECR spectrum measured by all the arrays is compatible within errors, if the energy estimations are calibrated.Namely, the observed position of the ankle and energy threshold of GZK suppression ¶ are in satisfactory agreement.
A difference in the energy of EAS primary particle estimated basing on the data of the two subsets of the Yakutsk array detectors -scintillators, Y Sc , and air Cherenkov light detectors, Y Ch , originates, presumably, in the systematic error of the relation used between the charged particle density at 600 m from the shower core and the light intensity at 150 m from the core.
The UHECR propagation modeling results are illustrated by the four examples in figure 20 in comparison with the Yakutsk array data.
Bahcall and Waxman gave two-component (Galactic+extra-galactic) model with the shape insensitive to the choice of absolute energy scale [44].
Another approach was used by Berezinsky et al. assuming UHECRs as extragalactic protons from uniformly distributed sources [45].The electron-positron pair production in collisions of protons with relic photons results in the energy spectrum of extragalactic CRs with the 'dip' feature in this model.
De Marko and Stanev' model fits the UHECR spectra measured by AGASA and HiRes with different injection spectra at CR sources that are uniformly and homogeneously distributed in the Universe [46].The best fit they found assuming that cosmic rays (E > 10 19 eV) are protons, varying the index, emissivity and cosmological evolution parameters of the injection spectrum, is given in the figure .Wibig and Wolfendale focus on the ankle in the primary energy spectrum attributing it to the rapid transition from Galactic to Extragalactic component of CRs [47].The sum of a smoothly falling Galactic spectrum and a power-low EG spectrum + fitted to the Yakutsk array scintillator data is shown.¶ except AGASA data + index -2.37 All the models used demonstrate GZK suppression * and the ankle features in agreement with the data from arrays.It is not surprising in view of the fact that the source emissivity and the ankle position in the energy scale are free parameters, and GZK effect is embedded in models.
The lack of EAS events above 10 19 eV observed with Cherenkov light detectors hinders in choosing a model of the better fit.The index of the energy spectrum observed in the interval (10 19 , 3 × 10 19 ) eV is d lg J/d lg E = −2.1 for Cherenkov light detectors data and is −2.7 for charged particle detectors data.If we use these values as the confidence bounds then model indices are within the interval.The model of Bahcall and Waxman is closest to our Cherenkov light detectors data and can be considered as preferable among equally matched.

Conclusions
The total flux of air Cherenkov light with subsidiary data on the electron and muon sizes at the ground level is used to estimate the energy of the primary CR particle initiating EAS.The relation between the total flux and ionization loss of electrons in the atmosphere is derived which depends on the extinction of light and X max parameters; the latter is the only parameter to accumulate the interaction model dependence of EAS development in this case.
The independent measurement technique based on the Cherenkov light detectors of the Yakutsk array enabled us to observe the cosmic ray energy spectrum in the range from E ∼ 10 15 eV to 6×10 19 eV.Two spectra measured with different detectors of the Yakutsk array -scintillators and Cherenkov light detectors, exhibit an ankle feature below E = 10 19 eV.The autonomous sub-array data confirm the previous observations of the knee at E ∼ 3 × 10 15 eV.A comparison of our results with the data of other EAS arrays shows the compatibility of spectra if the energy estimations are corrected.
The energy spectra predicted for several models of galactic CRs demonstrate agreement within experimental errors with our data in the knee region.However, only the single source model describes the fine structure of the spectrum observed attributing it to the contribution of nuclei from a recent nearby supernova.We have chosen this model as the best fit for our data below 10 18 eV.
At the highest energies extragalactic CRs forming a dip, or transition between galactic and extragalactic components above 10 18 eV are thought to be responsible for an ankle detected in the energy spectrum.A comparison of our data with models of these types shows a satisfactory agreement with all of them (the best fit with Bahcall & Waxman's model), so we cannot distinguish between different scenarios of the ankle formation basing on the energy spectrum measurement alone.There is a need for additional data, presumably, concerning the mass composition of UHECRs, in order to elucidate the origin of an ankle in the energy spectrum.Russian Academy of Sciences.This work is supported by RFBR grants #06-02-16973, #05-08-50045.

Appendix A. Energy balance of EAS components
The energy fractions of the EAS primary particle transferred to the shower components can be described on the basis of hadron transport equations.If E k , (k = N, π, µν, eγ) is the energy transferred to nucleons, charged pions, muons+neutrinos, electrons+photons, a few cascade parameters determine the ratios between E k -i.e. the energy balance in the shower.For instance, the transport equation for the charged pions density π(x, E) at depth x is: where the interaction mean free paths λ π , λ N are assumed to be constant; w ππ (E, U ), w πN (E, U ) are the spectra of charged pions produced in pion-air and nucleon-air interactions, can be transformed (integrating over E 2 ) to: In the energy range E B π the only parameters to define the solution are K N /λ N and λ π .It means that the energy transferred to charged pions is independent of the spectra of pions produced in nuclear interactions.Hence, in the general case, we can use simple δ-model with the production spectrum w ik (E, U ) = n s δ(E − U/n s ), where n s is the multiplicity of secondaries, to balance the components energy in a shower.Because of the net value of E µ+ν /E 0 ≤ 0.1, the uncertainty due to the simplified model should be of second order of magnitude.
To summarize, the main model parameters governing the energy balance in the shower are the average inelasticity coefficients, mean free paths, multiplicity of secondaries and the fragmentation rate of the primary nucleus.The influence of other model characteristics such as 'the form of the rapidity distribution of constituent quarks' in collisions is weak.
The exact analytic solution of equations in the case of constant λ π , λ N , K N and B π = 0 is shown in figure A1 together with δ-model results.The letter is demonstrating the influence of rising cross-sections and the real decay rate of charged pions.An asymptotic (x = ∞) estimation of E eγ with the CORSIKA/QGSJET code at E 0 = 10 18 eV [23] is shown in the figure as well.

Appendix B. Conversion of the measured cosmic ray intensity to the spectrum of the primary beam
There is a correction to be applied to the measured intensity of CRs before any comparison of the energy spectra observed by different EAS arrays.
The quantity Ê = 'primary particle energy' that has been estimated after a shower detection, and the actual energy of the particle, E, which has initiated the EAS, are different values, connected with each other by a relation to be found.Estimated energy has distribution around the given mean value g( Ê, E), formed by the instrumental errors and fluctuations of the shower parameters with a RMS deviation, σ.Energy fluctuation is small in comparison with instrumental errors.Our aim here is to calculate the exact difference between the observed intensity of cosmic rays, J( Ê)d Ê, and the original one J(E)dE in the case of a rapidly falling power law spectrum.
The problem was solved, in general, by Zatsepin and Kalmykov in the previous century.The measured number of EAS particles, the so-called shower size, N e , is connected with the primary energy.The function g(N e , E) depending on E/N e has been found by Zatsepin [48].Then Kalmykov has calculated the measured intensity in the case of a lognormal distribution of N e [49]: The resultant initial-to-observed intensity conversion factor is ).
The necessary conditions are a constant index and RMS error, or at least both changing only slowly with energy.A distinct feature of the break in the spectrum is its shift along the energy scale.Due to the smeared transition between the conversion factors below and above the break, its position moves upward in energy.
To illustrate this we have simulated the primary energy spectrum with ankle and knee: random primary energies have been generated according to a broken power law with the index γ = 3.3 and 2.9 below and above 10 18 eV in the case of ankle, and γ = 2.9 and 5 for the knee.The energy estimation is modelled by adding randomly a Gaussian error to ln E with σ = 0.32. Figure B1 shows the 'measured' and primary spectra.The different shifts in the logarithm of intensity on both sides of the break are clearly seen, as well as that the position of the observed break moves to the right.

Figure 1 .
Figure 1.The detector arrangement of the Yakutsk array.Charged particle detectors (open circles), Cherenkov light detectors of the C 1 subset (filled circles) and the C 2 subset (filled triangles), and the muon detectors (squares) are shown.

Figure 2 .
Figure 2. Local area network of the Yakutsk array.

Figure 5 .Figure 6 .
Figure5.The half width of the Cherenkov light pulse as a function of the shower core distance.Experiment: points[3]; theory: solid line with 1σ errors (dashed lines)[7].

Figure 7 .
Figure 7.A setup to calibrate the air Cherenkov light detector.

Figure 8 .
Figure 8. Atmospheric extinction of the Cherenkov light as a function of Xmax.θ = 20 0 (circles).Approximation 3 is shown by the solid line.

Figure 11 .
Figure 11.Air Cherenkov light radial distribution.Our experimental data are given by the black and white points alternately, not to confuse the adjacent data.The curves are approximations according to Equation (4) with Q 150 fitted to the particular data.An axis distance R = 150 m is indicated by the dotted line.

- 2 )Figure 13 .
Figure13.The correlation between charged particle density at 600 m from the core, S 600 , and air Cherenkov light intensity at 150 m, Q 150 , measured in the same showers with θ < 15 0 .

2 )Figure 15 .
Figure 15.Acceptance area of air Cherenkov light detector subsets of the Yakutsk array.The solid curve is for autonomous, C 2 , and the dotted one for the medium, C 1 , sub-array area.

1 )Figure
Figure Cosmic ray flux measured by arrays equipped with air Cherenkov light detectors.Vertical bars indicate statistical errors.

1 )Figure 17 .
Figure17.The energy spectrum of cosmic rays around the knee region.Observational data are from arrays equipped with Cherenkov light and scintillation detectors.Two panels show the spectra before (left) and after (right) energy corrections applied to the data, as described in the text.
σ N is RMS deviation of ln N e ; κ is spectrum index; a N = E Ne dNe dE .Here we also assume the lognormal distribution of y = ln Ê with the average value equal to ln E. The observed intensity of cosmic rays is then given by the convolution of the primary spectrum, J(z) = J 0 exp(−κz), and the distribution of instrumental errors

Figure
Figure Simulation of the spectrum measurement with ankle and knee features.

Table 2 .
The primary energy fractions of EAS components: ionization loss of electrons in the atmosphere, E i , and in the ground, Eg; energy of muons at the ground level, Eµ; energy of EAS components unobservable at the Yakutsk array, E unobs .E 0 = E i + Eg + Eµ + E unobs .θ = 0 0 .
0 Figure 12.The electromagnetic component energy estimation from the Yakutsk array data: C 1 subset (squares), C 2 subset (rhombuses); solid curve is the result of CORSIKA/QGSJET based estimation.

Table 4 .
185±0.02,Figure14.S 300 as a function of x = 1020/ cos θ for different CR intensities.Open symbols (S i ) are equi-intensity method results, filled ones (Q i ) are derived fixing Q 400 .The highest energy EAS events detected with the Yakutsk array

Table 5 .
Energy estimation errors, σ, and intensity conversion factors, R J , for EAS arrays.

Table 6 .
Correction factors to energy scales for the pairs of EAS arrays/detectors, R E , averaged in the region E > 10 18 eV.