Small-scale Cosmic-Ray Anisotropy Observed by the GRAPES-3 Experiment at TeV Energies

GRAPES-3 is a mid-altitude (2200 m) and near-equatorial (11.°4N) air shower array, overlapping in its field of view for cosmic-ray observations with experiments that are located in the Northern and Southern Hemispheres. We analyze a sample of 3.7 × 109 cosmic-ray events collected by the GRAPES-3 experiment between 2013 January 1 and 2016 December 31 with a median energy of ∼16 TeV for study of small-scale (<60°) angular-scale anisotropies. We observed two structures, labeled A and B, that deviate from the expected isotropic distribution of cosmic rays in a statistically significant manner. Structure A spans 50°–80° in R.A. and from −15° to 30° in decl. The relative excess observed in structure A is at the level of (6.5 ± 1.3) × 10−4 with a statistical significance of 6.8 standard deviations. Structure B is observed in the R.A. range 110°–140° and at decl. from −10° to 30°. The relative excess observed in this region is at the level of (4.9 ± 1.4) × 10−4 with a statistical significance of 4.7 standard deviations. These structures are consistent with those reported by Milagro, ARGO-YBJ, and HAWC. These observations could provide a better understanding of the sources of cosmic rays, their propagation, and the magnetic structures in our Galaxy.


INTRODUCTION
Cosmic rays (CRs) are charged particles because of which they undergo deflections by the randomized magnetic field in the interstellar medium through which they Corresponding author: P.K. Mohanty pkm@tifr.res.inpropagate.This leads to an isotropic distribution of CRs as observed on the Earth.However, in the past decade, several experiments located at different latitudes have observed the CRs anisotropy on different angular scales and energy ranges with an amplitude of ∼ 10 −4 to 10 −3 .Large-scale structures with a dominant dipolar component have been observed in TeV-PeV energy range by several ground based extensive air shower (EAS) experiments located in the Northern hemisphere such as Tibet-ASγ (M.Amenomori et al. 2006Amenomori et al. , 2017)), ARGO-YBJ (B.Bartoli et al. 2018), Milagro (A. A. Abdo et al. 2009), HAWC (A.U. Abeysekara et al. 2019), EAS-TOP (M.Aglietta et al. 2009), KASCADE-Grande (M.Ahlers 2019) and in Southern hemisphere such as IceCube (M.G. Aartsen et al. 2016).An excess has been observed within the right ascension (α) of 30 • to 120 • and a large deficit region has been observed within 150 • ≤ α ≤ 250 • .Anisotropies at ultra high energies (≥ 10 18 eV) where the interstellar magnetic field plays a relatively insignificant role for bending CRs, has been reported by the Pierre Auger Observatory (A.Aab et al. 2018) and the Telescope Array (R. U. Abbasi et al. 2020).The large-scale anisotropy with a strength of ∼ 10 −3 has been qualitatively explained by standard diffusion models describing propagation of CRs through the Galaxy following acceleration by supernova remnants (P.Blasi and E. Amato 2012;G. Giacinti and G. Sigl 2012;P. Mertsch and S. Funk 2015).
Small-scale anisotropy with angular width less than 60 • was first reported by the Milagro experiment showing excesses in two regions, namely region A and B with an excess at the level of 6.0×10 −4 and 4.0×10 −4 respectively (A. A. Abdo et al. 2008).Located at 36 • N latitude, the Milagro could observe the upper half of region A centered at α ≈ 69.4 • , δ ≈ 13.8 • , ranging in 66 • < α < 76 • and 10 • < δ < 20 • .The region B was reported to be a structure lying within 15 • < δ < 50 • with right ascension within 117 • < α < 141 • .The ARGO-YBJ also observed the same regions along with two additional regions namely '3' and '4' (B.Bartoli et al. 2013).The region 3 stretches in the Northern hemisphere within 234 • ≤ α ≤ 282 • .The region '4' is observed with 200 • ≤ α ≤ 216 • and 24 • ≤ δ ≤ 34 • and is similar to the the region 'C' observed by HAWC (A.U. Abeysekara et al. 2014).IceCube has also observed smallscale anisotropic structures in the Southern hemisphere by subtracting dipole and quadrapolar terms from the large-scale maps, and the region B appears to have a continuity in the Southern hemisphere (M.G. Aartsen et al. 2016).A full sky analysis of HAWC and IceCube combined also shows the two primary anisotropy structures in region 'A' and 'B' in which the region 'B' can be seen to extend throughout the entire declination range (A.U. Abeysekara et al. 2019).GRAPES-3 covers a declination range from −23.8 • ≤ δ ≤ 46.6 • with a 56% sky coverage and an overlapping field of view with the above discussed experiments.This paper reports the small-scale anisotropy structures and their features observed by the GRAPES-3 experiment.

