Study of the Nearby Source Based on Cosmic Ray Nuclear Spectra and Anisotropy

1. Model Description

1.1. Spatially dependent diffusion

The region where CRs diffuse in the Galaxy is called a magnetic halo, which is usually approximated as a cylinder with its radial boundary equal to the Galactic radius, i.e. $R=20$ kpc and its half thickness $z_h$ is about a few kpc. $z_h$ is usually determined by fitting the B/C ratio along with diffusion coefficient [55]. Both CR sources and the interstellar medium are usually assumed to be concentrated near the galactic disk, whose average thickness $z_s$ is roughly 200 pc.

HAWC experiment obtained by observing B/C that the diffusion coefficient of CRs near the source of the galactic disk was two orders of magnitude lower than that away from the source. In the past few years, the SDP model of CRs with different propagation coefficients in the near and away galactic disk has been proposed and widely applied [21,56,57]. In the SDP model, the galactic diffusion halo is divided into two regions, i.e inner halo (IH) and outer halo (OH). The galactic disk and its surrounding region is called the IH, while the diffusion region outside the IH is called the OH. In IH region, where are more sources, the activity of supernova explosion will lead to more intense turbulence. Therefore, the diffusion of CRs will be slow, and the diffusion coefficient will be less dependent on the rigidity. Whereas in OH region, the diffusion of CRs is less affected by stellar activity, and diffusion coefficient is consistent with the traditional propagation model and only depends on rigidity.

In this work, we adopt SDP model and the diffusion coefficient is parameterized as [58,59]

$$\begin{array}{cc}\hfill D_{xx}(r,z,\mathcal{R})& =D_{0}F(r,z){\left({\displaystyle \frac{\mathcal{R}}{{\mathcal{R}}_0}}\right)}^{{\delta }_{0}F(r,z)}\hfill \end{array}$$

(1)

