Two-mediator dark matter models and cosmic electron excess

The cosmic electron energy spectrum recently observed by the DAMPE experiment exhibits two interesting features, including a break around 0.9 TeV and a sharp resonance near 1.4 TeV. In this analysis, we propose a dark matter explanation to both exotic features seen by DAMPE. In our model, dark matter annihilates in the galaxy via two different channels that lead to both a narrow resonance spectrum near 1.4 TeV and electron excess events over an extended energy range thus generating the break structure around TeV. The two annihilation channels are mediated by two gauge bosons that interact both with dark matter and with the standard model fermions. Dark matter annihilations through the s-channel process mediated by the heavier boson produce monoenergetic electron-positron pairs leading to the resonance excess. The lighter boson has a mass smaller than the dark matter such that they can be on-shell produced in dark matter annihilations in the galaxy; the lighter bosons in the final state subsequently decay to generate the extended excess events due to the smeared electron energy spectrum in this process. We further analyze constraints from various experiments, including HESS, Fermi, AMS, and LHC, to the parameter space of the model where both excess events can be accounted for. In order to interpret the two new features in the DAMPE data, dark matter annihilation cross sections in the current galaxy are typically much larger than the canonical thermal cross section needed for the correct dark matter relic abundance. This discrepancy, however, is remedied by the nonperturbative Sommerfeld enhancement because of the existence of a lighter mediator in the model.


Introduction
Recently, DAMPE collaboration reported new measurements of the cosmic electron flux, which exhibit two exotic features in the energy spectrum, including a so-called break structure around 0.9 TeV and a sharp resonance near 1.4 TeV [1]. The morphology of the excess electron events near 1.4 TeV in the DAMPE data hints a nearby cosmic ray source; a number of papers have appeared to interpret the DAMPE narrow resonance near 1.4 TeV, including astrophysical sources [2][3][4][5] and dark matter (DM) sources . So far, most papers interpret only the 1.4 TeV excess as signals of a local source. In this analysis, we argue that both the 0.9 TeV break and the 1.4 TeV resonance in the DAMPE electron spectrum could have a common origin.
We propose a dark matter explanation to both the break and the resonance in the DAMPE electron data. A dark matter subhalo (SH) is assumed to exist in the vicinity of the solar system motivated by the morphology of the 1.4 TeV resonance excess. We assume a simple cosmic electron background, a single power-law form with only two parameters. Thus, both the break and the sharp resonance in the DAMPE data have emerged as results of excess electrons in dark matter annihilations in the nearby subhalo.
We propose a two-mediator dark matter model (2MDM) to interpret the excess electrons seen by DAMPE. In our model, dark matter can annihilate in the galaxy via two different annihilation channels due to the two mediators that interact both with dark matter and with the standard model (SM) sector. The two annihilation channels produce distinct signatures in cosmic electron flux because of different mass hierarchies between DM and mediators. One of the two mediators, denoted as V 1 , has a mass nearly twice of the DM particle; thus the dominated annihilation channel mediated by the V 1 boson is the χχ → V 1 → e − e + process which produces cosmic electrons (and positrons) with energy equal to the DM mass. This then leads to a sharp resonance in the energy spectrum, when DM annihilates in a nearby subhalo. Because the sharp resonance in the DAMPE data occurs around 1.4 TeV, we take the DM mass to be 1.5 TeV.
The other mediator, denoted by V 2 , is much lighter than DM such that the pairproduction of on-shell V 2 bosons (χχ → V 2 V 2 ) becomes the primary DM annihilation channel among the processes mediated by the V 2 boson. The lighter mediator V 2 in the final state further decays to produce SM fermions. If the V 2 boson can directly decay into a pair of electron and positron, the electrons (and positrons) have a box-shape energy spectrum which is centered at one-half of the DM mass and has a width determined by the mass ratio between V 2 and DM. The box-shape electron energy spectrum is further altered during the propagation between the source (the DM subhalo) and the observation point (the DAMPE satellite) to generate an extended excess in the electron energy spectrum. This then gives rise to a "break" structure roughly at one-half of the DM mass (∼ 750 GeV in our case) in the electron energy spectrum observed by DAMPE. In our model, because the V 2 boson is L µ − L τ gauged so that the electrons originating from V 2 decays only carry a fraction of the total energy, the electron energy spectrum is further smeared.
Thus, both the break and the sharp resonance observed by DAMPE arise due to the electrons coming from DM annihilations in the 2MDM model. We begin in section 2 by presenting the two-mediator DM models in which DM interacts with SM via two different gauge bosons. In section 3, we provide the cosmic electron background used in the analysis. In section 4, we describe the method to compute electron flux from DM annihilations both in the DM subhalo and in the Milky Way (MW) halo, as well as the method to compare our calculations with DAMPE data. In section 5, we compute the DAMPE signals expected in the two-mediator dark matter models. In section 6, we analyze HESS constraints on our DM model. In section 7, we study the constraints from the Fermi isotropic gamma ray background measurements. In section 8, we investigate AMS constraints on our DM model. In section 9, we calculate ATLAS constraints on our DM model. In section 10, we discuss Sommerfeld enhancement in our model and the impacts on DM relic abundance. In section 11, we summarize our findings.