THE GRAPES-3 EXPERIMENT AND DATA SET
The GRAPES-3 is a near equatorial and mid-altitude EAS experiment, located in Ooty, India (11.4 • N, 76.7 • E and 2200 m a.s.l.), designed for cosmic ray and gamma ray observations at TeV-PeV energies.It consists of an array of 400 plastic scintillator detectors spread over 25,000 m 2 area (S.K. Gupta et al. 2005; P. K. Mohanty et al. 2009) and a tracking muon detector of 560 m 2 area (Y.Hayashi et al. 2005) as shown schematically in Figure 1.The scintillator detectors are of 1 m 2 area each, placed in a hexagonal configuration with 8 m interdetector separation.The density of secondary particles and their relative arrival times in each EAS are recorded by the scintillator detectors.The EAS parameters such as core location, size and age were obtained by fitting Nishimura-Kamata-Greisen (NKG) function to the lateral distribution of particle densities (K.Kamata and J. Nishimura 1958;K. Greisen 1960).The direction of individual EAS in terms of zenith and azimuth angle (θ, ϕ) were obtained by a fitting a plane front to the relative arrival time data recorded on each scintillator detector after correcting the curvature in the shower front which is observed to be dependent on both size and age parameter (V.B.Jhansi et al. 2020;D. Pattanaik et al. 2022).About 20% of the EAS events mainly of low energies could not be reconstructed properly by the NKG function.Hence, the curvature corrections could not be performed for these events as the size and age parameters are not reliable.However, the θ and ϕ for these events were obtained only using a plane fit.This is not expected to impact the study of anisotropy since the angular scale of anisotropies relevant to this analysis are much larger than the improvement in resolution due to curvature corrections.For this analysis, we used 4 years of EAS data recorded from 1 January 2013 to 31 December 2016 which comprised 3.9 × 10 9 EAS events with a recorded live time of 1273.1 days.The EASs with zenith angle less than 60 • were selected.A total of 3.7 × 10 9 EASs which passed the selection criteria were used for the analysis.The median energy of the EAS was calculated to be 16.2 TeV based on energy-size relation using CORSIKA simulation as described in (D.Pattanaik et al. 2022).

