Large-scale distribution of cosmic rays in right ascension as observed by the Yakutsk array at energies above 10 18 eV
LLarge-scale distribution of cosmic rays in right ascensionas observed by the Yakutsk array at energies above 10 eV A.A. Ivanov, A.D. Krasilnikov, M.I. Pravdin, A.V. Sabourov
Shafer Institute for Cosmophysical Research and Aeronomy, 31 Lenin Avenue, Yakutsk 677980, Russia
Abstract
We present the results of searches for anisotropy in the right ascension ( RA ) distribution of arrival directions ofcosmic rays (CRs) detected with the Yakutsk array during the 1974–2008 observational period in the energyrange above 10 eV. Two methods of analysis are applied to two sub-samples of the data. Particularly,estimations of the first and second harmonic amplitudes are given, as well as the first harmonic phase inadjacent energy intervals. Analysis of variance demonstrates a significant contraction of the minimal widthof the RA distribution in the energy bin (10 , . × ) eV with respect to the isotropic distribution,which may be attributed to a possible source of CRs within the interval RA ∈ (15 , ). Keywords: cosmic ray, extensive air shower, harmonic analysis, analysis of variance
PACS:
1. Introduction
A conventional approach to shed light on the origin of cosmic rays (CRs) is to search for anisotropy in thearrival directions distribution. A number of attempts have been made to find excess fluxes of CRs correlatedwith large-scale structures in the nearby Universe, resulting in several indications of possible anisotropiceffects, but none of these effects have been confirmed independently (some examples can be found in Refs.[1, 2, 3]).By restricting the data under analysis exclusively to ground-based arrays, an essentially uniform exposurein the right ascension can be obtained. Under this restriction, in the present paper we examine the dataobtained from scintillation counters–surface detectors of the Yakutsk array. The results of previous analysisof the RA distribution of extensive air shower (EAS) primaries detected within the time period from January1974 to May 2000 were published in Refs. [4, 5, 6]. We have now analyzed the extended dataset up to June2008. Additional data include, specifically, 7598 EAS events above the threshold energy 1 EeV (=10 eV)[7]. The main aim of this analysis was to extend the time series of the data by using an observational time thatwas as long as possible. These efforts were aimed at testing the previous indications of possible anisotropyin the arrival directions of CRs in the extended dataset using different methods.
2. The Yakutsk array experiment and data selection for analysis
The Yakutsk array site is located near Oktyomtsy, the satellite village of Yakutsk, at geographical co-ordinates 61 . N, . E and at a mean altitude of 105 m above sea level. At present, it consists of 58ground-based and 4 underground scintillation counters to measure charged particles (electrons and muons),and 48 detectors of the air Cherenkov light. During its 40-year lifetime, the array has been reconfiguredseveral times. Before 1990, the total area covered by detectors was at its maximum ( S ∼
17 km ); now, it is8 . .EAS events were selected from the background using a two-level trigger of detector signals: The firstlevel is a coincidence of signals from two scintillation counters in a station within 2 µ s; the second level is acoincidence of signals from at least three nearby stations (not lined up) within 40 µ s. Stations spaced ∼ ∼ Email address: [email protected] (A.A. Ivanov)
Preprint submitted to Elsevier December 6, 2018 a r X i v : . [ a s t r o - ph . H E ] J u l X , m Y , m Figure 1: Arrangement of the Yakutsk array detectors (stations forming trigger-500 and trigger-1000; circles) and the showercores selected before (crosses) and after 1990 (triangles). period, more than 10 showers of primary energy above 10 eV were selected. Fig. 1 shows the positions ofdetectors of the Yakutsk array and the core positions of EASs with E > ∼
30 m for trigger-500 and ∼
50 m for trigger-1000 events.Arrival angles of the EAS primary particles were calculated in the plane shower front approximation usingdetection times at the stations. A clock pulse transmitter at the center of the array provided pulse timing to100 ns accuracy. Errors in arrival angles depend on the primary energy decreasing from ∼ at E = 1 EeVto ∼ above E = 10 EeV. More detailed information can be found in Refs. [8, 9, 10].In this work, a sample of the analyzed dataset consisted of EAS events detected in the period January1974 – June 2008 within the array area, with energy above 1 EeV, at zenith angles θ < .The energy estimation method was based on the total flux measurement of the air Cherenkov light andthe number of electrons and muons at observation level [8, 11, 12]. To unify the energies of showers detectedin a variety of years, we used the same S -to-energy relationship and attenuation length as in Ref. [5]. Thenumber of EAS events selected was 43710 at E > E v e n t r a t e Declination, degE>1 EeV E v e n t r a t e Declination, degE>5 EeV E v e n t r a t e Declination, degE>10 EeV E v e n t r a t e Right ascension, degE>1 EeV E v e n t r a t e Right ascension, degE>5 EeV E v e n t r a t e Right ascension, degE>10 EeV
Figure 2: Arrival directions in equatorial coordinates. Extensive air showers were selected with three threshold energies (his-tograms). Isotropic distributions are shown as dotted curves. A , % lg(E, eV) θ <60 θ <45 θ <30 A , % lg(E, eV) θ <60 θ <45 θ <30 Figure 3: The first (left panel) and second (right panel) harmonic amplitudes in the right ascension distribution of cosmic raysdetected in the period 1974–2008 at energies above 1 EeV. Amplitudes of the expected isotropic distribution ( θ < ) areindicated as dashed curves. Vertical bars are statistical errors; horizontal bars indicate energy bins. [13, 14].Fig. 2 shows the distribution of arrival directions of selected showers in the equatorial system comparedwith expected isotropic distributions.
3. Harmonic analysis
We used the Rayleigh formalism in harmonic analysis of the RA distribution of CR arrival directions[15]. If the large-scale anisotropy is present in the arrival directions of CRs, it can be expressed using theanisotropy coefficient for the intensity variation: A = I max − I min I max + I min , where I max and I min are the maximum and minimum of CR intensity, respectively, as a function of the arrivalangle, α . In the simplest cases with one source at α (for example, I max from the Galactic center, I min fromthe anticenter), the distribution can be described by a cosine wave I ( α ) = 1 + A cos( α − α ), or with the twoopposing sources as I ( α ) = 1 + A cos(2 α − α ) (for example, I max from the Galactic arm “in” and “out”directions, while I min is from perpendicular directions). In the first case, the coefficient is denoted A , whichis the amplitude of the first harmonic, and in the second case, the coefficient is denoted A , which is theamplitude of the second harmonic. In the general case of harmonic expansion of an arbitrary distribution: I ( α ) = 1 + ∞ (cid:88) k =1 A k cos( kα − kα ) . (1)If the observed distribution of arrival angles, α i , is given by a sum of delta functions f ( α ) = (cid:80) Ni =1 δ ( α − α i ),then A k = (cid:113) a k + b k , α = arctan( b k /a k ) , (2)where a k = N (cid:80) i cos( kα i ) , b k = N (cid:80) i sin( kα i ).Average values of the amplitudes are not zero even in the case of isotropic distribution, and converge tozero with N −→ ∞ (see Appendix A): A k = (cid:114) πN , σ ( A k ) = (cid:114) − πN . (3)The probability that isotropic amplitudes will be larger than the given A k is P ( > A k ) = exp( − N A k . (4)Earth’s rotation enables celestial regions to be represented as RA distributions. Scintillation countersof the Yakutsk array have a 24-hour duty cycle, resulting in almost uniform directional exposure. Small3 able 1: Yakutsk array observed amplitudes in energy bins, A k , and the probability that expected isotropic amplitudes willexceed observed amplitudes, P ( > A k ). θ ∈ (0 , ). Energy bins, A , P ( > A ), A , P ( > A ),lg( E, eV ) % % % %18.0-18.5 0.95 41.42 0.71 60.1018.5-19.0 2.08 64.03 0.99 90.3419.0-19.5 28.54 0.04 11.78 26.7519.5-20.0 16.89 74.10 38.13 21.73deviations from the uniformity are caused by maintenance of the detectors (mainly during working days) andsummertime shutdowns of the array due to thunderstorms.Diurnal and seasonal variations of the array exposure and atmospheric conditions result in spuriousamplitude (0 . ± . A k (cid:39) .
85% (Eq. 3). Therefore, in the energy region
E > RA distribution is not affected, because of diurnal smearing by the Earth’s rotation [18]. As a result, weignored the geomagnetic field effect at E > A k , k >
2, because we considered these to be insufficiently large at the angularscale. Observed amplitudes increased with energy, for both the first and the second harmonics, but theeffect was statistically insignificant against the background of “isotropic” amplitudes increasing because ofthe number of events decreasing with energy. Table 1 shows the probability that the amplitude of isotropicarrival directions will be larger than the observed amplitude by chance.A 3 . σ excess of the first harmonic amplitude over isotropic expectation was found in the energy bin(10 , .
8) EeV. Variation of zenith angle threshold did not eliminate the effect (Fig. 3) as well as doubling ofthe bin width (Table 1). It is interesting to note that the amplitude and phase of the first harmonic in thisenergy interval measured before 2000 and throughout the period had practically the same value [4, 6].Time variation of A (10 < E < .
6) observed by the Yakutsk array in six-year intervals is shown in Fig.4, left panel. Anisotropic amplitudes were observed before 1985 and after 1995. The variation of amplitudecan be attributed to the fluctuations of the EAS event number in reduced intervals of years.The statistical significance of the excess flux found a posteriori in a particular energy bin should beestimated using the penalty factor that is at least equal to the number of independent bins. In this case, we A / A i s o t r o p i c Years
E in (10,31.6) EeV A / A i s o t r o p i c lg(E WG , eV) PAO, 2004-2012Ykt, 1974-2008TA, 2008-2013
Figure 4: Ratio of observed and expected amplitudes of the first harmonic. Left panel: in six-year intervals of observation atYakutsk. Right panel: comparison of the Yakutsk array (Ykt), Pierre Auger Observatory (PAO), and Telescope Array (TA)results. Energy, E WG , is scaled using Working Group factors: 0.56 for Yakutsk, 0.91 for TA, and 1.1 for PAO energy estimations[14]. R A , d e g lg(E WG , eV) Ykt, 1974-2008PAO, 2004-2012TA, 2008-2013 -180-9009018017.5 18 18.5 19 19.5 20 R A , d e g lg(E WG , eV) Figure 5: The first harmonic phase as a function of energy. Left panel: comparison of data observed by the Yakutsk array(Ykt), Telescope Array (TA) [20], and by the Pierre Auger Observatory (PAO) [19]. Right panel: comparison of observationsat Yakutsk before and after 2000. consider an excess to be equally probable by chance in any of the bins.An alternative approach is to divide the observational data into two sub-samples: EAS events detectedbefore and after 2000. In this case, we are basing on the analysis of the Yakutsk array data consisting oftwo independent parts – first results (1974-2000) are published in 2001 [4], and the second part (2000-2008)is published in 2013 [7]. The first group was used to find the energy interval where the amplitude exceedsisotropic expectation, i.e., E ∈ (10 , .
6) EeV, where the first harmonic amplitude is (26 . ± . A = (41 . ± P ( > A ) = 0 .
118 in the same energy interval, indicating no possibility of anisotropyas suggested by the first harmonic amplitude in the RA distribution.No statistically significant deviations of harmonic amplitudes were found in the data from the PierreAuger Observatory (PAO) [19] and Telescope Array (TA) [20]. The data of PAO and TA are comparedwith our results in Fig. 4, right panel. The TA collaboration found no deviation from isotropy with energythresholds 10 EeV and 40 EeV. However, in the highest energy bin E >
57 EeV, they claim observation of ahotspot, with a statistical significance of 5.1 σ , centered at RA = 146 . , Dec = 43 . [20].Another item of interest can be found in the considerations of the PAO collaboration [19]. It was notedthat their phase measurements in adjacent energy intervals did not appear to be randomly distributed, butrather indicated a smooth transition between a common phase 270 consistent with a Galactic center regionbelow 1 EeV and another phase consistent with the RA of the anticenter above 5 EeV. This is potentiallyinteresting, because with a real underlying anisotropy, a consistency of the phase measurements in orderedenergy intervals is indeed expected to be revealed with a smaller number of events than required to detectthe amplitude with high statistical significance [19].We compared phases of the first harmonic in energy bins before and after 2000 (Fig. 5, right panel), and ourmeasurements in these two periods appear to be qualitatively consistent with the PAO conclusion. Althoughthe phase uncertainties in our data are relatively large, the phase of the first harmonic is not randomlydistributed, at least above 10 EeV, but appears to increase gradually in the direction of the anticenter (Fig.5, left panel).Additionally, analysis of the combined datasets of the PAO, TA, and Yakutsk array by the AnisotropyWorking Group for the CERN Symposium [21] also hinted at the same intriguing regularity in the phase ofthe dipole anisotropy. One conclusion was that larger statistics are needed to investigate these observations.Our present contribution in strengthening this hint is that the regularity is observed consistently indifferent time intervals and arrays.
4. Analysis of variance
We also tested for the probability of anisotropy in RA distribution using another method, i.e. analysisof variance. This method is able to determine systematic differences between the results of measurementscarried out under specific varying conditions [22].The distinctive feature of the isotropic distribution in RA , our null hypothesis, H , is that it has nomean value and the variance s i = π / s m = min (cid:113) ( α − α ) , where the trial mean α scans the whole α ∈ (0 , π ) interval. The5 able 2: Minimal width of the right ascension distribution of cosmic rays detected with the Yakutsk array in the two periods.Column heads: event numbers, N; observed, s m , and expected isotropic widths, s im ; chance probability that isotropic width isless than or equal to observed width, P . Energy bins, 1974-2000 2000-2008lg(
E, eV ) N s m , s im , P , N s m , s im , Pdeg deg % deg deg %18.00-18.25 17063 103.66 103.43 86.06 12081 103.20 103.34 28.0318.25-18.50 6765 103.14 103.14 45.31 3258 102.62 102.80 32.4318.50-18.75 2410 102.03 102.60 16.26 772 100.48 101.59 14.7318.75-19.00 761 102.94 101.57 93.99 178 102.21 99.02 96.8519.00-19.25 233 95.27 99.66 2.38 52 84.97 94.77 2.2419.25-19.50 72 96.88 96.17 52.31 23 100.92 89.95 99.6919.50-20.00 36 96.25 92.87 71.40 6 82.73 74.88 66.43distance between two points can be calculated directly, | α − β | , or around a circle, 2 π − | α − β | . We chosethe minimal of the two for all pairs. So, under the null hypothesis, the width in degrees is s im = 103 . . Onthe other hand, if there is a single source, SS, of CRs with width s (cid:28) s im in the isotropic background (ouralternative hypothesis, H ), then the aggregate distribution width may be sufficiently narrow depending onthe relative luminosity of the source.For example, if the SS fraction of the total flux is 50%, then the resulting width is s m = 73 . . This canbe easily identified by the analysis of variance. Even the flux fraction 10% from the single source resultingin width s m = 98 . can be distinguished with a sufficient number of EAS events.It is convenient to calculate the arrival directions and probabilities of the distribution widths under nulland other hypotheses by using the Monte Carlo method. In this case, we can apply the same procedures toform distributions and to calculate variances of the experimental data and random points in the uniform RA distribution. An example of the program is given in Appendix B.The data of the Yakutsk array were sampled in seven energy bins separated by lg E i = 0 . i, i = 0 , .., E i is in units of EeV. For each sample, the minimal width of the RA distribution was found (Table 2)compared with the expected width, s im , for an equal number of isotropic events. To avoid the penalty factorin the probability, we divided the dataset into two independent sub-samples: data observed before and after2000 (the terminal date is 05/31/2000). The former was used to find a bin with minimal distribution widthobserved; the latter was used to estimate the statistical significance in the energy bin fixed a priori fromthe first. The minimal width of RA distribution of 52 EAS events observed in the period 2000–2008 withenergies E ∈ (10 , .
8) EeV was s m = 85 . The chance probability that the isotropic distribution width ofthe 52 events would be less than or equal to the observed width was P ( ≤ s m ) = 2 . .
8% confidence level.In Fig. 6, our results covering the whole period 1974-2008 are shown together with that derived from thepublished data from PAO [23] and TA [20] above ∼
55 EeV. Consistent with the results of harmonic analysis,there was an energy interval (10, 17.8) EeV where the observed minimal width of the Yakutsk array datawas distinctly less than the isotropic expectation.Our alternative hypothesis, H , has two parameters to fit the experimental data: the source position in RA , and the fraction of the total flux produced by the source luminosity. We estimated the most probable H parameters fitting the data from Table 2 in the energy range (10 , .
8) EeV. The results showed that theSS position is α ∈ (15 , ), and the ratio of CR flux from SS to the total flux is I SS /I tot = 0 . ± .
5. Conclusions
We used two methods to examine the RA distribution of CR arrival directions measured with the Yakutskarray: harmonic analysis and analysis of variance. Resultant first and second harmonic amplitudes increasedwith energy but were consistent with expected amplitudes of the isotropic distribution. In the energy range(10 , .
6) EeV, our data observed in the period 1974–2000 exhibited an excess flux with the first harmonicamplitude of 26 . . A = 41 .
4% with a chance probability of ∼ s m - s m i , d e g lg(E WG , eV) TA 2008-2013Ykt 1974-2008PAO 2004-2009
Figure 6: Difference in observed and expected isotropic distribution widths in right ascension. Statistical errors for observedevent numbers are shown by the vertical bars.
Analysis of variance demonstrated the prominent excess in the same energy bin, namely, a significantcontraction of the minimal width of the RA distribution of CR arrival directions with respect to the isotropicdistribution. Downsizing in width was found in both independent parts of the Yakutsk array data, i.e.,data observed before and after 2000. The null hypothesis was rejected at the significance level ∼ , .
8) EeV with the source position RA ∈ (15 , ) and the ratioof CR flux from SS to the total flux I SS /I tot = 0 . ± . − , ), aswould be expected in the isotropic case, but exhibited a gradual increase with energy in RA , at least in theenergy interval above ∼
10 EeV. This behavior is inherent in both Yakutsk observation periods, i.e., beforeand after 2000, and is in agreement with the possible regularity in the phase of the dipole anisotropy observedby the PAO and TA collaborations.
Acknowledgments
We are grateful to the Yakutsk array staff for the data acquisition and analysis. The work is supportedin part by the Russian Academy of Sciences (Program 10.2) and RFBR (grants 11-02-00158, 13-02-12036).
Appendix A. The first harmonic amplitude of the isotropic right ascension distribution
Here, we apply the Rayleigh formalism, founded in [15], to illustrate a calculation of the first harmonicamplitude, A , for N isotropic points in the interval α ∈ (0 , π ).Treating the amplitude as a sum of N vectors of length | (cid:126)r i | = 2 /N with angle α i , we have the totallength, A = | (cid:126)R N | , which accumulates N equal random steps. The inductive argument about a relationshipbetween R i − and R i in a triangle of vectors (cid:126)R i = (cid:126)R i − + (cid:126)r i consists in R i = R i − + r i due to cos α i = 0.Consequently, A = 4 /N .Asymptotically, the distribution of the amplitude is circular Gaussian according to the central limittheorem: P ( (cid:126)R N ) = πσ exp( − A σ ), where N (cid:29)
1. The mean amplitude is A = 1 σ (cid:90) ∞ A exp( − A σ ) dA = √ σ Γ(1 .
5) = (cid:114) πσ , while the mean square of amplitude is A = 1 σ (cid:90) ∞ A exp( − A σ ) dA = 2 σ Γ(2) = 2 σ. A comparison with the vector length gives the following: σ = 2 /N , the mean amplitude A = (cid:112) π/N , andthe variance A − A = (4 − π ) /N . 7he probability of obtaining an amplitude greater than or equal to A is P ( ≥ A ) = 1 σ (cid:90) ∞ A A exp( − A σ ) dA = exp ( − N A . Appendix B. A Monte Carlo program to find the minimal width of RA distribution
The Fortran-90 program below illustrates the minimal width calculation for a distribution of N randompoints in the right ascension (RA) circle. The function Var computes the minimum variance with the trialmean scanning a circle (0 , ). The mean value and standard deviation of the distribution minimal widthare calculated with a sample of size M in the main program W.Program W;real Q(10000);N=52;M=100000 ! N points sampled M timesav=0;d=0;do k=1,Mdo i=1,N;Q(i)=RAN(ir)*360.0;enddo;rms=sqrt(Var(Q(1:N),N)) ! N random RA pointsav=av+rms/M;d=d+rms**2/M;enddo;d=sqrt(d-av**2) ! mean & deviation, degreesprint *,N,M,av,d;end programfunction Var(Q,N);integer N;real Q(N) ! Minimal variance of Q(N) points in a circle (0,360)smin=1e36;do k=1,36;A=k*10.0 ! A=trial means=0;do i=1,N;d=abs(Q(i)-A);d=min(d,360.-d);s=s+d**2/N;enddo ! var of Q(N) for Aif(s < =smin)smin=s;enddo;Var=smin;end function. References [1] A.M. Hillas, Astropart. Phys. 32 (2009) 160.[2] P. Sommers, S. Westerhoff, New J. Phys. 11 (2009) 055004.[3] A.A. Ivanov, Nucl. Phys. B (Proc. Suppl.) 190 (2009) 204.[4] A.D. Krasilnikov, A.A. Ivanov, M.I. Pravdin, in: Proc. of 27 th ICRC, Hamburg, 1 (2001) 398.[5] V.P. Egorova et al., J. Phys. Soc. Japan, B 70 (2001) 9.[6] A.A. Ivanov, A.D. Krasilnikov, M.I. Pravdin, JETP Lett. 78 (2003) 695.[7] A.A. Ivanov et al., in: Proc. of 33 d ICRC, Rio de Janeiro, icrc2013-0270, 2013.[8] A.A. Ivanov, S.P. Knurenko, I.Ye. Sleptsov, New J. Phys. 11 (2009) 065008.[9] A.A. Ivanov, S.P. Knurenko, M.I. Pravdin, I.Ye. Sleptsov, Moscow Univ. Phys. Bull. 65 (2010) 292.[10] Website of the Yakutsk array is at http://eas.ysn.ru.[11] A.A. Ivanov, S.P. Knurenko, I.Ye. Sleptsov JETP 104 (2007) 872.[12] M.I. Pravdin et al., Bull. Russian Acad. Sci: Phys. 71 (2007) 445.[13] A.A. Ivanov, Astrophys. J. 712 (2010) 746.[14] B.R. Dawson et al., EPJ Web of Conf. 53 (2013) 01005.[15] J.W. Strutt (Lord Rayleigh), Phil. Mag. 10 (1880) 73.[16] M.I. Pravdin et al., JETP 92 (2001) 766.[17] A.A. Ivanov et al., JETP Lett. 69 (1999) 288.[18] V.P. Egorova et al., in: Proc. of 26 th ICRC, Salt Lake City, 1 (1999) 403.[19] I. Sidelnik et al., in: Proc. of 33 d ICRC, Rio de Janeiro, icrc2013-0739, 2013.[20] M. Fukushima et al., in: Proc. of 33 d ICRC, Rio de Janeiro, icrc2013-1033, 2013; arXiv:1404.5890.[21] O. Deligny et al., EPJ Web of Conferences, 53 (2013) 01008.[22] A. Hald,