The two-mediator dark matter model
We consider DM models in which DM is a Dirac fermion and charged under two U (1) fields; the corresponding gauge bosons are V 1 and V 2 . Hinted by the sharp resonance near 1.4 TeV in the DAMPE data, we fix the DM mass at 1.5 TeV. The mass of the V 1 boson can be near 3 TeV; The mass of the V 2 is lighter than the DM mass such that V 2 can be on-shell produced in DM annihilations in the DM halo. Thus the relevant DM annihilation channels are The Feynman diagrams for the two annihilation channels are shown in Fig. (1). In our analysis, we consider two cases: (1) V 1 is electrophilic and V 2 is the L µ − L τ gauged; (2) V 1 is hidden and kinetically mixes with the SM hypercharge and V 2 is L µ − L τ gauged. Below we provide one concrete model in the latter case where a Dirac fermion dark matter particle χ couples to two spin-1 mediators V 1 and V 2 : where is the kinetic mixing parameter between gauge boson V 1µ and the SM

Cosmic electron background
Because DAMPE does not distinguish the charge of electron/positron events, we assume the following single power-law background (BG) for the total flux of electron and positron, where C and γ are free parameters to be determined by data. In our analysis, we use the first eight points and the last eight points in the DAMPE data [1] to fit this single power-law background, and obtain the following best-fit parameters: C = 458 (GeV m 2 s sr) −1 and γ= 3.25. We use these background parameters throughout our study. Since DAMPE is unable to discriminate electrons from positrons, we will use the word "electron" in this paper to collectively denote both electron and positron when there is no confusion.

Electron flux from DM annihilations
The sharp resonance of the excess events near 1.4 TeV in the DAMPE data hits a nearby electron/positron source. To fit the spectrum of the DAMPE data, we introduce a nearby DM subhalo with an NFW density profile [40] ρ(r) = ρ s (r/r s ) −γ (1 + r/r s ) 3−γ . (4.1) We use the following parameters for the subhalo: γ = 0.5, ρ s = 100 GeV/cm 3 , and r s = 0.1 kpc [6]. The distance between the subhalo and us (denoted by d s ) is also crucial to the cosmic ray spectrum. We find that the above subhalo with d s = 0.3 kpc can fit the DAMPE data well. The above values of the four parameters are assumed for the subhalo in our analysis if not specified otherwise. The electron/positron flux can originate from dark matter annihilations both in the Milky Way dark matter halo and in a nearby subhalo. The electron flux from DM annihilations in the MW halo (denoted by Φ χ−MW ) is computed via PPPC4DMID [41]. For the electron flux arising from DM annihilations in a nearby subhalo, we use the Green's function method [6,[42][43][44] where G(x, E; x s , E s ) is the Green's function, Q(x s , E s ) is the source term due to DM annihilation, v e is the electron velocity, and Φ χ−SH (x, E) is the electron flux due to DM annihilations in the subhalo, which has the unit of (GeV −1 m −2 s −1 sr −1 ). Here the subscript s indicates the quantities associated with the DM source. The Green's function can be calculated via GeV/s and the diffusion coefficient D(E) = D 0 (E/GeV) δ with D 0 = 11 pc 2 /kyr, and δ = 0.7 [45]. The source function due to DM annihilations is where m χ is the DM mass, ρ χ (x s ) is the DM mass density, σv is the velocity-averaged DM annihilation cross section, dN/dE s is the electron energy spectrum per DM annihilation. Thus, in our analysis, the total electron flux is given by Φ th = Φ BG + Φ χ−MW + Φ χ−SH where we consider three major contributions: the cosmic ray background, DM annihilations in the MW halo, and DM annihilations in the nearby subhalo. To compare our calculations with the DAMPE data, we further take into account the "bin effects" by performing the following computation In this section, we compute the electron flux expected in DAMPE from DM annihilations in the two-mediator DM model. DM annihilations both in the subhalo and in the MW halo are considered in our analysis. We use σv 1 to denote the velocity averaged DM annihilation cross section for the χχ → V 1 → ff process, which is mediated by the heavier gauge boson V 1 ; we use σv 2 to denote the velocity averaged DM annihilation cross section for the χχ → V 2 V 2 process where the lighter gauge boson V 2 is on-shell produced. The annihilation cross sections σv 1 and σv 2 are mainly responsible for the resonance and the break excess events in the DAMPE electron spectrum respectively.