ANALYSIS AND RESULTS
In this analysis, the reconstructed local direction (θ, ϕ) of each EAS event was converted into equatorial coordinate in terms of α and δ.The events were binned into a map in the celestial sphere using HEALPix package (K.M. Gorski et al. 2005).We set the parameter nside=64 which divides the sky into 45192 pixels of equal area of about (0.92 • ) 2 in α and δ.The EAS trigger rate exhibits variations caused by atmospheric pressure and temperature on a daily time scale at the level of 3-5% (M.Zuberi et al. 2017).Further, there are gaps in the data produced either due to scheduled maintenance activity or failure of electronics or DAQ system.Fluctuation in the rates due to abnormality in the detectors or electronics is another effect.It is difficult to identify and remove all these effects completely from the data.They can cause non-uniform exposure of events in the Celestial sphere, thereby producing spurious anisotropies which could be more than a order of magnitude larger than the genuine anisotropy of astrophysical origin that needs to be investigated.Therefore, the search for anisotropies in this map cannot be performed directly.Our search strategy is based on the estimation of a reference map which retains the anisotropies originated from instrumental or atmospheric effects while being insensitive to the astrophysical anisotropy.The reference map was estimated from the data itself using a method called time-scrambling as described in (D.E. Alexandreas et al. 1993).This method has been employed by several air shower experiments for successful extraction of anisotropy (B.Bartoli et al. 2013; M. G. Aartsen et al. 2016).In this method, the true recorded time of an EAS event is forgotten while assigning another time to it which is randomly picked up from a sample of recorded events within a time window of ∆t.The (θ, ϕ) for the event is kept same.Thus it changes the α of the event while keeping δ unchanged.Thus, the anisotropic structures of astrophysical origin smaller than an angular scale of ∆t × 15 • /1 hr are removed.Each recorded event is assigned with 20 random time values from the recorded time sample within the scrambling time window, and then it is filled into the reference map (N ref ′

i
).This essentially means that the reference map has 20 times more number of events than the data map.This is done to reduce the statistical fluctuation in the reference map.An average of all these 20 maps is then calculated to estimate the reference map as, /20.To re-state it, both the reference map and data map will contain the anisotropies produced by atmospheric and instrumental effects whereas the reference map will not have the anisotropies of astrophysical origin.The relative intensity for each pixel is calculated using, δI =  First we apply the time-scrambling method on data using a scrambling window of ∆t = 24 hrs.The stability in zenith and azimuth over long periods of time allows such a choice.Since in 24 hrs, there will be full coverage of the celestial sphere, anisotropies up to the largest angular scale present in the data should show up as they will be destroyed in the reference map.The relative intensity map does not show any discernible anisotropy due to high level of statistical fluctuation in the pixels.To improve the capability of identifying directional features, we employ a smoothing method similar to the method followed by (A.U. Abeysekara et al. 2014).This approach includes combining the event counts within individual pixels and adding counts from neighboring pixels positioned within a chosen smoothing radius of 10 • in this particular study.This approach serves to decrease statistical fluctuations in nearby pixels, thereby revealing any localized surpluses.However, it also establishes correlations in event statistics among neighboring pixels.Essentially, this is analogous to applying a smoothing operation to the map using a top-hat function of 10 • radius.The structures that we are trying to observe span over few tens of degrees, so a choice of a 10 • smoothing radius has ensured that these structures are not integrated over.The relative intensity map after smoothing is shown in Figure 2. The pixel-wise significance was The variation of events as a function of right ascension, obtained from the unsmoothed maps by summing events over declination bands spanning within −18.6 • to 41.4 • , is shown in Figure 4.The range is taken within ± 30 • from the latitude of the GRAPES-3 experiment.The plot clearly illustrates that the variation in the event count, resulting from non-uniform sky exposure due to detector and atmospheric effects, is mir-rored in the background estimation obtained through the time-scrambling algorithm.The systematic excesses for regions A and B are also seen in the relative intensity plot.A zoomed view of the two regions are shown in Figure 5 and Figure 6.The details are discussed below.
Region A is observed in the right ascension range of ∼50 • to 80 • and declination range of ∼-15 • to 30 • .The maximum relative intensity of the anisotropy in this region is (8.9 ± 2.1 ± 0.2) × 10 −4 at the pixel centered at (α = 63.9 • , δ = −7.2• ), where the first error is statistical and the second error is systematic.The calculation of the systematic error involves conducting an analysis using "anti-sidereal" time (K.Nagashima et al. 1998;A. A. Abdo et al. 2008).While sidereal time corresponds to the sky fixed frame, anti-sidereal time represents a nonphysical frame.Analyzing data based on sidereal time allows us to determine the presence of anisotropy, while analyzing data based on anti-sidereal time provides insight into the systematic error on sidereal anisotropy resulting from non-physical effects (R. Abbasi et al. 2010;F. J. M. Farley and J. R. Storey 1954).The same timescrambling analysis with a time window of 4 hrs was performed but with anti-sidereal time in order to estimate the systematics and the maximum strength observed in regions A and B are quoted as systematic errors.The results of anti-sidereal time have been shown in Figure 7 and no characteristic signal regions can be seen in this case implying that the systematics caused by non-physical effects are insignificant.The statistical error dominates as the systematic error is much lesser than statistical error.
Region A appears to be a circular structure with a tail like projection.The maximum pixel significance observed in this region is 5.8σ.To calculate the total significance of this region, the unsmoothed data and reference maps were used in order to avoid the correlations introduced by smoothing between the pixels.First, the structure was defined by selecting those pixels which have a significance of more than 2σ as shown in Figure 5 in bottom.The total number of events in data and reference maps from the region was obtained by summing up the pixel wise events in the selected area.The Li-Ma prescription was used to obtain the total significance which is 6.8σ.The relative excess number of events in this region is (6.5 ± 1.3) × 10 −4 .
Region B is an elongated structure observed within ∼110 • to 140 • of right ascension and almost throughout the full declination range (see Figure 6).The maximum relative intensity observed for this region is (5.6 ± 1.8 ± 0.1) × 10 −4 at the pixel centered at (α = 124.5 • , δ = 3.4 • ) and the significance of the pixel is 4.4σ.Similar to the criteria applied for region A, those pixels which have a significance of more than 2σ were selected to define region B. The overall relative intensity and significance of the region B is (4.9±1.4)×10−4 and 4.7σ, respectively.
The deficit seen around regions A and B are also consistent with the observations by Milagro, ARGO-YBJ and HAWC.The deficit observed between regions A and B is the most significant.The deficit structure has a significance of 3.7σ and a relative intensity of −(4.6 ± 1.7) × 10 −4 .

