Galactic Positron Excess from Selectively Enhanced Dark Matter Annihilation?

Precision measurements of the positron flux in cosmic ray have revealed an unexplained bump in the spectrum around $E\simeq 300\,\mathrm{GeV}$, not clearly attributable to known astrophysical processes. We propose annihilation of dark matter of mass $m_\chi = 780\,\mathrm{GeV}$ with a late-time cross section $\sigma v = 4.63\times 10^{-24}\,\mathrm{cm^3\,s^{-1}}$ as a possible source. The nonmonotonic dependence of the annihilation rate on dark matter velocity, owing to a selective $p$-wave Sommerfeld enhancement, allows such a large signal from the Milky Way without violating corresponding constraints from CMB and dwarf galaxy observations. We briefly explore other signatures of this scenario, and outline avenues to test it in future experiments.

Precision measurements of the positron flux in cosmic ray have revealed an unexplained bump in the spectrum around E 300 GeV, not clearly attributable to known astrophysical processes. We propose annihilation of dark matter of mass mχ = 780 GeV with a late-time cross section σv = 4.63 × 10 −24 cm 3 s −1 as a possible source. The nonmonotonic dependence of the annihilation rate on dark matter velocity, owing to a selective p-wave Sommerfeld enhancement, allows such a large signal from the Milky Way without violating corresponding constraints from CMB and dwarf galaxy observations. We briefly explore other signatures of this scenario, and outline avenues to test it in future experiments.

I. INTRODUCTION
Observations of the positron flux in cosmic rays have seen tremendous progress and excitement in recent years.
PAMELA [1], Fermi-LAT [2], MASS [3], Wizard/CAPRICE [4], AMS-01 [5], and HEAT [6] have reported the positron flux in the cosmic ray, with PAMELA and Fermi-LAT providing evidence for a rising positron flux around 100 GeV. The AMS-02 collaboration, in its latest published results with three times more statistics than previous measurements, has presented the positron flux up to energy 1 TeV, and the existence of a bump in the flux spectrum at E 300 GeV followed by a cutoff at E 810 GeV has been confirmed [7].
The positron flux spectrum shows several distinct features. The flattening at E 8 GeV seen in Fig. 1, is understood to be mainly caused by the diffusion of the positrons produced from the decay of the scattering products of cosmic ray with the interstellar gas [8,9], as described in Ref. [7]. However, the shape of the spectrum in the region E > 100 GeV, a rise followed by a sharp drop, is not yet fully understood.
DM annihilation/decay and particle acceleration by pulsars have been shown to be able to reproduce the positron spectrum but several concerns remain [10][11][12][13][14][15][16]. DM of mass (1.5 − 3) TeV and annihilation cross section (6 − 23) × 10 −24 cm 3 s −1 [11,12], such that it first annihilates into light intermediate states which subsequently decay to charged leptons, can produce a good fit to the earlier AMS data [17]. However, this large annihilation cross section into leptons is disallowed by the Fermi-LAT observations of the dwarf galaxies [18], which show no excess gamma ray flux. Such large DM annihilation cross section into visible sector particles during the recombination era is also highly constrained by CMB observation [19][20][21]. Moreover, note that this cross section is approximately two orders of magnitude larger than * anirbandas@theory.tifr.res.in † bdasgupta@theory.tifr.res.in ‡ anupam.ray@theory.tifr.res.in the canonical value, (2 − 3) × 10 −26 cm 3 s −1 , needed to produce the correct relic abundance of thermal DM [22]. The possibility of a boost in the annihilation due to the substructures in the DM halo was also explored [23], but such large boosts are unlikely due to substructures alone. Alternatively, the pulsars in the Milky Way (MW) have also been used to explain the positron spectrum with a possibility that the nearby young pulsars (e. g., Geminga and B0656+14) contributing the most [10][11][12]. However recently, Ref.
[24] used the HAWC and Fermi-LAT data to show that the diffusion of the positrons from the pulsars Geminga and B0656+14 is less efficient in the energy range (50 − 500) GeV in a two-zone diffusion model of cosmic ray, but can contribute a substantial fraction Overlaid is a white band of the region of parameter space that yields correct thermal relic abundance. The thick white line denotes the contour of σv relic = 2.5 × 10 −26 cm 3 s −1 , bordered by the dashed contours of σv relic = 3 × 10 −26 cm 3 s −1 (above) and 2 × 10 −26 cm 3 s −1 (below). The points marked with asterisks can provide large enough annihilation cross section to explain the positron flux excess. In this work, we use the point marked with orange asterisk.
of the positron flux. Therefore, the pulsar origin of the positron flux is also debated. All things considered, the shape of the positron flux spectrum in the high energy region still remains poorly understood.
In this work, we show that it is possible to explain the positron excess using selectively enhanced, velocitydependent p-wave annihilation of Majorana DM particles, following our earlier work [25]. An angular momentum and spin-dependent selection rule in Sommerfeld effect enhances the p-wave DM annihilation. The p-wave nature implies that the rate is a nonmonotonic function of DM velocity. With a suitable choice of parameters, it can be maximum in the Milky Way-like galaxies, and less in smaller dwarf galaxies. The DM annihilation rate in the MW is shown in the contour plot in Fig. 2. In this parameter space, there are a few islands, marked with yellow asterisks, that yield such large galactic annihilation signal, and satisfy the thermal relic constraint, i. e., σv relic (2 − 3) × 10 −26 cm 3 s −1 [22], the dwarf galaxy bound from Fermi-LAT [18], and DM annihilation constraint from CMB observation [19], all at the same time. A representative positron flux spectrum is shown in Fig. 1. In this scenario, Majorana DM particles of mass m χ = 780 GeV annihilate into light dark sector scalars which promptly decay into electrons and muons. The annihilation rate in the Milky Way is σv 4.63 × 10 −24 cm 3 s −1 . Our purpose, in the remaining part of the paper, is to provide some details of this result and the related technical background, and to outline other phenomenological consequences of this model. The paper is structured as follows: In Section II, we describe the dark sector of the model and its connection to the Standard Model (SM), followed by Section III discussing about the cosmological thermal history. In Section IV, we compute and explain the selective Sommerfeld enhancement in annihilation, and in Section V outline the relevant phenomenology. We summarize and conclude in Section VI.

