Cosmic Ray Energy Spectrum derived from the Data of EAS Cherenkov Light Arrays in the Tunka Valley

The extensive air shower Cherenkov light array Tunka-133 co lle ted data during 7 winter seasons from 2009 to 2017. From 2175 hours of data taking, we derived t h differential energy spectrum of cosmic rays in the energy range 6 · 1015 2 ·1018 eV. The TAIGA-HiSCORE array is in the process of continuous e xpansion and modernization. Here we present the results obtained wit h 28 stations of the first HiSCORE stage from 35 clear moonless nights in the winter of 2017-2018. The combin ed spectrum of two arrays covers a range of 2 · 1014 − 2 · 1018 eV.

Measurement of the energy spectrum and mass composition by the data of these arrays is important in order to understand the acceleration limit of the Galactic CR sources and the transition from Galactic to extragalactic CR.
Tunka-133 array contains 175 optical detectors, spaced over an area of 3 km 2 . Each detector contains a single * e-mail: v-prosin@yandex.ru PMT with a hemispherical photocathode of 20 cm diameter.
The next step of the experiments in the Tunka Valley is gamma-ray observatory providing search and study of individual cosmic ray sources. The observatory is called TAIGA (Tunka Advanced Instrument for cosmic rays and Gamma Astronomy). The TAIGA complex is designed to study gamma radiation and the charged cosmic rays in the energy range of 10 13 −10 18 eV [6]. The final version of the observatory will include a network of wide field of view (0.6 sr) timing Cherenkov light stations, called TAIGA-HiSCORE (High Sensitivity Cosmic Ray Explorer) [7,8], and up to 16 imaging atmospheric Cherenkov telescopes (IACT) (Fig. 2), covering an area of 10 km 2 . The capabilities of these Cherenkov arrays will be enhanced by muon The prototype of TAIGA consists now (2018) of 54 TAIGA-HiSCORE stations and one IACT. The layout of the TAIGA-HiSCORE stations is shown at the Fig. 3. The optical stations are distributed in a regular grid over an area of 0.5 km 2 with an inter-station spacing of 106 m. All stations are tilted into the southern direction by 25 • to increase the time for observarion of the first test objectthe Crab Nebula. Each optical station contains four PMTs with 20 or 25 cm diameter, namely EMI ET9352KB, or Hamamatsu R5912 and R7081. Each PMT is supplied with the Winston cone of 0.4 m diameter and a 30 • viewing angle (field of view is ∼ 0.6 sr). A detailed description of DAQ and synchronization systems is given in [13] The Tunka-133 array collected data during the 7 winter seasons of 2009 -2014 and 2015 -2017, accumulating information for 350 clear moonless nights. The total data set time is 2175 hours. The TAIGA-HiSCORE array is in the process of continuous expansion and modernization. Here we present the data obtained using the 28 stations, shown as the filled squares in the Fig. 3 for 35 clear moon- hours. The experimental data are processed using codes in which all approximating and recalculating functions are obtained from the analysis of artificial events generated by the CORSIKA code for the energy range 10 14 -10 18 eV [2,10,11]. For each shower, the arrival direction, the core position on the observation plane, and the energy of the primary particle are reconstructed. As a result, the combined differential all-particle energy spectrum in the energy range of 2 · 10 14 − 2 · 10 18 eV was obtained.

