A simple and natural interpretations of the DAMPE cosmic-ray electron/positron spectrum within two sigma deviations

Abstract

The DArk Matter Particle Explorer (DAMPE) experiment has recently announced the first results for the measurement of total electron plus positron fluxes between 25 GeV and 4.6 TeV. A spectral break at about 0.9 TeV and a tentative peak excess around 1.4 TeV have been found. However, it is very difficult to reproduce both the peak signal and the smooth background including spectral break simultaneously. We point out that the numbers of events in the two energy ranges (bins) close to the 1.4 TeV excess have $1\sigma$ deficits. With the basic physics principles such as simplicity and naturalness, we consider the $-2\sigma$, $+2\sigma$, and $-1\sigma$ deviations due to statistical fluctuations for the 1229.3 GeV bin, 1411.4 GeV bin, and 1620.5 GeV bin. Interestingly, we show that all the DAMPE data can be explained consistently via both the continuous distributed pulsar and dark matter interpretations, which have $\chi ^{2} \simeq 17.2$ and $\chi ^{2} \simeq 13.9$ (for all the 38 points in DAMPE electron/positron spectrum with 3 of them revised), respectively. These results are different from the previous analysis by neglecting the 1.4 TeV excess. At the same time, we do a similar global fitting on the newly released CALET lepton data, which could also be interpreted by such configurations. Moreover, we present a $U(1)_D$ dark matter model with the Breit–Wigner mechanism, which can provide the proper dark matter annihilation cross section and escape the CMB constraint. Furthermore, we suggest a few ways to test our proposal.

1 Introduction

Because of the strong radiative cooling via synchrotron and inverse Compton scattering (ICS) processes, the TeV electrons can only travel a short distance of about a few kpc in the Milky Way. Therefore, the nearby Cosmic Ray (CR) sources such as pulsars [1,2,3,4,5] and dark matter (DM) [6,7,8] can be probed via the high energy electrons and positrons. The spectra of the cosmic-ray electrons and positrons (CREs) have been measured up to TeV energy scales by the ground-based and space-borne experiments, for example, HESS [9, 10], VERITAS [11, 12], FermiLAT [13, 14], AMS-02 [15,16,17,18], and CALET [19, 20]. In particular, the excesses of the electrons [9, 21,22,23] and positrons [24,25,26,27] have been discovered too.

Table 1 The original and revised numbers of events and fluxes, $\varDelta N$, and $\varDelta N_{2\sigma }$ for the 1229.3 GeV bin, 1411.4 GeV bin, and 1620.5 GeV bin. Here, $\varDelta N$ and $\varDelta N_{2\sigma }$ are the adjusted numbers of events and the numbers of events for $2\sigma$ deviations from statistical fluctuations. Thus, we should require $|\varDelta N| \le |\varDelta N_{2\sigma }|$

Recently, the DArk Matter Particle Explorer (DAMPE), which is a new generation space-borne experiment to measure CRs, and which was launched in December 2015, has announced the first results of high energy CR electron plus positron ($e^{-}+e^{+}$) flux from 25 GeV to 4.6 TeV with unprecedentedly high quality [28]. The energy resolution of the DAMPE is better than 1.5% at TeV energies, and the hadron rejection power is about $10^5$. Thus, DAMPE is able to reveal (fine) structures of the electron and positron fluxes. The main DAMPE spectrum can be fitted by a smoothly broken power-law model with a spectral break around 0.9 TeV, which confirms the previous results by HESS experiment [9, 10]. And there exists a tentative peak-like flux excess around 1.4 TeV. Thus, the DAMPE results have stimulated the extensive studies [29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59]. The spectral break can be explained by the broad distributed pulsars, pulsar wind nebulae (PWNe), supernova remnants (SNRs) [30, 32], and by the dark matter annihilation and decay in the galaxy halo [32, 40, 43, 44]. Also, the tentative peak is always interpreted by local pulsars, PWNe, and SNRs [30, 32, 58], and by the DM sub-halos, clumps, and mini-spikes [29, 31, 33, 38,39,40, 44, 45, 51, 53]. Another important interpretation ascribes observed CR spectrum puzzling features to a nearly 2–3 Myr Super Nova [60, 61], which could naturally explain not only the proton to helium ratio, and the positron and anti-proton fluxes, but also the plateau in the cosmic-ray dipole anisotropy.

However, one can easily show that it is impossible to explain the spectral break and the tentative peak simultaneously [30, 32, 44, 45, 58]. In addition, we have 74, 93, and 33 events for three continuous bins or energy ranges [1148.2, 1318.3] GeV, [1318.3, 1513.6] GeV, and [1513.6, 1737.8] GeV, respectively, which for simplicity we shall call 1229.3 GeV bin, 1411.4 GeV bin, and 1620.5 GeV bin [28]. The number of events and fluxes for these bins are given in Table 1. From Fig. 2 of the DAMPE’s paper [28], it is obvious that the 1411.4 GeV bin has a little bit more than $3\sigma$ excess, while the 1229.3 GeV bin and 1620.5 GeV bin have about $1\sigma$ deficits. Therefore, it is very difficult to explain the events in these three bins, especially the first two, no matter by the pulsar or dark matter interpretations.