Electrophilic and gauged L µ − L τ
Here we first consider the two-mediator model in which the heavier mediator V 1 is electrophilic and the lighter mediator V 2 is the L µ − L τ gauge boson. In this case, for the annihilation process mediated by the V 1 , only χχ → V 1 → e + e − can occur. In the annihilation processes where V 2 is on-shell produced, V 2 further decays into µµ, τ τ , and νν final states with branching ratio BR= 1/3 for each final state. The energy spectrum of the µµ and τ τ final states exhibits a box-like distribution that is centered at m χ /2 (see e.g. [46][47][48][49][50][51] for early studies in the context of cosmic rays). In our analysis, because we assume a simple power-law background, there is a wide range of electron excess events, extending from about 50 GeV to almost over 1 TeV, as shown in the left panel figure of Fig. (2). To generate such an extended electron excess events, the mass of the V 2 boson has to be sufficiently small since the width of the box-shape energy spectrum is given by In addition, the V 2 boson in our study is also required to decay into the τ τ final state. Thus we take 10 GeV as the benchmark point for the V 2 boson mass, which is assumed throughout our analysis.
The  comes from the χχ → V 1 → e − e + annihilation channel; the 0.9 TeV break shown in Fig. (2) is primarily due to the χχ → V 2 V 2 annihilation channel. The formation of the break structure in the electron energy spectrum here is due to several aspects of the problem here. Because the mass of V 2 in our analysis is taken to be 10 GeV, the box-shape energy spectrum of µµ and τ τ have a wide energy range (extending almost from zero to 2m χ ). The energy loss in charged lepton decays and cosmic ray propagation in addition shapes the excess electrons. Thus one obtains an extended distribution of the excess electrons with a power-law break around TeV.  Table 1. Best-fitted cross sections of the χχ → e + e − and χχ → V 2 V 2 processes for different d s values. Here V 2 is the L µ − L τ gauge boson.
We further vary the distance between the subhalo and us on the right panel figure of Fig. (2), while keeping the rest of the parameters fixed for the subhalo. For each d s value, we find the best-fit DM annihilation cross sections for both channels in fitting the DAMPE data, which are shown in Table (1). We will adopt the case in which d s = 0.3 kpc as the benchmark model for our analysis, in which the DM annihilation cross sections for the two channels take σv 1 = 4.9 × 10 −26 cm 3 /s and σv 2 = 2.0 × 10 −24 cm 3 /s. Taking into account the Sommerfeld enhancement effects (see section (10) for the detailed discussions), we find that one should have g 2 = 0.68 in order to obtain σv(χχ → V 2 V 2 ) = 2.0 × 10 −24 cm 3 /s.

