Boosted Dark Matter in IceCube and at the Galactic Center

We show that the event excess observed by the IceCube collaboration at TeV--PeV energies, usually interpreted as evidence for astrophysical neutrinos, can be explained alternatively by the scattering of highly boosted dark matter particles. Specifically, we consider a scenario where a $\sim 4$ PeV scalar dark matter particle $\phi$ can decay to a much lighter dark fermion $\chi$, which in turn scatters off nuclei in the IceCube detector. Besides these events, which are exclusively shower-like, the model also predicts a secondary population of events at $\mathcal{O}(100 \text{TeV})$ originating from the 3-body decay $\phi \to \chi \bar\chi a$, where $a$ is a pseudoscalar which mediates dark matter--Standard Model interactions and whose decay products include neutrinos. This secondary population also includes track-like events, and both populations together provide an excellent fit to the IceCube data. We then argue that a relic abundance of light Dark Matter particles $\chi$, which may constitute a subdominant component of the Dark Matter in the Universe, can have exactly the right properties to explain the observed excess in GeV gamma rays from the galactic center region. Our boosted Dark Matter scenario also predicts fluxes of $\mathcal{O}(10)$ TeV positrons and $\mathcal{O}(100 \text{TeV})$ photons from 3-body cascade decays of the heavy Dark Matter particle $\phi$, and we show how these can be used to constrain parts of the viable parameter space of the model. Direct detection limits are weak due to the pseudoscalar couplings of $\chi$. Accelerator constraints on the pseudoscalar mediator $a$ lead to the conclusion that the preferred mass of $a$ is $\gtrsim 10$ GeV and that large coupling to $b$ quarks but suppressed or vanishing coupling to leptons are preferred.


I. INTRODUCTION
The IceCube experiment at the South Pole has recently made international headlines by discovering an excess of events in the energy range from 30 TeV to 2 PeV [1][2][3]. These events are usually interpreted as evidence for a flux of astrophysical neutrinos with a power-law spectrum ∼ E −2 ν , originating from the production and subsequent decay of charged pions, kaons, muons and neutrons produced in collisions of ultra-high energy charged cosmic rays with protons or photons in astrophysical sources.
In this paper, we explore another alternative idea, namely that IceCube may be observing dark matter (DM) particles with PeV energy directly (as opposed to observing only neutrinos from their annihilation or decay). The idea, which has first been put forward in [16], is the following: a heavy O(PeV) DM species φ, which makes up a substantial fraction of the dark matter in the Universe, decays to a much lighter species χ. The resulting flux of highly boosted χ particles scatters on nuclei in the IceCube detector and leads to the observed energy deposits E dep up to few PeV. An upper cutoff on E dep is naturally provided by the mass of φ, explaining the absence of events above a few PeV.
The idea of direct detection of boosted DM in large volume terrestrial experiments was introduced in ref. [17] in the context of new light (O(1) GeV) particles produced in the annihilation of O(100) GeV DM particles in the galactic halo. The authors focused on electron recoil signatures in Super-Kamiokande [18], Hyper-Kamiokande [19] and PINGU [20]. In a subsequent paper [21], also the possibility of detecting boosted particles from the annihilation of heavy DM captured in the Sun has been considered. Such signals can be enhanced if the heavy DM particles are self-interacting, so that their capture rate in the Sun is increased [22]. Also a model with "dark nucleosynthesis" could lead to mildly boosted dark sector particles emerging from the Sun [23]. The detection of boosted DM annihilation products at much lower energies GeV in direct DM detection experiments is discussed in [24]. The recoil energy spectrum and the annual modulation signal in this case are very distinct from those expected from scattering of ordinary non-relativistic DM. Compared to these previous works which focus on boosted DM with energies of 100 GeV, we study signals at even higher energies up to O(PeV), and we consider not only direct DM searches, but also indirect signatures which may be very relevant in our model.
In the context of the conventional neutrino interpretation of the IceCube events, there is an ongoing debate about the neutrino flavor ratios required to explain the data. The generic expectation for neutrino production from pion decay is that the flavor composition of the astrophysical flux at the source (S) is (ν e : ν µ : ν τ ) S = (1 : 2 : 0) S . After propagation and oscillation, the final flux at Earth (E) would have a composition of (ν e : ν µ : ν τ ) E ≈ (1 : 1 : 1) E . The IceCube events are categorized as track events and shower events, where the former are mostly from ν µ charged current (CC) interactions with nucleons in which the produced high energy muon leaves a track in the detector. The shower events are attributed either to neutral current (NC) interactions of neutrinos or to charged current interactions of ν e and ν τ . In several analyses, the flavor ratios of the IceCube events have been studied [25][26][27][28][29][30][31], and while in general, the data appears consistent with a (1 : 1 : 1) E flavor ratio, a mild lack of ν µ has been found.
A unique feature of the boosted DM scenario is that only shower events are predicted at PeV energy, while at the lower energies, the ratio of track events to shower events is similar to what is expected in the canonical interpretation of the data in terms of astrophysical neutrinos. The reason our model predicts also track events at low energy is that, in addition to the dominant flux of boosted χ particles, also a secondary flux of DM-induced neutrinos is expected. It arises when the particle that mediates DM-SM interactions-taken to be a pseudoscalar a here-is directly produced as final state radiation in the heavy DM decay, φ → χχa, and subsequently decays to SM particles. While the primary contribution to the IceCube data from χ scattering peaks at PeV energies but drops at lower energies due to the properties of the pseudoscalar interaction, the secondary neutrino flux peaks at O(100 TeV) energies. Thus, our scenario is also able to explain not only the observed ratio of shower to track events, but also the mild (though not yet significant) deficit of events in the intermediate energy range of few × 100 TeV.
Note that the IceCube collaboration has recently published a new analysis [32], the results of which are given separately for events coming from above, i.e. from the southern sky, and from below, i.e. from the northern sky. This analysis features a notable, but not yet statistically significant, bump in the event spectrum from the southern sky at around 80 TeV. Since the galactic center is located in the southern hemisphere, a decaying DM scenario like ours predicts a larger contribution from the southern sky than from the northern sky. Thus, this bump could be potentially interpreted as being due to the secondary neutrino flux discussed in the previous paragraph, which peaks at around 100 TeV.
In addition to the new window to the high energy Universe opened by IceCube, also observations at lower energies ∼ GeV have caused a stir recently. Namely, an excess of gamma rays from the vicinity of the galactic center was found in Fermi-LAT data, which could be explained by DM annihilation [33][34][35]. A good fit to the Fermi-LAT data is obtained for instance for a 30-40 GeV DM particle annihilating to bb with a thermally averaged cross-section of about σv rel ∼ 10 −26 cm 3 sec −1 , similar to the annihilation cross section expected for a thermal relic. In our scenario, a subdominant primordial population of the light DM species χ can naturally provide such signal by annihilation through s-channel exchange of the pseudoscalar mediator a. We will demonstrate that there is a viable region of parameter space which can explain the Fermi-LAT gamma ray signal and the IceCube signal simultaneously.
In the following, we first introduce our toy model of boosted DM in sec. II and then discuss the expected IceCube signals in sec. III. In particular, we show which regions of parameter space could explain the recently observed high-energy events. In sec. IV we review mechanisms for explaining the observed DM relic density [36] in the boosted DM scenario, and in sec. V we discuss the possibility that the galactic center gamma ray excess is explained by χχ annihilation along with the IceCube PeV events. We then discuss other constraints on the model in sec. VI, in particular limits from measurements of the cosmic positron and electron spectrum [37][38][39][40][41][42], from isotropic diffuse gamma rays [43,44], from direct detection experiments and from searches for the pseudoscalar mediator a in flavor physics experiments and at high energy colliders. We summarize and conclude in sec. VII.