II. THE MODEL
We extend the SM by a dark sector that has a complex scalar Φ and a Dirac fermion χ, with charges −2 and +1, respectively, under an approximate global U (1) symmetry. The right-handed leptons in the Standard Model, l R , are also charged +2 under this symmetry. The SM Lagrangian is therefore extended by L BSM ⊃ L dark +L portal , where The dark complex scalar Φ and the dark fermion χ are coupled to each other through a Majorana-type interaction [26] and the 5-dimensional nonrenormalizable operator, involving the SM Higgs doublet H, acts as a portal between the dark and the visible sectors. After the spontaneous breaking of the U (1) and the EW symmetries, it would lead to decay of the dark scalars into electrons and muons. Note that the U (1) is also explicitly broken at the treelevel by the electron and muon masses. It is also possible to include explicit U (1)-breaking terms in L dark itself. Such terms are crucial for determining the mass of the pseudo-Goldstone mode of Φ, but we will not model-build that aspect and treat the pseudo-Goldstone mass as a free parameter.
The scalar potential induces nonzero vacuum expectation value v Φ and splits Φ into a radial mode ρ and a pseudo-Goldstone mode η: Φ → (v Φ + ρ + iη) / √ 2, spontaneously breaking the U (1) symmetry. Due to the Majorana-type coupling with the scalar, χ also splits into two Majorana particles, χ 1 = (χ − χ c )/ √ 2 and χ 2 = (χ + χ c )/ √ 2, with masses m χ and m χ + ∆ where m χ ≡ M − f v Φ and ∆ ≡ 2f v Φ . After dark symmetrybreaking, the resultant interactions between the fermions and the scalars are With these interactions, {χ 1 , χ 2 } constitute a two-level inelastic self-interacting dark matter (SIDM). While the ρ couples similar DM states, η provides an off-diagonal interaction and couples different DM states. Additionally, the η-interaction allows χ 2 to decay to the lighter state, i. e., χ 2 → χ 1 η. The symmetry-broken scalar potential reads where m ρ = √ 2λv Φ and the mass of η arises from the U (1)-breaking terms and its interaction with other particles in the thermal bath.
The higher-dimensional part of the Lagrangian yields These terms allow the dark scalars ρ and η to decay into a pair of leptons. Further, the Higgs boson h develops new decay channels h → ρ l + l − and h → η l + l − and a new contribution to its decay into lepton pairs. The masses of the e and µ get a contribution from the dark symmetry breaking as well. At low energy, the DM annihilation phenomenology is controlled by the five free parameters {m χ , m ρ , m η , α, ∆}. The decay of the scalars is determined by c e , c µ , and Λ. We take the value of the SM vev as v H = 256 GeV.

