Explaining the DAMPE data with scalar dark matter and gauged $U(1)_{L_e-L_\mu}$ interaction

Inspired by the peak structure observed by recent DAMPE experiment in $e^+e^-$ cosmic-ray spectrum, we consider a scalar dark matter (DM) model with gauged $U(1)_{L_e-L_\mu}$ symmetry, which is the most economical anomaly-free theory to potentially explain the peak by DM annihilation in nearby subhalo. We utilize the process $\chi \chi \to Z^\prime Z^\prime \to l \bar{l} l^\prime \bar{l}^\prime$, where $\chi$, $Z^\prime$, $l^{(\prime)}$ denote the scalar DM, the new gauge boson and $l^{(\prime)} =e, \mu$, respectively, to generate the $e^+e^-$ spectrum. By fitting the predicted spectrum to the experimental data, we obtain the favored DM mass range $m_\chi \simeq 3060^{+80}_{-100} \, {\rm GeV}$ and $\Delta m \equiv m_\chi - m_{Z^\prime} \lesssim 14 \, {\rm GeV}$ at $68\%$ Confidence Level (C.L.). Furthermore, we determine the parameter space of the model which can explain the peak and meanwhile satisfy the constraints from DM relic abundance, DM direct detection and the collider bounds. We conclude that the model we consider can account for the peak, although there exists a tension with the constraints from the LEP-II bound on $m_{Z^\prime}$ arising from the cross section measurement of $e^+e^- \to Z^{\prime\ast} \to e^+ e^-$.


Introduction
One of the most important questions in current particle physics and cosmology is to understand the nature of cosmic dark matter (DM). Although many popular theories can predict viable DM candidates, no DM particle has been discovered by collider or direct detection (DD) experiments so far. So alternative ways, e.g. indirect detection of DM by seeking for its annihilation or decay products, become valuable in understanding the nature of DM.
Recently, the DArk Matter Particle Explorer (DAMPE) experiment released the new measurement of the total cosmic e + + e − flux between 25 GeV and 4.6 TeV and reported a hint of an excess in the e + e − spectrum at around 1.4 TeV [1,2]. Although such an excess may originate from certain new unobserved astrophysical sources, it may also be explained by DM annihilation [3]. Relevant discussion on this subject can be found in [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20], and one notable conclusion is that if one assumes the e + e − cosmic-ray spectrum to be generated directly from DM annihilation in a nearby clump halo, the best fit values for the DM particle mass, the DM clump mass and the annihilation luminosity are around 1.5 TeV, 10 6−8 M and 10 64−66 GeV 2 cm −3 , respectively, if the subhalo is about 0.1 ∼ 0.3 kpc away from the earth [3].
In this work, we attempt to construct a new theory to explain the excess by DM annihilation. For this end, one preliminary requirement on the theory is that the DM annihilation product should be rich in e + e − states. Other requirements include • I-ID: current DM annihilation cross section σv 0 1 × 10 −26 cm 3 /s with v ∼ 10 −3 c in the halo, and meanwhile satisfying the constraints from other DM indirect search experiments.
• IV-Collider: collider constraints from LHC and LEP-II measurements.
The requirements I-ID and II-RD prefer the annihilation process dominated by s-wave and at the same time without chiral suppression by light fermion masses [25,26]. Therefore a natural realization to explain the DAMPE excess is the Dirac DM scenario with certain lepton-specific gauge symmetry, where the DM annihilates directly to e + e − through schannel mediation of the new gauge boson Z [4,10,12,14,15]. However, as indicated by recent studies in [4,12,14,15], even though Z couples only with leptons, DM-nucleon scattering may still proceed by the t-change exchange of the Z − γ transition which is induced by lepton loops. As a result, DM direct detection bounds have tightly limited the parameter space into the resonant annihilation region, m Z 2m DM , so that it becomes less attractive, especially given the fact that the bound will be improved greatly in near future. On the other hand, the scenario with scalar DM pair annihilation into vector bosons χχ → Z Z → 2(e + e − ) can also satisfy the requirements I-ID and II-RD. In our previous work [13] we studied the scenario with the gauge boson as the mediator between the DM sector and SM sector corresponding to a gauged family symmetry U (1) B−3(Le−Lu−Lτ ) . We found that the minimal model can not explain the DAMPE results due to the tight constraint from the LHC search for new gauge boson, but a slightly extended model can do this without conflicting with the constraint [13].
In principle, the scenario of the scalar DM annihilation with vector portal may be transported to other gauged family symmetries. In this direction, we are particularly interested in the lepton-specific family symmetry U (1) Le−Lµ [27,28] due to the following considerations. Firstly, the model is the most economical anomaly-free theory which may explain the peak, so its capability in this aspect should not be neglected. Moreover, we note that because the gauge boson in the model decays democratically into e + e − and µ + µ − , instead of into single e + e − discussed in [15] or equally into three generation lepton pairs discussed in [13], its prediction on the e + e − spectrum should be different from the previous studies. Secondly, since quarks in the model are uncharged under the U (1) Le−Lµ symmetry, the constraints from LHC experiments and DM direct detection experiments can be greatly relaxed. This facilitates the model in interpreting the DAMPE excess. Finally, we note that the model may address the neutrino mass problem [4] and the discrepancy between theoretical prediction and E821 experimental measurement on muon anomalous magnetic momentum [29]. In this work, we only focus on the features of the model in explaining the DAMPE excess.
The paper is organized as follows. In Section 2, we show the favored values of m χ and m Z by fitting the e + e − spectrum predicted by the DM annihilation χχ → Z Z → lll l to the DAMPE results. In Section 3, we introduce the scalar DM model with U (1) Le−Lµ gauge boson as the portal to leptons. In Section 4, we present the interpretation of the DAMPE excess by the model while satisfying various experimental constraints. Finally, in Section 5 we draw our conclusion.