Kinetic mixing and gauged
Here we consider the two-mediator model in which the heavier mediator V 1 mixes with the SM hypercharge gauge boson via the kinetic mixing (KM) term and the lighter mediator The Lagrangian of this model is given in Eq. (2.3). Unlike the previous case in which the V 1 is electrophilic, the annihilation process mediated by the heavier mediator χχ → V 1 → ff now produces all SM fermions. The analysis regarding the V 2 boson is similar to that in the previous section.
Regarding the V 1 boson in the KM case, there usually are four free parameters in the calculation: the KM parameter , the gauge coupling g 1 , the DM mass m χ which is now fixed at 1.5 TeV, and the mediator mass m V 1 which is typically near 2m χ to provide a sufficient annihilation rate for the DM relic abundance. However, in our case, there is another lighter mediator which can significantly change the annihilation cross section mediated by V 1 boson, via the Sommerfeld enhancement mechanism. Thus, to correctly compute the DM annihilation cross section in the halo for the χχ → V 1 → ff process, one has to multiply the annihilation cross section due to V 1 (see e.g. Ref. [49]) and the Sommerfeld enhancement factor due to the lighter V 2 mediator (see section (10) for the detailed discussions) in the model.
In Fig. (3), we compute the electron flux arising from the 2MDM model in which the V 1 boson kinetically mixes with the SM hypercharge gauge boson, and the V 2 is the L µ − L τ gauged boson. Here the heavier V 1 boson can decay into various SM fermions and the branching ratios are determined primarily by the hypercharge quantum numbers of the SM fermions. Since the right-handed charge lepton has a relatively large hypercharge, the total branching ratio of the V 1 boson into the three generation charge leptons is rather large, =e,µ,τ BR(V 1 → + − ) 37%. We find that the DM annihilation cross section σv 1 = 3.9× 10 −25 cm 3 /s (for all SM final states in the χχ → V 1 → ff process), and σv 2 = 2.0× 10 −24 cm 3 /s (for the χχ → V 2 V 2 process) provide the best fit to the DAMPE data in this model. In Fig. (3), the peak comes from the contributions of the χχ → V 1 → ff processes (mainly due to the e + e − final state), whereas the break are primarily due to processes mediated by the L µ − L τ boson.  Figure 3. DAMPE electron energy spectrum. Overlaid is the sum of cosmic ray background and the electron flux from DM annihilations from the MW halo and from the nearby subhalo. Here V 1 kinetically mixes with the SM hypercharge boson and V 2 is the gauge L µ − L τ boson. The DM annihilation cross sections σv 1 = 3.9 × 10 −25 cm 3 /s (for all SM final states in the χχ → V 1 → ff process) and σv 2 = 2.0 × 10 −24 cm 3 /s (for the χχ → V 2 V 2 process) are used here.
We note that the best-fit cross section for σv 2 is about the same as in the previous model so that g 2 0.68. Taking into account the Sommerfeld enhancement factor, we find the model point ( , g 1 , m χ , M V 1 ) = (0.01, 0.1, 1500 GeV, 2994.2 GeV) in the parameter space can give rise to σv 1 = 3.9 × 10 −25 cm 3 /s. Here the mass of the V 1 boson is smaller than 2m χ so that the invisible decay V 1 → χχ cannot occur. In addition, the DM annihilation cross section at the early universe receives another suppression factor relative to that in the DM halo today, because the larger kinetic energy of the DM particles at the early universe moves the characteristic √ s of the DM annihilation process further away from the Breit-Wigner resonance relative to today [52]. Because the invisible decay of the V 1 boson is kinetically disallowed here and the branching ratios into charged leptons are rather significant, the discovery potential of LHC for such V 1 boson is high. The discussions on LHC constraints on this model are given in section 9.

HESS constraints
The gamma ray flux produced by DM annihilations can be calculated as follows where m χ is the DM mass, σv i is the velocity-averaged DM annihilation cross section for channel i, (dN γ /dE γ ) i is the gamma ray energy spectrum per annihilation for channel i, and J(∆Ω) is the J-factor for the region-of-interest (ROI). The differential flux dΦ γ /dE γ has unit of (GeV cm 2 s) −1 . The J-factor is computed via where ∆Ω is the solid angle of the ROI, ρ χ is the DM density, and s is the distance along the line of light. HESS searched for very high energy γ-rays in the inner region of the Milky Way halo, which is a circular region of 1 • radius excluding a ±0.3 • band in the Galactic latitude [53] [54]. With the 254-hour data accumulated [54], stringent upper bounds can be set on the DM annihilation cross sections for various SM final states. In a recent study [55], the HESS constraints on dark matter annihilations into on-shell mediators for various SM final states are analyzed. Our analysis here is similar to that in Ref. [55], but in our case, on-shell mediators annihilate into a collection of SM final states with branching ratios given in the 2MDM model.