III. COSMIC HISTORY
After the initial production of the dark sector particles through some mechanism like reheating, they decouple from the SM sector at a temperature T * . We assume this to occur at a scale much higher relative to other energy scales in the theory. After decoupling, the SM and the dark sectors evolve independently with separate temperatures T (SM) and T d (dark sector) with a ratio ξ d ≡ T d /T . After the dark sector symmetry breaking at T d v Φ , the dark sector consists of two Majorana fermions χ 1 , χ 2 and two scalars ρ, η. The comoving entropy of this sector changes only during the decays of χ 2 and ρ. In general ξ d can be written as Here g d (T ) and g SM (T ) are the relativistic degrees of freedom in the dark and the SM sectors, respectively, at temperature T . From Fig. 3, it is clear that the ratio ξ d is never too far away from unity. After ρ and η go out of chemical equilibrium, they decay. Note that ρ could decay into η if m η < m ρ /2. However, as we will show later, this not the case in the region of the parameter space we are interested in. If the dark sector decouples from the visible sector early enough then ρ, η do not affect the predictions for N eff and the Hubble parameter much during BBN, as was shown in [27]. Dark sector temperature T d rises after the DM decouples from the thermal bath. Later during the QCD phase transition at T ∼ 200 MeV, a large amount of entropy is dumped into the SM thermal bath, heating it up. This history depends on the relative values of M and T * (see Fig. 3 in Ref. [27]). In this paper, we shall only consider the case that the DM mass is above the QCD condensation scale. We shall always assume that M < v Φ , so that the symmetry breaking occurs earlier than the chemical freezeout of DM. The ensuing phenomenology was discussed in Ref. [26][27][28]. During relic annihilation, the interactions in Eq.(2) provide both annihilation and coannihilation, Even though in this case the DM is lighter than the symmetry breaking scale v Φ , the scalar particle can be much lighter than χ 1,2 if the hierarchy v Φ > M > µ Φ is maintained.
As χ 1 , χ 2 are both Majorana particles, their annihilation into two scalars does not have any s-wave ( = 0) component, and hence is p-wave ( = 1) suppressed. But coannihilation is not subject to such suppression. Even The Feynman graphs for DM annihilation and coannihilation.
after the DM is not in chemical equilibrium with the other particles in the dark sector, it exchanges kinetic energy with them as kinetic decoupling happens later than the chemical decoupling. In this case, the elastic scattering between χ 1 and η helps keep the DM in kinetic equilibrium with η. It was shown in Ref. [27] that a small mass split between the two DM states enhances this scattering cross section through a resonance. This delays the DM kinetic decoupling and may ameliorate the missing satellite problem [29,30].