From the theoretical physics point of view, we would like to explain nature with basic principles such as simplicity and naturalness, or say truth and beauty! In the words of Sir Isaac Newton, “Truth is ever to be found in the simplicity, and not in the multiplicity and confusion of things.” Therefore, to explain all the DAMPE data in a simple and natural way, we propose that the excess in the 1411.4 GeV bin and the deficits in the 1229.3 GeV bin and 1620.5 GeV bin arise from the $+2\sigma$, $-2\sigma$, and $-1\sigma$ deviations due to statistical fluctuations, which has happened frequently in collider experiments. Remarkably, we can indeed explain all the DAMPE data consistently via the pulsar and dark matter interpretations, which have $\chi ^{2} \simeq 17.2$ and $\chi ^{2} \simeq 13.9$ (for all the 38 points in DAMPE electron/positron spectrum with 3 of them revised), respectively. Our results are different from the previous analysis by neglecting the 1.4 TeV excess [43]. As a comparison, the newly released CALET lepton data is used to do a similar global fitting, which could give us some more supports on the origin of the lepton excess. In addition, we present a $U(1)_D$ dark matter model with Breit–Wigner mechanism, which can provide the proper dark matter annihilation cross section and escape the CMB constraint. Furthermore, we suggest a few ways to test our proposal and the 1.4 TeV excess.

2 Statistical fluctuations

In the DAMPE’s paper [28], the numbers of events and the CRE fluxes with $1\sigma$ statistical and systematic errors have been given in its Table 1. To evaluate the uncertainties for numbers of the events, we need to understand their relations. The relation between the number of events and fluxes in each energy bin is [15, 28]

$$\begin{aligned} \varPhi (e^{-}+e^{+}) = \frac{N(E) \cdot (1 - \varepsilon _{\mathrm {bg}}(E))}{A_{\mathrm {eff}}(E) \cdot T \cdot \varDelta E} \cdot \varepsilon _{\mathrm {other}}(E), \end{aligned}$$

(1)

where N is the number of ($e^{-}+e^{+}$) events, $A_{\mathrm {eff}}$ is the effective detector acceptance, T is the operating time, $\varDelta E$ is the energy range of the bin, $\varepsilon _{\mathrm {bg}}$ is the background fraction of the events, and $\varepsilon _{\mathrm {other}}$ represents the effects caused by other mechanisms which were not given in Table 1 of Ref. [28].

Taking $T = 530$ days and $\varepsilon _{\mathrm {other}}= 1.3$, we can reproduce the corresponding results in the 1229.3 GeV bin, 1411.4 GeV bin, and 1620.5 GeV bin within the uncertainty $< 0.1 \%$. Consequently, we use the formula

$$\begin{aligned} \varPhi (e^{-}+e^{+}) = \frac{N(E) \cdot (1 - \varepsilon _{\mathrm {bg}}(E))}{A_{\mathrm {eff}}(E) \cdot T \cdot \varDelta E} \cdot 1.3 \end{aligned}$$

(2)

in this letter to calculate the fluxes in these bins.

We calculate the $2\sigma$ deviations for the number of events ($\varDelta N_{2\sigma }$) from the flux statistical fluctuations as follows:

$$\begin{aligned} \varDelta N_{2\sigma } = \frac{\varDelta \varPhi (e^{-}+e^{+})_{2 \sigma _{\mathrm {stat}}}}{\varPhi (e^{-}+e^{+})} \cdot N. \end{aligned}$$

(3)

Thus, for the 1229.3 GeV bin, 1411.4 GeV bin, and 1620.5 GeV bin, we obtain $\varDelta N_{2\sigma } = \pm 18,~\pm 20,~\pm 12$, respectively. Assume $-2\sigma$, $+2\sigma$, and $-1\sigma$ deviations for these bins from statistical fluctuations, we have $\varDelta N = + 18,~- 20,~+ 6$, respectively. Therefore, the revised numbers of events for the 1229.3 GeV bin, 1411.4 GeV bin, and 1620.5 GeV bin, are 92, 73, and 39, respectively.

Furthermore, we reestimate the statistical uncertainties in these bins based on the revised numbers of events via the formula

$$\begin{aligned} \varDelta N_{1\sigma } \simeq \frac{1}{\sqrt{N}}, \end{aligned}$$

(4)

and then calculate the corresponding fluxes and their statistical uncertainties. The systematical uncertainties are assumed to be invariant. All the detailed pieces of information for these three bins are given in Table 1. By the way, as a cross check, with Eq. (4), we have reproduced similar $1\sigma$ statistical uncertainties of the original fluxes in DAMPE’s paper [28].