HESS constraints on DM annihilations in the Galactic center
In the following, we calculate the upper limit on the DM annihilation cross section σv χχ→V 2 V 2 from the HESS data, where V 2 is the gauge L µ − L τ boson. The method we use here is to rescale the limits calculated in Ref. [54] which analyzed 254-hour data recorded by HESS. The details of the method can be found in Appendix (B).
We first analyze the HESS limits on DM annihilations in the center of the galaxy. Because the DM distribution is not known to a good precision in the center of the galaxy and Here only gamma rays from the MW halo are considered. The limits are computed based on the exclusion limits for the χχ → µ + µ − and χχ → τ + τ − processes given in Ref. [54]. A light V 2 mass is assumed in this analysis.
the gamma rays are very sensitive to the DM density distributions in the Galactic center, several DM profiles are considered in the HESS analyses [53] [54]. We provide a comparison of the J-factors from different DM profiles in Appendix (A). Here we consider three different DM profiles, NFW, Isothermal, and Einasto, to interpret the HESS constraints. Fig. (4) shows the 95% CL limits on DM annihilation cross section σv χχ→V 2 V 2 where V 2 is the gauged L µ − L τ boson. For the 1.5 TeV DM annihilating into sufficiently light V 2 bosons, the HESS constraints are σv 2 1.1 × 10 −25 (4 × 10 −24 ) cm 3 /s for the NFW (Isothermal) profile. Thus the DM annihilation cross section σv 2 = 2.0 × 10 −24 cm 3 /s which is responsible for generating the break in the DAMPE data, is excluded if one considers the NFW or Einasto profile, but is still allowed if the isothermal profile is assumed.
6.2 HESS constraints on the location of the subhalo DM annihilations in the subhalo also contribute to the gamma ray flux observed by the HESS experiment. Because the HESS search region is 1 • around Galactic center, the gamma ray flux observed by HESS from the subhalo is a function of l SH that is the angle between the Galactic center and the center of the subhalo. We compute the J factor of the subhalo inside the HESS search region in the left panel figure of Fig. 5 for different d s values. The subhalo J factor increases when the subhalo moves towards either the Galactic center or us.
We further determine the minimum l SH value by saturating the HESS constraints on DM annihilations. To determine the minimum l SH value, we use σv DAMPE × (J iso MW + J SH (l min SH )) = σv HESS ×J iso MW , where σv DAMPE is the cross section needed for the DAMPE electron excess events, as given in Table (1), J iso MW is the J factor inside the HESS search region for the MW halo with the isothermal DM density profile which is 7.23 × 10 19 GeV 2 cm −5 , σv HESS is the HESS 95% CL upper bound on the DM annihilation cross section with the isothermal profile (which is 4 × 10 −24 cm 3 /s as given by the isothermal curve on Fig. (4)), J SH is the J factor inside the HESS search region for the subhalo. Because the gamma ray flux produced by the process χχ → V 1 → e + e − is much smaller than χχ → V 2 V 2 , we take σv DAMPE σv (χχ → V 2 V 2 ) in the calculation here. The right panel figure of Fig. 5 shows the lower bound on the l SH angle. When d s = 0.3 kpc, the subhalo has to be > 21 • away from the Galactic center to avoid HESS constraints.

HESS limits for both DM annihilation channels
Here we analyze the HESS constraints for the model in which V 1 kinetically mixes with the SM hypercharge and V 2 is L µ − L τ gauged. To take both channels into consideration, we use Φ γ ( σv 1 , m χ ) + Φ γ ( σv 2 , m χ ) = Φ 95 γ (m χ ), where Φ 95 γ is the 95% CL upper bound from the 254-h HESS data on the total gamma ray flux (in unit of cm −2 s −1 ) integrated over the energy range 160 GeV < E γ < m χ . Here Φ γ ( σv 1 , m χ ) and Φ γ ( σv 2 , m χ ) are the gamma rays from the two annihilation channels respectively. Fig. (6) shows the HESS limits for both σv 1 and σv 2 for the case where m χ = 1.5 TeV, where only contributions from the MW halo are considered. Our model is excluded if the NFW profile is used, but allowed if the isothermal profile is used for the DM distribution in the MW halo.

Fermi constraints
Similar to the gamma ray flux measured by HESS, the gamma ray flux observed by Fermi due to DM annihilations is calculated as follows, the Fermi isotropic gamma ray background region can thus be computed as follows where ρ χ is the DM density, b is the galactic latitude, is the galactic longitude, s is the distance between the point where DM annihilates and us. In this study, we take into account both the MW halo and the DM subhalo when calculating the J-factor. In this section, we consider the same isothermal DM profile for the MW halo as in the HESS analysis.

