Modeling of TeV Galactic Cosmic-ray Anisotropy based on Intensity Mapping in an MHD Model Heliosphere

Cosmic rays do not arrive at the Earth uniformly. Recent experiments have consistently observed small anisotropies with amplitudes of ∼0.2% in the arrival directions of cosmic rays at TeV energies. We perform the modeling of the cosmic-ray anisotropy at TeV energies using the intensity-mapping method based on Liouville’s theorem. This work indicates unrealistically small-scale ansiotropies with <∼10° angular scales in the distribution of the cosmic-ray relative intensity at the outer boundary of the heliosphere. We would possibly need to resolve an issue that the experimental data covers ten years while the Magnetohydrodynamics (MHD) model heliosphere used in this work is only a single snapshot at a certain moment. Performing the intensity-mapping for multiple snapshots of the MHD model heliosphere and taking the average of the results would improve the results of this work.


Introduction
Galactic cosmic rays (CRs) do not arrive at the Earth uniformly.Anisotropies with amplitudes of roughly 0.2% have been observed at tera-electronvolt (TeV) energies by recent cosmic-ray experiments (e.g.[1,2,3,4,5,6]), such as a large-scale deficit region named 'Loss-Cone' and a large-scale excess region named 'Tail-In'.The origins of these anisotropic features have yet to be revealed, although the anisotropy reflects the propagation of CRs in the magnetic fields of the heliosphere and the surrounding interstellar medium.
Previous papers [7,8] utilize the 'intensity-mapping' method, in which a model heliosphere is reconstructed by Magnetohydrodynamics (MHD) simulations, CR trajectories in the magnetic fields of the MHD model heliosphere are calculated, and Liouville's theorem is employed to map the distribution of the CR intensity at the Earth to the outer boundary of the heliosphere.Then, the obtained anisotropies of the CR intensity at the outer boundary is interpreted in physical terms by modeling the intensity distribution based on a superposition of CR flows.It was indicated in the previous paper [7] that, dominant in the interstellar medium outside the heliosphere is a dipole flow of CRs along the interstellar magnetic field (B ISM ), and that a density gradient of CRs also exists in the direction of the Vela supernova remnant.In order to make the modeling conclusive, however, two technical issues must be resolved.Firstly, the previous work [7] employed protons with monochromatic energy of 4 TeV for the intensity-mapping, although the energy change of those protons due to the acceleration and deceleration in the MHD model heliosphere were taken into account.Since CRs arriving at the Earth are composed of different atomic nuclei with different energies, the energy spectrum and the composition of observed CRs must be taken into account in the intensity-mapping process.Secondly, the χ 2 /d.o.f. of the fitting between the experimental data and the best-fit model anisotropy was 4.5 in the previous paper, which was not good enough.The χ 2 /d.o.f.value must be reduced to ∼ 1 by improving the modeling of the distribution of the CR intensity at the heliospheric outer boundary.