IV. ENHANCED ANNIHILATION
When the annihilating DM particles interact with each other through a long-range potential V (r), then their initial plane-waves are modified. The Sommerfeld effect is the change in the annihilation rate resulting from this wavefunction modification. It is a nonperturbative effect. It is quantified by the Sommerfeld factor S and is written as, Here σ and σ 0 are the total and the perturbative crosssections, respectively, for the orbital angular momentum = 0, 1, 2, and so on. The process is enhanced if S > 1, suppressed if S < 1.
This effect is most pronounced when the DM particles are nonrelativistic. In this limit, it is useful to use the formalism of nonrelativistic effective theory (NREFT). In NREFT, the wavefunction Ψ of the heavy particles obey the Schrödinger equation with the interaction potential V (r). Relative to the range of V (r), the actual annihilation process is short-range and is usually taken to be happening at the origin. Hence the Sommerfeld factor S is defined as the change in the amplitude of the wavefunction Ψ (r) at the origin due to the potential, Here, Ψ (r) free is the wavefunction in the absence of any potential, i. e., the free particle plane wave. Depending on the nature of the mediator and the charges of the interacting particles, V (r) can be either attractive or repulsive, causing enhancement or suppression, respectively, to the annihilation rate. The possible two-particle states in the present case are the following.
The η-interaction leads to mixing between the two states in the annihilation space. But the annihilation and coannihilation subspaces remain decouple as neither ρ nor η interaction can couple one to the other. As mentioned before, the annihilation subspace does not have any swave process, but coannihilation has both s and p-wave contributions. Although the two states for coannihilation consist of the same particles, we write them separately because it is easier to understand the transition from one to the other due to the η-exchange. We shall first compute the perturbative cross sections for these three processes using NREFT, and then calculate the Sommerfeld factors by solving the matrix Schrödinger equations numerically.
A. Annihilation matrices DM annihilation typically happens at a length scale 1/m χ . Therefore typical momentum exchange in such a process is ∼ O(m χ ). In NREFT, the annihilation process cannot be described using tree level graphs. However, like in any EFT, all information regarding these 'high energy' processes is contained the series of higher dimensional operators, Here c ij,kl (m χ ) are the Wilson coefficients computed at the scale m χ , and χ i χ j χ k χ l are four-Fermi operatorsχ 1 χ 1 χ 1 χ 1 , χ 1 χ 2 χ 1 χ 1 , χ 2 χ 2 χ 1 χ 1 , and so on. The c ij,kl coefficients can be calculated by matching a four-point amplitudes in the full theory with the corresponding four-Fermi operator in the effective Lagrangian.
To classify the effective operators according to the spin and angular momentum of the two-body states, we write the Dirac spinors of χ i using the Pauli two-component With these expressions for the spinors, we expand the operators in Eq.(10) in powers of |p|/m χ to find the effective operators for each spin and angular momentum where O ij,kl ( 2s+1 J ) are the effective operators consisting of the spinors ξ i , η i corresponding to spin s and angular momentum of the associated two-body states [31] (see Appendix A for details about these operators).
Once the Wilson coefficients are known, the annihilation cross sections into the light scalars can be found by using the Cutkosky theorem: the annihilation cross section of a process In addition to |p|/m χ , we also expand the operators in powers of other dimensionless parameters, such as m ρ /m χ , m η /m χ , and ∆/m χ . We computed the Wilson coefficients for both annihilation and coannihilation using the FeynCalc package in Mathematica 1 and classified them according to their spin and angular momentum [33]. The procedure is described in detail in Appendix A.
To the leading order in the dimensionless parameters mentioned above, we get the following annihilation matrices.
1 Mathematica notebooks are available at this https URL.
The = 0 coannihilation matrix has opposite signs for the off-diagonal entries relative to the diagonal elements. However, all elements of the = 1 matrices have similar sign. This fact is related to the particle exchange symmetry, and will be discussed in more detail when we describe our results.

