Upper limit on the fraction of the cosmic ray flux from a separable source in the energy interval (10 17 − 10 19 eV ) , estimated using data from the Yakutsk array

. Arrival directions of ultra-high energy cosmic rays (UHECRs) in the equatorial system, detected with the Yakutsk array in the energy range ( 10 17 − 10 19 eV), are re-analyzed using a new approximation for the zenith angle distribution of the event rate of extensive air showers (EAS). While the null hypothesis cannot be rejected based on the data used here, an upper limit on the fraction of cosmic rays from a separable source in the uniform background is derived as a function of declination and energy.


Introduction
Previous studies of cosmic ray (CR) arrival directions with the Yakutsk array have mainly been performed at the highest energies, where full trigger efficiency is reached. Our aim in this study is to examine the distribution in a wide energy range, 10 17 to 10 19 eV, with energy-dependent array exposure due to the absorption of showers in the atmosphere. This makes it possible to use lower-energy data in the analysis of CR arrival directions that could not previously be used.
For this purpose, we use the measured and expected-for-isotropy zenith angle distributions described in [1]. These distributions are then transformed to equatorial coordinates in order to test the null and alternative hypotheses for the distribution of the arrival directions of CRs (Fig. 1). The Yakutsk array data are shown by points in energy intervals, and expected-for-isotropy distributions are shown by curves derived from the uniform and zenith angle distributions.
To simplify the treatment of the shower attenuation effects, we use a threshold for particle density at 600 m from the shower core, ρ 600 , of 0.1 m −2 . This is chosen to be well above the intrinsic instrumental thresholds of the array. One benefit of using this technique is the a posteriori selection of showers almost independently of the shower core position within the array area.
Another assumption used to derive the expected zenith angle distribution of the EAS event rate involves fluctuations in the particle density. More specifically, it has previously been shown that fluctuations in some shower parameters within a narrow * e-mail: ivanov@ikfia.ysn.ru energy bin, for example the shower sizes N e , N µ and particle density, can be approximated at sea level by a log-normal distribution. In particular, it has been demonstrated with experimental data and a COR-SIKA simulation of the scintillation counter signal that y = ln(ρ 600 ) in EAS events can be approximated by a Gaussian. Assuming an isotropic flux of CRs in the energy range (0.1 − 10 EeV), we can derive an analytical expression to describe the zenith angle distribution of showers that survive after cutting at the particle density threshold and that reach detectors at sea level.
The estimation of the energy of the primary particle initiating an EAS is important [2][3][4][5]. In this paper, we use the most recent energy evaluation method proposed in [6]. The primary energy is estimated in the same way as for the SIBYLL-2.1 model [7]: where ρ 600 (0) is the particle density (m −2 ) at 600 m from the shower core, measured in a vertical shower. In inclined showers, ρ 600 (θ) is estimated using the constant intensity cuts method [8].
The selected sample of the Yakutsk array data consists of EAS events detected during the period January 1974 to June 2008, with axes within the stage II array area, energies between (10 17 − 10 19 eV), and zenith angles θ ∈ (0 0 , 60 0 ) [9].