DISCUSSION
By analyzing 3.7×10 9 EAS events collected over a period of 4 years, GRAPES-3 could significantly observe two excess regions namely A and B. The region A shows a tail like structure in the Northern hemisphere (δ > 0 • ).The shape of the structure is similar to the "region A" observed by Milagro (A. A. Abdo et al. 2008) and HAWC (A.U. Abeysekara et al. 2014), and "region 1" reported by ARGO-YBJ (B.Bartoli et al. 2013).Milagro (at 36 • N) observes the part of this structure in the Northern sky and the observations are continued to the Southern hemisphere by ARGO-YBJ (at 30 • N), HAWC (at 19 • N) and GRAPES-3 (at 11.4 • N).GRAPES-3 and HAWC lying closer to the Equator have the advantage of covering the southern part of region A. Region B is also observed by Milagro, ARGO-YBJ (referred to as "region 2") and HAWC as a continuous structure running almost throughout the entire declination band, similar to observations by GRAPES-3.The full sky analysis by HAWC and IceCube show that region B continues to the Southern hemisphere as well (A.U. Abeysekara et al. 2019) running through the entire declination band.The deficit regions seen around these excesses by GRAPES-3 are also coincident with observations by Milagro, ARGO-YBJ and HAWC.
ARGO-YBJ (B.Bartoli et al. 2013) and HAWC (A.U. Abeysekara et al. 2014) have also performed an analysis based on energy dependence by partitioning their datasets into multiple segments, some of which overlap with the median energy range of GRAPES-3 at about 16 TeV.The relative intensity of region A observed by GRAPES-3 is (6.5 ± 1.3) × 10 −4 .ARGO-YBJ divided their dataset into five segments, with the last two seg-   ergy.When considering HAWC's analysis, their dataset was divided into seven segments, and the relative intensity of region A was assessed for each segment around the peak.It was reported by HAWC that the relative intensity of region A increases with energy.When comparing our findings, we use the segment of data having a median energy of 14.7 +28.7 −9.9 TeV, which is the closest approximation to GRAPES-3's median energy and a relative intensity of ∼ (22.0±5.0)×10−4 was reported.For region B, the relative intensity of the observed structure is (4.9 ± 1.4) × 10 −4 by GRAPES-3.The average relative intensity of region B measured by ARGO-YBJ is 3.5 × 10 −4 .It varies in the range of (3.0 − 5.0) × 10 −4 for the last two energy bins that cover GRAPES-3's median energy.
There are different explanations for the origin of the small-scale anisotropic structures.Several models propose that the origin of the hotspots might be linked to the proximity of supernova explosions, events known for generating cosmic rays (A.D. Erlykin and A.W. Wolfendale (2006)).According to (M.Salvati and B. Sacco (2008)), regions A and B could result from a phenomenon associated with a supernova explosion that gave rise to the Geminga pulsar.Situated between the observed structures of regions A and B is the Geminga pulsar (α, δ=98.47 • , 17.77 • ).Some models suggest that these structures are a consequence of the turbulent magnetic field within the cosmic ray scattering length (G.Giacinti and G. Sigl 2012;M. Ahlers and P. Mertsch 2015) or CR scattering by Alfven waves created by turbulent cascades in local field direction (Malkov et al. 2010).Another set of models suggest that they are generated due to magnetic reconnections in the heilosphere (P.Desiati and A. Lazarian 2012).The model presented in L. Drury and F. Aharonian (2008) investigates the hypothesis that region A could potentially be explained by the production of secondary neutrons in the concentrated tail of interstellar material that forms downstream of the Sun's movement through the local interstellar medium (ISM).Some exotic models also suggest that decay of quark matter in pulsars or the acceleration of strangelets near molecular clouds are causes for the anisotropy (K.Kotera et al. 2013;M. Ángeles Pérez-Garcí et al. 2014).Thus, the examination of CR anisotropy assumes importance in comprehending the magnetic field structure, propagation, acceleration mechanisms of CRs, and their sources.