Fermi isotropic gamma ray background constraints
Here we compare the gamma ray flux produced by dark matter annihilations in the subhalo as well as in the MW halo, with the isotropic background measured by Fermi-LAT [56] to obtain constraints on our DM model. Because the galactic plane is masked in the Fermi IGRB analysis [56], the constraints from Fermi IGRB are minimized when the subhalo sits on the galactic plane. We use b SH to denote the galactic latitude of the subhalo center. Thus we will set b SH = 0 for our analysis unless specified otherwise. The left panel figure in Fig. (7) shows the Fermi IGRB data [56] and the gamma rays from DM annihilations for the case in which the heavier V 1 boson is electrophilic and the lighter V 2 boson is L µ − L τ gauged. The DM annihilation cross sections for the two annihilation channels are ( σv 1 , σv 2 ) = (4.9 × 10 −26 cm 3 /s, 2.0 × 10 −24 cm 3 /s), which are the same as those in the left panel figure of Fig. (2). Here the gamma ray flux arising from the χχ → e + e − process is only about 3% of that due to χχ → V 2 V 2 in this case. We find that the J-factor of the subhalo is about the same as the J-factor of the MW halo in the Fermi IGRB search region, J SH J M W 6 × 10 21 GeV 2 /cm 5 . We have plotted the gamma rays from the MW halo on the left panel figure of Fig. (7), as well as the gamma rays from both the MW halo and the subhalo. The predicted total gamma rays in our DM model do not exceed the current Fermi IGRB bound. Eγ (GeV) Figure 7. Left panel: Fermi IGRB data [56] and gamma rays from DM where V 1 is electrophilic and V 2 is L µ − L τ gauged. Right panel: Fermi IGRB data and gamma rays from DM where V 1 is kinetically mixed with SM hypercharge and V 2 is L µ − L τ gauged and has a mass of 10 GeV. DM annihilation cross sections are the same as those in Fig. (3). We use the default subhalo parameters for both figures. For the DM model in which the heavier V 1 boson kinetically mixes with the SM hypercharge gauge boson and the lighter V 2 boson is L µ − L τ gauged, the predicted gamma rays are shown on the right panel figure in Fig. (7). We use the following DM annihilation cross sections ( σv 1 , σv 2 ) = (3.9 × 10 −25 cm 3 /s, 2.0 × 10 −24 cm 3 /s) which are the same as those in Fig. (3). Unlike the DM model presented on the left panel figure of Fig. (3), the annihilation process mediated by the V 1 boson on the right panel figure of Fig. (3) has a larger cross section and various SM final states. We plotted the gamma rays from both annihilation channels on the right panel figure of Fig. (3). We find that the isotropic gamma ray measurements are beginning to probe this DM model at the high energy bins in the Fermi IGRB data.

Fermi constraints on the subhalo
Here we study the effects on the Fermi IGRB data by changing various parameters for the DM subhalo. The gamma ray flux is very sensitive to the distance between the subhalo and us. We compute the gamma rays expected at Fermi using different d s values on the left panel figure of Fig. (8). Different d s values not only lead to different J-factors in the Fermi search region, but also lead to different DM annihilation cross sections which are provided in Table (1), since one has to fit the DAMPE data. The predicted gamma rays become larger when the subhalo moves towards us. In order to evade the Fermi IGRB constraints, the subhalo has to be at least 0.3 kpc away from us. We also compute the gamma rays from the subhalo when it moves away from the Galactic plane. The gamma ray flux expected in Fermi is shown on the right panel figure of Fig. (8) for several different b SH values. If the subhalo moves away from the Galactic plane for more than 10 • , the gamma rays produced in the Fermi IGRB search region become significant above the current measurements.
We further study the gamma rays by changing the subhalo profile parameters (r s , ρ s ), in the Fig. (9) where d s = 0.3 kpc and γ = 0.5 are fixed. Two sets of parameters in addition to the default values for the subhalo are used here. For each case, the DM annihilation  Table (1). cross sections for the two different channels are chosen such that one obtains the least χ 2 fit to the DAMPE data. As shown in Fig. (9), the Fermi constraints can be significantly alleviated if the DM subhalo becomes smaller and denser.