3 Fitting procedure

In CR theory, the CR electrons are considered to be accelerated during the acceleration of CR nuclei at the sources, e.g. SNRs. On the other hand, the CR positrons are produced as secondary particles from CR nuclei interaction with the interstellar medium (ISM) [17, 24, 62, 63]. From the observed spectra of positrons and electrons [15,16,17,18], we can conclude that there should exist some extra sources producing electron–positron pairs. As we stated in the first section, these extra sources could be astrophysical sources or DM annihilation or decay. As a result, the CREs data contains: (i) the primary electrons; (ii) the secondary electrons; (iii) the secondary positrons; (iv) the extra source of electron–positron pairs. If we want to study the properties of the extra source, we should deduce the primary electrons and secondary electrons/positrons first.

The primary electrons are always assumed to have a power-law injection and the secondary electrons/positrons are determined mainly by the CR proton and helium nucleus interact with ISM. Consequently, we should do global fitting to all these related data simultaneously which can avoid the bias of choosing the lepton background parameters.

The public code dragon¹ [64] was used to do numerical calculations. Some custom modifications are performed in the original code, such as the possibility to use specie-dependent injection spectra, which is not allowed by default in dragon.

In view of some discrepancies when fitting the new data [65], we use a factor $c_{\,\mathrm {He}}$ to re-scale the helium-4 abundance (which has a default value of $7.199 \times 10^4$), which helps us to get a better global fitting.

3.1 Background

In this work, we use the widely used diffusion–reacceleration model, which can give consistent fitting results to the AMS-02 nuclei data (see, e.g., [66,67,68,69]). In the whole propagation region, a uniform diffusion coefficient ($D_{xx} = D_{yy}= D_{zz} = D_0\beta \left( R/R_0 \right) ^{\delta }$) is employed to describe the propagation.

The hardening of the nuclei spectra at $\sim 300 \,\mathrm {GeV}$ (which has been observed by ATIC-2 [70], CREAM [71], PAMELA [72], and AMS-02 [73, 74]) is considered by adding breaks in the primary source injections. At the same time, considering the observed significant difference in the slopes of proton and helium (of about $\sim 0.1$ [72,73,74]), we use separate primary source spectra settings for proton and helium. In summary, for the nuclei’s primary source injections, each of them has two breaks at rigidity $R_{\,\mathrm {A}1}$ and $R_{\,\mathrm {A}2}$. The corresponding slopes are $\nu _{\,\mathrm {A}1}$ ($R \le R_{\,\mathrm {A}1}$), $\nu _{\,\mathrm {A}2}$ ($R_{\,\mathrm {A}1} < R \le R_{\,\mathrm {A}2}$) and $\nu _{\,\mathrm {A}3}$ ($R > R_{\,\mathrm {A}3}$).

For the CR electrons’ primary source, we use one break $R_{e}$ for the electron primary source, and the corresponding slopes are $\nu _{e1}$ ($R \le R_{e}$) and $\nu _{e2}$ (($R > R_{e}$)).

In order to take into account the uncertainties when calculating the secondary CR particles’ fluxes, we employ parameters $c_{\,\bar{\mathrm {p}}}$ and $c_{e^{+}}$ to re-scale the calculated secondary flux to fit the data [75,76,77,78,79]. Note that the above-mentioned uncertainties may not be simply represented by a constant factor, but most probably are energy dependent [80, 81]. Here we expect that a constant factor could be a simple assumption.

The force-field approximation [82] is used to describe the effects of solar modulation effects. $\phi _{\mathrm {nuc}}$, $\phi _{\,\bar{\mathrm {p}}}$ and $\phi _{e^{+}}$ are used to modulate the local interstellar spectra of nuclei (proton and helium), anti-proton and positrons, respectively, which was based on the charge-sign dependent solar modulation. On the other hand, because the DAMPE lepton data $\gtrsim 20 \,\mathrm {GeV}$, we did not consider the modulation effects on the electrons (or leptons).

Fig. 1

The global fitting results of the DAMPE and CALET lepton flux for pulsar scenario. The $2\sigma$ (deep color) and $3\sigma$ (light color) bounds of total fitted results (red), contribution from background (blue) and pulsar (green) are also shown in the figure. We have $\chi ^{2} \simeq 17.2$ (DAMPE) and $\chi ^{2} \simeq 16.1$ (CALET). (Data sources: DAMPE [28], CALET [20], Fermi [94], HESS [9, 10] and AMS-02 [15].)

Fig. 2

The global fitting results of the DAMPE and CALET lepton flux for DM scenario. The $2\sigma$ (deep color) and $3\sigma$ (light color) bounds of the total fitted results (red), and the contribution from background (blue) and DM (green) are also shown in the figure. We have $\chi ^{2} \simeq 13.9$ (DAMPE) and $\chi ^{2} \simeq 25.7$ (CALET)