Harmonic analysis of CR arrival directions in a set of α-rings
The uniform right ascension distribution determines the application of harmonic analysis. All harmonic   [15]); −90 0 < δ < 45 0 (PAO 2017 [17]); 0 0 < δ < 90 0 (Yakutsk 2002 [16]). The results of this work (from Table, smoothed by curves) are denoted by the δ-intervals to which they belong. Right panel: Upper limit on SS fraction in declination bins. An upper limit based on Telescope Array's data for the fraction of EeV protons of galactic origin (TA, [19]) and Wibig & Wolfendale's upper limit on the fraction of galactic light nuclei (WW, [20]) in the CR beam are given for comparison.
amplitudes are zero in the case of an isotropic distribution of arrival directions, assumed as a null hypothesis, H 0 , and the first harmonic amplitude becomes non-zero if there is a source of CRs [10]. The phase of the first harmonic indicates the source position. Other ways to analyze the arrival directions of CRs have been studied, for instance, in [11,12].
The main disadvantage of harmonic analysis is its constraints on the one-dimensional uni-or bi-modal distribution in right ascension. A series expansion of the declination distribution consists of a set of nonzero amplitudes, all of which must be taken into account. Consequently, spherical harmonics cannot be straightforwardly applied to equatorial coordinates in a 2D map of arrival directions.
One suitable approach is to use one-dimensional analysis in a sliced map, that is, right ascension rings in declination bins, for example α ∈ (0, 360 0 ), 15 0 i ≤ δ < 15 0 (i + 1), where i = 0, .., 5. In this case, the declination distribution determines the relative number of CRs falling into rings and the resolving power of the array in a declination bin. In particu-lar, the Yakutsk array is blind to declinations below δ = 1.7 0 if showers are detected at θ < 60 0 .
The formalism of harmonic analysis in α-rings is the same as in previous studies (e.g. [9] and references therein). Namely, the amplitude and phase of the k-th harmonic are calculated by formulating the right ascension distribution as a sum of delta functions where a k = 2 N Σ i cos(kα i ); b k = 2 N Σ i sin(kα i ). The only difference lies in the sorting of arrival directions into declination bins.

Testing for uniformity of CR arrival directions
The resultant first harmonic amplitude in units of the isotropic amplitude, A iso 1 = π/N [9] is shown in the left panel of Fig. 2, in energy and declination bins. The null hypothesis cannot be rejected in any bin above 10 18 eV since the chance probability is greater than the critical value of 1%. However, in the declination bin (45 0 , 60 0 ) at energies in the range (10 17 − 10 18 eV), there is a significant deviation from isotropy, which is in qualitative agreement with the χ 2 test results. As a consequence, conventional harmonic analysis in right ascension with integrated declinations also exhibits a deviation from isotropy. A differential method has previously been proposed for reconstructing a genuine large-scale pattern that is obscured by instrumental and meteorological variations in the array exposure, known as the 'East-West method' [13,14]. Effects of experimental origin, which are independent of the incoming direction, are removed by subtracting the counting rates for showers coming from the East and West sectors. The disadvantages of this method are its reduced sensitivity with respect to the conventional Rayleigh test, and that it overlooks the details of CR arrival directions, apart from the hemisphere involved.
This method was successfully applied to PAO data below 1 EeV in order to derive the first harmonic amplitude and phase (−90 0 < δ < 25 0 ), without applying a correction for effects of instrumental or atmospheric origin [15]. In our case, however, the zenith and azimuth angles of CRs are measured and transformed to equatorial coordinates in order to set upper limits on the fraction of the CR flux from a separable source vs the uniform background in declination bins. Consequently, the East-West method can only be used for a distribution that is integrated over the observable declination range, to compare results with previous data.
Instead, as a way of distinguishing the instrumental and seasonal origins of an effect from astrophysical sources, the data set can be divided into seasonal subsets and the phase stability tested [16]. The results are presented in Table 1 for the declination bin in which the deviation from isotropy is located. The seasonal variation of the phase is evident, meaning that the observed anisotropy cannot be attributed to an extraterrestrial source with fixed angular coordinates.
The previous results of harmonic analysis in right ascension were derived with CR arrival directions integrated over the whole observable declination range due to the unknown zenith angle/declination distribution of the absorbed showers. The first harmonic amplitudes available in the energy range (0.1 − 10 EeV) are compared in Fig. 2.
Pravdin et al. [16] set an upper limit on A 1 using an anti-sidereal time vector caused by seasonal variations in the EAS event rate. Our results are in agreement with these limits, except for the declination interval (45 0 , 60 0 ) discussed above.
In spite of the different regions of observable sky, the data from PAO and the Yakutsk array demonstrate the first harmonic amplitudes, which do not significantly exceed the isotropic expectation at energies below 8 EeV [15]. At higher energies, however, the PAO collaboration reveals a large-scale anisotropy in the arrival directions of CRs, A 1 /A iso 1 = 4.8, which is more than the 5.2σ level of significance [17]. Unfortunately, the anisotropy dipole points to the declination δ = −24 +24 −13 degrees, which is invisible to the Yakutsk array.
In general, there is no statistically significant deviation from isotropy in the arrival directions of CRs exceeding instrumental and seasonal effects in any of the energy and declination bins observed by the Yakutsk array, and the null hypothesis is therefore not rejected based on the data analyzed here.