B. Sommerfeld factors
The nonrelativistic wavefunction of the DM particles obey the matrix Schrödinger equation, Here u (r) is a 2 × 2 matrix with the radial parts of the full wavefunction Ψ (r), and D ab = − 1 2µ d 2 dr 2 + ( +1) 2µr 2 − E δ ab is the differential operator. The elements of u (r) are the transition amplitudes u (r) ab ≡ b|a of going from state |a to |b , with |a , |b being one of the three possible two-body states in Eq. (6). As the particles are in scattering state before annihilation, the total energy of the a th two-body state is where v is the velocity of individual particles in the center-of-mass frame, and k a = µ a v rel is the momentum of the a th two-body state. The potential V (r) arises from the exchange of the scalars. It is given by with V 11 = −α e −mρr /r, V 12 = V 21 = −α e −mηr /r, and V 22 = −α e −mρr /r + 2∆ for annihilation and V 22 = −α e −mρr /r for coannihilation. The extra term 2∆ in V 22 for annihilation comes from the mass gap between the two two-body states |χ 1 χ 1 and |χ 2 χ 2 . No such mass gap exists for coannihilation. We follow Ref. [34] to compute the Sommerfeld factors for the multilevel DM. We quote the final expression for the Sommerfeld factor, S a 2s+1 J , that is given as No summation is implied over the repeated index a in the denominator. The matrix T is defined as the complex conjugate of the amplitude of the irregular solutions v (r) of Eq.(15) in the large r limit, The detailed procedure to compute the T matrix is given in the Appendix B. We solved the Schrödinger equation using two methods -i) directly using NDSolve in Mathematica, and ii) using the variable phase method [35,36]. The direct method works well as long as the kinetic energy of the incoming particles is above the threshold for the |χ 2 χ 2 state, i. e., When the incoming particles are below threshold, the wavefunctions for |χ 2 χ 2 final state are not scattering state anymore. In this case, the direct method fails as it simultaneously tries to solve for scattering and bound state solutions within a single system. Ref. [36] proposed that the variable phase method can be used in such cases. Instead of solving for the full solution, this method solves for the modification in the wavefunction from free particle solution because of the long range potentials. We have verified the results from the variable phase method with direct method solutions for above threshold parameters. We found excellent match between the results from two methods.
To reduce the number of parameters, we fix m η = 0.7m ρ for all the results shown in this paper. In Fig. 5 and 6, we show the Sommerfeld factors and the annihilation rates, respectively, as functions of DM velocity v. For the figures in the next two subsections, we choose a representative set of values for the parameters: α = 0.0373, m ρ /m χ = 0.0023, ∆/m χ = 0.001, v = 100 km s −1 which approximately corresponds to the point marked with orange asterisk in Fig. 2.
Both p-wave annihilation and coannihilation show large enhancement in the small velocity limit in Fig. 5. The enhancement factor has a 1/v 3 -dependence in the~υ intermediate velocity regime. The difference between the annihilation and coannihilation factors is due to the mass gap between the two states. However, the s-wave Sommerfeld factor is less than one. The p-wave annihilation and the s-wave coannihilation rates are shown in Fig. 6. The ∼ v 2 -scaling of the bare p-wave cross section yields the ∼ v 2 and ∼ v −1 behavior of σv p in the small and intermediate velocity regimes, respectively. We note that the relic annihilation cross section is dominated by annihilation because of a residual Sommerfeld enhancement factor S p 2 even for v relic 0.2. Therefore, p-wave annihilation process contributes more by a factor of ∼ 2 than coannihilation in the early Universe. This more carefully computed result disagrees somewhat with the conclusions made in previous papers [25,27,28]. However, we must note that at this large velocity v ∼ 0.2, the partial wave expansion of the cross section may not be very accurate as higher partial wave contributions become important. Fig. 7 shows the m ρ -dependence of the factors. In addition to a large overall enhancement for p-wave, a resonance feature is also present for certain values of m ρ . We shall explain the origin of these features in the rest of this section.

C. Particle exchange symmetry
We shall use the particle exchange symmetry following the Ref. [25]. Suppose, A and B are two fermions. They can form two two-body states, namely |AB and |BA of total angular momentum and spin s. However, these two states consist of the same set of particles and are related to each other through This factor has three components: (−1) from relative angular momentum, (−1) s+1 from their spins, and (−1) due to Wick exchange of two fermions. We can directly apply this to the coannihilation channel in the present model implying that two states |χ 1 χ 2 and |χ 2 χ 1 in the matrix Schrödinger equation are related to each other, It further implies that we can combine the two equations into a single equation with an effective potential which is a linear combination of the diagonal and off-diagonal potentials, Therefore V eff has the following forms for the two cases, We note that in the = 0 case, the difference between the two potentials is acting as the effective 1D potential. When r 1/m ρ and 1/m η , V eff saturates at the s-wave p-wave value (m ρ − m η ). For larger r, it gradually decreases to zero never becoming negative. The nature of the net potential thus becomes repulsive (See Fig. 8). However for = 1, the diagonal and off-diagonal potentials are added with the same sign, and hence yield an even stronger attractive potential. This explains the behaviour of the different Sommerfeld factors in Figs. 5 and 7. The effective repulsive nature of V eff leads to the suppression in S co-ann s .