3.2 Extra sources

We consider both pulsar and DM scenarios to generate the CRE excesses in the observed spectrum by the DAMPE experiment. For the pulsar scenario, a continuous distributed pulsar background was used [43, 79]. The injection spectrum of such sources is assumed to be a power law with an exponential cut-off,

$$\begin{aligned} q_e^{\,\mathrm {psr}}(p) = N_{\,\mathrm {psr}}(R/\mathrm {10\,\mathrm {GeV}})^{-\nu _{\,\mathrm {psr}}} \exp {(-R/R_\mathrm {c})}, \end{aligned}$$

(5)

where $N_{\,\mathrm {psr}}$ is the normalization factor, $\nu _{\,\mathrm {psr}}$ is the spectral index, and $R_\mathrm {c}$ is the cut-off rigidity. For the DM scenario, we employ the Einasto profile [83,84,85,86],

$$\begin{aligned} \rho (r)=\rho _\odot \exp \left[ -\left( \frac{2}{\alpha }\right) \left( \frac{r^{\alpha }-r_\odot ^{\alpha }}{r_{s}^{\alpha }} \right) \right] , \end{aligned}$$

(6)

with $\alpha \approx 0.17$, $r_{s}\approx 20 \,\mathrm {kpc}$, and $\rho _\odot \approx 0.39 \,\mathrm {GeV}\,\mathrm {cm}^{-3}$ is the local DM relic density [87,88,89,90,91]. The source term, which we use to add the CRE particles from the annihilations of the Majorana DM particles, is

$$\begin{aligned} Q({r},p)=\frac{\rho ({r})^2}{2 m_{\chi }^2}\langle \sigma v \rangle \sum _{f} \eta _{f} \frac{\mathrm{d}N^{(f)}}{\mathrm{d}p} , \end{aligned}$$

(7)

where $\langle \sigma v \rangle$ is the velocity-averaged DM annihilation cross section multiplied by DM relative velocity (referred as cross section), $\rho ({r})$ is the DM density distribution, and d$N^{(f)}/\mathrm{d}p$ is the injection energy spectrum of CREs from DM annihilating into the Standard Model (SM) final states via leptonic channels $f{\bar{f}}$ ($e^{-}e^{+}$, $\mu \bar{\mu }$, and $\tau \bar{\tau }$) with $\eta _{f}$ ($\eta _{e}$, $\eta _{\mu }$, and $\eta _{\tau }$) the corresponding branching fractions. Here, we normalized $\eta _{f}$ as $\eta _{e}+ \eta _{\mu }+ \eta _{\tau }= 1$.

The parameters related to the extra source of the leptons for the pulsar scenario are $(N_{\,\mathrm {psr}}, \nu _{\,\mathrm {psr}}, R_{c})$, and for the DM scenario they are $(m_{\chi }, \langle \sigma v \rangle , \eta _{e}, \eta _{\mu }, \eta _{\tau })$.

3.3 Data sets and parameters

As in Ref. [43], we perform a global fitting on the data set including the proton fluxes from AMS-02 and CREAM [71, 73] helium flux from AMS-02 and CREAM² [71, 74], $\,\bar{\mathrm {p}}/\mathrm {p}$ ratio from AMS-02 [92], positrons flux from AMS-02 [16], and CRE flux from DAMPE [28], which could account for the primary electrons, the secondary leptons, and the extra leptons in a self-consistent way.³ Moreover, the employed AMS-02 positron flux is used to calibrate the positron contribution in the DAMPE CRE flux in energy region $\lesssim 300 \,\mathrm {GeV}$.⁴ The framework of the fitting procedure is the same as our previous work [43, 66], where the details can be found.

Considering the systematics between different CRE spectra observed by different experiments (see Figs. 1 and 2), we take the newly released CRE spectrum from CALET [20] as a comparison to do global fitting as that on DAMPE CRE spectrum. Because both of these experiments are implemented in space and have a similar ability to obtain CREs data, a reasonable explanation on lepton excess should explain them simultaneously. In this case, DAMPE and CALET could be considered as the maximum and minimum CREs flux examples.

Altogether, the data set in our global fitting is

$$\begin{aligned} D =&\{D^{\text {AMS-02}}_{\,\mathrm {p}}, D^{\text {AMS-02}}_{\,\mathrm {He}}, D^{\text {AMS-02}}_{\,\bar{\mathrm {p}}/\mathrm {p}}, D^{\text {CREAM}}_{\,\mathrm {p}}, \\&D^{\text {CREAM}}_{\,\mathrm {He}}, D^{\text {AMS-02}}_{e^{+}}, D^{\text {DAMPE}}_{e^{-}+e^{+}} / D^{\text {CALET}}_{e^{-}+e^{+}} \}~. \end{aligned}$$