Upper limit on the fraction of CR flux from a separable source
We use an alternative hypothesis, H 1 , in which a separable source of CRs, SS, gives a fraction f SS , while all other sources form the uniform background, providing a fraction 1 − f SS of the total flux. The mass composition of astroparticles from the hypothesized SS is constrained by deflections in the galactic magnetic field; these are thought to be neutral particles (such as photons or neutrinos) if a source is not nearby. At energies above 1 EeV, however, the magnetic dispersion of protons is less than the declination bin width of 15 0 , and protons may therefore be supposed to be emitted.
If the alternative hypothesis is actually implemented, and the amplitude of H 1 is realized for the entire CR population, then two consequences are possible: a) the null hypothesis can be rejected, where exp(−NA 2 H1 /4) < 0.01 (estimating the statistical power of the Rayleigh test [18]); and b) an upper limit on the SS fraction can be set, where the observed amplitude, A exp 1 , is significantly less than the expected amplitude, A H1 .
In this paper, the second of these approaches is employed. The probability P (A 1 ≤ A exp 1 ) for the alternative hypothesis is calculated using the Monte Carlo method. Isotropic (1 − f SS )N points and f SS N in α SS , where N is the number of CRs detected in a particular α-ring, are sampled M = 10 5 times to calculate A 1 . An upper limit f thr SS is determined, where the probability is equal to P crit = 0.01, and thus the SS fraction above the threshold can be rejected.
The resultant upper limit on the SS fraction as a function of energy is shown in the right panel of Fig. 2. The energy dependence can be explained by the decreasing number of EAS events detected as the energy increases. The declination dependence is due to the irregular array exposure of SS that is assumed to lurk in a particular δ-bin.
As a reference, upper limits are shown for the fraction of EeV protons of galactic origin derived by the Telescope Array collaboration [19] and the fraction of galactic light nuclei in the CR beam set by Wibig and Wolfendale [20]. In general, these limits are comparable to conventional harmonic analysis results (marked '0 0 < δ < 90 0 ') and the SS fraction limit in the declination interval 0 0 < δ < 15 0 . In the case of conventional analysis, a separable source may be anywhere in the range δ ∈ (0 0 , 90 0 ), and thus the SS fraction limit is larger for the source luminosity.

Conclusion
The arrival directions of UHECRs detected with the Yakutsk array in an equatorial system are analyzed in the energy range (10 17 − 10 19 eV), based on the zenith angle distribution of EAS event rate measured using a ρ 600 threshold rather then numerous different detector density thresholds. A comparison of the measured and expected distributions is carried out in order to check for a statistically significant deviation from isotropy in the arrival directions of CRs.
A harmonic analysis in a set of α-rings bordered within δ bins is applied on a 2D equatorial map. Above 10 18 eV, the null hypothesis cannot be rejected, and this is in alignment with previous analysis of the Yakutsk array data. Below this threshold, in the energy interval (10 17 − 10 18 eV), there is significant deviation from uniformity in the declination bin δ ∈ (45 0 , 60 0 ), but this effect is a consequence of instrumental and seasonal variations in the array exposure, and cannot be associated with a CR source located somewhere in the equatorial region.
By making use of an alternative hypothesis containing a separable source in an otherwise isotropic set of CR sources, a stringent upper limit on the fraction of the total CR flux from such a source is established; at least, in declination intervals above δ = 15 0 .