Favored masses by the DAMPE excess
In this section, we determine the mass ranges of χ and Z by fitting the e + e − spectrum generated by the process χχ → Z Z → lll l to the DAMPE data. For this end, we first adopt the parameterization of the cosmic ray background according to the formulae in [8] and use the package LikeDM [30] to calculate the propagation of the e + e − flux in the background. The needed data to determine the parameters for the background are taken as those of the AMS-02 e + fraction and e + e − flux, and the DAMPE e + e − flux [8]. Then we add the contribution of the local subhalo directly by noting that such a component only affects the energy bin ∼ 1.5 TeV. We adopt a Navarro-Frenk-White profile [31] with a truncation at the tidal radius [32] for DM density distribution inside the subhalo, and use the subhalo mass M halo as an input to determine the profile with the method introduced in the appendix of [3]. The propagation of the nearby e + e − can be calculated analytically under the assumption of spherically symmetric geometry and infinite boundary conditions [33]. More details of the procedure were introduced in [4,10,13,15]. Finally, we construct a likelihood function by comparing the predicted spectrum with the AMS-02 data and the m DM (GeV)  Total fit DM Background AMS data DAMPE data Figure 1. Fitting the e + e − spectrum generated by the process χχ → Z Z with Z → e + e − , µ + µ − to the AMS-02 and DAMPE data with the background parameterization taken from [8]. Left panel: χ 2 map projected on ∆m − m χ plane with ∆m ≡ m χ − m Z and the color bar denoting the χ 2 values. The best fit point locates at (7 GeV, 3060 GeV), and the contour of χ 2 = χ 2 best + 2.3 is also plotted. Right panel: The cosmic e + e − spectrum of the best fit point generated by the DM annihilation process in comparison with the AMS-02 and DAMPE data. DAMPE data, which are distributed from 0.5GeV to 24GeV in 36 energy bins for the former [34] and from 24GeV to 4.57TeV in 40 energy bins for the latter [1]. Obviously, this function depends on a series of parameters such as m χ , m Z , σv 0 , M halo , the distance of the subhalo way from the earth d as well as the parameters appearing in the background. Since DAMPE experiment has collected more than thirty electron/positron events in each energy bin around 1.4TeV [1], we can determine the favored ranges of these parameters by the binned likelihood fit adopted in this work. For the flux in the ith bin, we get its theoretical value by averaging the flux over the width of the energy bin, i.e. Φ i = 1 ∆E Φ(E)dE. Before we proceed, let's illustrate two key features of the spectrum. One is that the height of the spectrum is roughly decided by the product σv 0 × (M halo /m ) 0.76 [35], which implies that one may fix σv 0 or M halo while varying the other one in performing the fit. Since M halo is also needed to determine the DM profile, in practice we fix M halo and vary σv 0 to simplify the calculation. We note that so far M halo is still unknown and may vary from 10 6 M to 10 8 M [3], thus σv 0 can be chosen at the order of 10 −(24∼26) cm 3 /s in the analysis, which is widely adopted in many papers. The other feature is that the e + e − spectrum generated by the two-step annihilation χχ → Z Z → lll l exhibits boxshaped features, and the peak-like structure is significant only when the mass splitting ∆m ≡ m χ − m Z is small. Since the location and the shape of the spectrum are decided by the parameters m χ and ∆m, the excess should impose non-trivial constraints on their values.
In our analysis, we fix d = 0.1 kpc, M halo = 1.9 × 10 7 m and scan the masses m χ and ∆m. For each set of m χ and ∆m, we generate the e + e − spectrum arising from both the background and the signal χχ → Z Z with Z → e + e − , µ + µ − , and vary the parameters  Total fit DM Background DAMPE data Figure 2. Best-fit of the e + e − spectrum generated by the process χχ → Z Z with Z → e + e − , µ + µ − to the DAMPE data with the background parameterization taken from [16]. The left panel is based on background-only hypothesis with corresponding χ 2 = 25.5, while the right panel corresponds to background + signal fit with χ 2 = 17.0.
for the background and σv 0 to get the maximum value of the constructed likelihood function. Then we obtain the χ 2 of the fit by the formula χ 2 = −2 ln L max . In the Left panel of Fig.1, we show the map of the χ 2 on the ∆m − m χ plane with the color bar denoting the values of the χ 2 . We find that the best fit point locates at (7 GeV, 3060 GeV) with σv 0 = 2.34 × 10 −26 cm 3 /s and its prediction on the χ 2 is smaller than that obtained from the background-only spectrum by about 5.8. In this panel, we also plot the constant contour of χ 2 = χ 2 best + 2.3. The region bounded by this contour is interpreted as the best region of the two-step DM annihilation process to explain the DAMPE excess at 1σ level. In the Right panel of Fig.1, we present the e + e − spectrum predicted by the best fit point. In obtaining this panel, we take the parameters for the background from the maximization of the likelihood in calculating the χ 2 . This panel indicates that by choosing appropriate (∆m, m χ ), the process χχ → Z Z → lll l is indeed able to re-produce the DAMPE peak.
We checked that if we adopt the parameterization of the background by the formulae in [16], we can roughly reproduce the results of [16]. In this case, due to the differences in background parameterization and data set 1 , the best point locates at (8 GeV, 3061 GeV) with σv 0 = 2.8 × 10 −26 cm 3 /s, and the χ 2 of the best point improves that of the backgroundonly hypothesis by about 8.5. In Fig.2, we plot the best-fit spectrum with the background considered in [16], and the left and right panels correspond to background-only hypothesis and background + signal hypothesis, respectively. This figure indicates that adopting the parameterizations of the background in [16] can also fit the data well, and does not result in a significant difference in underlying physics.
In this section, we introduce the key features of the scalar DM model which extends the SM gauge group G SM ≡ SU (3) C × SU (2) L × U (1) Y by the U (1) Le−Lµ gauge symmetry. The particle contents of the model and their charges under the gauged groups are presented in Table 1, where we introduce a complex scalar φ χ as DM field and a complex scalar φ s to be responsible for the broken of the U (1) Le−Lµ symmetry as well as for the generation of the Z mass by its vacuum expectation value (vev) v s [27,28]. We impose an odd Z 2 parity for φ χ and an even parity for the other fields to guarantee the stability of the DM candidate. The Lagrangian relevant to our discussion includes In the above expressions, D µ = ∂ µ − ig Y Y φ Z µ with Z denoting the gauge field of the new symmetry, F µν = ∂ µ Z ν − ∂ ν Z µ , and m i , λ i and λ ij with i, j = H, χ, s are all real free parameters. In Eq.(3.1) we neglect the kinematic mixing terms of the two U (1) gauge fields. Note that the λ χs term must vanish if Y χ = Y s . Also note that in Eq.(3.2) we do not consider the special cases of Y s = 0 and 2Y χ ± Y s = 0. As a result, cubic operators such as φ ( * ) s φ χ φ χ are absent. Moreover, due to the assigned Z 2 symmetry, bilinear operators such as φ * χ φ s are also absent. In our discussion, we set λ χH = λ sH = 0 and choose proper m 2 H and λ H so that the properties of the H-dominated scalar are the same as those of the Standard Model (SM) Higgs boson. The minimization conditions of the potential for the fields φ χ and φ s are then given by where v χ denotes potential vev of the field φ χ . The first condition indicates that for λ χ ≥ 0, the field φ χ will not develop a vev if m 2 χ +(λ χs +2λ χs )v 2 s > 0. In this case, m 2 s = −2λ s v 2 s and one can replace m 2 s by v s as a theoretical input parameter. After the symmetry breaking, the terms induced by λ χs and λ χs in the potential can be rewritten as where we define the imaginary and real parts of the field φ χ as χ ≡ φ χ,I and χ ≡ φ χ,R , respectively, and neglect the field φ s,I since it becomes the longitudinal component of the gauge particle Z . This expression indicates that λ χs can induce a mass splitting between χ and χ , ∆m 2 χ χ ≡ m 2 χ − m 2 χ = 8λ χs v 2 s , and that there is no mixing between χ and χ . Since the parameters λ χ and λ χs are unimportant to our discussion, we hereafter treat them as zero .
The gauge charges of the fields φ χ and φ s play an important role in DM physics, and some important features about them are listed below.
• Because in most cases the DM annihilation χχ → Z Z proceeds mainly via the quartic χχZ Z coupling, the DM relic density measured by Planck experiment [21,22] and σv ∝ Y 4 χ g 4 Y require Y χ g Y 1.2 for m χ ∼ m Z ∼ 3 TeV [13]. This implies that the larger value of Y χ one takes, the smaller g Y should be chosen. In this case, the LEP-II constraints on Z mass from e + e − → Z * → e + e − becomes relaxed by the formula | g Y Ye m Z | 2.02 × 10 −4 GeV −1 presented in [36], where Y e = 1 in our model.
• If Y χ = Y s , the λ χs term in Eq.(3.2) must vanish due to gauge invariance. In this case, χ and χ degenerate in mass and the complex field φ χ acts as the DM candidate. On the other hand, since the Z φ * χ φ χ coupling is non-zero, the DM-nucleon scattering can proceed by the t-channel exchange of the Z − Z transition induced by e/µ loops [37][38][39]. Consequently, the spin-independent (SI) DM-proton scattering cross section can be approximated by [26,39] This expression indicates that with Y χ g Y 1.2 for m Z 3 TeV, Y χ must be larger than about 2 in order to survive the tight bound from the DM direct detection Table 2. Status of the model G SM × U (1) Le−Lµ confronting the four conditions.