The parameter sets for the pulsar scenario is

$$\begin{aligned} {\theta }_{\,\mathrm {psr}}&= \{ D_{0}, \delta , z_{h}, v_{A}, | N_{\,\mathrm {p}}, R_{\,\mathrm {p}1}, R_{\,\mathrm {p}2}, \nu _{\,\mathrm {p}1}, \nu _{\,\mathrm {p}2}, \nu _{\,\mathrm {p}3}, \\&\quad R_{\,\mathrm {He}1}, R_{\,\mathrm {He}2}, \nu _{\,\mathrm {He}1}, \nu _{\,\mathrm {He}2}, \nu _{\,\mathrm {He}3}, | \phi _{\mathrm {nuc}}, \phi _{\,\bar{\mathrm {p}}}, c_{\,\mathrm {He}}, c_{\,\bar{\mathrm {p}}}, | \\&\quad N_{\,\mathrm {e}}, R_{\,\mathrm {e}1}, \nu _{\,\mathrm {e}1}, \nu _{\,\mathrm {e}2}, | \\&\quad N_{\,\mathrm {psr}}, \nu _{\,\mathrm {psr}}, R_{c}, | \\&\quad \phi _{e^{+}}, c_{e^{+}} \}, \end{aligned}$$

for the DM scenario it is

$$\begin{aligned} {\theta }_{\,\mathrm {DM}}&= \{ D_{0}, \delta , z_{h}, v_{A}, | N_{\,\mathrm {p}}, R_{\,\mathrm {p}1}, R_{\,\mathrm {p}2}, \nu _{\,\mathrm {p}1}, \nu _{\,\mathrm {p}2}, \nu _{\,\mathrm {p}3}, \\&\quad R_{\,\mathrm {He}1}, R_{\,\mathrm {He}2}, \nu _{\,\mathrm {He}1}, \nu _{\,\mathrm {He}2}, \nu _{\,\mathrm {He}3}, | \phi _{\mathrm {nuc}}, \phi _{\,\bar{\mathrm {p}}}, c_{\,\mathrm {He}}, c_{\,\bar{\mathrm {p}}}, | \\&\quad N_{\,\mathrm {e}}, R_{\,\mathrm {e}1}, \nu _{\,\mathrm {e}1}, \nu _{\,\mathrm {e}2}, | \\&\quad m_{\chi }, \langle \sigma v \rangle , \eta _{e}, \eta _{\mu }, \eta _{\tau }, | \\&\quad \phi _{e^{+}}, c_{e^{+}} \}. \end{aligned}$$

Note that most of these two scenarios’ parameters in the set ${\theta }_{\,\mathrm {psr}}$ and ${\theta }_{\,\mathrm {DM}}$ are the same, except for those which account for the extra sources of leptons.

4 Results

When the Markov chains have reached their equilibrium state, we take the samples of the parameters as their posterior probability distribution functions. The best-fit values, statistical mean values, standard deviations and allowed intervals at $95 \%$ CL for parameters in the sets ${\theta }_{\,\mathrm {psr}}$ and ${\theta }_{\,\mathrm {DM}}$ for DAMPE and CALET are shown in the appendix, Table 2 (DAMPE, pulsar scenario), Table 3 (CALET, pulsar scenario), Table 4 (DAMPE, DM scenario), and Table 5 (CALET, DM scenario), respectively. For the best-fit results of the global fitting, we got $\chi ^{2}/\mathrm{d.o.f} = 243.13/299$ (DAMPE, pulsar scenario), $\chi ^{2}/\mathrm{d.o.f} = 229.98/301$ (CALET, pulsar scenario), $\chi ^{2}/\mathrm{d.o.f} = 262.94/297$ (DAMPE, DM scenario), and $\chi ^{2}/\mathrm{d.o.f} = 265.03/299$ (CALET, DM scenario).

4.1 Background

The best-fitting results and the corresponding residuals of the proton flux, helium flux and $\,\bar{\mathrm {p}}/\mathrm {p}$ ratio for the pulsar scenario are showed in the appendix, Fig. 5, and for the DM scenario they are showed in the appendix, Fig. 6. In these figures, we can see that the nuclei data is perfectly reproduced, which would provide a good precondition for the fitting on the lepton data.

Although the main purpose of this work does not focus on the CR propagation models, we would like to emphasize some points here: (i) As shown in the appendix, Tables 2, 3, 4, and 5,we got larger best-fit values of $D_{0}$ and $z_{h}$ than previous work [66]. This is so mainly because the newly released AMS-02 nuclei spectra favor large values of $D_{0}$ and $z_{h}$.⁵ Moreover, the employed two breaks in the nuclei’s primary source injection strengthen the classical degeneracy between $D_{0}$ and $z_{h}$ based on the data set we used in this work (without B/C); these both got larger best-fit values in this work. (ii) The employed two breaks in the nuclei primary source injection accounted for the observed hardening in the observed spectra, rather than use only one break and let $\delta$ compromise the different slopes in high energy regions, which leads to a smaller value of $\delta$ and fitting uncertainties on $\delta$ ($\lessapprox 0.01$) than previous work.