Intensity Mapping
In this work, we use the data obtained from November 1999 to May 2010 by the Tibet ASγ experiment.Figure 1 shows the obtained CR relative intensity distribution in equatorial coordinates.We pixelize the sky in our field of view (−20 • < decl.< 80 • ) in 2056 pixels, using the HEALPix algorithm [9] with N side = 16.Each pixel has an approximate size of 3.7 • × 3.7 • .We estimate the rigidity distribution of observed CRs using detailed Monte Carlo (MC) simulations of air-shower generation and detector response.We generate air showers using CORSIKA v7.4000 [11] with EPOS LHC [12] for the high-energy hadronic interaction model and FLUKA v2011.2b[13,14] for the low-energy hadronic interaction model, in the energy range from 0.3 TeV to 10 PeV based on a model of the CR energy spectrum and chemical composition derived from direct measurements [10].We feed the generated air showers into the detector response simulation developed by GEANT v4.10.00 [15] and analyze them in the same way as the experimental data.The obtained rigidity distribution of CRs is shown in Figure 2 for three declination bands.
We use the fourth-order Runge-Kutta method to calculate the trajectories of CRs in the MHD model heliosphere.To smooth out possible seasonal effects, we place the Earth at four positions on the ecliptic plane at a distance of 1 AU from the Sun.Using the HEALPix algorithm [9] we pixelize the sky with N side = 32 in the declination range from −20 • to 80 • , and shoot CR anti-partiles from each pixel center into the MHD model heliosphere, sampling their rigidity  from Figure 2. We use the same MHD model heliosphere as the one employed in [7].We track the CRs in the MHD model heliosphere and record their momentum directions at the outer boundary, and then the CR intensity at the outer boundary in the CR momentum direction (I ISM ) is regarded as equal to that at the Earth (I E ).We evaluate the distribution of the relative intensity of CRs at three boundaries at distances of 630 AU, 1580 AU and 3980 AU from the Sun, with the nose direction truncated at the surface where the direction (strength) of the magnetic field becomes different from that outside the heliosphere by < 0.1 • (0.1%).
The 'declination bias' exists in the experimental data --the average of the CR intensity has been normalized to unity in each declination band since the detection efficiency of the experiment has not been absolutely calibrated along declination.The following steps, therefore, need to be taken to estimate the I ISM distribution at the outer boundary: 1) construct I ISM , a model of the CR intensity distribution at the outer boundary, 2) map I ISM to that at the Earth I E based on the CR trajectory calculation in the MHD model heliosphere, 3) normalize the average of I E to unity in each declination band, and 4) obtain χ 2 between the experimental data and the normalized I E .We repeat these four steps until the χ 2 is minimized, and obtain the best-fit model of I ISM .We express the distribution of the CR relative intensity at the outer boundary in spherical harmonics Y lm at step 1): where f lm 's are free fitting parameters.We increase the maximum order l max of the spherical harmonics until we get the χ 2 /d.o.f.close to unity.

Results and Discussions
The results are summarized in Figure 3. Panels (a), (b) and (c) are the best-fit model distributions of the CR relative intensity at the outer boundaries at distances of 630 AU, 1580 AU and 3980 AU from the Sun, respectively.The maximum order l max of the spherical harmonics Y lm in Eq.( 1) is l max = 4 for Panel (a), l max = 8 for Panel (b), and l max = 20 for Panel (c), and the obtained χ 2 /d.o.f.values are 0.96 (probability 89%), 0.98 (71%) and 0.94 (95%), respectively.From Figure 3, one can notice that l max decreases as the distance of the boundary from the Sun decreases.In terms of physics, it should be a natural assumption that only large-scale anisotropies exist in the distribution of the CR intensity at the helisopheric outer boundary, and that small-scale anisotropies are added to the intensity distribution by the modulation of CR trajectories in the heliospheric magnetic field.The tendency seen in panels (a) − (c) in Figure 3, however, disagrees with this assumption.Considering that the MHD model heliosphere employed in this work has been refined already to give reasonable agreement with experimental results (for details see [7]), we consider that our intensity mapping needs to be improved further.One possible issue would be that, the CR scattering with magnetic irregularities in the heliosphere is not taken into account in the intensity mapping process when we calculate CR trajectories in a single snapshot of the magnetic field structure of the MHD model heliosphere, while the experimental data used in this work, covering a period of ∼ten years, contains this CR scattering effect.This discrepancy would be resolved by performing the intensity-mapping for multiple snapshots of the MHD model heliosphere and taking the average of the results, which would lead to remove the unnatural tendency in Figure 3 and enable us to derive CR flows in the interstellar medium outside the heiosphere.

Figure 1 .
Figure 1.CR relative intensity ditribution in equatorial cooridnates obtained from the data of the Tibet ASγ experiment taken from November 1999 to May 2010.

Figure 2 .
Figure 2. Distribution of CR rigidity as seen by the Tibet ASγ experiment for three declination bands, reproduced with detailed MC simulations.

Figure 3 .
Figure 3. Best-fit model distributions of the CR relative intensity at the heliospheric outer boundaries at distances of 630 AU (a), 1580 AU (b) and 3980 AU (c) from the Sun, respectively.