D. Velocity dependence & resonances
The 1/v 3 -dependence of S p in Fig. 5 can be understood by taking the limit m ρ /m χ 1 and α/v 1. In this Coulomb limit, an analytic expression for S for general angular momentum is given by [37] where S 0 is the s-wave factor in the Coulomb approximation Therefore for small velocity, S ∼ 1/v 2 +1 . However for finite m ρ , the Sommerfeld factor does not grow indefinitely for smaller velocity. Instead, it saturates to a constant value when the de Broglie wavelength of the DM particles gets much longer than the range of the potential, i. e., approximately when the condition µ 1 v/m ρ 1 is met. All of these features are evident in Fig. 5.
Resonances from virtual bound states arise and cause large Sommerfeld enhancement as can be seen in Fig. 7. This occurs whenever the range of the potential matches a multiple n of the Bohr radius 1/(αm χ ) of the DM particles, i. e., [37] αm χ κm ρ = n 2 , n = 1, 2, 3, . . . .
Here κ π 2 /6. Note that this resonance condition is only approximately true in the present model as we have two mediators with slightly different masses. Also the resonance positions are shifted by the mass gap ∆ [34].
A consequence of the nonmonotonic velocity dependence of the p-wave Sommerfeld factor is that it predicts different annihilation rates of DM in astrophysical objects of different size. In the net p-wave annihilation rate, the perturbative cross section provides a velocity scaling ∼ v 2 . Therefore, S p (σv) p increases as ∼ v 2 for small velocity, reaches a maximum, and then falls as ∼ 1/v as shown in Fig. 6. The position of the maximum annihilation depends on the ratio m ρ /m χ . Hence this model naturally predicts large DM annihilation signal in the galaxies, but very small signal from the dwarf galaxies.

A. Galactic positron excess
The AMS-02 and several other experiments have observed an excess in the positron spectrum from the galaxy, as mentioned in the introduction. DM annihilation in our model can be used to explain such observation. To this end, we show the variation of the p-wave annihilation rate (σv) p in the Milky Way with m ρ /m χ and α in Fig. 2. The cross section is enhanced for large α and small m ρ /m χ ratio, with a resonance feature as discussed above. The points within the overlaid white band yield a relic annihilation cross section within (2 − 3) × 10 −26 cm 3 s −1 to satisfy the relic abundance constraint [22].
The points marked with asterisks yield DM annihilation cross section σv > 10 −24 cm 3 s −1 , sufficient to explain the AMS-02 positron flux excess, without running into any problems with thermal relic or dwarf galaxy constraints. The typical annihilation cross section in dwarf galaxies is 6 × 10 −28 cm 3 s −1 , which is well below the Fermi-LAT bound in the e + e − and µ + µ − channels [18]. Because of the suppression in the s-wave channel as shown in Fig. 6, DM annihilation rate during recombination era is also predicted to be much smaller than the current experimental bound [19]. We show the positron flux from such a representative point in Fig. 1. We used the publicly available code PPPC4DMID to compute the positron flux spectrum in the χ 1 χ 1 → V V → 4e, 4µ channel, where the scalar V represents ρ and η, from the galaxy after diffusion [38]. 9. The model parameter space in the ce − cµ plane for Λ = 10 13 GeV. The gray line is the contour of BR(e + e − ) = 38% and BR(µ + µ − ) = 62%. In the blue-shaded region, ρ and η decay times are larger than τc. See text for definition. The orange-shaded region is excluded by lepton mass measurement.
An astrophysical background of positron flux was assumed, following Ref. [7], whereÊ = E +φ e + is the positron energy in the interstellar space. The dark sector and background parameters   used in Fig. 1 are listed in Table I. For m χ = 780 GeV, the required mediator mass is about m ρ 1.5 GeV (see Fig. 2). Therefore ρ and η can decay to electron and muon. Note that we are not using the background parameters quoted in Ref. [7] as this is an independent model. The fit shown in Fig. 1 corresponding to the parameter values quoted in Table I were found by a crude scan over the parameter space. Although the fit of the theoretical curve to the data appears acceptable, the goodness of the fit is rather modest. The annihilation products ρ, η couple to the visible sector through the charged leptons e and µ with branching ratios 38% and 62%, respectively, to get an acceptable fit. One may ask if it is possible to get these branching ratios within our model. The ΦHl L l R term, shown in Eq.(1), leads to decay of ρ and η generated in the DM annihilations. We compare their decay times with τ c ( 47 yrs) which is the typical time required by them to escape the MW without any scattering. The decay widths are parametrized by c e , c µ , and Λ. For c e , c µ 10 −5 and Λ = 10 13 GeV, their lifetimes τ ρ,η are less than τ c , as shown in Fig. 9 and 10. For our purpose, their decay can be assumed to be prompt.