SUMMARY
We have investigated small-scale anisotropy in CR distribution by analyzing 3.7×10 9 events at median energy of ∼ 16 TeV, collected by the GRAPES-3 experiment from 1 January 2013 to 31 December 2016 with a live time of 1273.1 days.Time scrambling method was employed to estimate the background map.Two smallscale structures namely regions A and B have shown excesses in the CR flux at the level of (6.5 ± 1.3) × 10 −4 and (4.9±1.4)×10−4 , respectively with statistical significance of 6.8σ and 4.7σ, respectively.These structures are in agreement with previous results reported by Milagro, ARGO-YBJ, and HAWC.

Figure 1 .
Figure 1.Schematic diagram of the GRAPES-3 experiment with the scintillator detectors represented by (■) and the muon detector modules represented by (□).
of events in the i-th pixel of data and reference map respectively.

Figure 2 .
Figure 2. Anisotropy (top) and significance (bottom) obtained after time-scrambling with ∆t = 24 hrs and smoothing radius of 10 • .The large-scale deficit can be seen while small scale structures A and B are prominently seen.

Figure 3 .
Figure 3. Anisotropy (top) and significance (bottom) observed with a scrambling window of 4 hrs and a smoothing radius of 10 •

Figure 4 .
Figure 4. 1-d projection of anisotropy showing the variation in number of events with right ascension

Figure 5 .
Figure 5. Anisotropy (top) and significance (bottom) of Region A as observed with a scrambling window of 4 hrs.The pixels with higher than 2σ significance have been shown in the bottom plot and the region defined has been used to calculate the significance of the entire structure.

Figure 6 .
Figure 6.Anisotropy (top) and significance (bottom) of Region B as observed with a scrambling window of 4 hrs, and plotted similar to region A

Figure 7 .
Figure 7. Anisotropy (top) and significance (bottom) obtained after performing the same analysis with anti-sidereal time.No characteristic signal region can be observed.

Table 1 .
The peak intensities of regions A and B as reported by ARGO-YBJ, HAWC and GRAPES-3