Data processing and EAS parameters reconstruction
Data processing for the Tunka-133 array is described in [2,3]. The position of the EAS core is reconstructed by fitting the measured values of pulse amplitudes with the amplitude-distance function (ADF) [2]. An example of CORSIKA simulated ADF is shown in Fig. 4. The shape of ADF is rather complicated with change of a function type at the distances R kn , 200 m and 400 m. But we succeeded in [2] to describe all the variable detailes with a single variable so called steepness b A . The shower arrival direction, characterized by the zenith (θ) and azimuthal (φ) angles of the axis, is reconstructed by fitting the measured time delays with a special curved front [12].
The shower energy is reconstructed from the the Cherenkov light flux density at a distance of 200 m from the core (Q 200 ). For interpolation to 200 m from the measured values of Q i , the special lateral distribution function (LDF) described in [11] is used. This function in contrast with ADF connects pulse area or the light flux density with EAS core distance. The relationship of energy with Q 200 is also obtained from the CORSIKA simulation [2]. The new version of this simulation is shown in Fig. 5. During this simulation the modern model of particle interaction QGSJET-II-04 was used. The simulation is made for 2 zenith angles 30 • and 45 • , 2 sorts of the primary particles (p and iron) and the raw of 9 energies from 3 · 10 14 to 10 18 eV.
It is essential to emphasise that EAS Cerenkov light flux reflects the intgral of the EAS cascade curve and thus don't depend on the interaction model assumed for the simulation. Figure 4 shows that the conversion from Q 200 to the primary energy almost don't depend on the primary mass and the zenith angle of the EAS. The new simulation confirms the previous one: The main EAS parameters according to the TAIGA-HiSCORE array data are reconstructed using the same algorithms and fitting functions as for the Tunka-133 array. In particular, for showers with energies greater than 10 15 eV, the shower energy is reconstructed from the Cherenkov light flux density at the core distance of 200 m Q 200 .
The effective area for the selection of events using this method of core reconstruction is the area of the ellipse with the semiaxes of 300 and 225 m shown in the Fig. 3 with the black curve.
For the energy range less than 10 15 eV, not all showers have a measurement of the light flux at the core distance of 200 m. Therefore, another method was developed for reconstructing the energy from the closest to the core detector readings. The EAS core position in this case is found as the gravity center of the measured amplitudes at the 4 stations closest to the core. The calculation shows that with existing geometry, the light flux density is mea- sured for this case on average at the core distance of 70 m. The correlations of Q 70 with the shower zenith angle and the primary energy were found from experimental data for the energy range 10 15 -3 · 10 15 eV, in which both Q 70 and Q 200 can be measured for each shower.
So derived recalculation of the Q 70 from the measured zenith angle to the vertical direction: log 10 (Q 70 (0)) = log 10 (Q 70 (θ)) + 1.13 · (sec(θ) − 1) The conversion from Q 70 (0) to the primary energy derived from the experimental correlation is as follows:  The method of reconstruction of the EAS core position as the amplitude gravity center leads to the large errors at the edge of the location of the array stations. To obtain the undistorted spectrum, a strip of 50 m width at the edge of the array is excluded from the effective area. Thus the area of the ellipse with semiaxes of 250 and 175 m is used as the effective area. This reduced ellipse is shown in Fig. 3 with a green curve.
The constants in the equations for Q 200 and Q 70 are obtained from the absolute energy calibration. The absolute primary energy calibration is carried out, as for all Tunka Valley Cherenkov light experiments, by normalizing the obtained every night integral spectra to the integral spectrum obtained in the Tunka-25 experiment [14], normalized in its turn to the absolute intensity of cosmic rays, obtained in the QUEST experiment [15].
The variations of these coefficients are treated as the reflection of the atmospheric transparency variations. The atmospheric transparency monitoring is provided by measuring the count rate of EAS (4 stations coincidence inside the time gate of 10 µs) during every 100 sec. As the integral spectrum index is about 2, the relative transparency coefficient is: The periods of stable count rate are selected for data processing. We use the data of clear nights only. Wile data processing the transparency coefficients are estimated for the total time of the processing period. The maximum single night deflection from the mean value of the transparency coefficient is less than 7% only.
It is essential to emphasise that our energy measurements are relative. Connection with absolute energy scale is made by normalizing to the known integral cosmic ray flux at the fixed energy only. This method let us avoid such problems as measurement of absolute sensitivity of Cerenkov light detectors, the study and monitoring of the detailed parameters of atmosphere and deconvolution of spectra taiking into account the energy resolution. The normalization point is close to the knee where many different experiments give the integral cosmic ray flux very close to that measured at the QUEST experiment [15].