B. Other BSM effects
The ΦHl L l R operator also modifies e and µ masses, as well as their coupling with the SM Higgs h. The electron and muon masses are predicted to change by an amount ∼ c l v Φ v H /Λ. The present experimental preci- sions of their mass measurement are very high, at the level of ∼ 10 −8 [39]. Specifically, electron mass measurement constrains our parameter space significantly. The exclusion regions due to this are shown as orange-shaded regions in Figs. 9 and 10. The branching ratio of Higgs decay into muons is also affected. However, the resulting change in the Higgs signal strength in the µ + µ − channel is very small relative to the current measurement, hence we do not consider it any further [40,41]. In Figs. 9 and 10, we show the contour of BR(e + e − ) = 38% and BR(µ + µ − ) = 62% with a gray solid line. Everywhere on the gray lines, the model produces a good fit to data, though the blue shaded region in the lower left is disfavored because the scalars decay too slowly (τ > τ c ) and the orange shaded region in the top right is disfavored because the fractional shift to lepton masses is larger than 10 −8 .

C. Self-scattering
The long range attractive potential would enhance DM self-scattering cross section as well. Large self-scattering cross section, say σ/m χ 1−10 cm 2 g −1 , may affect DM halo formation and produce core at the center [42,43]. At the same time, the Bullet cluster observation puts an upper limit σ/m χ 0.1 cm 2 g −1 on this cross section [44].
In Fig. 11  in Fig. 2, and has σ T 4 × 10 −5 cm 2 g −1 at dwarf-scale velocity v 10 −4 which is not large enough to produce cores in the DM halos. This is not surprising because indirect detection and CMB experiments put strong constraints on SIDM models with light mediators, and it is difficult to satisfy these constraints along with large selfscattering cross sections to produce core in DM halo [45].
However, multilevel SIDM models can generically have interesting scattering phenomenology originating from the inelasticity in the system. DM particles in the ground state can upscatter to the excited state and quickly deexcite to the ground state by emitting a light mediator particle that escapes the halo. As a result, the halo dissipates energy at a certain rate and cools [46,47]. This cooling mechanism can significantly affect the halo dynamics and profile in other parts of the parameter space where self-scattering is large [48,49].

D. Kinetic decoupling of dark matter
The scattering between DM χ 1 and η keeps the DM in kinetic equilibrium with the relativistic η particles longer than the usual. The kinetic decoupling happens roughly when the momentum transfer rate Γ ηχ = (T d /m χ )n χ σ ηχ becomes smaller than the Hubble expansion rate, Γ ηχ H(T ). A t-channel resonance due to the mass gap ∆ between χ 1 and χ 2 enhances the cross section, and delays the kinetic decoupling of DM further. Such late kinetic decoupling of DM suppresses structure formation at low halo mass scale, and can explain the apparent under-abundance of satellite galaxies around the MW [30]. It was already pointed out in Ref. [27] that for m χ 1 GeV, ∆/m χ 10 −4.5 and massless η (acting as dark radiation), the kinetic decoupling temperature is ∼ 0.5 keV which is essential to solve the missing satellite problem. This temperature decreases further by a factor of a few when a more exact analytic expression for the cross section is used and the thermal distributions of χ 1 and η are taken into account, as shown in Fig. 12.
However, note that explaining the positron excess requires the DM mass to be m χ ≈ 800 GeV which would increase the decoupling temperature by several orders of magnitude from the keV scale. Also in the present scenario, the mediator mass m η ∼ 1 GeV which ceases to be relativistic much earlier. Thus in the model we present, there is no delayed kinetic decoupling. Of course, one could consider DM with M ≈ 800 GeV and a nearly massless η, but with much smaller ∆/m χ to keep the decoupling temperature at keV scale. The Sommerfeld effect phenomenology would change significantly in such a scenario and requires a separate study.