4.2 Extra sources

The fitting results of the pulsar and DM scenario on DAMPE and CALET CRE spectrum are given in Figs. 1 and 2, respectively, which also show some CREs’ spectrum from other experiments. Cleaner fitting results are shown in Figs. 3 and 4.⁶ From these figures, we can conclude that the two scenarios could provide excellent fittings to the DAMPE CRE spectrum within $3\sigma$ fitting deviation, which do not need to employ extra local sources. At the same time, the CALET CRE spectrum could also be fitted by the two scenarios, although the fitting result is not that good, because of a suspected bump at about 0.9–1 $\,\mathrm {TeV}$.

For the best-fit result on the DAMPE and CALET CRE spectrum, we got $\chi ^{2} \simeq 17.2$ (DAMPE, pulsar scenario), $\chi ^{2} \simeq 16.1$ (CALET, pulsar scenario), $\chi ^{2} \simeq 13.9$ (DAMPE, DM scenario) and $\chi ^{2} \simeq 25.7$ (CALET, DM scenario). This indicates that the CALET CREs data disfavors the DM scenario because of the defective fit on the suspected bump at about 0.9–1 $\,\mathrm {TeV}$. This needs more events accumulated in the future.

The detailed results of the constraints on the parameters can be found in the appendix, Tables 2, 3, 4, and 5. As regards most of the parameters there are slight differences between DAMPE and CALET fitting because of the systematics between them. One point we want to mention is that in the fitting results of CALET CREs data, the DM scenario, $\eta _{e}\simeq 0.930$, $\eta _{\mu }\simeq 0.042$, and $\eta _{\tau }\simeq 0.028$, which is largely different from the DAMPE results.

In the following, we will focus on the analysis of the DAMPE results, whose method could be extended to dealing with the CALET results without difficulties.

Fig. 3

The global fitting results and the corresponding residuals to the lepton (DAMPE and CALET) and positron (AMS-02) flux for pulsar scenario. The $2\sigma$ (deep color) and $3\sigma$ (light color) bounds of the total fitted results (red), and the contribution from background (blue) and pulsar (green), are also shown in the figure. Different from Figs. 1 and 2, all the experimental data which did not participate in the global fitting has not been not represented

Fig. 4

The same as Fig. 3, but for the DM scenario

For the pulsar scenario, the fitting results give $\nu _{\,\mathrm {psr}} \simeq 0.62$, which is obviously different from the fitting results in previous work (see, e.g., [95]). In standard pulsar models, the injection spectrum indices of CREs from pulsars are always in the range $\nu _{\,\mathrm {psr}} \in [1.0,2.4]$ [96,97,98]. As a result, more attention should be paid in future research. This may indicate: (i) there is something wrong or there is an inaccuracy with the classical pulsar CRE injection model; (ii) the CRE excess is not contributed primarily by pulsars. Moreover, the cut-off is $R_c \simeq 692$ GV. In the previous work [43], where the 1.4 TeV peak excess was neglected, we found that the spectral index of the injection is $\nu _{\,\mathrm {psr}} \simeq 0.65$ and the cut-off is $R_c \simeq 650$ GV. Thus, there exist about $+5\%$ and $-5\%$ deviations for $\nu _{\,\mathrm {psr}}$ and $R_c$, respectively.

For the DM scenario, we obtain $\langle \sigma v \rangle \simeq 4.07 \times 10^{-23} \,\mathrm {cm}^{2} \,\mathrm {s}^{-1}$ and $m_{\chi }\simeq 1884 \,\mathrm {GeV}$. The value of $\langle \sigma v \rangle$ is about 3 orders larger than that of thermal DM [99]. To explain this discrepancy, we will present a concrete model in the next section. Moreover, we have $\eta _{e}\simeq 0.465$, $\eta _{\mu }\simeq 0.510$, and $\eta _{\tau }\simeq 0.025$. So the DM annihilation into $\tau \bar{\tau }$ is highly suppressed, which provides some hints to constructing an appropriate DM model. In our previous work [43], where the 1.4 TeV peak excess was neglected, we have $\langle \sigma v \rangle \simeq 1.48 \times 10^{-23} \,\mathrm {cm}^{2} \,\mathrm {s}^{-1}$, $m_{\chi }\simeq 1208 \,\mathrm {GeV}$, $\eta _{e}\simeq \eta _{\mu }\simeq 0.5$, while $\eta _{\tau }$ is highly suppressed. Thus, we have similar results on the branching fractions, but different DM masses and annihilation cross sections.