where r and z are cylindrical coordinates, $\mathcal{R}$ is particle’s rigidity and $D_0$ is a constant. The total half thickness of the propagation halo is $z_0$, and the half-thickness of the IH is $\xi z_0$. The parameterization of $F(r,z)$ can be parameterized as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill F(r,z)=\left\{\begin{array}{cc}g(r,z)+[1-g(r,z)]{\left({\displaystyle \frac{z}{\xi z_0}}\right)}^n,\hfill & \left|z\right|\le \xi z_0\hfill \\ 1,\hfill & \left|z\right|>\xi z_0\hfill \end{array}\right.\end{array}\end{array}$$

(2)

where $g(r,z)=N_m/[1+f(r,z)]$, and $f(r,z)$ is the source density distribution. In this work, we adopt numerical package DRAGON to solve the transport equation [60].

1.2. Background sources

The injection spectrum of background sources is assumed to be a power-law of rigidity with a high-energy exponential cutoff, $q\left(\mathcal{R}\right)\propto {\mathcal{R}}^{-\nu }exp(-\mathcal{R}/{\mathcal{R}}_{\mathrm{c}})$. The cutoff rigidity of each element could be either Z- or A-dependent.

The spatial distribution of sources takes the form of SNR distribution [61]:

$$f(r,z)\propto {(r/r_{\odot })}^{1.69}exp[-3.33(r-r_{\odot })/r_{\odot }]exp(-|z|/z_s),$$

(3)

where $r_{\odot }=8.5$ kpc and $z_s=0.2$ kpc.

1.3. Nearby source

We solve the time-varying propagation equation of CRs from the nearby sources assuming a spherical geometry with infinite boundary conditions and using Green’s function method [62,63].

The CR density of nearby source as a function of location,time and rigity is calculated by

$$\varphi (r,\mathcal{R},t)=\frac{q_{\mathrm{inj}}\left(\mathcal{R}\right)}{{\left(\sqrt{2\pi }\sigma \right)}^3}exp\left(-\frac{r^2}{2{\sigma }^2}\right),$$

(4)

where $q_{\mathrm{inj}}\left(\mathcal{R}\right)\delta \left(t\right)\delta \left(\mathbf{r}\right)$ is the instantaneous injection spectrum of a point source, $\sigma (\mathcal{R},t)=\sqrt{2D\left(\mathcal{R}\right)t}$ is the effective diffusion length within time t, $D\left(\mathcal{R}\right)$ is the diffusion coefficient which is adopted as the value nearby the solar system. The injection spectrum is also parameterized as a cutoff power-law form, $q_{\mathrm{inj}}\left(\mathcal{R}\right)=q_0{\mathcal{R}}^{-\alpha }exp(-\mathcal{R}/{\mathcal{R}}_{\mathrm{c}}^{\prime })$. The normalization $q_0$ is obtained through fitting to the CR energy spectra.

In this work, five nearby sources, i.e. Geminga, Monogem, Vela, Loop I and Cygnus are analyzed and their location and age are shown in Table 1.

1.4. Anisotropic Diffusion and Large-Scale Anisotropy

The amplitude of the dipole anisotropy is proportional to the spatial gradient of the CR density and the diffusion coefficient. The anisotropy can be written as [44,64]

$$\delta ={\displaystyle \frac{3D}{v}}{\displaystyle \frac{\nabla \psi }{\psi }}.$$

(5)

By observing neutral particles passing through the heliosphere boundary, the IBEX experiment revealed that the LIMF, in the range of 20 pc, follows $(l,b={210.5}^{\circ },-{57.1}^{\circ })$, which coincides with the anisotropic phase below the 100 TeV region [65]. Some work has also revealed that TeV cosmic ray anisotropy is related to LIMF [44,66]. CRs diffuse anisotropically in magnetic field, and the diffusion tensor $D_{ij}$ associated with the magnetic field is written as

$$D_{ij}\equiv D_{\perp }{\delta }_{ij}+\left(D_{\parallel }-D_{\perp }\right)b_i{b}_j,b_i={\displaystyle \frac{B_i}{|\overrightarrow{B}|}}$$

(6)

Where $D_{\Vert }$ and $D_{\perp }$ are the diffusion coefficients aligned parallel and perpendicular to the ordered magnetic field, $b_i$ is the i-th component of the unit vector [67], respectively. The values $D_{\Vert }$ and $D_{\perp }$ of rigidity in this work refer to the form [42,68] and are shown as follows,

$$\begin{array}{cc}\hfill D_{\Vert }& =D_{0\Vert }{\left(\frac{\mathcal{R}}{{\mathcal{R}}_0}\right)}^{{\delta }_{\parallel }},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill D_{\perp }& =D_{0\perp }{\left(\frac{\mathcal{R}}{{\mathcal{R}}_0}\right)}^{{\delta }_{\perp }}\equiv \epsilon D_{0\Vert }{\left(\frac{\mathcal{R}}{{\mathcal{R}}_0}\right)}^{{\delta }_{\perp }},\hfill \end{array}$$

(8)

where $\epsilon ={\displaystyle \frac{D_{0\perp }}{D_{0\Vert }}}$ is the ratio between perpendicular and parallel diffusion coefficient at the reference rigidity ${\mathcal{R}}_0$.

Under the anisotropic diffusion model, the form of formula 5 can be written as,

$$\delta ={\displaystyle \frac{3}{v\psi }}{D}_{ij}{\displaystyle \frac{\partial \psi }{\partial x_j}}.$$

(9)

2. Results and Discussion

2.1. $\mathrm{B}/\mathrm{C}$ ratio

We obtained the propagation parameters by fitting the B/C ratio. The comparison of the B/C ratio between the model prediction and the observation data of AMS-02 is shown in Figure 1, which indicates that the relevant parameters are reasonable. The corresponding propagation parameters are respectively $D_0=4.87\times {10}^{28}{\mathrm{cm}}^2$, ${\delta }_0=0.58$, $N_{\mathrm{m}}=0.62$, $\xi =0.1$, $n=4$. The Alfvénic velocity is $v_A=6\mathrm{km}·{\mathrm{s}}^{-1}$, and thhe half thickness of the propagation halo is $z_h=5\mathrm{kpc}$.

Figure 1. Fitting to B/C ratio with the Model prediction. The B/C data points are taken from AMS-02 experiment [69].

Figure 1. Fitting to B/C ratio with the Model prediction. The B/C data points are taken from AMS-02 experiment [69].

2.2. Proton and Helium spectra of five nearby sources

First, we calculate the proton and helium spectra for different five nearby sources as listed in Table 1. Z-dependent cutoff is applied to the injection spectra of background sources with a high-energy exponential cutoff. We obtained the cutoff rigidity by fitting the proton and helium spectra observed by KASCADE. Similarly, the injection spectra of the nearby sources are also based on Z-dependent cutoff. In order to better match the CR spectra, the injection spectra of the background sources are fine-tuned, when we calculate the contribution of different nearby sources to CR spectra. The corresponding injection parameters of background and nearby sources are shown in Table 2. Only the injection spectra for Geminga and Monogem sources are listed here.

Figure 2 presents the spectral results of proton (left) and helium (right), where the solid gray line is the contribution of the background sources, the dashdotted lines in different colors represent the contributions from different single nearby sources, and the solid lines in corresponding colors display the sum of single nearby sources and background sources. It can be seen that the contribution of Geminga, Monogem and Loop I SNR can account for spectral hardening at ∼ 200 GeV and softening features at ∼ 10 TeV, but Vela and Cygnus can not. Because Vela and Cygnus source are younger than the others, low-energy CRs produced by them are difficult to reach the solar system, so the flux of CRs below 100 TeV is lower.

Figure 2. Energy spectra of protons (left) and helium nuclei (right) from the different single nearby sources. The data points are taken from DAMPE [70,71], AMS-02 [11,72], CREAM-III [73], NUCLEON [74], KASCADE [75] and KASCADE-Grande [76] respectively. The grey solid lines (BKG) represent the fluxes of background sources, and the dashdotted lines in different colors are the fluxes from different single nearby SNR sources respectively, and the solid lines of corresponding color represent the sum contributions of the background and nearby sources.

Figure 2. Energy spectra of protons (left) and helium nuclei (right) from the different single nearby sources. The data points are taken from DAMPE [70,71], AMS-02 [11,72], CREAM-III [73], NUCLEON [74], KASCADE [75] and KASCADE-Grande [76] respectively. The grey solid lines (BKG) represent the fluxes of background sources, and the dashdotted lines in different colors are the fluxes from different single nearby SNR sources respectively, and the solid lines of corresponding color represent the sum contributions of the background and nearby sources.

2.3. Anisotropy of Geminga, Monogem and Loop I

Since only the energy spectra of Geminga, Monogem and Loop I sources are consistent with the experimental data, we will only analyze the anisotropy of these three sources next.

The anisotropic propagation diffusion coefficient of CRs is a tensor related to the magnetic field. In this work, we adopted the LIMF observed by IBEX. The parameters of parallel diffusion coefficient $D_{\Vert }$ are set as those in Section 2.1. CRs from TeV to PeV energy region are thought to travel faster parallel to the magnetic field than perpendicular to it, therefore we set $D_{\Vert }>D_{\perp }$, $\epsilon =0.01$ and the difference between ${\delta }_{\perp }$ and ${\delta }_{\Vert }$ is 0.32.

The amplitude and phase of CR anisotropy with Geminga SNR are shown in Figure 3. It is obvious that both phase and amplitude agree well with experimental data, which indicates that Geminga SNR is a good candidate source. Below 100 TeV, the phase points in the direction of the LIMF, which indicates that Geminga source and the LIMF deflection dominate the anisotropic phase, although the nearby source flux is sub-dominant. Above 100 TeV, the phase points to Galactic Center(GC) indicates background sources dominate, since galactic CR sources are more abundant in the inner galaxy

Figure 3. The amplitude (left) and phase (right) of anisotropy with the contribution from nearby Geminga SNR source. The data points are taken from underground muon detectors: Norikura [77], Ottawa [78], London [79], Bolivia [2], Budapest [80], Hobart [80], London [80], Misato [80], Socorro [80], Yakutsk [80], Banksan [81], Hong Kong [82], Sakashita [83], Utah [84], Liapootah [85], Matsushiro [86], Poatina [87], Kamiokande [88], Marco [89],SuperKamiokande [27]; and air shower array experiments: PeakMusala [89], Baksan [91], Norikura [92], EAS-TOP [38,93,94], Baksan [95], Milagro [29], IceCube [30,32], Ice-Top [33], ARGO-YBJ [36], Tibet [46,96,97], HAWC [43], HAWC-IceCube [43].

Figure 3. The amplitude (left) and phase (right) of anisotropy with the contribution from nearby Geminga SNR source. The data points are taken from underground muon detectors: Norikura [77], Ottawa [78], London [79], Bolivia [2], Budapest [80], Hobart [80], London [80], Misato [80], Socorro [80], Yakutsk [80], Banksan [81], Hong Kong [82], Sakashita [83], Utah [84], Liapootah [85], Matsushiro [86], Poatina [87], Kamiokande [88], Marco [89],SuperKamiokande [27]; and air shower array experiments: PeakMusala [89], Baksan [91], Norikura [92], EAS-TOP [38,93,94], Baksan [95], Milagro [29], IceCube [30,32], Ice-Top [33], ARGO-YBJ [36], Tibet [46,96,97], HAWC [43], HAWC-IceCube [43].

Figure 4 shows the amplitude and phase of the anisotropy with Monogem SNR source. It can be seen from Figure 4 that the amplitude is consistent with the experimental data, while the phase clearly contradicts the observations. The anisotropic phase of Monogem is similar to that of Geminga except at about 100 TeV. At about 100 TeV, the anisotropy contribution of the background sources and the nearby source is close, so the phase points to the direction of their vector synthesis. The direction of Monogem and background source synthes is contrary to the experimental observation, because Monogem $(l,b={203.0}^{\circ },{12}^{\circ })$is located above the galactic disk.

Figure 4. The amplitude (left) and phase (right) of the dipole anisotropy with the contribution from nearby Monogem SNR source.

Figure 4. The amplitude (left) and phase (right) of the dipole anisotropy with the contribution from nearby Monogem SNR source.

Figure 5 shows the anisotropy with Loop I SNR source. It is clear that neither amplitude nor phase agree with the measurements of experiment. Loop I source located in the direction of the GC are in the same direction as the background sources located mainly in the inner Galaxy, so their amplitude of synthesis increases as the energy increases, meanwhile the phase always points toward the GC.

Figure 5. The amplitude (left) and phase (right) of anisotropy with the contribution from nearby Loop I SNR source.

Figure 5. The amplitude (left) and phase (right) of anisotropy with the contribution from nearby Loop I SNR source.

By studying the anisotropy of three sources namely Geminga, Monogem and Loop I, we found that the location of the nearby source and LIMF dominate the phase of anisotropy below 100 TeV. For below 100 TeV, the anisotropic phase of Geminga and Monogem located in the anti-GC direction and close to the galactic disk is consistent with the experimental observations, while that of Loop I located in the GC direction is not. At about 100 TeV, the anisotropy of the Geminga source located below the galactic disk closer to the LIMF direction agrees with the experimental observations because the contributions of the nearby source and background sources are close, while that of the Monogem source located above the galactic disk does not. For above 100 TeV, the background source dominates, so the phase points towards GC.

3. Summary

In view of the anomalous results of recent experiments on CR spectra and anisotropy, we analyze the contributions of different five nearby sources including Geminga, Monogem, Loop I, Vela and Cygnus to CR nuclear spectra and anisotropy, based on the anisotropic diffusion of CRs induced by the LIMF nearby the solar system and SDP model. CR spectra results indicate that only older nearby sources such as Geminga, Monogem and Loop I can explain the nuclear spectral hardening at $\sim 200$ GeV, while younger ones such as Vela and Cygnus can not. CR anisotropy results reveal that the location of the nearby source and LIMF dominate together the phase of anisotropy below 100 TeV. Nearby Geminga source, which is located at the anti-GC, below the galactic disk and near the LIMF direction, can simultaneously explain the proton and helium spectral hardening at $\sim 200$ GeV, softening at ∼ 10 TeV and the anisotropy from 100 GeV to PeV energy region well, while other sources can not. The future high-precision measurements of CR spectra and anisotropy by experiments such as LHAASO and HAWC will help to check our model.