VI. SUMMARY & OUTLOOK
In this work, we have explored the possibility of explaining the bumps in the galactic positron flux, seen by AMS-02 and several other experiments, using Sommerfeld-enhanced DM annihilation without overpredicting the DM annihilation signal from dwarf galaxies.
We computed the Sommerfeld-enhanced annihilation rate in a two-level SIDM model, following our previous work in Ref. [25]. We show that even though the treelevel annihilation process is p-wave suppressed due to the Majorana nature of the DM, the nonperturbative Sommerfeld effect enhances the rate by ∼ 6 orders of magnitude. At the same time, the s-wave coannihilation is suppressed. This selective Sommerfeld enhancement/suppression was explained to be arising due to a particle exchange symmetry, and occurs only in multilevel DM models.
The dominating p-wave annihilation rate has a unique dependence on DM velocity which leads to lower annihilation rate in smaller objects, like dwarf galaxies, but large rate σv ∼ 10 −24 cm 3 s −1 in MW-sized galaxies. This can naturally explain that perhaps the positron flux excess in the cosmic ray observed by AMS-02 and several other experiments is due to the decay of the DM annihilation products into the charged leptons through a higher-dimensional operator. We investigated this possibility and found reasonably good match between the theoretical and observed spectra. This is the first work, to the best of our knowledge, interpreting the positron excess as a result of DM annihilation without violating the bounds from dwarf galaxy and CMB observations. One way to test this model is to search for gamma ray lines in, e. g., H.E.S.S. and Fermi-LAT data, originating as final state radiation from the electron and muon. Bremsstrahlung radiation typically shows up as a peak in the gamma ray spectrum around DM mass. However, in the present scenario DM annihilation products ρ/η cannot emit photon. Only the e ± /µ ± , to which the scalars decay into, can emit bremsstrahlung photon. Hence the resulting gamma-ray spectral shape is expected to be different from the usual case. A search for such gamma-ray spectrum can be done from the H.E.S.S. and Fermi-LAT data [50,51]. Non-observation of such gamma ray emission from dwarf galaxies could serve as an evidence for this model. We intend to do such analysis in a future work. Moreover, the high energy charged leptons in the decay cascade can also inverse Compton scatter with photons to produce additional signature [52]. The self-scattering cross section in such a DM model turns out to be inadequate to explain the existence of the cored density profiles of the DM halos.
Even though the positron excess has been seen by both Fermi-LAT and AMS-02, there exists tension between their data in the region (20 − 200) GeV. The reason behind this is not known. Future independent cosmic ray experiments with improved statistics and systematics, and more precise directional flux measurements will pin down the source of the excess.

ACKNOWLEDGMENTS
We are grateful to Torsten Bringmann, Subhajit Ghosh, and Tuhin Roy for useful discussions. This work was partially funded through a Ramanujan Fellowship of the Dept. of Science and Technology, Government of India, and the Max-Planck-Partnergroup "Astroparticle Physics" of the Max-Planck-Gesellschaft awarded to B.D and has received partial support from the European Union's Horizon 2020 research and innovation programme under the Marie-Sklodowska-Curie grant agreement Nos. 674896 and 690575. AD thankfully acknowledges the hospitality provided by Koushik Dutta at SINP during the final stages of this work. We acknowledge use of the FeynCalc package [33].
To find the annihilation matrices in Eq.(14) for a particular process |a → X A X B where |a = {|χ 1 χ 1 , |χ 2 χ 2 } for annihilation, and {|χ 1 χ 2 , |χ 2 χ 1 } for coannihilation, we performed the following steps-3. The T matrix obtained by inverting the B matrix which is defined below This procedure works well when the heavier annihilation channel |χ 2 χ 2 is kinematically open. When that is not the case, the wavefunctions χ 2 χ 2 |χ 1 χ 1 or χ 2 χ 2 |χ 2 χ 2 are exponentially growing/decaying. The matrix inversion in Eq.(B13) with those solutions becomes difficult and gives rise to numerical instabilities. To mitigate this problem, we have followed the modified variable phase method as prescribed in Ref. [36].