AMS constraints
We do not attempt to explain the AMS positron excess. However, the two-mediator DM model cannot produce too many positrons so that they violate the AMS data on the positron fraction measurement [57].
To compute the AMS constraints on the DM model, we extrapolate our simple cosmic ray electron/positron background given by Eq. (3.1) down to low electron energy range. We further assume that the background of the positron fraction take the following simple . We use first 15 data points in the AMS positron fraction data [57] to find the best-fit parameters: C f = 11.2 and γ f = 0.31. The positron fraction including contributions both from the background and from DM annihilations is thus computed by where Φ χ is the cosmic flux including both electron and positron due to DM annihilations.
We use (f AMS  (10) shows the AMS constraints on the DM annihilation cross section mediated by the L µ − L τ gauge boson using the positron fraction data. The most stringent limit comes from the highest energy bin in the AMS data, which provides the 95% CL upper bound as σv 2 3 × 10 −24 cm 3 /s for the L µ − L τ gauge boson. The predicted positron fraction values at the AMS energy range in our model lie below the AMS measurements. We note in passing that the gap between our predicted positron fraction and the actual AMS data could be due to astrophysical sources.

LHC constraints
Here we study the LHC constraints on the V 1 boson that is kinetically mixed with the SM hypercharge. In this case, the V 1 boson couples to all SM fermions due to the kinetic mixing parameter , which is given in Eq. (2.3). Thus, the V 1 boson can be produced in the Drell-Yan process at the LHC and can be searched for by reconstructing the dilepton final states. Here we utilize the recent ATLAS data [58] to put constraints on the kinetic mixing parameter between V 1 boson and the SM hypercharge boson. Fig. (11) shows the ATLAS upper bound on the dilepton production cross section, using 36.1 fb −1 data at the 13 TeV colliding energy. Predicted dilepton signals arising  Fig. (11). The limit on will certainly improve when all data currently accumulated at the LHC are analyzed (about 150 fb −1 data have been collected by ATLAS and by CMS individually so far [59]). However, to reach the sensitivity of probing the model point considered in our analysis, = 0.01 for a 3 TeV Z boson, more data in future LHC runs are probably needed.

Sommerfeld enhancement
The cross section of the process χχ → V 2 V 2 is larger than the canonical thermal DM annihilation cross section by about two orders of magnitude, which would suppress the DM abundance significantly. However, we should take into account the Sommerfeld enhancement induced via V 2 exchanges between DM particles in the annihilation processes as illustrated in Fig. (12), since the mediator V 2 is light and the velocity of DM is low in the MW halo.  The Sommerfeld enhancement factor S can be approximated by [60] [61] [51] S = π v sinh X cosh X − cos (2π/¯ 2 ) − X 2 , (10.1) where¯ 2 = (π/12) 2 and X = v /¯ 2 , and 2 = m V 2 /(α 2 m χ ), v = v/α 2 with α 2 = g 2 2 /(4π). We take v = 10 −3 as the typical DM velocity in the halo. The left panel figure of Fig. (13) shows the Sommerfeld enhancement factor as a function of the gauge coupling g 2 where the mediator V 2 mass is 10 GeV and the dark matter mass is 1.5 TeV. For the process χχ → V 2 V 2 , the DM annihilation cross section is given by σv 2 = σv 0 2 × S(g 2 ) where σv 0 2 g 4 2 /(16πm 2 χ ) is the annihilation cross section without taking account the Sommerfeld enhancement effect. By equating this expression with 2.0 × 10 −24 cm 3 /s, the needed cross section to fit DAMPE, one obtains g 2 = 0.68. Thus one obtains the corresponding Sommerfeld enhancement factor S 93. We further plot σv 2 as a function of coupling g 2 in the right panel figure of Fig. (13).
For the process χχ → V 1 → ff , one also has to consider the same enhancement due to the V 2 mediator, so that the DM annihilation cross section should be computed via σv 1 = S × σv 0 1 , where the superscript 0 indicates the cross section without taking the Sommerfeld enhancement into account. Using S 93, we find that the model point ( , g 1 , m χ , M V 1 ) = (0.01, 0.1, 1500 GeV, 2994.2 GeV) in the parameter space of the KM model can give rise to σv 1 = 3.9 × 10 −25 cm 3 /s which is needed to fit the DAMPE data.
The DM relic abundance which is primarily determined by the DM annihilation cross section at the so-called freeze-out epoch at the early universe. Typical freeze-out occurs at the temperature T m χ /(20 − 25) such that the DM velocity is approximately v 1/4, where DM annihilation cross section no longer receives significant Sommerfeld enhancement that is present at the current galaxy. We compute the DM annihilation cross section for the processes χχ → V 1 → ff (KM) and χχ → V 2 V 2 at the freeze-out and find that σv 1 = 1.0 × 10 −28 cm 3 /s, and σv 2 = 2.2 × 10 −26 cm 3 /s, when T = m χ /25. Thus the total DM annihilation cross section is approximately 2.2 × 10 −26 cm 3 /s at freeze-out which is very close to the canonical thermal DM annihilation cross section needed to generate the right DM relic density in the universe. We note that there is a three orders of magnitude boost on σv 1 at current galaxy relative to the early universe, owing to both the Breit-Wigner enhancement and the Sommerfeld enhancement in this annihilation channel.

Conclusions
There are two exotic features present in the new cosmic electron spectrum observed by the DAMPE collaboration, including a break at 0.9 TeV and a peak at 1.4 TeV. We propose to simultaneously explain both features in the DAMPE data via annihilations from one DM species that interacts with SM via two different mediators. Thus two different DM annihilations channels via the two mediators generate the two new features in the cosmic electron energy spectrum near TeV. The annihilation process mediated by the heavier V 1 boson generates the 1.4 TeV peak; the annihilation process mediated by the lighter V 2 boson produces the extended break near 0.9 TeV.
In this work, we consider two concrete examples of the two-mediator DM models. In both cases the lighter V 2 boson is L µ − L τ gauged and has mass 10 GeV such that V 2 can be on-shell produced in annihilations of DM which is taken to be 1.5 TeV. We consider the heavier V 1 boson to be either electrophilic or kinetically mixed with the SM hypercharge.
We assume a single power-law cosmic electron background which contains only two parameters and a DM subhalo which is 0.3 kpc from us. Both electrophilic and KM V 1 bosons provide good fits to the 1.4 TeV excess, with the annihilation cross section 4.9×10 −26 cm 3 /s and 3.9 × 10 −25 cm 3 /s respectively; the L µ − L τ gauge boson V 2 provides a good fit to the break with the annihilation cross section 2.0 × 10 −24 cm 3 /s. Several experimental constraints on the DM models are analyzed, including HESS, Fermi IGBG, AMS positron fraction and LHC dilepton searches. Gamma rays expected at the HESS search region are mainly coming from annihilations via the V 2 boson due to the larger cross section. HESS constraints are very sensitive to the DM density profile for the MW halo. The needed cross section for the V 2 process is excluded if one assumes the NFW or Einasto profile for the MW halo, but still allowed if the isothermal profile is considered. In addition, a substantial amount of gamma rays also arise in DM annihilations via the kinetic-mixing V 1 boson; we find that the gamma rays from both annihilation channels are consistent with HESS data assuming the isothermal profile for the MW halo. We also find that the subhalo cannot be put at the Galactic center direction since it would contribute a significant amount of gamma rays to the HESS search region. Fermi isotropic gamma ray background constraints are sensitive to the distance between the subhalo and us. We find that our models do not violate the Fermi isotropic gamma ray background if the subhalo is placed at 0.3 kpc from us. We also note that one can begin to probe our model with more data accumulated at Fermi. DM annihilations in our model cannot provide satisfactory explanations to the AMS positron fraction excess. Nonetheless, one can use the AMS data to put the constraints on DM models by demanding that the predicted positron fraction in DM models not exceed the AMS measurement. We find that the highest energy bin in the AMS data gives the most stringent bound on the L µ − L τ gauge boson process, and will probe our model in the near future. LHC constraints on the KM V 1 boson are analyzed in the dilepton channel. For a 3 TeV V 1 boson, the upper bound on is about 0.1.
The DM annihilation cross sections needed to fit DAMPE data are much larger than the canonical thermal cross section. This discrepancy can be nicely explained by the Sommerfeld enhancement due to the light V 2 mediator in the models. Taking into account the non-perturbative Sommerfeld enhancement corrections present in the current galaxy, our model is consistent with the relic density requirement in the thermal DM framework.

B Rescaling method for HESS limits
Here we describe the "rescaling" method used in our analysis to obtain the HESS limits. We first digitize the 95% CL upper bound on DM annihilation cross section σv 95 (m χ ) from the HESS 254-h data analysis [54] for one specific annihilation channel, for instance the χχ → µ + µ − annihilation channel. The 95% CL upper bound on the total gamma ray flux Φ 95 γ (m χ ) can be obtained by integrating the differential flux given in Eq. (7.1) over the gamma ray energy range 160 GeV < E γ < m χ , where we used the Einasto J-factor J E given in Ref. [54]. The HESS constraints can then be calculated for different DM annihilation channels by integrating the differential flux given in Eq. (7.1) taking into account different halo profiles and different gamma ray energy spectra. Fig. (14) shows the HESS 95% CL bound on Φ 95 γ (m χ ) for the χχ → µ + µ − and χχ → τ + τ − annihilation channels. For the on-shell produced vector bosons in the L µ − L τ model, we use the average value of Φ 95 γ (m χ ) of the χχ → µ + µ − and χχ → τ + τ − annihilation channels to compute the HESS limits.  Figure 14. HESS 95% CL upper bound on the total gamma ray flux Φ 95 γ in the energy range 160 GeV < E γ < m χ [54]. The 2µ (2τ ) curve is obtained based on the HESS upper bound on the annihilation cross section curve in the χχ → µ + µ − (χχ → τ + τ − ) channel [54]. The "average" curve is the arithmetic mean of the 2µ and 2τ curves.