II. THE FRAMEWORK
While most of the qualitative results of this paper apply to any PeV-scale boosted DM model, we consider as a specific example a toy model featuring a dark sector that contains two DM particles: a heavy real scalar φ with mass m φ ∼ O(PeV) and a light Dirac fermion χ with mass m χ ∼ O(10) GeV. We denote the relic abundance of φ by f φ Ω DM and the relic abundance of χ by f χ Ω DM , where Ω DM 0.258 is the total dark matter density in the Universe [45]. We will discuss in sec. IV how f φ and f χ could be determined in the early Universe. We assume that there are no other dark relics besides φ and χ, i.e. we assume f φ + f χ = 1. The dark sector Lagrangian reads  Table I. Summary of our two benchmark points (BP), both of which can explain the IceCube event excess and the galactic center gamma ray excess. In both models, the pseudoscalar a is assumed to couple dominantly to b quarks. We also give the calculated values of the velocity-averaged χ annihilation cross section σv rel bb (relevant for the galactic center gamma ray excess), of the fractional abundance of the light DM species f χ = 1 − f φ and of the branching ratio for the radiative decay φ → χχa. Note that benchmark point 1 can be realized only in the Vector-like quark model since in the MSSM-like and Flipped scenarios, laboratory constraints on g Y b are too strong (see sec. VI).
Here, the coupling constant y φχ determines the φ → χχ decay rate. We assume y φχ to be tiny, so that the lifetime of φ is significantly longer than the age of the Universe. One possible way of explaining the smallness of y φχ could be to envision φ as a composite particle made up of superheavy constituents Q φ and held together by a new confining gauge interaction. When this new gauge symmetry is broken by a tiny amount, a correspondingly small mixing between the Q φ and χ could be generated. Note that we do not include quartic couplings of φ or Higgs portal couplings in eq. (1) since these interactions will not be relevant to our phenomenological discussion. A possibly problematic term could be an operator of the form φ(H † H), but we assume the mechanism that suppresses y φχ also forbids or suppresses this operator. The light DM χ interacts with SM particles through a pseudoscalar mediator a [46][47][48][49]. This pseudoscalar couples to light DM and Standard Model fermions through the Lagrangian where g χ and g Y f are real couplings of a to light DM χ and to Standard Model fermions f , respectively, m f are the SM fermion masses and v 246 GeV is the vacuum expectation value (vev) of the SM Higgs field. While generically all g Y f are free parameters, we will specifically consider natural scenarios in which the g Y f are generation-independent.
Throughout most of the paper, we will consider two benchmark points in the parameter space of the model, defined in table I. The heavy DM mass m φ and lifetime τ φ , the light DM mass m χ , and the couplings g χ and g Y b are chosen such that both the IceCube excess of high energy events as well as the galactic center gamma ray excess are explained. We assume the mass of a to satisfy m a 10 GeV since constraints are weak in this case (see sec. VI), thus allowing large couplings g Y f to fermions. This is important for the model to fit the IceCube data and is also interesting because it allows for a detectable indirect signal from the annihilation of non-relativistic relic χ particles.
Since the coupling of the pseudoscalar a to SM fermions in eq. (2) should be considered as an effective operator after the spontaneous breaking of electroweak symmetry, we need to discuss possible ultraviolet completions for such an operator. We consider here three interesting models which can provide such a coupling. MSSM-like model. In the first model, the pseudoscalar a mixes with an extended Higgs sector, for example with the pseudoscalar A 0 in a type-II Two Higgs Doublet Model (2HDM), by a term of the form iaH † 1 H 2 + h.c. [50]. In this case, the Higgs couplings to quarks and leptons are the same as in the Minimal Supersymmetric Standard Model (MSSM). We therefore denote this model as MSSM-like. The relations for the couplings between effective operator model and the complete renormalizable model read [50] where tan β = v 2 /v 1 is the ratio of the two Higgs vevs and sin θ is the mixing angle between the pseudoscalar a and the A 0 boson of the 2HDM. g Y d , g Y and g Yu are the generation-independent normalization factors of the Yukawa-like couplings for down-type quarks, leptons and up-type quarks, respectively. Since the pseudoscalar a couples to SM fermions only through its mixing with A 0 , all of these couplings are suppressed by sin θ. As mentioned in the Introduction, we are interested in particular in scenarios with large coupling between the pseudoscalar a and bottom quarks to optimally fit the galactic center gamma ray excess. This requires large tan β to lift up the coupling to down-type quarks. Already at this stage, we can see that the MSSM-like model will be constrained by experiments sensitive to anomalous couplings of the charged leptons (which are also tan β-enhanced) and by searches for an extended Higgs sector. As we will see in sec. VI, these constraints lead to the conclusion that the IceCube events and the Fermi gamma ray excess can be simultaneously explained in the MSSM-like model only when the pseudoscalar a is heavy (m a m h /2). Flipped model. The second model, which we call Flipped is a flipped Two Higgs Doublet Model [51][52][53][54][55][56]. This means that one Higgs doublet couples to up quarks and leptons, while the other couples to down quarks. The difference between this model and the MSSM-like model is that the coupling to leptons in the Flipped scenario is proportional to cot β rather than tan β and is thus suppressed rather than enhanced in the large tan β region. Therefore, limits from the lepton sector will be significantly weaker. The couplings to up-type quarks and down-type quarks are the same as in the MSSM-like model.
Vector-like quark model. The third model has no extended Higgs sector, and the pseudoscalar mediator a does not directly couple to SM quarks. Instead, it couples to new, heavy vector-like quarks, which in turn mix with the SM quarks [57]. Since a has no couplings to leptons in this model and since there is no extended Higgs sector, we expect constraints to be weaker than in the other two scenarios. However, the mass of the heavy vector-like quark should be large to avoid LHC limits. Highly boosted χ particles from the DM decay process φ → χχ can scatter on atomic nuclei in the IceCube detector through their coupling to the pseudoscalar mediator a (see fig. 1 (a)). At the high energies we are interested in, the scattering is deep inelastic. Phenomenologically, this process is very similar to neutral current scattering of neutrinos, hence its characteristic signature is a shower-like event topology. The deposited (or visible) energy E dep in this case is the energy of the recoil nucleus or its fragments.
The total number of shower events from χ scattering in a given E dep bin [E min dep , E max dep ] is given by [26] N sh,NC Figure 1. The Feynman diagrams for (a) the scattering of light DM particle χ on nucleons and (b) the 3body decay φ → χχa, which produces a flux of high energy pseudoscalars whose decay products contribute to astrophysical neutrino, gamma ray and positron fluxes.
Here, T is the observation time, m N is the nucleon mass, dσ p(n) /dE dep is the differential scattering cross section on protons (neutrons The galactic contribution is given by [58], Here, m φ and τ φ are the mass and lifetime of the heavy DM particle φ, respectively, ρ halo is the DM density distribution in the Milky Way, r(s, ψ) is the position vector relative to the origin at the galactic center, s is the distance along the line of sight and ψ is its angular direction. We integrate the flux over the solid angle Ω ψ and integrate along the line of the sight s. The energy spectrum of boosted χ particles is simply dN χ /dE χ = δ(E χ − m φ /2). The spectrum of antiparticles, dNχ/dEχ is the same. The extragalactic contribution to the flux of χ particles is [58] dΦ EG In this expression, H(z) H 0 Ω Λ + Ω m (1 + z) 3 is the Hubble expansion rate as a function of redshift z. It depends on the Hubble constant H 0 = H(0), the dark energy density Ω Λ ∼ 0.692 and the matter density Ω m ∼ 0.308. The cold dark matter density Ω DM is 0.258, and the critical density of the Universe ρ c is given by ρ c 4.9 × 10 −6 GeV/cm 3 [45]. Note that we do not account here for attenuation of the χ flux due to scattering on the interstellar and intergalactic medium. This attenuation is already small for neutrinos [58], and the χ scattering cross section is even smaller than the neutrino scattering cross section.
The differential cross section for χ scattering on a proton p (neutron n) of mass m N is where x is the Bjorken scale variable, s = m 2 χ + x 2 m 2 N + 2xm N E χ is the center of mass energy, and Q 2 = 2xm N E dep is the momentum transfer in the scattering. E χ is the energy of the incoming particle χ and the nucleon is assumed to be at rest initially. E dep is the energy transferred to the hadronic system in the lab frame during the scattering. We are interested in events with a large deposited energy E dep 10 TeV in this analysis due to the IceCube energy threshold.
The factor f p(n) q (x) in eq. (9) is the parton distribution function (PDF) for protons (neutrons) and quark flavor q. We use the PDFs from NNPDF3.0 [59], which are valid in the range x ∈ [10 −9 , 1] and Q 2 ∈ [2 GeV 2 , 10 8 GeV 2 ] and contain the most recent deep inelastic scattering data. 1 In the calculation, we set the PDFs equal to 0 when Q 2 is smaller than 2 GeV 2 . Because the cross section is proportional to 1/(Q 2 + m 2 a ) 2 , it becomes large when Q 2 is small. Our cutoff at low Q 2 would therefore affect the results for m a GeV. In the following, however, we focus on the mass range m a > 10 GeV, and we have checked that in this case the contribution of the Q 2 < 2 GeV 2 region to the cross section is negligible. If one is interested in an extremely light t-channel mediator with m 2 a 2GeV 2 , then the Q 2 2GeV 2 and x 1 region, corresponding to exchange of a nearly on-shell a, is important. In this region, the PDF description breaks down and one should instead calculate the cross section for a * absorption by protons along the lines of the equivalent photon approximation in deep-inelastic scattering of electrons on protons. We have used the central values of the PDFs, but have checked that varying them within the error band changes the total cross section by only O(10%). The impact on the energy dependence of the differential cross section is also negligible.

III.2. Secondary signal: neutrino flux from 3-body decays of heavy DM
As mentioned in the Introduction, the boosted DM scenario predicts not only a population of high energy events from the scattering of boosted χ particles, but also a contribution at lower energy from neutrinos produced in the 3-body decay φ → χχa (see fig. 1 (b)), followed for instance by a → bb. Since boosted DM can only explain the IceCube events if χ particles can scatter on nucleons, a mediator particle like a is always needed and the existence of the 3-body decay process is thus very generic. Making a heavy does not significantly influence the 3-body decay rate unless m a becomes comparable to m φ . The differential decay width of the 3-body decay is, in the limit m a → 0 and at leading order in m χ , where E a is the energy of a in the rest frame of φ. The branching ratio is where is the rate of the dominant 2-body decay. In the second line of eq. (13), we have assumed that g χ is small so that the 3-body decay width is much smaller than the 2-body decay width. We can see from Table I that this assumption is satisfied at our benchmark points. We plot the energy spectrum of a particles from 3-body decay of φ in fig. 2. The decay of a to light quarks or b quarks produces neutrinos after parton showering, hadronization and hadron decay. We take the spectra of the secondary neutrinos from each a decay in the a rest frame from [61] and boost them into the laboratory frame by folding with the E a distribution from fig. 2 [62]. Multiplying by BR 3 (φ → χχa) gives us the number dN ν /dE ν of neutrinos per energy interval dE ν per φ decay. The flux of secondary neutrinos is then obtained from equations very similar to eqs. (7) and (8) by simply replacing the factor dN χ /dE χ by dN ν /dE ν . The strength of the indirect signal is proportional to g 2 χ f φ /τ φ once the masses m φ , m χ and m a are fixed. In principle, one might also include a factor of the form exp[−Abs(E χ , z)] in the expression for the extragalactic flux to account for the absorption of neutrinos in interactions with the cosmological relic neutrino background and with the intergalactic medium [58]. However, these effects are negligible in our analysis and we therefore do not include such an attenuation factor. Moreover, the high energy neutrino flux reaching the detector from below is affected by neutrino interactions during passage through the Earth. In particular, at energies above ∼ 100 TeV, the neutrino-nucleon interaction cross section is so large that the Earth can attenuate the neutrino flux. On the other hand, electron and muon neutrinos can be regenerated in the decay of tau leptons produced in ν τ CC interactions. The net effect of both absorption and regeneration is a reduction of the neutrino   Figure 3. The galactic and extragalactic neutrino fluxes from 3-body decay of heavy DM, φ → χχa, followed by a → bb. We have added up the neutrino and antineutrino fluxes and have also summed over neutrino flavors. The horizontal dashed line shows the generic flux expected from astrophysical sources, E −2 ν , normalized such that optimum agreement with the IceCube data is achieved [3]. The model parameters are set to the benchmark values given in the plot. flux by about 15% at neutrino energies ∼ 100 TeV [29], and we therefore neglect this small effect in our calculation.
We plot the expected contributions to the neutrino flux from galactic and extragalactic φ → χχ + (a → bb) decays in fig. 3 for our two benchmark points. Since the neutrinos originate mostly from meson decays after hadronization of the b quarks, their flavor ratio after propagation is naturally (1 : 1 : 1) E . Therefore, we have summed the different flavors, as well as the neutrino and antineutrino fluxes, in fig. 3. We see that the secondary neutrinos are softer by about one order of magnitude compared to the boosted DM particles χ. The extragalactic flux is in general softer than the galactic one due to redshift.
Note that, besides the secondary neutrino flux, there is also a population of boosted DM events from scattering of the χ particles produced in 3-body decay φ → χχa. We neglect these events for the following reasons: first, the 3-body branching ratio is almost two orders of magnitude smaller than the 2-body branching ratio. Second, the spectrum of χ particles from 3-body decays is softer than the one from 2-body decays and would therefore contribute only in a regime with larger expected backgrounds. Third, a 3-body decay φ → χχa produces only two χ particles, but typically more than two neutrinos [61]. Thus, the flux of χ particles from 3-body decay is subdominant compared to the secondary neutrino flux. Fourth, the χ scattering cross section on nucleons is usually smaller than the neutrino charged current cross section.

III.3. Fitting procedure
To determine the preferred parameter regions for the boosted DM scenario, we use the log likelihood ratio (LLR) method. The LLR is defined as follows: , B i and O i are the predicted signal event rate, the predicted background event rate, and the observed event rate in the i-th energy bin, respectively. ∆B i is the 1σ error on the background prediction. When the nuisance parameter x is 1 (−1), the error x ∆B i (x) describes the upper (lower) limits of the error band, and when x = 0 the background takes its central value. The term f Gauss (x) corresponds to a normal distribution in x and is the Gaussian pull term for the nuisance parameter x. By using only one nuisance parameter, we effectively assume that the background uncertainty is correlated between bins.

III.4. Results
We show the results of our fit in fig. 4 and compare the best fit points to the IceCube data in fig. 5. For the mediator mass m a = 12 GeV (80 GeV), the three panels of fig. 4 give the best fit points (black (red) "+" signs) and preferred parameter regions (black unshaded contours (red shaded contours)) at 1, 2, 3σ confidence level. For m a = 80 GeV, the best fit point, marked by a red "×" sign, corresponds to one of our benchmark points from table I, while for m a = 12 GeV, the benchmark point (indicated by the black "×" sign) is slightly shifted compared to the best fit in order to be consistent also with the galactic center excess and with all constraints. The larger value of m a is particularly interesting for the MSSM-like and Flipped models, where it helps to evade important constraints from B s → µ + µ − decays and from h → aa decays. (see sec. VI.4). Note that we parameterize the parameter space in fig. 4 in terms of three parameters: the heavy DM mass m φ ; the combination g 2 Y b g 2 χ f φ /τ φ of the a coupling constants, the cosmological abundance f φ of the heavy DM particle φ and its lifetime τ φ , to which the χ scattering rate is proportional; and the ratio g 2 χ f φ /τ φ to which the interaction rate of secondary neutrinos is proportional. In the upper left hand plot, we also show constraints from the diffuse γ ray flux (see sec. VI.2) as thick black (red) lines. We always fix the mass of the light DM particle at m χ = 30 GeV, as motivated by the galactic center gamma ray excess, see sec. V. As expected, the best fit point is always around m φ ∼ 4 PeV due to the lack of IceCube events above 2 PeV. In fig. 5, we compare the IceCube data from ref. [3] to our predictions at the benchmark points. We also show the individual contributions to the spectrum separately: the atmospheric ("ATM") neutrino background (red dotted), the galactic (brown dashed) and extragalactic (black dot-dashed) fluxes of boosted χ particles, and the flux of secondary neutrinos from φ → χχ + (a → bb) decay (purple dashed).
We see that both the galactic and extragalactic χ fluxes contribute at PeV energies, with the latter being somewhat softer due to redshift. Actually, the integrated fluxes of the two components are comparable, but since the scattering cross-section is higher when the energy of the incoming χ particle is larger, the softer component is subleading experimentally. Below 1 PeV, the boosted DM event rates drop because of the Q 2 dependence of the scattering matrix element, eq. (10). In their place, the secondary neutrino flux takes over below ∼ 500 TeV, so that a good fit to the IceCube data is obtained at all energies. Note that the normalization of the secondary neutrino flux is set by the parameter combination g 2 χ f φ /τ φ and is thus not directly correlated with the boosted DM scattering rate, which is proportional to g 2 Y b g 2 χ f φ /τ φ . Comparing our two benchmark values of m a (shaded vs. unshaded contours in fig. 4, left vs. right panel in fig. 5), we observe that the choice of m a has a small influence on the spectral shape of the DM contributions, but its main impact is on the overall rate. Therefore, at larger m a , the best   Figure 4. Preferred parameter regions for the boosted DM scenario from our fit to IceCube high energy data [3]. The three panels show 2-dimensional projections of the 3-dimensional parameter space spanned by the heavy DM mass m φ , the product g 2 χ g 2 Y b f φ /τ φ to which the scattering rate of boosted χ particles is proportional, and the combination g 2 χ f φ /τ φ to which the flux of secondary neutrinos is proportional. (Here, g χ and g Y b are coupling constants, f φ is the cosmological abundance of φ, and τ φ is its lifetime.) Solid black unshaded (red dashed shaded) contours show the preferred parameter regions at 1, 2, 3σ for m a = 12 GeV (m a = 80 GeV) and the black (red) "+" signs indicate the best fit points. At m a = 80 GeV, the best fit point is identical to one of our benchmark points (red "×" sign) from table I, while for m a = 12 GeV we define our benchmark point (black "×" sign) slightly away from the best fit. This way, both benchmark points can also explain the galactic center gamma ray excess and evade all constraints. In the upper left hand plot we also show as a thick black (thick red) curve the strongest exclusion limits on the m a = 12 GeV (m a = 80 GeV) benchmark model, coming from diffuse γ ray searches (see sec. VI.2). We use m χ = 30 GeV for the mass of the light, boosted, DM particle here, motivated by the galactic center gamma ray excess, but note that m χ does not affect the IceCube event rate as long as m χ m φ .  Figure 5. Comparison of IceCube high energy data [3] to the prediction at our two benchmark points (see Table I. We plot the signals from galactic (brown dashed) and extragalactic (black dot-dashed) φ → χχ decays, as well as the contribution from secondary neutrinos produced in φ → χχ+(a → bb) (purple dashed) separately. The red dotted lines show the atmospheric neutrino background ("ATM"), the blue bars depict the background uncertainty and the solid blue lines show the total expected event rate. We have taken the mass of the pseudoscalar mediator m a to be 12 GeV (80 GeV) in the left panel (right panel ). We always use m χ = 30 GeV for the mass of the light (boosted) DM particle here, motivated by the galactic center gamma ray excess, but note that m χ does not affect the IceCube event rate as long as m χ m φ .
fit value of g 2 Y b g 2 χ f φ /τ φ is significantly larger than at smaller m a . When a is heavy, one either needs large g Y b g χ coupling to keep the scattering cross section of the boosted DM particle χ on nucleons unchanged, or the flux of χ particles must be enhanced by decreasing the heavy DM lifetime τ φ . Note that the two benchmark models shown in figs.4 and 5 explain not only the IceCube data, but also the galactic center gamma ray excess (see sec. V) and are consistent with all constraints (see sec. VI).
An interesting aspect of our boosted DM scenario is that a dip in the event spectrum is predicted between recoil energies of ∼ 400 TeV and 1 PeV. This dip is more pronounced at larger m a , see right panel of fig. 5. This is in excellent agreement with the current data, which does not feature any events in this energy range. Therefore, if this lack of events should become statistically significant in the future, the boosted DM scenario would provide one possible explanation of it. Another interesting aspect of our scenario is that, at low energies, where the flux is dominated by neutrinos, the expected flavor ratio is (1 : 1 : 1) E after propagation for most decay modes of a. Thus the ratio of shower and track events is predicted to be the same as for the canonical astrophysical neutrino interpretation at E dep few × 100 TeV. On the other hand, at E dep ∼ 1 PeV, the predicted event rate is entirely dominated by the DM contribution, which only provides shower events. This is a unique feature of this model and can be tested with future data.
Let us also remark that a recent IceCube analysis [32] which separates events from the northern sky and from the southern sky, exhibits a noticeable, but not yet statistically significant, bump at energy deposits around 80 TeV in the southern sky. If this bump should become significant in the future, it could be interpreted as being due to a relatively large secondary neutrino flux in the boosted DM scenario. Since the galactic center, from where most of these secondary neutrinos are expected to come, is located in the southern sky, and because neutrinos from the northern hemisphere suffer some attenuation in the Earth, our model could explain why a similar bump is not observed in the northern sky.
Note that, without the neutrinos from the 3-body decay φ → χχ + (a → bb), the IceCube fit of our boosted DM scenario becomes much worse because the prediction would fall short of the observed number of events at energies ∼ 100 TeV. This could be avoided if a mediator with scalar rather than pseudoscalar couplings to fermions, or a vector boson mediator is considered. In this case, the boosted DM scattering cross section would not be proportional to (Q 2 ) 2 , and scattering of χ particles could explain the IceCube event excess across the spectrum. However, as we will argue in sec. VI.3, direct detection constraints in this case may be prohibitive. Ways to avoid these constraints include models with inelastic DM scattering or with a very small m χ 3 GeV, below the direct detection threshold. The second possibility would preclude a simultaneous explanation of the IceCube events and the galactic center gamma ray excess.
Let us finally discuss the morphology of the IceCube signal from boosted DM. While the extragalactic flux dΦ EG χ /(dE χ dΩ ψ ) is isotropic, the galactic component dΦ GC χ /(dE χ dΩ ψ ) peaks in the galactic center region. (Here ψ denotes the direction of sight.) The angular resolution in IceCube is about 10 • -20 • for shower events [1]. With this resolution and more statistics, a morphology study of the high energy events would provide an important consistency check of the boosted DM hypothesis.

IV. DARK MATTER RELIC DENSITY
An important problem of the boosted DM scenario which we have not addressed yet is how a particle with a mass of order PeV can account for the observed DM density in the Universe. For instance, thermal freeze-out is not a possibility at masses above few hundred TeV due to unitarity constraints [63]. A long-lived dark matter particle with a mass of O(PeV) can nevertheless have the correct abundance in the Universe [64][65][66][67][68].
Non-thermal production mechanisms for PeV DM include [67]: (1) production in cascade decays of the inflaton. In this mechanism, the DM abundance depends on the number density of inflatons and on the branching ratio of inflaton decay to DM. (2) production through inelastic scattering between high energy particles from inflaton decay and the hot plasma. When high-energy daughter particles scatter on the thermalized plasma, DM can be produced until the daughter particles' energy become less than E th = m 2 φ /(4T ). (3) For low reheating temperature, DM could be thermally produced with the correct relic abundance even when the maximum temperature of the Universe during reheating, T max , is larger than m φ , as long as the reheating temperature (defined as the temperature at which the inflaton energy density equals the radiation energy density) is smaller than m φ . The reason is that the continuing decays of the inflaton produce entropy after DM freeze-out, diluting the DM abundance. The authors of ref. [67] show that these mechanisms can account for the abundance of DM with O(PeV) mass. Mechanism (2) can achieve this even if the inflaton does not decay to DM and is thus highly model independent. PeV DM φ produced through this mechanism can for instance account for the observed abundance of DM in the Universe if the reheating temperature of order 10 GeV and the mass of inflaton is of order 10 15 GeV. [67].
In addition to the non-thermally produced relic abundance of heavy DM particles φ, there could also be a thermally produced population of the light DM species χ if the thermally averaged cross section σv rel for χχ annihilation through s-channel exchange of the mediator a is not too large. This is naturally realized in our scenario. σv rel receives contributions from two classes of processes, shown in fig. 6: annihilation to ff and, if m a < m χ , also annihilation to aa. The thermally averaged annihilation cross sections read [48] σv rel ff where m f are the SM fermion masses, the sum runs over all SM fermions f , Γ a is the total decay width of a, the color factor N f c is 3 if f is a quark and 1 if f is a lepton, and T is the temperature. The thermally averaged cross section for annihilation to leptons is completely analogous to eq. (16) except for the color factor. Note that eqs. (16) and (17) are approximate results, with only the leading terms in the relative velocity v rel kept. The proportionality to T in eq. (17) arises because the process χχ → aa is p-wave suppressed. When evaluating σv rel aa for calculating the relic density of χ, we set T to its typical value at freeze-out: T F m χ /20 [69]. Due to the temperature dependence, annihilation to aa can be important in determining the thermal relic abundance of χ, but does not lead to observable indirect signals today, where the relic population of χ is nonrelativistic. χχ → ff , on the other hand, is an s-wave process and is therefore relevant both today and in the early Universe.
At our first benchmark point from table I (m a = 12 GeV), it is indeed the interplay of the annihilation processes χχ → aa and χχ → bb that sets the relic density of χ, f χ 0.6. At the second benchmark point (m a = 80 GeV), annihilation to aa is kinematically forbidden at freeze-out, therefore χχ → bb accounts for the relic density f χ 0.33 alone.
In fact, the thermal production of χ has some subtlety to it if the abundance of the heavy species φ is explained by a low reheating temperature T RH . The freeze-out temperature T F of χ is of order T F ∼ m χ /20 ∼ 1.5 GeV at our benchmark points. If T RH T F , the relic abundance Ω χ of χ will be smaller than predicted from the naive estimate for Dirac fermions, Ω χ h 2 ∼ 6 × 10 27 cm 3 sec −1 / σv rel . If T RH T F , the thermal production of χ is not affected. This is possible with a ∼ 10 15 GeV inflaton field with T RH ∼ 10 GeV that could provide the correct relic abundance for φ [67]. For simplicity, we assume in the following that this second case is realized. We moreover assume in the following that φ and χ have comparable relic density, and that together they account for all the DM in the Universe (i.e. f φ + f χ = 1).

V. THE GALACTIC CENTER GAMMA RAY EXCESS
The fact that the light DM species χ in the boosted DM scenario can have a non-negligible relic abundance and a relatively large annihilation cross section to SM fermions in the present day Universe indicates that there may be interesting indirect signatures, in addition to the primary signal from highly boosted χ particles from φ decay.
In particular, the boosted DM scenario can fit the excess of gamma rays which has been observed from the direction of the galactic center at energies of few GeV [33][34][35]. It has been argued that, if the dominant DM annihilation channel is χχ → bb, as in our boosted DM scenario, a 30-40 GeV DM particle with σv rel bb in the range 1.4-2.0 × 10 −26 cm 3 sec −1 provides a good fit to the data. Since in our scenario the light DM species χ constitutes only a fraction f χ of the total DM relic density, its annihilation cross section today has to be correspondingly larger by 1/f 2 χ . At our benchmark points from table I, the predicted annihilation cross sections are σv rel bb ∼ 2.8 × 10 −26 (18 × 10 −26 ) cm 3 /sec. Here the first number stands for the benchmark point with m a = 12 GeV, while the second one (in parenthesis) is for the benchmark point with m a = 80 GeV. With f χ = 0.6 (0.33) (see sec. IV), and taking into account that we chose m χ ∼ 30 GeV at the benchmark points, we thus see that the galactic center gamma ray excess could be explained by our boosted DM scenario. Note that for the special case m a ∼ 2m χ , this could be achieved even for much smaller couplings g Y b and g χ because the annihilation would be resonantly enhanced.

VI. CONSTRAINTS
Constraints on the boosted DM scenario arise on the one hand from indirect DM searches sensitive to high-energy particles from the 3-body decay φ → χχa, followed by decay of the mediator a into SM particles including positrons and gamma rays. We will discuss these possibilities in secs. VI.1 and VI.2, respectively. On the other hand, direct DM searches could hope to directly observe the relic population of light DM particles χ, see sec. VI.3. Finally, the mediator a could be directly produced in accelerator experiments, leading to constraints as well (see sec. VI.4).

VI.1. Positron flux from 3-body decay φ → χχa
The e ± flux at any given point x in the galaxy is given by [61] dΦ e ± (E e , x) where ρ(x) gives the DM density distribution in the galaxy, Γ 3 (φ → χχa) is the 3-body decay rate from eq. (11), E S e is the e ± energy at production, and dN f e ± (E S e )/dE S e is the e ± spectrum at production for a decay to ff . We obtain dN f e ± (E S e )/dE S e in analogy to the secondary neutrino spectrum discussed in sec. III.2 by folding the e ± spectrum in the a rest frame (taken from [61]) with the energy distribution of a particles from φ → χχa (see eq. (11) and fig. 2). The sum in eq. (18) runs over all final states of a decay, and BR(a → ff ) are the corresponding branching ratios. The factor b(E e , x) describes energy loss during propagation [61]. Finally, I(E e , E S e , x) is the generalized halo function, which can be understood as a Green's function of the diffusion-loss equation, describing the probability for an e ± with initial energy E S e to be detected with energy E e . We take the halo function from ref. [61], assuming a Navarro-Frenk-White (NFW) DM density  profile [70] and the MED propagation model [71]. The dependence of our results on the DM density profile is quite small because the dark matter decay rate only depends linearly on the DM density. The uncertainty from the propagation model could change our constraints, but we have checked that even for the propagation model MAX from ref. [71], the predicted flux is at most a factor of 2 larger than for the MED model. In fig. 7, we have plotted the positron flux at Earth from the φ → χχa decay, where a dominantly decays into bb. We fix the mass parameters at our benchmark values m φ = 4.5 PeV (3.9 PeV), m χ = 30 GeV and m a = 12 GeV (80 GeV) in the left panel (right panel). Once the masses are fixed, dΦ e ± /dE e depends on the model parameters through the ratio g 2 χ f φ /τ φ . The background model for the e + flux is taken from refs. [72,73], while the background model for the combined e + + e − flux is taken as a fitting function from ref. [41]. We compare to the AMS-02 e + flux data [39] as well as the Fermi-LAT [40] and H.E.S.S. [41,42] e + + e − flux data to provide a constraint on this decay. Note that when comparing to Fermi-LAT and H.E.S.S. data, which includes both e + and e − , the signal flux is twice the e + signal flux. The error bars in the H.E.S.S. data do not contain systematic uncertainties, while those in the Fermi-LAT and AMS-02 data do. By requiring that the signal flux should be outside the 1σ error bar for any of these data points, we find constraints on the coupling g χ , the relative abundance of the heavy DM f φ , and its lifetime τ φ : We see from table I that our two benchmark points easily satisfy these constraints. The cosmic electron background is complicated and model dependent. The background model from ref. [41] has a lot of parametric freedom regarding in particular the overall normalization, which could alleviate the constraints. Our constraints should therefore be considered as very conservative. Even for the most conservative assumption of zero background, we would still obtain a constraint on g 2 χ f φ /τ φ by requiring that the predicted signal does not significantly overshoot the data. The dominant constraint in this case would come from the last two bins of H.E.S.S. data, and the constraint would be weaker by a factor of ∼ 5 compared to eq. (19). Also including the systematic error of the H.E.S.S. data would make the constraint even weaker.

VI.2. Gamma ray flux from 3-body decay φ → χχa
The secondary gamma ray flux from the decay φ → χχa may contribute to gamma ray searches, in particular to gamma ray searches in the galactic center region and in measurements of the diffuse isotropic gamma ray flux, i.e. the residual flux obtained after subtracting the contribution from known astrophysical sources. We focus here on the diffuse flux because we will see that the strongest limits are coming from air shower detectors located in the northern hemisphere and thus unable to observe the galactic center [74]. The only exception is a γ ray search carried out by the IceCube collaboration using the IceTop array [75]. This search, however, is only sensitive at energies above 1 PeV, where the secondary γ ray flux from decay of ∼ 4 PeV DM particles is already negligible. Moreover, it is worth emphasizing that searching for signals of decaying DM in DM-rich, but also foreground-rich, regions like the galactic center is much less promising than searching for annihilating DM in these regions. The reason is that the DM decay rate depends linearly on the DM density ρ(x), while the annihilation rate scales as ρ(x) 2 .
The procedure for calculating the diffuse gamma ray flux is similar to the one for the secondary neutrino fluxes described in sec. III and for the e ± fluxes described in sec. VI.1. In particular, we can use eqs. (7) and (8) after replacing E χ by the γ energy E γ and the DM spectrum dN χ /dE χ by the gamma ray spectrum at production dN γ /dE γ . Note that dN γ /dE γ must be normalized such that its integral over E γ gives the average number of photons produced in each a decay, accounting for two body decays without photon emission and for three body decays that lead to the radiation of photons. We obtain dN γ /dE γ by boosting the γ ray spectra in the a rest frame (taken from [61]) into the lab frame according to the energy spectrum of a particles given by eq. (11) and fig. 2 and multiplying by BR 3 (φ → χχa) from eq. (13). For the gamma ray flux, also an absorption factor of the form must be included in eq. (8) to describe the attenuation of extragalactic gamma rays on their way from the source to us. We take this factor from ref. [61].
We then obtain the diffuse gamma ray flux conservatively according to the formula [76] Here, dΦ GC /(dE γ dΩ)| minimum denotes the minimum of the differential galactic flux over solid angles, which we take to be the flux from the direction opposite to the galactic center [76]. We have checked that using instead the average of the differential flux over a cone with opening angle 90 • , centered around the direction opposite to the galactic center, would change dΦ GC /(dE γ dΩ)| minimum by O(20%).
We plot the galactic and extragalactic contributions to the diffuse gamma ray flux in fig. 8. We see that the contribution from φ decay in the galaxy dominates over the extragalactic flux due to the attenuation factor eq. (20), which suppresses the extragalactic gamma ray flux.
Note that we neglect the low energy contribution from inverse Compton scattering (ICS) of high-energy e ± from the decay of heavy DM φ on CMB photons, starlight, and light rescattered   Figure 8. The diffuse galactic (solid purple) and extragalactic (dashed brown) gamma ray fluxes from φ → χχa decay, followed by a →bb. The galactic flux is assumed to have in every direction the magnitude it has in the direction opposite to the galactic center [76], evaluated assuming a Navarro-Frenk-White DM density profile [70]. We include only prompt gamma rays, neglecting the low energy contribution from inverse Compton scattering because we have checked that the limit is dominated by the prompt signal. We compare to the Fermi-LAT measurement of the diffuse gamma ray flux from ref. [43], using foreground model C defined in this reference, and to the limits from air shower detectors [77][78][79]. The model parameters are fixed at the values given by our first (second) benchmark point from on dust. We estimate [80,81] that the energy spectrum of ICS photons induced by φ decay peaks at 1-100 GeV. Following [9], we have then estimated that the energy density in ICS gamma rays predicted at our benchmark points is at least one order of magnitude lower than the energy density measured by Fermi-LAT at 1-100 GeV [9,43]. Similarly, also the contribution from bremsstrahlung of e ± on dust is negligible.
To set limits on the parameter space of boosted DM, we compare to the diffuse gamma ray spectra from Fermi-LAT [43] and to the flux limits from the air shower detectors KASCADE [77], GRAPES-3 [78] and GAMMA [79], see also [74]. From fig. 8, we see that the constraint will come mostly from the air shower detectors and the last bin of Fermi-LAT data. By requiring that the predicted signal is smaller than the limit from the air shower detectors, we obtain the constraints We see that both of our benchmark points from table I satisfy these constraint.

VI.3. Direct detection
In the boosted DM scenario, conventional DM direct detection experiments can only constrain the thermally produced population of light DM particles χ, not the population of heavy DM particles φ. The density of φ particles and thus also the flux of boosted χ particles from φ decay are too small to be observed in these detectors. Therefore our discussion of direct detection will focus on the non-relativistic population of the light DM species χ. The cross section for χ-nucleus scattering is [48] dσ dE r = m T 32π where E r is the nuclear recoil energy, v is the DM velocity, Q 2 = 2m T E r ∼ 100 MeV 2 is the 4-momentum transfer squared, m T is the mass of the target nucleus and m N is the nucleon mass. The quantities F N,N Σ are the pseudoscalar form factors of the target nucleus (see e.g. [82]), and the effective nucleon couplings g N , g N depend on the g Y f (see also ref. [49]). For our choice m a 10 GeV, we have m 2 a Q 2 , so that Q 2 is negligible in the denominator. The factor (Q 2 ) 2 in the numerator arises because, in the non-relativistic limit,χγ 5 χ ∝ Q 2 . Direct detection constraints are in general very weak in our boosted DM model due to the (Q 2 ) 2 suppression unless the mediator mass m a is extremely small. The resulting limit on g χ g Y f is therefore much weaker than the value needed by the thermal relic density [48].
Departing for a moment from our toy model with a pseudoscalar mediator, we note that in general, boosted DM models with interaction cross sections strong enough to explain the IceCube events would also lead to a large signal in direct detection experiments. From a model building point of view, there are several ways of circumventing this, other than using a pseudoscalar coupling as in our toy model. (1) Construct a model in which the scattering of the light DM particles χ on nuclei is inelastic [83]. If the mass splitting δm between the ground state of χ and the excited state χ * which is produced in the scattering is sufficiently large, it will lead to vanishing event rates in direct searches, but will have no influence on boosted DM collisions as long as δm is small compared to the energy of the boosted DM particles. (2) Assume the relic abundance of the light DM species is sufficiently low to avoid direct detection limits. This would of course preclude a simultaneous explanation of the IceCube events and the galactic center gamma ray excess. (3) Choose the light DM mass smaller than ∼ 3 GeV, below the energy threshold for direct detection. This would also preclude an explanation of the galactic center gamma ray excess.

VI.4. Constraints from flavor physics experiments and from collider searches
In the following, we discuss constraints on our boosted DM scenario from experiments at flavor factories and at high energy colliders and indicate for each constraint to which of the three renormalizable models from sec. II it applies.
A large number of constraints arises from Kaon and B meson decays [48]. Searches are sensitive to the production of the pseudoscalar a in decays of these mesons if a subsequently decays to leptons, photons or invisible particles. Since we are considering the case m a 10 GeV, those constraints are, however, significantly weakened by the fact that a would have to be off-shell.
B s → µ + µ − is the only search channel sensitive to an off-shell pseudoscalar. If we consider a renormalizable model for the pseudoscalar a in the framework of a Two Higgs Doublet Model, as in the MSSM-like and Flipped models from sec. II, a couples to the SM by mixing with the heavy pseudoscalar A 0 . The mixing angle is denoted by θ. The branching ratio for B s → µ + µ − in the MSSM-like model is given in ref. [50,84]. The contribution from a to the amplitude is proportional to tan 2 β sin 2 θ. The constraint for m a ∼ 10 GeV is tan β sin θ = 2g Y d g Yµ 0.4 (0.51) for charged Higgs boson masses of m H ± ∼ 800 (400) GeV, while the constraint for m a ∼ 80 GeV is about tan β sin θ = 2g Y d g Yµ 3.8 (4.8) [50]. For the Flipped model, where lepton couplings are proportional to cot β, the amplitude from a exchange is proportional to tan β cot β sin 2 θ = sin 2 θ. Therefore, the constraint is sin θ = 2g Y d g Yµ 0.4 (0.51) for charged Higgs boson masses of m H ± ∼ 800 (400) GeV when m a ∼ 10 GeV. For m a ∼ 80 GeV, there is no constraint on sin θ.
Note that the B s → µ + µ − constraint does not apply to the Vector-quark model because a does not couple to leptons in this model.
An additional constraint, which is independent of the couplings of the pseudoscalar a to fermions, arises from the exotic decay h → aa. In the context of the MSSM-like and Flipped models, the branching ratio for this decay is constrained by [50,85,86] BR(h → aa) 0.02 m A 800 GeV where m A is the mass of the heavy pseudoscalar. If m A 800 GeV, sin θ has to be smaller than 0.02. If m a becomes comparable to m h /2, the above constraint is weakened, and for m a > m h /2 it is completely absent. It is also absent in the Vector-quark model. We should also consider constraints from the LEP experiments, which have searched for e + e − → hA 0 , where A 0 is the pseudoscalar Higgs boson appearing in the MSSM [87]. While these searches exclude A 0 masses below 90 GeV, they do not apply to models with an extra pseudoscalar a, like the scenarios we are considering here [88].
If a is heavy enough to decay to χχ, ref. [57] shows that searches for b jets and missing energy can provide an excellent constraint on the pseudoscalar a. The dominant processes are gg → bba with a decaying to χχ subsequently. The current CMS and ATLAS searches [89,90], which are optimized for final states with two b quarks, lead to the constraint for m a ∼ 100-250 GeV and assuming [57]. If g χ is significantly larger than g Y b √ 2m b /v, the limit will become somewhat weaker since the probability for radiating an on-shell a particle changes [57].
In the intermediate mass region 20-80 GeV, ref. [88] also discusses the processes gg → bba and ( -) b g → ( -) b a, but considering the subsequent decays a → µ + µ − and a → τ + τ − . By looking for these leptonic final states, the high luminosity LHC can be sensitive to g Y b ∼ 7 with 100 fb −1 of integrated luminosity, assuming that g Y f is universal for down type quarks and charged leptons (as in the MSSM-like model). Since this assumption is not satisfied in the Flipped model, which has suppressed couplings of a to leptons, and in the Vector-quark model, in which a does not couple to leptons at tree level, the constraint would be significantly weaker or completely absent in these models.
In the light mass region m a ∼ 5.5-14 GeV, CMS has searched for a → µ + µ − in the context of the Next-to-Minimal Supersymmetric Standard Model [91]. The upper limit on the cross section for the process pp → a → µ + µ − is around 2-4 pb. This translates into a constraint of g Y d ∼ 2 in the MSSM-like model, where g Y = g Y d [88]. The Flipped and Vector-quark models are not restricted by this constraint due to the smallness or complete absence of leptonic couplings of a.
Let us summarize the most stringent constraints for the three models defined in sec. II (see also the last column in table II below). For the MSSM-like model, the most stringent limit comes from B s → µ + µ − . It rules out the MSSM-like model as a UV-completion for our m a = 12 GeV benchmark point, while for the m a = 80 GeV benchmark point, it is a viable possibility. For the Flipped model, the coupling between leptons and the pseudoscalar a is suppressed once we are in the large tan β region. But the constraint from h → aa still implies that the mixing angle sin θ between a and the heavy pseudoscalar A 0 should be very small. If we require that tan β 50, this disfavored also the Flipped model as a UV completion for our m a = 12 GeV benchmark point. At m a = 80 GeV, the h → aa constraint is absent because the decay is kinematically forbidden. For the Vector-quark model, only the perturbativity of the Yukawa couplings involving the heavy quarks, together with the LHC limits on their mass, imposes a very weak constraint g Y b 20 [57]. IceCube galactic center e ± diffuse γ Lab Vector-quark 20 Table II. Summary of constraints on the boosted DM scenario for two different benchmark values for the mass m a of the pseudoscalar that mediates interactions between the light DM species χ and the SM. Since IceCube sees both the scattering of highly boosted χ particles and secondary neutrinos from φ → χχ + (a → bb) decay, the experiment constrains two independent combinations of the pseudoscalar couplings to DM (g χ ) and b quarks (g Y b ), the lifetime of the heavy DM particle, τ φ , and its fractional abundance in the Universe, f φ . Note that we always assume here that a couplings to SM fermions other than the b quark are negligible.
Requiring that the galactic center gamma ray excess can be explained constrains the light DM mass m χ and an additional combination of coupling constants. Further constraints come from secondary e ± and γ rays from φ → χχ + (a → bb) and from laboratory searches for the pseudoscalar mediator a.

VII. SUMMARY AND CONCLUSION
In summary, we have discussed the possibility that the high energy event excess observed by the IceCube collaboration is explained by the scattering of highly boosted DM particles on atomic nuclei in the detector. We have constructed a simple toy model in which a DM particle φ with a mass of order PeV can decay into a much lighter DM species χ. The χ particles, in turn, interact with atomic nuclei through a t-channel mediator a, thus explaining the IceCube signal at PeV energies.
The experimental constraints on this toy model are summarized in table II for two different benchmark values of the pseudoscalar mass m a . At both benchmark points, we have assumed that the mediator a has significant coupling to b quarks, while its couplings to light quarks and to leptons are suppressed. This is naturally realized in UV-complete models with either an extended Higgs sector or with the introduction of vector-like quarks (see sec. II). The highest energy events in IceCube set the scale for the heavy DM mass m φ and the normalization of the scattering cross section. At lower energy, IceCube is sensitive to the secondary neutrino flux from the 3-body decay φ → χχ + (a → bb), see fig. 5. This provides a constraint on the branching ratio for this decay. Since the same decay also leads to secondary electron/positron and gamma ray fluxes, e ± and γ ray data from AMS-02, Fermi-LAT, HESS and several air shower arrays provide a constraint on its branching ratio as well. Moreover, the boosted DM scenario is constrained by searches for the new pseudoscalar particle a in flavor physics experiments and at high energy colliders.
We have shown that, besides explaining a population of high energy events in IceCube, the boosted DM scenario can simultaneously also account for the gamma ray excess observed in Fermi-LAT data from the direction of the galactic center. This is possible because the light DM species χ can have a non-negligible thermally produced relic abundance, and can annihilate in the Milky Way today. Fermi-LAT data then identifies a preferred range for the light DM mass m χ and its couplings to ordinary matter.
The boosted DM scenario shares some features with interpretations of the IceCube data in terms of DM decay directly to SM particles, including neutrinos. First, the morphology of the signal is similar in the two cases, with a mild peak expected in the galactic center region. Moreover, it is worth mentioning that the most recent IceCube data [32] provides a mild hint at a bump-like feature at ∼ 80 TeV from the southern sky. Since this is where the galactic center is located, such a bump could be explained by the secondary neutrino flux in the boosted DM scenario. The second common feature between boosted DM and more conventional decaying DM explanations of the IceCube data is the rapid drop of the signal at energies larger than half of the heavy DM mass. With more statistics collected, these features can help to distinguish DM interpretations of the IceCube data from an interpretation in terms of isotropic astrophysical neutrino emission.
A unique feature of the boosted DM scenario is the prediction that, at PeV energies, where the IceCube signal is explained by scattering of boosted DM particles, only shower-like events should be observed. At lower energies ∼ 100 TeV, however, where the secondary neutrino flux from the 3-body decay φ → χχ + (a → bb) contributes, both shower and track-like events are predicted, with a ratio very similar to the one expected from astrophysical neutrino sources. Between the two populations of events, a mild dip in the energy spectrum is predicted. These features distinguish the boosted DM scenario from astrophysical explanations of the IceCube data and from interpretations in terms of neutrinos from DM decay.  Figure 9. Effective target mass for neutrino interactions in IceCube as a function of the incoming neutrino energy E ν for CC ν e (red solid) and NC (gray solid) interactions [1]. The dashed brown line shows our prediction for the effective target mass in the NC case from eq. (A2), which is in excellent agreement with the results from [1].