Condition
Result Details I-ID √ χχ → Z Z proceeds mainly by the quartic χχZ Z interaction with Br(Z → e + e − ) = 40% for m χ ∼ m Z ∼ 2 × 1.5 (TeV) in our model.

II-RD
√ Same as I-ID since the annihilation is a s-wave dominated process.

III-DD √
The splitting between m χ and m χ forbids the inelastic scattering χN → χ N by kinematics, while the elastic scattering χN → χN proceeds at two-loop level.

IV-Collider √
Y χ = 2 corresponds to g Y 0.6 for m Z 3 TeV, which is on the edge of being excluded by the LEP constraint on m Z . A larger Y χ can alleviate the tension. Since Z does not couple with quarks, LHC experiments do not provide any constraints. experiments such as XENON-1T [23] and PandaX-II [24]. Considering that the bound will be further improved greatly in near future, we conclude that the case is disfavored if no DM signal is confirmed.
• If Y χ = Y s , the non-vanishing λ χs term can induce a sizable splitting between m χ and m χ , especially m χ > m χ for a positive λ χs [13]. In this case, χ acts as the DM candidate. More importantly, there is no Z χχ or Z χ χ interaction, but only Z χ χ interaction [13]. The remarkable implication of the mass splitting is that if it is larger than the DM kinetic energy today ∼ 1 MeV, the inelastic scattering process χN → χ N with N denoting nucleon is kinematically forbidden, while the elastic scattering process χN → χN only starts from 2-loop level since Z does not couple directly with quarks [13]. As a result, the DM direct detection experiments will no longer limit the model.
• In principle, our framework may also account for non-zero neutrino masses and mixing angles if one further introduces additional scalars [4]. For this purpose, Y s is naturally chosen to be 2 so that it can couple to the right-handed neutrinos to generate their masses. For such a model, only Z among the gauge bosons couples with the righthanded neutrinos, which has important implications at colliders [40,41].
Based on the above arguments, we take Y χ = Y s = 2 in the following as an example to study the implication of the model in DM physics.