Energy Spectrum
To construct the spectrum, the results of processing the Tunka-133 array were used to record events with zenith angles θ ≤ 45 • and the core position in a circle with a radius of R c ≤ 450 m for energies E 0 ≤ 10 17 eV and in a circle of radius of R c ≤ 800 m for showers with energy E 0 ≥ 10 17 eV. The event selection efficiency reaches ≥ 95% for energies E 0 ≥ 6 · 10 15 eV for a circle with a radius of 450 m and for energies slightly less than 10 17 eV for a circle with a radius of 800 m. Thus, about 375,000 events were used to construct the spectrum. Among them about 4200 events have the energy more than 10 17 eV.
To construct the spectrum, events with a zenith angle θ ≤ 45 • were selected based on the results of processing the TAIGA-HiSCORE array. The spectrum contains more than 170,000 events with energies greater than 10 15 eV and about 700,000 events in the energy range of 3 ·10 14 −10 15 eV. Points in the range of 2 ·10 14 −3 ·10 14 eV are built according to one night data with a uniquely good transparency (10/28/2018) and contain 29000 events. The resulting combined differential energy spectrum is shown in Fig. 3 together with the spectrum of the Tunka-25 array [14].
The latter experiment has a much lower energy threshold and can be used for experimental estimation of the Tunka-133 efficiency. The differential efficiency evaluated as the ratio of the flux recorded with Tunka-133 to the flux recorded with Tunka-25 as a function of primary energy is shown in Fig. 4. Only the statistical errors of this ratio are shown in Fig. 4. The experimental estimation is compared to the simulated efficiency. The MC simulation was made for some fixed energies assuming the mean apparatus parameters. The main simulation assumption is two or more hit clusters, because a single cluster event does not provide the Cherenkov light flux measurement at a core distance 200 m, used for the energy measurement. Figure 9 shows that a 95% efficiency of EAS registration is reached for E 0 ≥ 6 · 10 15 eV for the effective area of 450 m radius. The same high efficiency is reached for an energy somewhat less than E 0 = 10 17 eV for all events with R c ≤ 800 m. Expanding the energy range we used only events with small R c (R c ≤ 450 m) to reconstruct the spectrum for E 0 < 10 17 eV), and all events to reconstruct the spectrum for E 0 ≥ 10 17 eV. To estimate the intensity for the most energetic points of the spectrum we doubled the bin size from 0.1 to 0.2 in log E 0 for energies larger than 3 · 10 17 eV.
The initial part of the spectrum obtained from TAIGA-HISCORE data in the energy range (2 · 10 14 − 3 · 10 15 eV) can be approximated by a power law with an index of 2.73±0.01. In addition to the statistical error, there may be a systematic error associated with a possible inaccuracy in the index of the recalculation formula from Q 70 to energy.  At high energies, the spectrum exhibits a number of deviations from a power law. In the range 3 · 10 15 − 6 · 10 15 eV, a gradual steepening of the spectrum occurs. The subsequent points up to the energy of 2 · 10 16 eV can be approximated by a power law with the index 3.28 ± 0.01. Further, the spectrum becomes essentially flatter, and in the range of 2·10 16 ≤ E 0 ≤ 3·10 17 eV, in general, does not contradict to the power law with the index γ = 2.99 ± 0.01. At high energies, the index sharply increases to γ = 3.34 ± 0.09 Figure 10. Comparison of energy spectra of different experiments in the wide energy range 10 14 -10 20 eV (the second "knee"). In Fig. 10 spectrum is compared with a number of other works. At the left edge, our spectrum is in agreement with the spectra of all particles obtained in direct experiments with balloon ATIK-2 [16] and satellite NUCLEON (CLEM) [17]. The most recent and the most abundant statistics in this field of energy has been obtained by the ground-based HAWC experiment [18] in the Mexican mountains. This spectrum coincides perfectly both with direct experiments and with our results. In the energy range of 10 16 -10 17 eV, there is agreement between the result of our work and the spectra of the KASCADE-Grande [19] and IceTop [20] arrays. Noticeable in Fig. 10 difference between these spectra and the spectrum of the Tunka-133 array can be eliminated by increasing the energy estimation by 3% for KASCADE-Grande and by the same decrease in the energy estimation for IceTop. Such shifts are significantly less than the absolute accuracy of these experiments. At energies that are extremely large for the Tunka-133 experiment, our spectrum does not contradict in the frame of statictic errors to the results of the Telescope Array (TA) [14] and PAO [14] experiments.

Conclusion
Thus, the joint spectrum of the Cherenkov arrays Tunka-133 and TAIGA-HiSCORE covers 4 orders of magnitude in energy using a single technique and demonstrates excellent agreement between the results of direct satellite and balloon experiments with the results of giant ground-based arrays.
The first seasons of operation of the TAIGA-HiSCORE and first TAIGA-IACT demonstrated good performance of the installation and showed yet preliminary but interesting results. During the winter season 2018-2019 the TAIGA configuration will include 54 operational wide angle stations arranged over an area of 0.5 km 2 , and one IACT. During the next year it is planned to finish deployment of the first stage of TAIGA with 110 TAIGA-HiSCORE stations and 3 IACTs on an area of 1 km 2 .