5 Model building

Because we have $\eta _{e}\sim 0.465$, $\eta _{\mu }\sim 0.510$, and $\eta _{\tau }\sim 0.025$, the constraints from the Fermi-LAT observations of dwarf spheroidal galaxies [100,101,102,103,104,105] can be avoided [32]. To escape the constraints from the Planck observations of CMB anisotropies [106], we employ the Breit–Wigner mechanism [107,108,109,110,111,112,113,114]. We consider the dark $U(1)_D$ model where the SM fermions and Higgs fields are neutral under it. We introduce one SM singlet Higgs field S, one chiral fermionic dark matter particle $\chi$, and three pairs of the vector-like particles (${\widehat{XE}}_i, {\widehat{XE}}_i^c)$, whose quantum numbers under $SU(3)_C\times SU(2)_L \times U(1)_Y \times U(1)_D$ are

$$\begin{aligned}&S: (\mathbf {1}, \mathbf {1}, \mathbf {0}, \mathbf {2}),\, \chi : (\mathbf {1}, \mathbf {1}, \mathbf {0}, \mathbf {-1}) \nonumber \\&{\widehat{XE}}_i: (\mathbf {1}, \mathbf {1}, \mathbf {-1}, \mathbf {-2}),\, {\widehat{XE}}^c_i: (\mathbf {1}, \mathbf {1}, \mathbf {1}, \mathbf {2}). \end{aligned}$$

(8)

The relevant Lagrangian is

$$\begin{aligned} -\mathcal{L}= & {} -m_S^2 |S|^2 + \frac{\lambda }{2} |S|^4 + (M^V_{ij} {\widehat{XE}}^c_i {\widehat{XE}}_j \nonumber \\&+ y_{ij} S {\widehat{E}}_i^c {\widehat{XE}}_j +y S \chi \chi + \mathrm{H.C.}), \end{aligned}$$

(9)

where ${\widehat{E}}_i^c$ are the right-handed charged leptons. For simplicity, we choose $M^V_{ij} = M^V_{i} \delta _{ij}$ and $y_{ij} = y_i \delta _{ij}$. After S acquires a vacuum expectation value (VEV), the $U(1)_D$ gauge symmetry is broken down to a $Z_2$ symmetry under which $\chi$ is odd. Thus, $\chi$ is a DM matter candidate. For simplicity, we assume that the mass of the $U(1)_D$ gauge boson is about twice of the $\chi$ mass, i.e., $M_{Z'} \simeq 2 m_{\chi }$, while the Higgs field S and vector-like particles are heavier than $M_{Z'}$. Moreover, ${\widehat{E}}_i^c$ and ${\widehat{XE}}_i^c$ will be mixed due to the $M^V_{i}\,\,\,\, {\widehat{XE}\,\,}^c_i \,\,\,\,{\widehat{XE}}_i$ and $y_{i} S {\widehat{E}}_i^c {\widehat{XE}}_i$ terms, and we obtain the mass eigenstates $E_i^c$ and $XE_i^c$ by neglecting the tiny charged lepton masses,

$$\begin{aligned} \left( \begin{array}{c} E_i^c \\ XE_i^c \end{array} \right) = \left( \begin{array}{cc} \cos \theta _i &{} \sin \theta _i \\ -\sin \theta _i &{} \cos \theta _i \end{array} \right) \left( \begin{array}{c} {\widehat{E}}_i^{c} \\ {\widehat{XE}}_i^{c \prime } \end{array} \right) , \end{aligned}$$

(10)

where $\tan \theta _i = -y\langle S \rangle /M^V_i$.

Neglecting the charged lepton masses again, we obtain

$$\begin{aligned} \sigma v = \sum _{i=1}^3\frac{g'^4\sin ^2\theta _i}{6\pi }\frac{s-m_\chi ^2}{(s-m_{Z'}^2)^2+(m_{Z'}\varGamma _{Z'})^2}~,~\, \end{aligned}$$

(11)

where $m_\chi = y \langle S \rangle$, and $g'$ and $M_{Z'}$ are the gauge coupling and gauge boson mass for the $U(1)_D$ gauge symmetry.