Numerical results
In the numerical calculations, we use the public package SARAH [42] to implement the model, SPheno [43,44] to obtain the mass spectrum and micrOMEGAs [45,46] to calculate the DM relic abundance in which the threshold effects may be important when m χ ∼ m Z [47]. We scan the following parameter space by the package Easyscan-HEP [48] which is based on Markov Chain Monte Carlo (MCMC) sampling technique [49], 0 ≤ λ s , g Y , λ χs ≤ 1, 2.9 TeV ≤ m χ ≤ 3.2 TeV, 1 TeV ≤ v s ≤ 5 TeV. (4.1) In the scan, we impose the constraints (see the four conditions listed in Section I) where both of the two physical scalars φ s,R and χ are required to be heavier than 3.5 TeV so that they do not affect the DM annihilation. Moreover, we require that m χ and ∆m are in the region bounded by the constant χ 2 contour in Fig.1.
The numerical results are provided in Fig.3 on the plane of g Y versus m Z , where the upper left region of the LEP-II bound, i.e. g Y m Z 2.02 × 10 −4 GeV −1 [36], has been excluded. This figure indicates that the LEP-II bound has tightly limited the ability of the model to explain the DAMPE excess, nevertheless there still exist surviving samples that are capable of the explanation. As we discussed in previous section, the tension can be  [50], which is different from our previous work [13]. In order to illustrate clearly the status of our explanation confronting the four conditions listed in Section I, we summarize the details in Table 2.
Before we end the discussion, we have the following comments.
• As was pointed out in [4,15], the DM explanation of the DAMPE peak is consistent with other DM indirect detection experimental results such as the H.E.S.S. data on the annihilation χχ → V V → 2(e + e − ) [51,52], the Fermi-LAT data in the direction of the dwarf spheroidal galaxies [53], the Planck CMB data which is sensitive to energy injection to the CMB from DM annihilations [54,55], and the IceCube data on DM annihilation into neutrinos [56]. It also survives the upper bounds from XENON-10 and XENON-100 experiments on the DM scattering off electron [57].
• The DAMPE explanation imposes non-trivial requirements on g Y and m Z which will be further tested at the future International Linear Collider (ILC). For example, it is expected that the sensitivity g Z m Z 2.2×10 −5 GeV −1 will be achieved at the ILC with √ s = 1 TeV and a luminosity of 500 fb −1 [58]. To gain a more complete picture, we also show the modified production cross section σ(e + e − → µ + µ − ) in Fig.4 for m Z = 3 TeV with two benchmarks g Y = 0.6, 0.5. This figure indicates that the gauge boson Z can indeed make sizable difference in the production rate at the ILC.
• In our scalar DM model, Z contributes to the muon anomalous magnetic momentum g − 2 by [59] For g Y 0.6 and m Z 3 TeV, this contribution is about 3.35 × 10 −11 which is too small to account for the measured value of muon g − 2.
• Throughout this work, we have chosen the scalar DM model with the U (1) Le−Lµ symmetry to explain the DAMPE excess. Alternatively, one may also choose a similar model but with U (1) Le−Lτ gauge symmetry. We find that roughly the same results for the e + e − spectrum and the favored (∆m, m χ ) region as those in Fig.1 can be obtained. Therefore we conclude that the excess may also be explained by the U (1) Le−Lτ symmetry.

