Northern sky Galactic Cosmic Ray anisotropy between 10-1000 TeV with the Tibet Air Shower Array

TL;DR

This study analyzes Galactic cosmic ray anisotropy between 10 and 1000 TeV using Tibet Air Shower Array data, revealing energy-dependent features and a new anisotropy component at high energies, contributing to cosmic ray origin understanding.

Contribution

The paper extends anisotropy analysis to higher energies and declination range, identifying a new anisotropy component above 100 TeV and combining northern and southern sky data for full-sky mapping.

Findings

Anisotropy maps differ above 100 TeV from lower energies.

A new anisotropy component appears around 100 TeV.

Amplitude of first-order anisotropy increases sharply above 100 TeV.

Abstract

We report the analysis of the $10-1000$ TeV large-scale sidereal anisotropy of Galactic cosmic rays (GCRs) with the data collected by the Tibet Air Shower Array from October, 1995 to February, 2010. In this analysis, we improve the energy estimate and extend the declination range down to $-30^{\circ}$. We find that the anisotropy maps above 100 TeV are distinct from that at multi-TeV band. The so-called "tail-in" and "loss-cone" features identified at low energies get less significant and a new component appears at $\sim100$ TeV. The spatial distribution of the GCR intensity with an excess (7.2$\sigma$ pre-trial, 5.2$\sigma$ post-trial) and a deficit ($-5.8\sigma$ pre-trial) are observed in the 300 TeV anisotropy map, in a good agreement with IceCube's results at 400 TeV. Combining the Tibet results in the northern sky with IceCube's results in the southern sky, we establish a full-sky picture of the anisotropy in hundreds of TeV band. We further find that the amplitude of the first order anisotropy increases sharply above $\sim100$ TeV, indicating a new component of the anisotropy. All these results may shed new light on understanding the origin and propagation of GCRs.

Tables

Table 1. Table 1 Fitting results of the first harmonic (Amplitude, Phase, and Reduced χ 2 superscript 𝜒 2 \chi^{2} ) in the sidereal (Columns 2-4), solar (Columns 5-7), and Antisidereal (Columns 8-10) Times. The number of events in each energy sample is given in Column 11.

Energy	$A_{sid}$	${\varphi }_{sid}$	${\chi }_{sid}^2/ndf$	$A_{sol}$	${\varphi }_{sol}$	${\chi }_{sol}^2/ndf$	$A_{asid}$	${\varphi }_{asid}$	${\chi }_{asid}^2/ndf$	Number of Event
TeV	${10}^{-4}$	[∘]		${10}^{-4}$	hr		${10}^{-4}$	hr
15	8.5$\pm$0.2	21.9$\pm$1.6	911./16	4.1$\pm$0.2	6.05$\pm$0.22	69.7/16.	0.53$\pm$0.24	22.2$\pm$1.7	24.4/16	$2.33\times {10}^{10}$
50	5.3$\pm$0.4	20.8$\pm$4.7	152.9/16	4.6$\pm$0.4	6.37$\pm$0.35	46.7/16	0.39$\pm$0.43	22.5$\pm$4.2	46.7/16	$3.97\times {10}^9$
100	2.7$\pm$0.6	326.8$\pm$12.0	67.6/16	4.0$\pm$0.6	5.91$\pm$0.53	14.2/16.	0.80$\pm$0.56	22.3$\pm$2.7	9.8/16	$1.96\times {10}^9$
300	6.0$\pm$1.4	267.1$\pm$13.5	27.0/16.	3.5$\pm$1.4	6.53$\pm$1.56	19.0/16.	1.7$\pm$1.4	5.3$\pm$3.2	11.6/16	$2.71\times {10}^8$
1000	13.0 $\pm$3.0	286.6$\pm$12.6	9.3/16	9.8$\pm$2.8	6.97$\pm$1.13	10.9/16.	1.0$\pm$2.9	14.8$\pm$11.0	6.3/16	$5.72\times {10}^7$

Equations

$$R(\alpha )=1+A_1\cos (\alpha -{\varphi }_1),$$

Full text

1. Introduction