For $m_{Z'}\simeq 2m_\chi$, $Z'$ decays dominantly into leptons, and the decay width is

$$\begin{aligned} \varGamma _{Z'}=\sum _{i=1}^3 \frac{g'^2\sin ^2\theta _i}{6\pi } m_{Z'}. \end{aligned}$$

(12)

To explain the DM best-fit results, we choose

$$\begin{aligned}&g' \simeq 0.028,~ m_\chi \simeq 1884~\text{ GeV },~ \frac{m_{Z'}-2m_\chi }{m_{Z'}} \simeq 3.0\times 10^{-6},~\nonumber \\&\sin \theta _e \simeq 0.21,~~~\sin \theta _\mu \simeq 0.22,~~~ \sin \theta _\tau \simeq 0.05. \end{aligned}$$

(13)

Then we obtain $\langle \sigma v\rangle \simeq 4.07\times 10^{-23} \text{ cm }^3\text{ s }^{-1}$, and $\eta _e :\eta _\mu : \eta _\tau \simeq 0.465 : 0.510 : 0.025$. Of course, there exists fine-tuning between $m_{Z'}$ and $m_\chi$, which deserves further study. For some solutions, see Ref. [113].

6 Discussions and conclusion

First, we would like to point out that if the numbers of events in the 1229.3 GeV bin and 1411.4 GeV bin are exchanged, we can also explain the DAMPE’s data similarly. Of course, the most important question is how to test our proposal that there exists statistical fluctuations in the 1229.3 GeV bin, 1411.4 GeV bin, and 1620.5 GeV bin. For the data analysis, we suggest that one chooses different energy ranges to study the data again. For example, we can shift the energy ranges by $\pm 50$ GeV and $\pm 100$ GeV for the high energy bins, and then study the corresponding events and fluxes. In the future, DAMPE will provide us more accurate spectrum data reaching up to $\sim 10 \,\mathrm {TeV}$, which can give us an unprecedented opportunity to study the origin and propagation of CREs. We predict that the CRE spectrum would be more continuous. In particular, the peak excess in the 1411.4 GeV bin and the deficits in the 1229.3 GeV bin and 1620.5 GeV bin will all decrease! Moreover, if the 1.4 TeV peak signal was proved to be correct, we do need a local source of high energy CREs. Another experiment is needed as a cross check if such a signal arises from DM annihilation; for example, our recent work [115] proposed a novel scenario to probe the interaction between DM particles and electrons for the DM mass range $5 \,\mathrm {GeV}\lesssim m_{\chi }\lesssim 10 \,\mathrm {TeV}$.

In summary, with the physics principles of simplicity and naturalness, we proposed that there exists $-2\sigma$, $+2\sigma$, and $-1\sigma$ deviations due to statistical fluctuations for the 1229.3 GeV bin, 1411.4 GeV bin, and 1620.5 GeV bin of the DAMPE data. Interestingly, we showed that all the DAMPE data can be explained consistently via both the pulsar and the dark matter interpretations, which have $\chi ^{2} \simeq 17.2$ and $\chi ^{2} \simeq 13.9$ (for all the 38 points in DAMPE electron/positron spectrum with 3 of them revised), respectively. These results are different from the previous analysis by neglecting the 1.4 TeV excess. At the same time, we employed the newly released CALET CREs spectrum to do a similar global fitting, which cold also be fitted by continuous distributed pulsar and DM scenarios. Moreover, we presented a $U(1)_D$ dark matter model with Breit–Wigner mechanism, which can provide the proper dark matter annihilation cross section and escape the CMB constraint. Furthermore, we suggested a few ways to test our proposal.

Appendix

See Tables 2, 3, 4, and 5.

Table 2 Constraints on the parameters of DAMPE CREs spectrum, pulsar scenario. The prior interval, best-fit value, statistic mean, standard deviation and the allowed range at $95\%$ CL are listed for the parameters. For the best-fit values, we have $\chi ^{2}/\mathrm{d.o.f.} = 243.13/299$

Table 3 Constraints on the parameters of CALET CRE spectrum, pulsar scenario. The prior interval, best-fit value, statistic mean, standard deviation and the allowed range at $95\%$ CL are listed for the parameters. For the best-fit values, we have $\chi ^{2}/\mathrm{d.o.f.} = 229.98/301$

Table 4 Constraints on the parameters of DAMPE CRE spectrum, DM scenario. The prior interval, best-fit value, statistic mean, standard deviation and the allowed range at $95\%$ CL are listed for the parameters. For the best-fit values, we have $\chi ^{2}/\mathrm{d.o.f.} = 262.94/297$

Table 5 Constraints on the parameters of CALET CRE spectrum, pulsar scenario. The prior interval, best-fit value, statistic mean, standard deviation and the allowed range at $95\%$ CL are listed for the parameters. For the best-fit values, we have $\chi ^{2}/\mathrm{d.o.f.} = 265.03/299$

Fig. 5

The global fitting results and the corresponding residuals to the proton flux, helium flux, $\,\bar{\mathrm {p}}/\mathrm {p}$ ratio, and positron flux for two scenarios. The $2\sigma$ (deep red) and $3\sigma$ (light red) bound are also showed in the figures

Fig. 6

The global fitting results and the corresponding residuals to the proton flux, helium flux, $\,\bar{\mathrm {p}}/\mathrm {p}$ ratio, and positron flux for two scenarios. The $2\sigma$ (deep red) and $3\sigma$ (light red) bound are also showed in the figures

See Figs. 5 and 6.