Conclusion
In this work, we discussed the feasibility of scalar DM annihilation to explain the tentative peak observed by the recent DAMPE experiment in the cosmic e + e − flux. Assuming that the two-step process χχ → Z Z → ¯ with ( ) = e, µ is responsible for the e + e − spectrum, we first determined by fitting the spectrum to the DAMPE data the favored masses of the DM χ and the mediator Z , which are presented on ∆m − m χ plane in Fig.1. Then we considered a scalar DM model with gauged U (1) Le−Lµ symmetry to realize the process and investigated the features of the model in explaining the excess. Our results indicate that the model can account for the excess without conflicting constraints from DM direct and indirect detection experiments as well as collider experiments. Depending on the assignment of DM charge Y χ under the U (1) Le−Lµ symmetry, one can alleviate the tension between the DAMPE interpretation and the LEP-II bound on m Z . These observations are presented in Fig.3 and also summarized in Table 2.
Before we end this work, we want to emphasize some facts. First, the scalar DM model with the U (1) Le−Lµ symmetry is one of the most economical anomaly-free theory which may be used to explain the DAMPE excess. So its features in this aspect should be investigated carefully, especially in light of the fact that another economical anomalyfree U (1) extension of the SM studied in [13] fails in doing this. As we have shown in this work, the theory can explain the excess without conflicting with any experimental measurements, and a large Y χ is favored to relax the tension with the LEP-II constraint. These conclusions are rather new. Second, as we discussed in Section II, the measured e + e − flux constrain nontrivially the parameters space of m χ and m Z . In getting the χ 2 map on ∆m−m χ plane in Fig.1, we solved the propagation equation of electron/positron in cosmic ray and sought for the maximum value of the likelihood function, which depends on seven variables, for each set of (m χ ,∆m). So the involved calculation is rather complex and time consuming. We remind that, as far as we know, such a χ 2 map was not obtained in previous literatures to explain the DAMPE excess. Finally, we'd like to clarify the relation of this work with our previous work [13], where we utilized another U (1) extension of the SM to explain the excess. In either work, a scalar DM candidate from a minimal and anomaly-free theory was considered to generate the e + e − flux by the two step annihilation process χχ → Z Z → lll l . As a result, the scalar potentials in the two works and also the relevant experimental constraints are same, and therefore we organize the discussions in a similar way. However, the underlying physics is quite different, which can be seen in the following aspects: • the decay products of the mediator Z , or equally speaking the final state in DM annihilation. For the theory considered in [13], Z decays democratically into e + e − , µ + µ − and τ + τ − to explain the excess. This was motivated by the original work [3] about the DAMPE experiment, where the authors pointed out by simulation that pure e + e − final state or equally mixed e + e − , µ + µ − and τ + τ − final state is capable of predicting the right shape and height of the spectrum. In this work, however, Z decays with an equal rate into e + e − and µ + µ − states, and whether the final states can explain the excess is still unknown before our study. We stress that different final states can result in significant difference in the e + e − spectrum, and consequently can favor different ranges of physical inputs. For example, as was illustrated in a later discussion in [20], the former state prefers a larger mass splitting and a larger annihilation cross section than the latter state.
• more important, the status of the theory confronting the constraints. Explicitly speaking, the anomaly-free condition for the former theory requires that the quark fields in the SM carry a charge of 1 3 under the new U (1) symmetry. This charge assignment can induce the tree-level process pp → Z → l + l − at the LHC, and the analysis by ATLAS collaboration on dilepton signal has pushed the lower mass bound of Z up to about 4 TeV, which implies that the minimal framework fails to account for the excess [13]. By contrast, the explanation in this work is free of such a problem.