The arrival directions of Galactic cosmic rays (GCRs) are nearly isotropic due to deflections in the Galactic magnetic field (GMF).Only weak anisotropy is expected from the diffusion and/or drift of GCRs in GMF. Observations of ground-based air shower arrays and underground muon detectors do show the existence of small anisotropies with relative amplitudes of the order of $10^{-4}$ to $10^{-3}$ at energies from 100 GeV to hundreds of TeV (see Figure 5). However, the variation of the amplitude with energy seems to be difficult to interpret in terms of the conventional GCR diffusion model in the Galaxy ($e.g.$ (Moskalenko et al., 2002; Ahlers & Mertsch, 2016)). The study of GCR anisotropy, therefore, is important to understand the origin and propagation of GCRs.

Only a few results of the anisotropy in the energy range from hundreds of TeV up to $\sim 10$ PeV have been reported, primarily due to the low fluxes of cosmic rays (CRs) in this energy range. EAS-TOP collaboration reported for the first time a detection of anisotropy at $\sim 200$ TeV (Aglietta et al., 1996). With the accumulation of data, they improved their result later and reported a sharp increase of the anisotropy amplitude at primary energies around $\sim$370 TeV (Aglietta et al., 2009). At the PeV energy region, the Akeno experiment reported an increase of the CR anisotropy amplitude in 1986 (Kifune et al., 1986). No hint of the anisotropy, on the other hand, has been found in the KASCADE data at higher energies between 0.7 and 6 PeV (Antoni et al., 2004). Recently, the IceCube collaboration reported the anisotropy observed in the southern sky, showing a new feature different from that obtained by EAS-TOP (Abbasi et al., 2012). A clear deficit with a post-trial significance of $-6.3\sigma$ at 400 TeV was detected, which was then confirmed by the result from Ice-Top (Aartsen et al., 2013). The Ice-Top data further revealed the existence of anisotropy at energies up to 1 PeV (Aartsen et al., 2013).

The Tibet Air Shower (AS) array collaboration presented the first two-dimensional anisotropy measurements in an energy region from several TeV to several hundred TeV. The anisotropy features, known as the “tail-in” and “loss-cone” features, were observed with very high significances (Amenomori et al., 2006). A new component anisotropy at multi-TeV energies from the Cygnus direction was also reported (Amenomori et al., 2006). It has been shown that the amplitude of the first order anisotropy decreases above a few hundred TeV, indicating the co-rotation of GCRs around the Galactic center. With more data accumulated, hints of $\sim 300$ TeV anisotropies have been revealed (Feng et al., 2009; Amenomori et al., 2013). The anisotropy feature was found to be different from those in lower energy regions and in agreement with IceCube’s result at 400 TeV (Abbasi et al., 2012). These analyses of Tibet AS array data cover decl. from $-15^{\circ}$ to $75^{\circ}$, yet leaving a gap to be connected with IceCube’s result in the southern sky. Here we extend these analyses to include events with zenith angle up to $60^{\circ}$, which corresponds to a coverage of decl. from $-30^{\circ}$ to $90^{\circ}$ (Amenomori et al., 2015). Combining with IceCube’s results, we present for the first time a full-sky anisotropy observed at hundreds of TeV. By improving the reconstruction of primary energy, we will also extend the analyzed energy range to two decades between 10 TeV and 1 PeV, which is also the widest coverage of such works.

2. Analysis

2.1. Experiment and Data reconstruction

The Tibet AS Array is located at Yangbajing in Tibet, China ($90.522^{\circ}$E, $30.102^{\circ}$N, 4300 m above sea level, 606 g/cm2 atmospheric depth). The detector array consists of plastic scintillation detectors with an area of 0.5 m2 each. The effective area of the Tibet AS array has been gradually enlarged, via adding the same-type detectors to the array. The Tibet I array was constructed in 1990, using 65 plastic scintillation detectors placed on grids with 15 m spacing. It was then upgraded to 221 detectors on 15 m grids, covering a total of 36,900 m2, known as the Tibet II array. It began operation in 1995 October, with a trigger rate of $\sim 230$ Hz. The Tibet II was then upgraded to the current Tibet III, a denser array with 7.5 m grids, in 1999 and 2003 (Amenomori et al., 2003). The trigger rate is $\sim$1700 Hz for the Tibet III array.

In order to maintain the uniformity of the array performance, we analyze the data keeping the same configuration of the Tibet II array throughout the observation period from 1995 October to 2010 February, so that the full data sample taken by Tibet II and Tibet III array can be used in the present analysis. The traditional shower reconstruction procedure is applied to get all the parameters of one shower, such as the core position, zenith, and azimuth angles ($\theta$, $\phi$) of the incident direction and shower size $\sum\rho_{FT}$ (the sum of the number of particles per m2 counted by all the fast-timing [FT] detector). The following three criteria are applied to select events for further analyses: (1) each AS event should fire four or more detectors, with each recording 1.25 or more particles; (2) the AS core position should be located inside the array; and (3) zenith angle $\theta<60^{\circ}$.

2.2. Estimation of the CR Energy

The ASs reaching the array with a large zenith angle $\theta$ travel through a larger slant atmospheric depth than the vertical ones. This leads to a zenith angle dependence of the relation between $\sum\rho_{FT}$ and the primary particle energy. In most of the previous works of the Tibet AS$\gamma$ Collaboration (Amenomori et al., 2003, 2005a, 2006; Feng et al., 2009; Amenomori et al., 2013), the shower size $\sum\rho_{FT}$ is solely adopted to infer the primary energy of an AS without considering the zenith angle dependence. This approximation works well for small zenith angles ($\theta<\sim 40^{\circ}$), considering the natural fluctuation in the development of the AS and the limited resolution of the primary energy.

In this work, we intend to explore the anisotropy with decl. down to $-30^{\circ}$ by including showers with zenith angles up to $60^{\circ}$. Therefore, the zenith angle dependence of the energy reconstruction must be taken into account. We develop a two-dimensional selection criterion in the $\sum\rho_{FT}-sec\theta$ plane for the energy reconstruction.

The uncertainty of the CR energy reconstruction has been estimated with a full Monte Carlo (MC) simulation of CR interactions in the atmosphere by CORSIKA (version 6.204; (Heck et al., 1998). The hadronic interaction model QGSJET01c and the detector response modeled by Epics (version 8.65; (Kasahara, 2006) are used. In this simulation, we adopt the composition and spectrum models of primary CRs given in (Hörandel, 2003) as inputs.

Figure 1 shows the simulated distribution in the primary CR energy as a function of $\sum\rho_{FT}$ and sec $\theta$. It shows that for a given range of $\sum\rho_{FT}$, small zenith angle events (at sec $\theta\sim 1$) are dominated by CRs with lower average energy, compared with large zenith angle events (at sec $\theta\sim 2$). We show regions of constant primary energies (15, 50, 100, 300, and 1000 TeV) in the plane of ($\sum\rho_{FT}$, sec $\theta$) as regions delimited by the dashed lines in Figure 1. This grouping enables us to select events in five energy samples with minimal overlappings.

Figure 2 shows the simulated primary energy distributions of the five energy samples, as indicated by the dashed lines in Figure 1. The uncertainty of the primary energy estimate is dominated by the fluctuation of the AS. Event numbers in five energy samples are $2.33\times 10^{10}$ (15 TeV), $3.97\times 10^{9}$ (50 TeV), $1.96\times 10^{9}$ (100 TeV), $2.71\times 10^{8}$ (300 TeV), and $5.72\times 10^{7}$ (1 PeV), as listed in Table 1.

2.3. Analysis of the first harmonics of the anisotropy

We employ the All-Distance Equi-Zenith Angle Method (Amenomori et al., 2005a, 2006), which has been shown to be sensitive to probe the large-scale anisotropy, to analyze the data. Details about this method can be found in (Amenomori et al., 2005a). The sky of the horizontal coordinate is divided into cells with a bin size of $1^{\circ}$ in both zenith (from $0^{\circ}$ to $60^{\circ}$) and azimuth (from $0^{\circ}$ to $360^{\circ}$). For the equatorial coordinate, the sky is divided into cells of $2^{\circ}\times 2^{\circ}$ between $0^{\circ}$ and $360^{\circ}$ in the R.A. ($\alpha$) and between $-30^{\circ}$ and $90^{\circ}$ in the decl. ($\delta$). The two-dimensional (2D) map in the equatorial coordinate is then smoothed in a window, changing the window widths between 5 and 30*∘*.

To quantify the magnitude of the anisotropy, we project the two-dimensional (2D) anisotropy map before smoothing onto the R.A. axis, through averaging the relative intensities in all declinations from $-30^{\circ}$ to $90^{\circ}$, to derive the one-dimensional (1D) profile of the anisotropy. The R.A. is binned into 18 bins for this 1D analysis, and the 1D profile of the anisotropy is fitted by the first order harmonic function in form of

[TABLE]

where $R(\alpha)$ denotes the relative intensity of CRs at R.A. $\alpha$, $A_{1}$ is the amplitude of the first harmonics, and $\phi_{1}$ is the phase at which $R(\alpha)$ reaches its maximum.

3. Results

3.1. Sidereal Anisotropy map at 300 TeV

Figure 3 shows the significance map and the relative intensity map of the sidereal anisotropy for the 300 TeV energy sample. A smoothing with an optimized window width of $30^{\circ}$ is applied in this figure. We combined the 300 TeV and 1 PeV samples together in this figure to increase the statistics. The total event number used in this figure is $3.28\times 10^{8}$, and the median energy is approximately 300 TeV.

From the significance map, we find that two regions are significant-that is, an excess centered at ($\alpha=263^{\circ}$, $\delta=11^{\circ}$) with a significance of $7.2\sigma$ and a deficit centered at ($\alpha=93^{\circ}$, $\delta=-25^{\circ}$) with a significance of $-5.8\sigma$. Note that the significance values are the pre-trial results. We conservatively estimate a trial factor by assuming that all scans give statistically independent results. Since the search for this excess is performed over about 60$\times$180 cells, and 26 different smoothing radii, the total trial factor is estimated to be about $2.81\times 10^{5}$. The post-trial significance of the excess is $\sim 5.2\sigma$. The deficit is no longer significant, fail to reach the 5$\sigma$ level, after the correction for the trials.

Because the acceptance of the detector decreases for larger zenith angles, the relative intensity map is similar but not completely the same as the significance map. An excess region centered at ($\alpha=269^{\circ}$, $\delta=-13^{\circ}$) with a maximum excess of +1.38$\times 10^{-3}$, and a deficit region centered at ($\alpha=87^{\circ}$, $\delta=-29^{\circ}$) with a maximum deficit of -1.80$\times 10^{-3}$ can be seen. Both the excess and deficit regions are consistent with the results of IceCube at 400 TeV in the southern hemisphere (Abbasi et al., 2012). Combining these results gives us a full-sky picture of the sidereal anisotropy of GCRs at hundreds of TeV.

The bottom panel in Figure 3 also shows the 1D projection of the relative intensity before the smoothing onto the R.A. axis. The correlation among different bins is carefully considered when calculating the statistical errors and fitting the data with the harmonic function in equation 1. If the correlation is not considered correctly, the errors of the fitting parameters would be underestimated. The blue curve shows the best-fitting result, with the fitting parameters indicated in the figure. The significance of non-zero amplitude is $5.6\sigma$, which shows that the obtained first harmonics are indeed significant. The reduced $\chi^{2}$ value is 26.7/16, which means that the first harmonic function can describe the 1D projected profile well.

One of the possible origins of the sidereal anisotropy is the Compton-Getting (CG) effect, due to the orbital motion of the solar system around the Galactic center (Compton & Getting, 1935). The relative intensity of this effect is carefully calculated by the MC method. Considering the location of the Tibet AS array ($90.522^{\circ}$E, $30.102^{\circ}$N), the velocity (220 km s*-1*) of the orbital motion of the solar system around the Galactic center and the spectrum index (2.7) of the CRs energy spectrum, the intensity of the sky that is divided into cells of $2^{\circ}\times 2^{\circ}$ between $0^{\circ}$ and $360^{\circ}$ in the R.A. ($\alpha$) and between $-30^{\circ}$ and $90^{\circ}$ in the decl. ($\delta$) is calculated. Then the identical analyses are performed to this MC data sample.

The R($\alpha$) expected from this CG effect, shown as the black dashed line in Figure 3, has a maximum at ($\alpha=315^{\circ}$, $\delta=0^{\circ}$) and a minimum at ($\alpha=135^{\circ}$, $\delta=0^{\circ}$). Neither the amplitude nor the phase of the large-scale anisotropy observed in this work can be described in terms of the CG effect.

3.2. Variation of CR sidereal anisotropy with the energy between 10 and 1000 TeV

Figure 4 shows the variation of the sidereal anisotropy with the energy between 10 TeV and 1 PeV. At 15 TeV and 50 TeV, the tail-in and loss-cone features (Amenomori et al., 2006) are observed with very high significances. An intensity excess in the Cygnus region can also be seen. However, these features become less significant above 100 TeV, being replaced with some new features. At 300 TeV and 1 PeV, the anisotropy maps are distinctly different from those in 15-50 TeV. We can clearly see the phase of the 1D projection changing with the primary energies, as seen in Table 1, showing the best-fit parameters.

Figure 5 compares the amplitude and phase obtained in this work with those reported so far from the deep underground muon experiments and extensive AS experiments. Our results are in close agreement with other results in similar energy regions in both the amplitude and the phase. It is interesting to note that a sharp increase of the amplitude above 100 TeV can be seen in the upper panel. The origins of this feature cannot be explained with the conventional diffusion scenario of GCRs, and may provide us with a new hint for understanding the origin and propagation of GCRs.

3.3. Anisotropy in Solar Time and Antisidereal Time

In order to confirm that the obtained anisotropy is not affected by the seasonal variation of the AS array performance, identical analyses are performed in the solar time and antisidereal time frames in five energy samples. Figure 6 shows the local solar time and antisidereal time daily variations measured by Tibet AS Array in five energy samples, and the best-fit parameters are also shown in Table 1. The amplitude and phase in the solar time frame are in good agreement with the expectation from the CG effect due to the terrestrial orbital motion around the Sun ($A_{\rm sol,CG}=0.047\%$ and $\phi_{\rm sol,CG}=6.0$ hr). In all five energy samples, no significant anisotropy is observed in the antisidereal time frame, ensuring that no additional correction is required for the seasonal effects. The observed results in the solar and antisidereal time frames support the reliability of the observed sidereal anisotropy.

4. Conclusion and Discussion

Fifteen years data recorded by the Tibet AS array have been analyzed to study the sidereal anisotropy of CRs. In this work, we improve the estimate of the primary CR energies through a 2D cut in the $\sum\rho_{FT}-$sec$\theta$ plane, to explore the anisotropy including larger zenith angle events. For the first time, we extend the analyzed decl. down to $-30^{\circ}$ to complete a full-sky coverage of the anisotropy at hundreds of TeV energies by combining with the IceCube’s results at the South Pole. The 2D anisotropy map at $\sim 300$ TeV obtained in this work is smoothly connected with IceCube’s results at 400 TeV. The energy dependence of the large-scale sidereal anisotropy has been derived between 10 TeV and 1 PeV. We measured the energy dependence of the first harmonics of the anisotropy above 100 TeV, which may be associated with local origins of GCRs.

The CG effect expected from the orbital motion of the solar system around the Galatic center is not observed at 300 TeV, as shown in Figure 3. The basic picture that GCRs are co-rotating with the local Galactic neighbors still holds at this energy (Amenomori et al., 2006). As pointed out earlier, the GCR rest frame may have a smaller relative velocity with a different direction from neighboring stars and the interstellar medium (Abbasi et al., 2012). This scenario is possibly responsible for the GCR anisotropy observed at hundreds of TeV.

The strongest excesses at hundreds of TeV are from the direction of the Galactic center, which may imply a Galactic center origin of GCRs at these energies (Guo et al., 2013). It is interesting to note that the highest-energy CR accelerators have been identified by the HESS telescope in the Galactic center (HESS Collaboration et al., 2016). However, the energy dependences of the amplitude and phase cannot be easily understood in a simple diffusion scenario with any types of GCR sources.

The sharp increase of the amplitude above 100 TeV may imply an evolution of propagation parameters, such as spatial parameters (Tomassetti, 2015; Guo et al., 2016). The knowledge of the propagation of GCRs needs to be further improved for our full understanding the properties of the anisotropy, especially in this high-energy region where the conventional diffusion/drift models may not work any more. Finally, we add to note that the measurements of the anisotropy above PeV, which is possibly associated with the knee of GCRs, are very important to advance our understanding of origin and propagation of GCRs.

5. Acknowledgments

The collaborative experiment of the Tibet Air Shower Arrays has been performed under the auspices of the Ministry of Science and Technology of China and the Ministry of Foreign Affairs of Japan. This work was supported in part by a Grant-in-Aid for Scientific Research on Priority Areas from the Ministry of Education, Culture, Sports, Science and Technology, by Grants-in-Aid for Science Research from the Japan Society for the Promotion of Science in Japan, and by the Grants from the National Natural Science Foundation of China and the Chinese Academy of Sciences. Zhaoyang Feng is supported by the Natural Sciences Foundation of China (Nos.11405182, Nos.1135010). C. Liu is supported by the Natural Sciences Foundation of China (Nos. 11405180).