Boosted Dark Matter and its implications for the features in IceCube HESE data

Boosted Dark Matter and its implications for the features in IceCube HESE data


We study the implications of the premise that any new, relativistic, highly energetic neutral particle that interacts with quarks and gluons would create cascade-like events in the IceCube (IC) detector. Such events would be observationally indistinguishable from neutral current deep-inelastic (DIS) scattering events due to neutrinos. Consequently, one reason for deviations, breaks or excesses in the expected astrophysical power-law neutrino spectrum could be the flux of such a particle. Motivated by features in the recent 1347-day IceCube high energy starting event (HESE) data, we focus on particular boosted dark matter () related realizations of this premise. Here, is assumed to be much lighter than, and the result of, the slow decay of a massive scalar () which constitutes a major fraction of the Universe’s dark matter (DM). We show that this hypothesis, coupled with a standard power-law astrophysical neutrino flux is capable of providing very good fits to the present data, along with a possible explanation of other features in the HESE sample. These features include a) the paucity of events beyond PeV b) a spectral feature resembling a dip or a spectral change in the 400 TeV–1 PeV region and c) an excess in the TeV region. We consider two different boosted DM scenarios, and determine the allowed mass ranges and couplings for four different types of mediators (scalar, pseudoscalar, vector and axial-vector) which could connect the standard and dark sectors.We consider constraints from gamma-ray observations and collider searches. We find that the gamma-ray observations provide the most restrictive constraints, disfavouring the allowed parameter space from IC fits, while still being consistent with the allowed region. We also test our proposal and its implications against the (statistically independent) sample of six year through-going muon track data recently released by IceCube.

IceCube, Dark matter, Ultra-high energy neutrinos


1 Introduction and Motivation

In this section, we shall begin with a summary of the -day IceCube (IC) high-energy starting event (HESE) neutrino data, focussing on events with deposited energies greater than around TeV, and discuss some of its features, especially those that are of particular interest for this study. We shall then introduce two possible scenarios of boosted dark matter, which, in combination with a power-law astrophysical flux, can provide a good fit to these features.

1.1 IceCube High Energy Starting Events (HESE) and features of the 1347-day data

The observation of 54 HESE (i.e., events with their interaction vertices inside the detector) Aartsen:2015zva (); 1311.5238 (), with deposited energies between 30 TeV to a maximum energy of 2.1 PeV by the IceCube experiment (IC) has opened an unprecedented window to the universe at high energies.1 The data constitute an approximately signal in favour of a non-atmospheric and extra-terrestrial origin of the events.2 It is generally believed, but not conclusively known, that the highest energy cosmic rays ( GeV), for which observations now extend to GeV, and ultra-high energy (UHE) neutrinos with energies greater than TeV, share common origins and are produced by the same cosmic accelerators. The specific nature of these accelerators, however, remains unknown, although over the years, anticipating their detection, several classes of highly energetic cosmic astrophysical sources have been studied as possible origins of these particles. For general discussions of this topic, we refer the reader to 0710.1557 (); astro-ph/0204527 (); 1202.0466 (); 1301.1703 (); 1305.4123 (); 1306.2309 (); 1306.5021 (); 1312.6587 (); 1410.3680 (); 1407.7536 (); 1405.5487 ().

Subsequently, based on the recent IC data, many authors have considered a host of source classes and possibilities for explaining both the origin and some emerging spectral features in the IC data. These efforts have been motivated, at least in part, by evidence that the data, to an extent, diverge from expectations. The considered candidate sources include gamma-ray bursts 1211.1974 (); 1212.1260 (); 1307.7596 (); 1306.2274 (); 1310.7061 (); 1405.2091 (); 1405.5487 (); 1411.7491 (); 1504.00107 (); 1512.08513 (); 1512.01559 (), star-burst galaxies 1306.3417 (); 1405.7648 (); 1404.1189 (); 1406.1099 (); 1509.00983 (), active galactic nuclei 1303.0300 (); 1305.7404 (); 1406.2633 (); 1406.0645 (); 1407.0907 (); 1408.3664 (); 1410.8124 (); 1411.2783 (); 1411.3588 (); 1501.07115 (), remnants of hyper-novae 1310.1263 () and of supernovae 1501.02615 (), slow-jet supernovae 1407.2985 (), microquasars 1410.0348 (), neutron star mergers 1306.3006 (), blackholes 1512.08596 (), cosmic-ray interactions 1307.2158 (); 1310.5123 (); 1311.0287 (); 1404.6237 (); 1405.3797 (); 1411.6457 (); 1412.8590 (), the galactic halo 1403.3206 (), galaxy clusters 1410.8697 (), dark matter decay 1303.7320 (); 1308.1105 (); 1312.3501 (); 1410.5979 (); 1403.1862 (); 1407.3280 (); 1408.1745 (); 1411.1071 (); 1503.02669 (); 1503.04663 (); 1505.06486 (); 1506.08788 (); Boucenna:2015tra (); 1508.02500 (); 1601.02934 (); 1606.04517 (); 1607.05283 (), and exotic particles, processes or possibilities 1305.6907 (); 1402.1681 (); 1402.6678 (); 1404.0622 (); 1404.2279 (); 1404.2288 (); 1404.2932 (); 1404.7025 (); 1407.2848 (); 1409.4180 (); 1409.5896 (); 1410.0408 (); 1410.3208 (); 1411.5318 (); 1411.5889 (); 1507.03015 (); 1507.03193 (); 1606.06238 (); 1606.07903 (); Dev:2016uxj ().

It is generally accepted, however, that the charged particles in a source which link the acceleration of cosmic-rays to the acceleration of astrophysical neutrinos attain their high energies via Fermi shock acceleration Fermi:1949ee (), and as a generic consequence, the neutrinos resulting from them are expected to follow a spectrum 0710.1557 (); astro-ph/0204527 (). Some variation from this general spectral behaviour may occur, however, depending on the details of the source, as discussed, for instance, in 2005PhRvL..95r1101K ().

IceCube is sensitive to high energy neutrinos via their electroweak charge and neutral current (CC and NC respectively) deep inelastic (DIS) interactions with nucleons in ice, which result in the deposition of detectable energy in the form of Cerenkov radiation. An event may thus be classified as3

  • a track, produced by  CC and a subset of  CC interactions (where a produced decays to a ) , characterized by a highly energetic charged lepton traversing a significant length of the detector, or

  • a cascade, produced by either i)  CC interactions, ii) a subset of  CC interactions or iii) NC interactions of all three flavours. Cascades are characterized by their light deposition originating from charged hadrons and leptons, distributed around the interaction vertex in an approximately spherically shaped signature.

Additionally, because neutrino production in astrophysical sources stems from photohadronic interactions producing light mesons, such as pions and kaons, and to a smaller degree, some heavier charmed mesons, including , and their subsequent decays, the flux ratio at source is expected to be . However, standard oscillations between the three flavours over cosmological distances renders this ratio close to 0505017 () by the time they arrive at earth. In this situation, cascade events are expected to constitute about 75–80% of the total observed sample Beacom:2004jb (). The background to the HESE events is provided by the rapidly falling atmospheric neutrino flux and the muons created in cosmic-ray showers in the atmosphere.

We now describe the significant features of the HESE data, some of which are fairly firm even at the present level of statistics, and others which, while interesting and suggestive, are emergent and need further confirmation via more observations before they can be considered as established. (We note that the energies quoted below refer to those deposited by the primary in IC.)

  • The data, to a high level of significance (about , as mentioned earlier), indicate that above a few tens of TeV, the sources of the events are primarily non-atmospheric and extra-terrestrial in nature.

  • Due to the lack of multi-PeV events, including those from the Glashow Resonance Glashow:1960zz (); Bhattacharya:2011qu (); Barger:2014iua () in the range 6-10 PeV, a single power-law fit to the flux underlying the observed events now disfavours the expected spectral index from Fermi shock acceleration considerations, , by more than . Indeed, for an assumed spectrum, and with the corresponding best-fit normalization to the flux, about 3 additional cascade events are expected between 2 PeV and 10 PeV, largely due to the expected presence of the Glashow resonance. However, in spite of IceCube’s high sensitivity at these energies, none have been observed thus far. The present best fit value of is consequently significantly steeper, being around Aartsen:2015zva (); halzen_pheno16 ().

  • The data, when subjected to directional analyses Aartsen:2015zva (); 1309.2756 (); 1309.4077 (); 1311.5864 (); 1311.7188 (); 1406.0376 (); 1406.2160 (); 1407.2243 (); 1410.5979 (); 1501.05158 (); 1507.05711 (); 1509.00517 (); 1509.03522 (); 1510.00048 (); 1511.09408 (); 1603.06733 (), at its present level of statistics, is compatible with an isotropic diffuse flux, although several studies among the ones cited above indicate the presence of a small galactic bias. The accumulation of more data will be able to ascertain whether the galactic bias is real, in which case it would imply important (and possibly new) underlying physics.

  • The three highest energy events Aartsen:2015zva (), with the estimated (central value) of the deposited energies of PeV, PeV and PeV are all cascade events from the southern hemisphere. At these energies, i.e.  PeV, the earth becomes opaque to neutrinos, thus filtering out neutrinos coming from the northern hemisphere.

  • Below PeV, there appears to be a dip in the spectrum, with no cascade events between roughly TeV and PeV.4

  • At lower energies, in the approximate range of TeV, there appears to be an excess, with a bump-like feature (compared to a simple power-law spectrum), which is primarily present in events from the southern hemisphere 1410.1749 (). The maximum local significance of this excess is about , which is obtained when the lowest estimates for the conventional atmospheric neutrino background is adopted, with the prompt component of the background assumed to be negligible MESE ().

  • Finally, and importantly, the data when interpreted as being due to a single astrophysical power-law neutrino flux, appears to require an unusually high normalization for this flux, which is at the level of the Waxman-Bahcall (WB) bound  Waxman:1998yy (); Bahcall:1999yr () for neutrino fluxes from optically thin sources of high energy cosmic rays and neutrinos. This is an aspect that is difficult to understand within the confines of the standard interpretive mechanism, which connects ultra-high energy neutrino fluxes to observations of the highest energy cosmic-rays 5.

1.2 Deep Inelastic Scattering of Boosted Dark Matter in IceCube

As proposed in 1407.3280 (), if there is a source of long-lived, highly relativistic and energetic neutral particles in the present Universe which can interact with quarks or gluons, the signal produced by them in IceCube would, in all likelihood, be indistinguishable from the NC DIS cascade of a neutrino primary. To the extent that the astrophysical neutrino flux is expected to follow a simple power-law behaviour, one could argue that features in the HESE data (as described in the previous subsection) which deviate from this, such as statistically significant excesses, spectral breaks or line-like features, could indicate the presence of such a particle6. Although there are strong constraints on the presence of additional relativistic degrees of freedom during the epochs of recombination and big bang nucleosynthesis, such particles might be injected at later times by the slow decay of a heavy particle, which, overall, is the approach we adopt here.

We consider the case where this heavy particle constitutes a significant part of the dark matter (DM) density of the Universe. Its late-time decay produces a highly energetic flux of light dark matter (LDM) particles, which can then give rise to a subset of the NC DIS events at IC. We note that this is different from the scenario where the heavy dark matter (HDM) particle directly decays to standard model particles, leading to a neutrino flux in IC, as discussed in, for instance 1303.7320 (); 1308.1105 (); 1312.3501 (); 1410.5979 (); 1403.1862 (); 1408.1745 (); 1411.1071 (); 1503.04663 (); 1505.06486 (); 1506.08788 (); 1507.01000 (); 1508.02500 (); 1601.02934 (); 1606.04517 (); 1607.05283 (). In the scenario(s) discussed here, in order to have NC DIS scattering with nuclei, the LDM particles need to couple to the SM quarks (or gluons) with appropriate strength. It is then possible that these interactions could keep them in chemical equilibrium with the SM sector in the early Universe. Thus, the standard thermal freeze-out mechanism will give rise to a relic density of the LDM particles as well in the present Universe, though the exact value of their present-day density would in general depend upon all the annihilation modes open and the corresponding annihilation rates. It is important to note that the couplings relevant for the IC analysis provide only a lower bound on the total annihilation rate. For our purpose, the precise relic density of LDM is not of particular relevance, and we simply need to ensure that it annihilates sufficiently fast in order not to overclose the Universe, while its relic abundance should not be too high, in order to allow for a sufficient HDM presence in the universe. The latter is needed to produce enough of the relativistic LDM flux from its late time decays. In other words, scenarios where the LDM abundance is small are preferred but not required. Similarly, for phenomenological analysis of the IC data, the production mechanism of the HDM particle does not play any essential role. Therefore, we abstain from discussing specific cosmological models for HDM production in this article, and instead refer the reader to possibilities discussed in Refs. 1402.2846 (); 1201.3696 (); hep-ph/0203118 (); hep-ph/0205246 (). We further note that general considerations of partial-wave unitarity of scattering amplitudes imply an upper bound on the mass of any DM particle that participates in standard thermal equilibrium production processes and then freezes out. Such a particle should be lighter than a few hundred TeV, as discussed in PhysRevLett.64.615 (). As we shall see, the HDM under consideration here is necessarily non-thermal due to this reason7.

In what follows, we pursue two specific realizations (labelled Scenario I and II below) of such a dark matter sector, which, in combination with a power-law astrophysical component, provide a good description to the features in the IC data described in the previous subsection. For each realization, we perform a likelihood analysis to fit the IC HESE data and its observed features, in terms of a combination of four distinct fluxes. These fluxes are:

  1. Flux-1: An underlying power-law flux of astrophysical neutrinos, , whose normalization () and index () are left free.

  2. Flux-2: A flux of boosted light dark matter (LDM) particles (), which results from the late-time decay of a heavy dark matter (HDM) particle (). When is much lighter than , its scattering in IC resembles the NC DIS scattering of an energetic neutrino, giving rise to cascade-like events.

  3. Flux-3: The flux of secondary neutrinos resulting from three-body decay of the HDM, where a mediator particle is radiated off a daughter LDM particle. The mediator then subsequently decays to SM particles, producing neutrinos down the decay chain. Since the NC DIS scattering that results from Flux-2 requires a mediator particle which couples to both the LDM and the SM quarks, such a secondary neutrino flux is always present.

  4. Flux-4: The conventional, fixed, and well-understood, atmospheric neutrino and muon background flux, which is adapted from IC analyses Aartsen:2015zva (); 1311.5238 ().

Scenario I : PeV events originating from DIS scattering of boosted LDM at IC

In Scenario I, the three highest energy PeV events, which are cascades characterized by energy depositions (central values) of PeV, PeV and PeV, are assumed to be due to Flux-2 above, requiring an HDM mass of PeV. Both Flux-1 and Flux-3 contribute to account for rest of the HESE events, including the small bump-like excess in the TeV range. This scenario, in a natural manner, allows for the presence of a gap, or break in the spectrum between 400 TeV to 1 PeV8.

A similar scenario has previously been studied in Refs. 1407.3280 (); 1503.02669 (), in which the 988-day HESE data were taken into account. While Ref. 1407.3280 () ascribed the events below a PeV upto tens of TeV entirely to the astrophysical flux (Flux-1), Ref. 1503.02669 (), ascribed these as being generated by the secondary neutrino flux from three-body HDM decay (Flux-3). In this study we do not make any assumption regarding the specific origin of these sub-PeV events, and allow any viable combination of Flux-1 and Flux-3 in the fitting procedure. As we shall see later, one of our main results from the fit to the HESE data within Scenario I is that with the current level of statistics, a broad range of combinations of Flux-1 and Flux-3 can fit the sub-PeV events, while the PeV events are explained by Flux-2. We note in passing that, in Ref. 1503.02669 () the DM model parameter space was guided by the requirement that the LDM annihilation in the present Universe explain the diffuse gamma ray excess observed from the Galactic centre region TheFermi-LAT:2015kwa () in the Fermi-LAT data. In the present study, the focus is entirely on satisfactorily fitting the IC events.

Scenario II : PeV events from an astrophysical flux and the TeV excess from LDM DIS scattering

In Scenario II, we relax the assumption made regarding the origin of the three PeV events in Scenario I, and perform a completely general fit to both the PeV and the sub-PeV HESE data, with all four of the flux components taken together. Essentially, this implies that the mass of the HDM particle is now kept floating in the fit as well. We find that both the best-fit scenario and the statistically favoured regions correspond to a case where the PeV events are explained by the astrophysical neutrino flux (Flux-1), while the excess in the TeV window primarily stems from the LDM scattering (Flux-2). Flux-3, which now populates the low 1–10 TeV bins becomes inconsequential to the fit, since the IC threshold for the HESE events is 30 TeV. Expectedly, in order for the astrophysical flux to account for the PeV events, the slope of the underlying power-law spectrum in Scenario II is significantly flatter compared to that in Scenario I.

In addition to performing general fits to the PeV and sub-PeV HESE data as described above, we also explore, for both Scenarios I and II, the extent to which different Lorentz structures of the LDM coupling with the SM quarks impact the results. While a vector mediator coupling to the SM quarks and the LDM was considered in Ref. 1407.3280 (), a pseudo-scalar mediator was employed in Ref. 1503.02669 (). Adopting a more general approach, we consider scalar, pseudo-scalar, vector and axial-vector mediators. However, we find (expectedly) that if the LDM relic density is appreciable, strong limits on the spin-independent coherent elastic scattering cross-section with nuclei of the relic LDM component come into play and restrict the available parameter space for scalar and vector mediators. There are also interesting differences between the pseudo-scalar and axial-vector scenarios insofar as fitting the IC data, as we shall show in later sections.

Finally, as emphasized in Ref. 1503.02669 (), the three-body decay of the HDM particles that gives rise to the secondary neutrino flux (Flux-3 above), also produces a flux of diffuse gamma-rays in a broad energy range, which is constrained from the measurements by the Fermi-LAT telescope Ackermann:2014usa () at lower energies, and by the cosmic ray air shower experiments (KASCADE Feng:2015dye () and GRAPES-3 Gupta:2009zz ()) at higher energies Ahlers:2013xia (). We find that the parameter space of the proposed dark matter scenarios that can fit the IC data is significantly constrained by the upper bounds on residual diffuse gamma ray fluxes9.

The rest of this paper is organized as follows: Sec. 2 examines the different ways the LDM particle can interact with SM quarks, and summarizes the current constraints on the effective couplings and the mass parameters, using gamma ray and collider data. We also discuss the general method used to calculate the contribution made by the HDM three-body decay to galactic and extra-galactic gamma-ray fluxes. Sec. 3 focusses on Scenario I and describes our procedure for deriving best-fits to the observed IC HESE data for it, and the results obtained for different choices of the mediator. The validity of these results is then examined in the light of various constraints. Similarly, Sec. 4 repeats this for Scenario II. Although the focus of this work is on understanding the HESE data, IC has recently released a statistically independent sample of high energy muon track events 1607.08006 () for the neutrino energies between TeV to PeV, where the interaction vertex is allowed to be outside the detector. Sec. 5, examines both the scenarios considered here in the light of this data sample. Finally, our findings are recapitulated and summarized in Sec. 6.

2 LDM interaction with quarks: simplified models and current constraints

This section provides further details on how we model the interaction of the LDM with the SM quarks. In what follows, we shall work with a representative model where the HDM () is described by a real scalar field, and the LDM () is a neutral Dirac fermion, both of which are singlets under the standard model gauge interactions. The interaction of the heavy dark matter particle with the LDMs is described by an Yukawa term of the form .

We further assume that the LDM particles are stabilized on the cosmological scale by imposing a symmetry, under which the LDM field is odd, and all other fields are even. The LDM can interact with the SM fermions (quarks in particular) via scalar, pseudo-scalar, vector, axial-vector or tensor effective interactions. To describe such effective interactions we introduce a simplified model, where the interactions are mediated by a even spin-0 or spin-1 particle. The LDM can also couple to SM fermions via a odd mediator, which carries the quantum numbers of the SM fermion it couples to. We do not consider the t-channel models or the tensor type interaction in this study.

In the following sub-sections we shall describe the simplified model setup and mention the generic constraints on the couplings of a spin-0 or spin-1 mediator to the LDM and the SM fermions. Such constraints on the coupling and mass parameters can be modified within the context of a specific UV complete scenario, especially if it necessarily involves other light degrees of freedom not included in the simplified model. However, since the primary focus of this study is to determine the combination of different fluxes which can fit the features observed in the IC data, the simplified models chosen are sufficient for this purpose. Our approach allows us to draw general conclusions regarding the possible contributions of LDM scattering and the secondary neutrino fluxes, while being broadly consistent with constraints from experiments and observations.

2.1 Spin-0 mediators

The parity-conserving effective interaction Lagrangian (after electroweak symmetry breaking) of the LDM with SM fermions , involving a scalar mediator or a pseudo-scalar mediator can be written as follows:


Here is the mass of the SM fermion , () represents the coupling of the LDM with the scalar (pseudoscalar) mediator, and GeV) stands for the vacuum expectation value of the SM Higgs doublet (in the presence of other sources of electroweak symmetry breaking the definition of will be appropriately modified). The sum over fermion flavours can in principle include all SM quarks and leptons, although for our current study, the quark couplings are more relevant. We shall take the coupling factors and , which appear in the coupling of fermion flavour with the scalar and the pseudo-scalar mediators respectively, to be independent of the quark flavour for simplicity.

A SM singlet spin-0 mediator cannot couple in a gauge-invariant way to SM fermion pairs via dimension-four operators. One way to introduce such a coupling is via mixing with the neutral SM-like Higgs boson after electroweak symmetry breaking. Such a mixing, if substantial, can however modify the SM-like Higgs properties leading to strong constraints from current LHC data. Other possible ways include introducing a two Higgs doublet model (and mixing of the singlet scalar with the additional neutral scalar boson(s)), or introducing new vector-like fermions to which the singlet scalar couples, and which in turn can mix with the SM fermions Izaguirre:2014vva (). In all such cases the couplings of the singlet-like scalar to SM fermions should be proportional to the fermion Yukawa couplings in order to be consistent with the assumption of minimal flavour violation, thus avoiding flavour-changing neutral current (FCNC) constraints D’Ambrosio:2002ex ().

2.2 Spin-1 mediators

The effective interaction Lagrangian involving a spin-1 mediator, , to SM fermions and the LDM can be written as follows:


Here the subscripts , and refer to vector, axial-vector, left-chiral and right-chiral couplings respectively. The left and right handed SM fermion currents are invariant under the SM gauge transformations. Therefore, in general, both vector and axial-vector interactions are present with coefficients and . In order to obtain only vector or axial-vector SM fermion currents at a low energy scale, we need to set or , respectively.

If the couples to charged leptons, there are strong upper bounds on its mass from collider searches for dilepton resonances from the LHC. In order to avoid them, we assume the leptonic couplings to be absent. In a minimal scenario with only the SM Higgs doublet giving mass to all the SM fermions, we encounter further relations from gauge invariance (here, is the gauge field corresponding to the gauge interaction) on the coupling coefficients to quarks and leptons Kahlhoefer:2015bea (). This is because if left and right handed SM fermions have different charges under the new gauge group, the SM Higgs doublet needs to be charged under as well. Thus, when a single Higgs doublet gives rise to the mass of both SM quarks and charged leptons, if the quarks are charged under , so would be the leptons. However, such constraints can be avoided in a non-minimal scenario, for example in a two Higgs doublet model, where different Higgs bosons are responsible for giving mass to quarks and leptons, thereby making their charges uncorrelated. We keep in view such considerations related to ultra-violet completion for this study, although we do not fully flesh out their consequences.

2.3 Constraints on the couplings and the mass parameters

Figure 1: The interactions corresponding to decay (left), mediator decay (centre) and scattering (right) involving a generic mediator, along with relevant coupling constants.

Figure 1 shows the main interaction vertices which are relevant for both Scenario I and II. represents the coupling between the HDM and LDM leading to the slow decay of the former, with lifetime . The other couplings shown correspond to the vertices of either (a) SM quarks or (b) the LDM interacting with a generic mediator, which can be a pseudo-scalar () or a scalar () or a spin-1 boson () which couples via vector and/or axial-vector couplings, as discussed in the previous section.

The rate of LDM DIS scattering at IC is proportional to , where and are the mediator-quark pair and mediator-LDM pair couplings, respectively. It is also proportional to , or, equivalently, inversely proportional to 10. Finally, the IC event rates are also proportional to the fractional contribution of the HDM to the total DM density, . Here, (with being the normalized Hubble constant) from recent PLANCK results1502.01589 ().

From the above considerations, the IC event rate from LDM DIS scattering, for a given choice of mediator mass , is determined by the quantity . It is useful to determine its maximum allowed value. In order to keep the couplings perturbative, we require . We also require the lifetime of the HDM to be longer than the age of the Universe seconds. And since , we obtain the upper bound, . If the value of exceeds this maximum, the couplings will not be perturbative, or the HDM would have decayed too quickly to have an appreciable density in the present Universe.

The secondary neutrino flux from the three-body decay of (Flux-3 in Sec. 1.2), is proportional to (again, in the limit where the two-body decay width is much larger than the three-body width). It is also inversely proportional to the life-time of the HDM, . In addition to the mass of the , is determined by when the two body decay to LDM pairs dominate. Thus, the parameters relevant for fitting the features in the IC data in our work are , mass of the mediator particle (), and . The results do not depend on , as long as it is significantly lower than .

It is useful to examine the ball-park numerical values of some of the quantities which are used to fit the IC events using DIS -nucleon scattering. The cross section depends essentially on and the mediator mass, . Hence, given a certain value of , and a value for the factor , one could obtain minimum value of the couplings needed to fit an observed number of cascade events. This is given by , assuming , and seconds. A typical value that occurs in the fits is, for instance, , and using this leads to a lower bound . Assuming, for simplicity, , each coupling should thus be greater than about .

As mentioned earlier, (in Sec. 1.2), the most restrictive constraint on the value of comes from the upper bound on the flux of diffuse gamma rays. We defer a detailed discussion of our computation of the gamma ray flux from the three-body decay of the HDM, and the resulting constraints to Sec. 2.3.2. Significant constraints also arise from collider experiments, where the mediator and the LDM particles can be directly produced, and we discuss these in the next sub-section.

The relic density of , which we denote as , is not of direct relevance to our study, which focusses on the IC events coming either from DIS scattering of the LDM in IC, and on the flux of secondary neutrinos from the three-body decay of the HDM. However, direct detection constraints can be important if there is a significant density of the LDM in the current Universe. It is well-known that if is significant, the spin-independent direct detection bounds on the scalar and vector interactions are very strong, and thus would force us to focus on either pseudo-scalar or axial-vector couplings (or relegate us to corners of values which are not yet probed by the direct detection experiments). For our purpose, we could either assume that this is the case, or, equivalently, that the density is indeed small. If the latter, within the simplified model setup discussed above, the relic density of can be diluted to very small values in two ways. The first is by increasing , and restricting to values of , such that the dominant annihilation mode of is to the mediator pair, which can then decay to the SM fermions even via a small . The second way (albeit fine-tuned), is by setting close to , thereby allowing for a resonant annihilation of LDM pairs to SM quarks. Since the IC event rates do not depend upon as long as it is significantly smaller than the HDM mass, both these approaches do not affect the IC event rates. Finally, there can always be additional annihilation modes of the LDM not described by the simplified models which do not affect the IC computations, but help make small.

With respect to the choices of mediators, we note that as far as the IC DIS scattering cross-sections are concerned, the exact Lorentz structure of the couplings is not important. However, as we shall see later, the two-body branching ratio of the HDM to LDM pair is sensitive to the Lorentz structure.

Collider constraints

The collider constraints are sensitive to the interplay of several couplings and mass parameters relevant to our study, specifically, and . A scalar or pseudo-scalar mediator particle which dominantly couples to heavy fermions can be produced in association with one or two b-quarks (involving the parton level processes and respectively). Such a final state may be accessible to LHC searches if the (pseudo-)scalar decays further to an LDM pair . However, in case, , the (pseudo-)scalar would decay back to the SM fermion pairs, thereby making the search considerably harder due to large SM backgrounds. On the other hand, off-shell production does lead to a cross-section in the one or two b-jet(s) and missing transverse momentum (MET) channel. Furthermore, an effective coupling of to gluon pairs is also generated by the top quark loop, and therefore, mono-jet and missing energy searches are also relevant. These bounds have been computed in, for example, Ref. Buckley:2014fba (). The current bounds from these searches are weaker than , across the range of and of our interest Buckley:2014fba (). As we shall see later, the coupling values required in our study are well within the current collider limits. For individual couplings, values of should be allowed, although the LHC bounds are very sensitive to the ratio , which determines the rate of events with MET.

In the case of a spin-1 mediator with either vector or axial vector couplings to SM quarks, the strongest collider constraints come from dijet resonance searches, where the mediator is produced on-shell, and decays back to the SM quarks. Depending upon the values of and , monojet and MET searches could also be important, especially if a) the mediator width is large, making the resonance searches harder, or if b) for a given value of the product , such that the branching ratio to LDM pairs dominates the on-shell mediator decay (when ). Bounds on couplings in the axial-vector case have been discussed in Ref. Chala:2015ama (), which combines the results of different experiments spanning a range of centre of mass energies, including (8 TeV) LHC (ATLAS and CMS), Tevatron and UA2. Similar considerations and bounds would apply to the vector mediator case. For values of , bounds from dijet searches cover masses in the range of GeV to TeV, depending upon the ratio , across the range of values. For a detailed discussion of these bounds for different values of and , we refer the reader to Ref. Chala:2015ama (). With the recent TeV LHC data, ATLAS limits on the coupling to quarks vary in the range of to , as is varied in the range to TeV, when the mediator decay to LDM pairs is absent ATLAS:2016lvi (). Thus, we conclude that the collider bounds on the spin-1 boson couplings are in the range of , and the values required to fit the IC event rates are very much allowed by collider constraints.

Contributions to Galactic and Extra-Galactic Gamma-Ray Fluxes from HDM Decay

The three-body decay of the HDM to a pair of LDMs and a mediator particle (where the mediator particle is radiated by an LDM in the final state), will necessarily contribute to a diffuse gamma ray flux spanning a wide range of energies. This sub-section describes the general method we use to calculate these contributions. The mediator particles lead to hadronic final states via their decays to quark pairs or to hadronically decaying tau pairs, with gamma rays originating from the decays of neutral pions produced in the cascade. Leptonic decays of the mediator can also give rise to high-energy photons via bremsstrahlung and inverse Compton scattering. In the computation of the gamma ray constraints, we only consider the hadronic decay modes of the mediator via quark final states, since the coupling of the mediator to quarks is essential to explaining the IC events in our scenario. For the case of a (pseudo)scalar mediator, the leptonic couplings are expected to be small due to the smaller Yukawa couplings of the charged leptons, while for the case of (axial-)vector mediators, as discussed in Sec. 2.2, consistency with dilepton resonance search constraints favour a setup in which the leptonic couplings are absent. We note in passing that the same three body decays would also lead to signatures in cosmic rays, and there can be additional constraints from measurements of positron and anti-proton fluxes. Due to the large uncertainties in diffusion and propagation models of cosmic rays, we do not include these constraints in our analysis.

The gamma ray flux, like the secondary neutrino flux which we calculate below in Section 3, has a galactic and an extra-galactic component Cirelli:2010xx ():


The extra-galactic flux is isotropic and diffuse (after subtracting out contributions from known astrophysical sources), while the minimum of the galactic flux is an irreducible isotropic contribution to the diffuse flux Cirelli:2010xx (). Since the most important constraints on very high-energy gamma-rays come from air-shower experiments, observations of which are confined to the direction opposite to the Galactic center, we take this minimum to be the flux from the anti-Galactic center, following Refs. Cirelli:2010xx (); Cirelli:2012ut ().

Unlike the neutrino flux, the extra-galactic gamma-ray component suffers significant attenuation due to pair creation processes, and consequently in the energy region of interest here, one finds the galactic component to be the dominant one from any given direction in the sky. This is given by


where, is the total decay width of the HDM, is its mass, and the line of sight integral over the DM halo density is performed along the direction of the anti-GC. We take the DM density profile in our galaxy to be described by a Navarro-Frenk-White distribution astro-ph/9508025 ():


with the standard parameter choices, and kpc. Here, represents the gamma-ray spectra per decay of the HDM in the HDM rest frame. We take the prompt gamma ray energy distribution in the rest frame of the mediator from PPPC4 1012.4515 (), and then subsequently fold it with the three-body differential energy distribution of the mediator obtained using CalcHEP 1207.6082 (), and finally boost the resulting gamma ray spectra to the rest frame of the decaying HDM.

The extra-galactic component of the flux is given by Cirelli:2010xx ()


where, the Hubble constant is given by , with being the present Hubble expansion rate, and and are the matter, DM and dark energy densities respectively, in terms of the present critical density, . We take the values of all relevant cosmological parameters from recent Planck best fits 1303.5076 (); 1502.01589 (). The attenuation factor describes the absorption of gamma rays described above, as a function of the redshift and observed gamma-ray energy , which we take from PPPC4 tables 1012.4515 ().

Having established the framework and general considerations for our study, and outlined the constraints to which it is subject, in the sections to follow we proceed with the specific calculations necessary to demonstrate how IC data may be understood in scenarios combining boosted dark matter and astrophysical neutrinos.

3 Scenario I: PeV events caused by LDM scattering on Ice and its implications

In this section we consider a scenario where boosted DM scattering off ice-nuclei leads to the three events at energies above a PeV seen in the 1347-day HESE sample. In the present data-set, these events are somewhat separated from the others, since there appear to be no HESE events in the region TeV PeV, providing some justification for considering them as disparate from the rest.

Both a) the details of the scattering cross-section of the LDM with ice-nuclei, and b) the three-body spectrum leading to the secondary neutrino flux in sub-PeV energies depend on the particle mediating the interaction. Thus we first examine different mediator candidates — pseudo-scalar, scalar, vector and axial vector — and determine how the corresponding fits and parameters change when a specific choice is made.

As discussed in Sec. 2.3, for the (dominant) two-body decay of the HDM () into a pair of LDM (), the corresponding event rate for scattering is proportional to . The observed rate of the PeV events in IC, along with their deposited energies, then determines a) the ratio of couplings and lifetime , and b) the mass () of the HDM (), using the usual two-body decay kinematics 1407.3280 (). Specifically, if the mean inelasticity of the interaction of the LDM with the ice-nuclei, mediated by a particle is given by , then we require the LDM flux from HDM decay to peak around energies , where represents an estimated average deposited energy at IC for such events.

In this scenario, events in the sub-PeV energy range are then explained by a combination of events from Flux-1 (an astrophysical power-law neutrino flux), Flux-3 (the secondary flux of neutrinos from three-body HDM decay) and Flux-4 (the standard atmospheric neutrino and muon flux), as outlined in Sec. 1.2. For Flux-4, we use the best-fit background estimates from the IC analysis. We determine the best-fit combination of Flux-1 and Flux-3, which, when folded in with the IC-determined best-fit Flux-4 will explain all the sub-PeV observed events in the 1347-days HESE sample. The parameters relevant to this sub-PeV best-fit are , , (the number of sub-PeV events from Flux-1), and (the power-law index for Flux-1).

The total number of shower events within each IC energy bin is given by 1404.0017 ():


Here is the inelasticity parameter, defined in the laboratory frame by , with being the energy deposited in the detector and denotes the energy of the incident dark matter, the runtime of the detector (1347 days) and is the Avogadro number. The limits of the integration are given by and . and are the minimum and maximum deposited energies for an IC energy-bin. is the energy dependent effective detector mass for neutral current interactions obtained from 1311.5238 (). is the differential scattering cross-section, which we quantify below.

The total flux is composed of two parts, the Galactic component and the red-shift () dependant extra-Galactic component . They are given byEsmaili:2012us (); 1308.1105 () :




For the two-body decay , the flux at source is given by :


where, denotes the incident energy at IC for each particle.

We next describe the computation of the secondary neutrino flux due to the three-body decay mode, where one of the daughters () is the mediating particle in scattering. The general procedure is the same as outlined in 1503.02669 (). In our representative calculation here, is assumed to decay to a pair, which by further hadronisation and decays leads to the secondary neutrino spectrum. It is straightforward to obtain the resulting neutrino flux in the rest frame of (see, e.g., 1012.4515 ()), using event generators that implement the necessary showering and hadronisation algorithms, such as PYTHIA8 0710.3820 (). This flux is then boosted to the lab-frame, which is, approximately, the rest frame.

This boosted flux in the rest frame is used in conjunction with Eq. 3 to get the final flux of the secondary neutrinos. The neutrino event rates from this source are determined by folding this flux with the effective area and the exposure time of the detector 1311.5238 ().

Having obtained the event rates for the secondary neutrinos, one defines the necessary to quantify our goodness of fit to the observed data:


Minimizing this  determines the best-fit point in the parameter space of , . It should be noted that the sub-PeV events in Scenario I are due both to the decay of the mediator and a uniform power-law spectrum typical of diffuse astrophysical sources, which is why the overall function is dependent on all the four parameters shown above.

We now turn to discussing the results for specific mediators.

3.1 Pseudoscalar mediator

When the mediator is a pseudo-scalar particle, the corresponding double differential cross-section is given by :


where is the Bjorken scaling parameter, are the masses of the nucleon, LDM, and the mediator respectively, and . is the parton distribution function (PDF) of the quark in the nucleon. We henceforth use the CT10 PDFs Lai:2010vv () throughout our work.

Eq. (12) allows us to compute the event rates (using Eq. (8)) and the mean inelasticity of the scattering process. In Fig. 2 we show the total deep inelastic cross section and the average inelasticity (), and compare them with the case Gandhi:1995tf (); Gandhi:1998ri (); 1106.3723 ().

Figure 2: Representative plots showing the relative behaviour of and neutral current cross sections (left). Average inelasticities are also plotted for both cases (right).

Fig. 3 shows the individual flux components that contribute to the PeV and the sub-PeV events in Scenario I. This is a representative plot, and the parameters that were used while calculating the fluxes are the best-fit values shown in Table. 1.

Figure 3: Relevant fluxes that contribute towards the PeV and the sub-PeV events in Scenario I. The galactic flux is not shown since it originates from the two body decay of , and is given by the simple form in Eq. 10, unlike the extra-galactic flux, which exhibits a dependance. The values of parameters used to calculate the fluxes are given in Table. 1.
Parameter [GeV] (all flavour)
12.0 0.32 2.57
5.3 0.50 2.61
Table 1: The best fit values of relevant parameters in case of a pseudoscalar mediator , when it dominantly decays to and respectively.  is given in units of.

As discussed previously, in Scenario I, the sub-PeV events depend on the mediator mass , the ratio and on the HDM mass . The three PeV events, on the other hand depend on , the ratio and as well as on . Treating the PeV events as arising from two-body decay of the to using gives us PeV. A major fraction of the sub-PeV events arise from the secondary neutrino flux, and for this we carry out calculations in two different kinematic regions: a) where the mediator mass lies above the production threshold, and b) where it lies below this threshold, making the main decay mode. The results for best fits to the data using events from all of the above fluxes, and considering both kinematic regions, are shown in Fig. 4. The solid red line represents the total of the contributions from the various fluxes, and we find that it provides a good description to the data across the energy range of the sample. The best fit values of the parameters are given in Table 1. The corresponding normalisation of the astrophysical flux is shown in terms of the flux at the 100 TeV bin .

We note the following features of Fig 4, which also conform to emergent features of IC data:

  • The secondary neutrino event spectrum has a shape that would allow it to account for a ‘bump’, or excess, such as presently seen in the vicinity of 30–100 TeV.

  • The astrophysical neutrino contribution, especially in the case, is not a major component. This is unlike the standard situation where only astrophysical neutrinos account for events beyond TeV, requiring a flux very close to the Waxman-Bahcall bound.

  • A dip in the region 400–1000 TeV occurs naturally due to the presence of fluxes of different origin in this region.

  • Over the present exposure period, no HESE events are expected in the region beyond 2–3 PeV, since the only contributing flux here is the astrophysical flux, which is significantly lower in this scenario as opposed to the IC best-fits. With more exposure, some astrophysical events can be expected to show up in this region.

Figure 4: Best-fit events (stacked bars) from a combination of secondary ’s, astrophysical ’s and background in the sub-PeV energies, with LDM events explaining the PeV+ events. The best-fit value of PeV. Left: Decays to . Right: The mediator mass limited to below production threshold, so that it can dominantly decay only to pairs.

Parameter correlation analyses

It is useful to examine the parameter space for Scenario I allowed by IC data. We use the case of a pseudo-scalar mediator as representative, and examine the correlations and degeneracies between the parameters. We give contour plots between pairs of parameters for each of the LDM decay scenarios considered above, i.e. for decay to and to . Noting that the sub-PeV events in the HESE sample that do not have their origin in the atmosphere are, in our scenario, either from the secondary neutrino flux or from the astrophysical (power-law) neutrino flux, we denote the total number (in the 1347-day sample) of the former by , and that of the latter by .

For each case we start with the best-fit values obtained in the previous section for each of the parameters in the set: . We note that is proportional to , whereas the primary DM component of the event spectrum, coming from scattering off ice nuclei at PeV energies is related to , and . For a fixed , specifying the is tantamount to specifying the overall astrophysical flux normalisation in the uniform power-law spectrum .

The total number of signal events observed in the 1347-day IC sample is 35 at its best-fit value, with a variation of 29–42 (20–57). This assumes the conventional atmospheric background is at the expected best-fit, and the prompt background is zero. Selecting two parameters for each analysis, we vary their values progressively from their best-fits, while marginalizing over the other parameters over their allowed () ranges. For each pair of the chosen two-parameter subset, we compute the where represent the value of the two chosen parameters in the iteration. With the resulting we plot and contours enclosing the allowed variation of these parameters (Fig. 5, Fig. 6 and Fig. 7). Due to the sparse statistics presently available, the 3 allowed regions in these plots permit the IC data to be fit well for a wide range of values of the chosen variables.

Figure 5: and allowed regions for parameters and (left) and and (right) for mediator decays to . The solid dot in each case represents the corresponding best-fit point in the parameter subspace.
Figure 6: and allowed regions for parameters and (left) and and (right) for mediator decays to . The solid dot in each case represents the corresponding best-fit point in the parameter subspace.
Figure 7: Plot showing allowed regions satisfying gamma ray constraints in the case when pseudoscalar mediator decays to (Left) and to (Right). Regions above the red line are constrained by observations of the diffuse gamma ray flux.
Figure 8: Diffuse gamma-ray flux for the best-fit parameter choice in the pseudo-scalar mediator scenario, where the mediator dominantly decays to (left) and (right). The current constraints from Fermi-LAT data Ackermann:2014usa () at lower energies, and cosmic ray air shower experiment (KASCADE Feng:2015dye () and GRAPES-3 Gupta:2009zz ()) data at higher energies are also shown.

Following the discussion in Sec. 2.3, the only major constraint on the parameters in the pseudo-scalar mediator scenario stems from the upper bound on diffuse gamma-ray fluxes, while the current collider constraints restrict the values of the couplings to values. The sum of the galactic and extragalactic gamma ray fluxes corresponding to the best fit parameter points are shown in Fig. 8. They are compared with both the Fermi-LAT data Ackermann:2014usa () at lower energies, and cosmic ray air shower experiment (KASCADE Feng:2015dye () and GRAPES-3 Gupta:2009zz ()) data at higher energies. These constraints significantly restrict the available parameter-space, and, indeed, our best-fit values for the  lie in a disfavoured region. We find, however, that a reasonable region of the allowed parameters-space is nonetheless consistent with these constraints, and that the allowed region for is larger than that for . Fig. 7 reflects these conclusions.

3.2 Scalar mediator

In this section we explore the case when the mediator in Scenario I is a scalar. The relevant double differential scattering cross-section in this case is given by:


where, the various quantities used are as before (Eq. (12)).

The parameter values at the best-fit point are shown in Table 2, and we show the corresponding event rates in Fig 9. It is interesting to note that, compared to the pseudo-scalar case, due to the additional terms contributing to the differential scattering cross-section (in particular, the term), the best fit value for turns out be smaller in the scalar case, while rest of the relevant parameters take similar values.

Best fit parameters [GeV] (all flavour)
5.3 0.29 2.63
Table 2: The best fit values of relevant parameters in the case of a scalar mediator , when it decays dominantly to . The best fit value of here is 5.3 PeV.  is given in terms of.
Figure 9: Same as Fig. 4, for the scalar mediator scenario, with the mediator dominantly decaying to .

The gamma-ray constraints on the scalar mediator case are found to be similar to the pseudo-scalar case, and as discussed in Sec. 2.3, the collider constraints on the coupling parameters are also of similar magnitude. As further explained in Sec. 2.3, although we restrict ourselves to regions of parameter space where is very small, for parameter values where becomes appreciable, there are additional constraints from relic density requirements as well as direct detection bounds. The spin-independent direct detection bounds in particular are very stringent in the scalar mediator scenario, unless the DM mass lies below , where the nuclear-recoil experiments lose sensitivity. Overall, stronger constraints notwithstanding, we find that the best-fit point lies in an allowed region of the parameter space, and provides an excellent fit to the data, with explanations for the observed features identical to those described in the last subsection.

3.3 Vector and axial-vector mediators

The double differential cross section in the case of a vector mediator is given by:


where, is the coupling of to the quark , and .

To evade the strong bounds particular to vector (and axial-vector) mediators coming from dijet resonance searches in collider experiments, as discussed in Sec. 2.3.1, we impose a penalty on the computation whenever the combination of the coupling constant and extends into a region disfavoured at more than 90% confidence level. Once we have thus determined the allowed region of the parameter space, we show the results (Fig. 10) corresponding to a benchmark point in this space, defined by the values in Table 3, that maximises the contribution from secondary neutrinos from DM decay (Flux-3), and correspondingly deems the astrophysical neutrino component insignificantly small (which we consequently do not show). An increased flux for the latter can be accommodated by a corresponding scaling down of the value of and so on.

Benchmark Values [GeV]
Table 3: Benchmark values of relevant parameters in the case of a vector mediator , when it decays to all possible pairs. The value of used here is 5.0 PeV. As noted in the text, we have chosen a benchmark point in the parameter space that maximises the secondary contribution from DM decay, and consequently deems the astrophysical flux negligible. The latter has therefore not been shown here.

As seen in Fig. 10, unlike the pseudo-scalar and the scalar cases, we note that the galactic and the extra galactic secondary flux events remain approximately flat with decreasing energy below PeV. This results in the absence of a dip or deficit in the region 400 TeV–1 PeV which is one of the features of the present IC data that we would like to reproduce in Scenario I. This can be mitigated by increasing the mass of the mediator (see Fig 11). A comparison with the pseudoscalar mediator event spectrum, where this problem is absent, is shown for a fixed mass, in the right panel Fig. 11.

We now turn to the relevant gamma-ray constraints, along the same lines we studied it for the case of a pseudo-scalar mediator. While the differential three-body decay width of the HDM follows somewhat different distributions for different choices of mediator spin and CP properties, the very large boost of the mediator particle washes out these differences to a large extent, and we arrive at a similar spectral shape as discussed for the spin-0 mediators above. We find that the corresponding constraints are not severe, but may have mild tension in some energy regions. As far as relic density and spin-independent direct detection bounds are concerned, similar considerations as in the scalar mediator case would also apply to the vector mediator scenario, and we refer the reader to the discussion in Sec. 3.2.

Figure 10: Event rates for the benchmark parameter values shown in Table 3. In keeping with the description in text, the correspondingly tiny number of events from the astrophysical flux have not been shown here.
Figure 11: Left: PeV events in the vector mediator scenario, with different choices for the mass. A larger value of the mass is more likely to explain the dip at around PeV. Right: PeV events in vector and pseudoscalar case with a mediator mass fixed to 20 GeV. The pseudoscalar scenario, as discussed earlier, explains the dip more accurately because of it’s sharply falling event rates, unlike in the vector scenario.

Even though the differential cross-section behaves similarly in the vector and axial-vector scenarios (in small and limit), there are additional important considerations particular to the axial-vector case that limit the available parameter space very stringently. As explained earlier, in order to accommodate the PeV events by DIS scattering, we require that the three body decay width of the HDM is much smaller than its two body decay width. However, as shown in Fig. 12, the three-body branching ratio starts to dominate for values as low as in the axial-vector case, whereas for scalar, pseudo-scalar or vector mediators, the three-body branching ratio becomes large only for . Thus, since the PeV event rate is proportional to , to obtain the required number of events in the PeV region, the value of needs to be pushed higher than its perturbative upper bound of . Ultimately, we find that it is not possible to fit both the PeV and the sub-PeV events while simultaneously satisfying the perturbativity requirement for an axial-vector mediator.

Figure 12: Variation of three body branching ratio with for the vector, axial-vector and the pseudoscalar mediators. The scalar mediator scenario shows a similar behaviour as the pseudo-scalar one.

4 Scenario II: Excess events in the TeV region caused by LDM scattering on Ice and its implications

As discussed in Sec. 1.2, in Scenario II, we relax the assumption made regarding the origin of the three PeV events in Scenario I, and perform a completely general fit to both the PeV and the sub-PeV HESE data, with all four of the flux components taken together. This essentially implies that the HDM mass is also left floating in the fit in the entire range . Therefore, the space of parameters now comprises the set and .

We find that doing this causes the best-fit HDM mass to float to a value TeV, so that the resulting LDM spectra from its decay are naturally able to explain the bump, or excess in the 50–100 TeV energy range that is seen in the IC data. At the present time, this feature has a statistical significance of about . An important consequence of this is that the flux of secondary neutrinos from mediator decay, which played an important role in Scenario I, now populates the low energy bins (between TeV to TeV) and falls outside the range relevant to our fit (the IC threshold for the HESE events is 30 TeV ). This flux is thus subsumed in the atmospheric background. At energies of around a TeV, where the secondary neutrino flux from three-body decays of HDMs in this scenario might have been otherwise important, the atmospheric neutrino flux is already about a 1000 times higher, and completely overwhelms it. Furthermore, the full-volume IceCube is only sensitive to contained events depositing at least about 10 TeV in the detector, hence this flux is also largely rendered unobservable because it lies outside the HESE sensitivity range.

Note that Scenario II also suggests that the other currently emergent features, the cluster of 3 events close to 1–2 PeV and the dip in the 400 TeV–1 PeV region, which were very important motivations for Scenario I, may not survive with time. Thus, at the current level of statistics, this fit gives primacy to the 50–100 TeV excess. In Scenario I, the PeV events, assumed to arise from the two-body decay of HDM, will (in the form of cascades resembling NC neutrino events) steadily increase in number and manifest themselves as an excess or bump, whereas in Scenario II they would just become part of the overall astrophysical power-law neutrino spectrum without a special origin. The related dip, or deficit, currently seen in the TeV to PeV region would gradually become prominent and significant in Scenario I, but would get smoothed over in Scenario II. Consequently,, in Scenario II the only relevant fluxes are the astrophysical flux and the flux originating from the two body decay of , in addition to, of course, the background atmospheric flux. We show the representative contributing fluxes in Fig. 13.

The best fit parameters for the fit in Scenario II are given in Table 4, and the corresponding results are shown in Fig. 14, for the pseudo-scalar mediator scenario (left column), and the axial-vector mediator scenario (right column). As in Scenario I, the scalar and pseudo-scalar mediators lead to similar fits. However, unlike in Scenario I, since the secondary neutrino flux lies outside the energy range under study, both vector and axial-vector mediators lead to similar results for Scenario II. Therefore, we have not shown the scalar and vector cases separately.

Parameter [GeV] [TeV] (all flavour)
Pseudoscalar 16.1 680 2.31
Axial-vector 470 2.30
Table 4: The best fit values of relevant parameters in case of a pseudoscalar and axial-vector mediator for Scenario II. is given in units of .

The similarity in the number of events originating from DM and from astrophysical neutrinos in the two cases is not surprising. In both cases, only the small excess in the vicinity of TeV is due to DM cascades, the remaining events conform to the expected astrophysical neutrino spectrum, which then sets the normalization and the index. Consequently, we also note an important difference between the astrophysical fluxes in Scenario II compared to Scenario I, i.e. in Scenario I this flux is usually sub-dominant to the secondary neutrino flux, whereas in Scenario II it accounts for all events except those comprising the excess in the range TeV. The difference in in the two cases is due to the variation in the values of for the two type of mediators.

Figure 13: Relevant fluxes for Scenario II. The corresponding parameters are given in Table 4. As before the monochromatic spike at due to the galactic flux is not shown here.
Figure 14: The total event rate is shown as the red solid curve. This comprises events from LDM scattering, astrophysical neutrinos and the atmospheric background. Events from the astrophysical power-law spectrum are shown as orange bars and stacked bars shaded in green show the LDM events over and above the astrophysical events. The other events over and above the green/yellow bars are due to atmospheric neutrinos and muons. The left hand side shows the pseudo-scalar case while the right hand side gives the case of an axial-vector type mediator.

4.1 Gamma-ray constraints on Scenario II

As for Scenario I, the diffuse gamma-ray constraints provide the most significant restrictions on our parameter space, and lead to upper bounds on . The behaviour of the differential -ray flux is sensitive to the mediator mass and the type of mediator under study, as shown in Fig. 15.

Figure 15: Diffuse gamma-ray flux for pseudo-scalar (left) and axial-vector case (right). The maximum allowed values of have been used for the flux computation here.

Using results on the diffuse gamma ray fluxes from Fermi-LAT, KASCADE and GRAPES3 data, we obtain upper bounds on for the pseudo-scalar and axial-vector cases, respectively, as follows :


The upper bounds on that result from the above are significantly more stringent for the axial-vector case, and rule out the best-fit case shown in Fig. 14 for this mediator. The best-fit shown for the pseudo-scalar case is broadly consistent with the current gamma-ray constraints.

5 Muon-Track events

Our discussion so far has been confined to the HESE events, whose starting vertices are, by definition, contained within the IC instrumented volume. More recently, however, a 6-year analysis of through-going muon track events at IC has been reported 1607.08006 (). The events in this data sample include those with interaction vertices outside this volume. There are events both in the PeV and the sub-PeV regions. When fit with a uniform astrophysical power-law flux, this sample prefers a stronger astrophysical spectrum, with . This is notably different from the conclusion from the HESE analysis, which suggests , whilst disfavouring a spectrum with at more than . This tension could, perhaps be a hint for additional flux components which cannot be accounted for in a simple power-law picture. Indeed, as pointed out in 1607.08006 (), a possible reason for the tension could be a flux component from galactic sources, which becomes sub-dominant as the energy increases. We note that the secondary neutrino flux from the galaxy, which dominates the sub-PeV contribution in Scenario I, is a possibility that conforms to this requirement.

While we have not attempted a full comparative study of this sample in the context of our scenarios here, we have tried to get an approximate idea of the track event predictions that Scenario I and II would give. In Scenario I, for example, contributions to these events would arise from the secondary neutrino and astrophysical fluxes. We can then compare the predicted event rates with those predicted by the IC best-fit astrophysical flux (with index 2.13, from 1607.08006 ()). We show the comparisons in Fig. 16 for the pseudoscalar mediator in Scenario I. The through-going track events span the energy range from 190 TeV to a few PeV 1607.08006 () . For both the cases when pseudoscalar and we have taken a value of which satisfies all constraints. For the astrophysical flux, the values of the index and the normalisations were however fixed to their best-fit values (Fig. 16).

We find good overlap with the IC prediction (i.e., the red and black curves) in the lower part of the energy range of interest, i.e. 190 TeV to TeV (where most of the observations lie); however, for higher energies the curves differ, and Scenario I predicts substantially less through-going muon track events. We note that statistics in higher energy region are sparse, making definitive conclusions difficult. In the multi-PeV region, for instance, the highest energy event in this 6-yr sample 1607.08006 (), has a deposited energy of PeV, and an estimated muon energy of about PeV. It is difficult to say if this is an unusually high energy event isolated in origin from the rest; for a detailed discussion of possibilities, see 1605.08781 ().

Similarly, we show the IC prediction along with the expectation for Scenario II in Fig 17. Although our Scenario II flux is somewhat lower than the IC fit, the agreement overall is reasonable (given the present level of statistics), since the astrophysical power-law flux is a dominant contributor in Scenario II, unlike in Scenario I. Further confirmation will have to await more data, especially in the high energy region ( PeV).

Figure 16: Muon track events for the pseudoscalar case in Scenario I and their comparison with the IC predicted best fit. The black line represents the IC power-law prediction and should be compared to our total prediction for throughgoing track events in the energy region 190 TeV to a few PeV (red line).
Figure 17: Muon track events in Scenario II. Shown for the case of pseudoscalar (left) and axial-vector type mediators (right). In Scenario II, the astrophysical flux is the main contributor to the track events. In our notation . Best fit values of and are used in the above plot.

6 Summary and Conclusions

By steadily accumulating high energy events over the last four years in the energy range TeV to PeV, IC has conclusively established the presence of a diffuse flux or fluxes which have a non-atmospheric origin and (at least partially) extra-galactic origin, the source(s) of which are at present largely unknown.

Standard expectations dictate that this signal is due to a flux of astrophysical neutrinos, primarily from sources outside of our galaxy, and that it should correspond to a uniform power-law flux, characteristic of Fermi shock acceleration, with index approximately . Features in the data seem to indicate that there are deviations from these expectations, which may signal the presence of one or more additional fluxes. These features include a) a lack of cascade events beyond 2.1 PeV, in spite of both IC’s sensitivity in this region, and the presence of the Glashow resonance around 6.3 PeV; b) a possible dip in the spectrum between 400 TeV–1 PeV; c) a low energy excess of around significance over and above the IC best-fit power-law spectrum in the energy range TeV. In addition, an overall puzzling feature of the flux is its unexpected proximity to the WB bound, since standard expectations would argue for a neutrino flux that is a factor of a few below this upper limit.

In this work, we have explored the idea that some of the events in IC which cause the overall signal to deviate from the standard power-law originate from the scattering of boosted DM on ice. We have considered two scenarios, both involving the incidence of such fermionic dark matter (LDM), which is produced (in the context of a minimal two-component dark matter sector) from the slow decay of its (significantly) heavier cousin (HDM). The LDM, upon scattering off the ice-nuclei inside IC, mimicks standard model neutrino-nucleon neutral current scattering, but, in general, with weaker interaction strengths. If the HDM has a mass 5–10 PeV, the LDM flux can be shown to peak in a cluster around the 1–2 PeV energies and, with the right parameters, can explain the IC PeV events. This forms the basis of Scenario I, which accounts for the rest of the events (at sub-PeV energies), by a combination of those from astrophysical sources and a secondary neutrino flux originating from the decay of the mediator involved in the LDM-nucleon scattering. It is interesting to note that the secondary neutrinos naturally provide a bump in the region TeV once the parameters for the three PeV events from LDM scattering are fixed.

On the other hand, in Scenario II, for lighter masses of the HDM 500–800 TeV, the LDM flux leads to scattering events in the sub-PeV 30–100 TeV energies and is helpful in explaining this low-energy excess over and above a harder (compared to Scenario I) astrophysical power-law flux. In both scenarios, in order to explain observations, our work incorporates the detection of boosted DM by IC, in addition to its detecting UHE neutrinos. This allows the standard astrophysical flux to stay appreciably below the WB bound for Scenario I, and, to a lesser extent, for Scenario II.

Four different mediators which connect the SM and DM sectors are considered, specifically, scalar, pseudo-scalar, vector and axial-vector. For Scenario I, we find excellent fits to the IC data in both spin-0 mediator cases — the LDM scattering explains the three PeV events with a hard cut-off set by the HDM mass. It has a soft astrophysical power-law flux that dies out around energies of 400 TeV, and a small but significant neutrino flux from the decay of the mediator that helps explain the small bump around 30-100 TeV, making the full spectrum a better match to the data than a power-law-only spectrum. However, for the pseudo-scalar, stringent constraints from -ray observations rule out the region of parameter space where the best-fit itself lies. The allowed parameter-space region around the best-fit is quite large, nevertheless, and we find that a significant portion of this is as yet allowed by the -ray bounds.

For spin-1 mediators, in Scenario I we find significantly increased tension between constraints and best-fit parameters in the vector mediator case, but are, nevertheless, able to fit the IC data well for specific values of the parameters within the allowed regions. The case for the axial-vector mediator is, unfortunately, more pessimistic: we find that perturbativity requirements on the coupling constants prevent a simultaneous fit to the observed PeV and sub-PeV data.

If, with future data, Scenario I were to sustain, we would expect to see a gradual statistical improvement in the evidence for a dip-like structural feature around 400-800 TeV, since this region marks the interface of fluxes of different origins. One would see a paucity of events beyond 2.1 PeV, due to a significantly lower astrophysical flux compared to current IC predictions. In addition, a PeV event spectrum predominantly from LDM scattering (due to HDM decay) predicts i) a significantly enhanced ratio of cascade-to-track events approximately in the (0.75-2.5 PeV) region, ii) a build-up in the number of such cascade events in this region as the HDM decay and LDM scattering proceed, and iii) a small but non-zero number of up-going cascades in this energy region over time from the northern hemisphere compared to the case where these events would have been due to a neutrino flux (because of the relatively lower -nucleon cross section and consequent reduced screening by the earth.)11. Finally, through-going muon track events beyond PeV are also expected to be lower in number in this scenario than what current IC power-law fit predictions suggest. The overall signal would also exhibit a gradual galactic bias with more statistics, since generically, in DM scenarios, the contributions from our galaxy and from extra-galactic DM are roughly of the same order. Such a directional bias is not expected in a genuinely isotropic flux12.These features would be in contrast to what one would expect to see if the standard astrophysical power-law flux explanation were indeed responsible for the observed events and will be discernable as statistics increase.

Scenario II, on the other hand, is designed to explain only the event excess at 50–100 TeV energies as being due to DM scattering on ice, with the other events, including those above 1 PeV, attributed to an astrophysical neutrino power-law spectrum. It thus assumes that the other features, including the 400-800 TeV dip and the existence of a cut-off beyond 2 PeV, which are part of Scenario I, would gradually disappear and smooth out over time. It is in good agreement with the HESE data, and because it requires a harder astrophysical spectrum to explain the highest energy events, it is also in better agreement with IC’s six-year through-going muon track data. Indeed, its predictions for both cascade and track events (both starting and through-going) are only slightly below those of the official IC fits. The secondary neutrinos produced from the decays of the mediator in this scenario peak at energies around a TeV and lie in a region dominated by the conventional atmospheric background. They are thus not consequential to our considerations here. With respect to the different types of mediators, this scenario is somewhat less constrained over-all compared to Scenario I. The best-fits we obtain for the vector and axial-vector cases are disallowed by gamma-ray observations; nevertheless, good fits in the region are possible. The scalar and pseudo-scalar mediators make for better agreement, with their best-fits being allowed.

To conclude, we have shown that present differences in the IC data in comparison to what is expected from standard astrophysical diffuse neutrino fluxes may be explained by assuming that the full spectrum is made up of multiple flux components, with one significant component being the flux of a boosted DM particle. Depending on the HDM mass, the LDM flux either peaks at PeV energies (Scenario I) and explains the PeV events in the 4-yr HESE sample, or at lower energies (Scenario II) and aids in explaining the 50–100 TeV excess. In Scenario I, the excess at 50-100 TeV is naturally accounted for by a secondary neutrino flux from HDM decay. In both cases, the different components conspire in ways that explain the IC data better than any single component flux can. This is in spite of strong constraints from -ray observations, which limit but do not completely exclude the available 3 parameter space around the corresponding best-fits. On this note, it is worth mentioning that our work skirts the recent strong constraints 1612.05638 () on masses and lifetimes of heavy DM decays as explanations of IC events, as they apply to scenarios in which such DM decays directly to SM particles. Finally, we have also discussed signatures that would, with future data, help distinguish each case under consideration from fits with a solely uniform power-law flux. More data over the next few years should be able to conclusively support or veto such multi-component explanations of high-energy observations at IC compared to other, more standard expectations.


RG thanks Carlos Arguelles, Ben Jones and Joachim Kopp for helpful discussions, and the Fermilab Theory Group and the Neutrino Division for hospitality while this work was in progress. AG and RG also acknowledge the DAE Neutrino Project Grant for continued support. AB is grateful to Jean-René Cudell, Maxim Laletin and Daniel Wegman Ostrosky for helpful discussions. AB is supported by the Fonds de la Recherche Scientifique-FNRS, Belgium, under grant No. 4.4501.15. SM is supported in part by the U.S. Department of Energy under grant No. DE-FG02-95ER40896 and in part by the PITT PACC.


  1. In addition to the analysis presented by the IceCube collaboration in Aartsen:2015zva (); 1311.5238 (), a recent analysis of the HESE data may be found in 1605.01556 ().
  2. The statistical significance is dependent upon the largely theoretically modelled upper limits of the prompt neutrino flux from heavy meson decays. The value corresponds to the scenario where the prompt flux is assumed to be absent. Nonetheless, even with the highest upper limits from present computations, the statistical significance of a new signal over and above the atmospheric background is well above .
  3. This classification allows us to categorize most events. There are other, potentially important types of events, however, which have not yet been observed; e.g. the double bang events signalling the CC production of a highly energetic lepton hep-ph/9405296 (), and the pure muon and contained lollipop Bhattacharya:2011qu () events which would unambiguously signal the detection of the Glashow resonance PhysRev.118.316 (); Berezinsky:1981bt (); Gandhi:1995tf ().
  4. A recent analysis 1611.07905 () statistically reinforces the presence of a break in the spectrum in the region 200 - 500 TeV, which could have a bearing on this feature.
  5. The WB bound is valid for sources which produce neutrinos as a result of or interactions. It assumes that they are optically thin to proton photo-meson and proton-nucleon interactions, allowing protons to escape. Such sources are characterized by an optical depth which is typically less than one. As explained in Bahcall:1999yr (), the bound is conservative by a factor of .
  6. Alternatively, such features could, of course, also indicate that the conventional neutrino astrophysical flux, while originating in standard physics, is much less understood than we believe, and may have more than one component.
  7. We note that a two-component WIMP-like DM scenario, with the lighter particle (of mass (1 GeV)) being boosted after production (via annihilation in the galactic halo of its heavier partner of mass GeV) and subsequently detected in neutrino experiments has been discussed in 1405.7370 (). Boosted thermal DM detection from the sun and the galactic center due to annihilation of a heavier counterpart at similar masses and energies has been discussed in 1410.2246 (); 1411.6632 (); 1611.09866 ().
  8. The statistical significance of such a break has now increased due to the recent release of six-year muon track data 1607.08006 (); see, for instance, the discussion in  1611.07905 (). Additionally, as we shall see below, by providing a significant fraction of the events directly (via Flux 2) or indirectly (via Flux 3) from DM, this scenario does not require the astrophysical neutrino flux to be pushed up uncomfortably close to the Waxman-Bahcall bound, unlike the standard single power-law interpretation.
  9. We note that stronger constraints, based on IC data and Fermi-LAT, as discussed recently in 1612.05638 () are evaded in our work since they are derived assuming the two-body decay of dark matter directly to SM particles, e.g..
  10. This assumes that the two-body decay to is the dominant mode.
  11. We note that IC has already observed an upgoing cascade in this energy region, with deposited energy PeV Ishihara ()
  12. We stress that in our scenario also, the events due to the astrophysical neutrino flux (Flux-1) will be isotropic in distribution. The directional bias will be exhibited by only those events that originate in DM( Flux-2 and Flux-3).


  1. IceCube collaboration, M. G. Aartsen et al., The IceCube Neutrino Observatory - Contributions to ICRC 2015 Part II: Atmospheric and Astrophysical Diffuse Neutrino Searches of All Flavors, in Proceedings, 34th International Cosmic Ray Conference (ICRC 2015), 2015. 1510.05223.
  2. IceCube collaboration, M. G. Aartsen et al., Evidence for High-Energy Extraterrestrial Neutrinos at the IceCube Detector, Science 342 (2013) 1242856, [1311.5238].
  3. A. C. Vincent, S. Palomares-Ruiz and O. Mena, Analysis of the 4-year IceCube HESE data, 1605.01556.
  4. J. K. Becker, High-energy neutrinos in the context of multimessenger physics, Phys. Rept. 458 (2008) 173–246, [0710.1557].
  5. F. Halzen and D. Hooper, High-energy neutrino astronomy: The Cosmic ray connection, Rept. Prog. Phys. 65 (2002) 1025–1078, [astro-ph/0204527].
  6. G. Sigl, High Energy Neutrinos and Cosmic Rays, Proc. Int. Sch. Phys. Fermi 182 (2012) 145–184, [1202.0466].
  7. M. D. Kistler, T. Stanev and H. Yüksel, Cosmic PeV Neutrinos and the Sources of Ultrahigh Energy Protons, Phys. Rev. D90 (2014) 123006, [1301.1703].
  8. N. Gupta, Galactic PeV Neutrinos, Astropart. Phys. 48 (2013) 75–77, [1305.4123].
  9. R. Laha, J. F. Beacom, B. Dasgupta, S. Horiuchi and K. Murase, Demystifying the PeV Cascades in IceCube: Less (Energy) is More (Events), Phys. Rev. D88 (2013) 043009, [1306.2309].
  10. L. A. Anchordoqui, H. Goldberg, M. H. Lynch, A. V. Olinto, T. C. Paul and T. J. Weiler, Pinning down the cosmic ray source mechanism with new IceCube data, Phys. Rev. D89 (2014) 083003, [1306.5021].
  11. L. A. Anchordoqui et al., Cosmic Neutrino Pevatrons: A Brand New Pathway to Astronomy, Astrophysics, and Particle Physics, JHEAp 1-2 (2014) 1–30, [1312.6587].
  12. K. Murase, On the Origin of High-Energy Cosmic Neutrinos, AIP Conf. Proc. 1666 (2015) 040006, [1410.3680].
  13. W. Winter, Describing the Observed Cosmic Neutrinos by Interactions of Nuclei with Matter, Phys. Rev. D90 (2014) 103003, [1407.7536].
  14. S. Dado and A. Dar, Origin of the High Energy Cosmic Neutrino Background, Phys. Rev. Lett. 113 (2014) 191102, [1405.5487].
  15. I. Cholis and D. Hooper, On The Origin of IceCube’s PeV Neutrinos, JCAP 1306 (2013) 030, [1211.1974].
  16. R.-Y. Liu and X.-Y. Wang, Diffuse PeV neutrinos from gamma-ray bursts, Astrophys. J. 766 (2013) 73, [1212.1260].
  17. S. Razzaque, Long-lived PeV–EeV neutrinos from gamma-ray burst blastwave, Phys. Rev. D88 (2013) 103003, [1307.7596].
  18. K. Murase and K. Ioka, TeV–PeV Neutrinos from Low-Power Gamma-Ray Burst Jets inside Stars, Phys. Rev. Lett. 111 (2013) 121102, [1306.2274].
  19. N. Fraija, GeV–PeV neutrino production and oscillation in hidden jets from gamma-ray bursts, Mon. Not. Roy. Astron. Soc. 437 (2014) 2187–2200, [1310.7061].
  20. M. Petropoulou, D. Giannios and S. Dimitrakoudis, Implications of a PeV neutrino spectral cutoff in GRB models, Mon. Not. Roy. Astron. Soc. 445 (2014) 570–580, [1405.2091].
  21. S. Razzaque and L. Yang, PeV-EeV neutrinos from GRB blast waves in IceCube and future neutrino telescopes, Phys. Rev. D91 (2015) 043003, [1411.7491].
  22. I. Tamborra and S. Ando, Diffuse emission of high-energy neutrinos from gamma-ray burst fireballs, JCAP 1509 (2015) 036, [1504.00107].
  23. N. Senno, K. Murase and P. Meszaros, Choked Jets and Low-Luminosity Gamma-Ray Bursts as Hidden Neutrino Sources, Phys. Rev. D93 (2016) 083003, [1512.08513].
  24. I. Tamborra and S. Ando, Inspecting the supernova–gamma-ray-burst connection with high-energy neutrinos, Phys. Rev. D93 (2016) 053010, [1512.01559].
  25. K. Murase, M. Ahlers and B. C. Lacki, Testing the Hadronuclear Origin of PeV Neutrinos Observed with IceCube, Phys. Rev. D88 (2013) 121301, [1306.3417].
  26. L. A. Anchordoqui, T. C. Paul, L. H. M. da Silva, D. F. Torres and B. J. Vlcek, What IceCube data tell us about neutrino emission from star-forming galaxies (so far), Phys. Rev. D89 (2014) 127304, [1405.7648].
  27. I. Tamborra, S. Ando and K. Murase, Star-forming galaxies as the origin of diffuse high-energy backgrounds: Gamma-ray and neutrino connections, and implications for starburst history, JCAP 1409 (2014) 043, [1404.1189].
  28. X.-C. Chang and X.-Y. Wang, The diffuse gamma-ray flux associated with sub-PeV/PeV neutrinos from starburst galaxies, Astrophys. J. 793 (2014) 131, [1406.1099].
  29. I. Bartos and S. Marka, Spectral Decline of PeV Neutrinos from Starburst Galaxies, 1509.00983.
  30. O. E. Kalashev, A. Kusenko and W. Essey, PeV neutrinos from intergalactic interactions of cosmic rays emitted by active galactic nuclei, Phys. Rev. Lett. 111 (2013) 041103, [1303.0300].
  31. F. W. Stecker, PeV neutrinos observed by IceCube from cores of active galactic nuclei, Phys. Rev. D88 (2013) 047301, [1305.7404].
  32. C. D. Dermer, K. Murase and Y. Inoue, Photopion Production in Black-Hole Jets and Flat-Spectrum Radio Quasars as PeV Neutrino Sources, JHEAp 3-4 (2014) 29–40, [1406.2633].
  33. F. Krauß et al., TANAMI Blazars in the IceCube PeV Neutrino Fields, Astron. Astrophys. 566 (2014) L7, [1406.0645].
  34. F. Tavecchio, G. Ghisellini and D. Guetta, Structured jets in BL Lac objects: efficient PeV neutrino factories?, Astrophys. J. 793 (2014) L18, [1407.0907].
  35. S. Sahu and L. S. Miranda, Some possible sources of IceCube TeV–PeV neutrino events, Eur. Phys. J. C75 (2015) 273, [1408.3664].
  36. O. Kalashev, D. Semikoz and I. Tkachev, Neutrinos in IceCube from active galactic nuclei, J. Exp. Theor. Phys. 120 (2015) 541–548, [1410.8124].
  37. F. Tavecchio and G. Ghisellini, High-energy cosmic neutrinos from spine-sheath BL Lac jets, Mon. Not. Roy. Astron. Soc. 451 (2015) 1502–1510, [1411.2783].
  38. S. S. Kimura, K. Murase and K. Toma, Neutrino and Cosmic-Ray Emission and Cumulative Background from Radiatively Inefficient Accretion Flows in Low-Luminosity Active Galactic Nuclei, Astrophys. J. 806 (2015) 159, [1411.3588].
  39. M. Petropoulou, S. Dimitrakoudis, P. Padovani, A. Mastichiadis and E. Resconi, Photohadronic origin of -ray BL Lac emission: implications for IceCube neutrinos, Mon. Not. Roy. Astron. Soc. 448 (2015) 2412–2429, [1501.07115].
  40. R.-Y. Liu, X.-Y. Wang, S. Inoue, R. Crocker and F. Aharonian, Diffuse PeV neutrinos from EeV cosmic ray sources: Semirelativistic hypernova remnants in star-forming galaxies, Phys. Rev. D89 (2014) 083004, [1310.1263].
  41. S. Chakraborty and I. Izaguirre, Diffuse neutrinos from extragalactic supernova remnants: Dominating the 100 TeV IceCube flux, Phys. Lett. B745 (2015) 35–39, [1501.02615].
  42. A. Bhattacharya, R. Enberg, M. H. Reno and I. Sarcevic, Charm decay in slow-jet supernovae as the origin of the IceCube ultra-high energy neutrino events, JCAP 1506 (2015) 034, [1407.2985].
  43. L. A. Anchordoqui, H. Goldberg, T. C. Paul, L. H. M. da Silva and B. J. Vlcek, Estimating the contribution of Galactic sources to the diffuse neutrino flux, Phys. Rev. D90 (2014) 123010, [1410.0348].
  44. H. Gao, B. Zhang, X.-F. Wu and Z.-G. Dai, Possible High-Energy Neutrino and Photon Signals from Gravitational Wave Bursts due to Double Neutron Star Mergers, Phys. Rev. D88 (2013) 043010, [1306.3006].
  45. X.-Y. Wang and R.-Y. Liu, Tidal disruption jets of supermassive black holes as hidden sources of cosmic rays: explaining the IceCube TeV-PeV neutrinos, Phys. Rev. D93 (2016) 083005, [1512.08596].
  46. A. Neronov, D. V. Semikoz and C. Tchernin, PeV neutrinos from interactions of cosmic rays with the interstellar medium in the Galaxy, Phys. Rev. D89 (2014) 103002, [1307.2158].
  47. J. C. Joshi, W. Winter and N. Gupta, How Many of the Observed Neutrino Events Can Be Described by Cosmic Ray Interactions in the Milky Way?, Mon. Not. Roy. Astron. Soc. 439 (2014) 3414–3419, [1310.5123].
  48. B. Katz, E. Waxman, T. Thompson and A. Loeb, The energy production rate density of cosmic rays in the local universe is at all particle energies, 1311.0287.
  49. K. Fang, T. Fujii, T. Linden and A. V. Olinto, Is the Ultra-High Energy Cosmic-Ray Excess Observed by the Telescope Array Correlated with IceCube Neutrinos?, Astrophys. J. 794 (2014) 126, [1404.6237].
  50. M. Kachelrieß and S. Ostapchenko, Neutrino yield from Galactic cosmic rays, Phys. Rev. D90 (2014) 083002, [1405.3797].
  51. L. A. Anchordoqui, Neutron -decay as the origin of IceCube’s PeV (anti)neutrinos, Phys. Rev. D91 (2015) 027301, [1411.6457].
  52. Y. Q. Guo, H. B. Hu and Z. Tian, On the Contribution of "Fresh" Cosmic Rays to the Excesses of Secondary Particles, 1412.8590.
  53. A. M. Taylor, S. Gabici and F. Aharonian, Galactic halo origin of the neutrinos detected by IceCube, Phys. Rev. D89 (2014) 103003, [1403.3206].
  54. F. Zandanel, I. Tamborra, S. Gabici and S. Ando, High-energy gamma-ray and neutrino backgrounds from clusters of galaxies and radio constraints, Astron. Astrophys. 578 (2015) A32, [1410.8697].
  55. B. Feldstein, A. Kusenko, S. Matsumoto and T. T. Yanagida, Neutrinos at IceCube from Heavy Decaying Dark Matter, Phys. Rev. D88 (2013) 015004, [1303.7320].
  56. A. Esmaili and P. D. Serpico, Are IceCube neutrinos unveiling PeV-scale decaying dark matter?, JCAP 1311 (2013) 054, [1308.1105].
  57. Y. Ema, R. Jinno and T. Moroi, Cosmic-Ray Neutrinos from the Decay of Long-Lived Particle and the Recent IceCube Result, Phys. Lett. B733 (2014) 120–125, [1312.3501].
  58. A. Esmaili, S. K. Kang and P. D. Serpico, IceCube events and decaying dark matter: hints and constraints, JCAP 1412 (2014) 054, [1410.5979].
  59. A. Bhattacharya, M. H. Reno and I. Sarcevic, Reconciling neutrino flux from heavy dark matter decay and recent events at IceCube, JHEP 06 (2014) 110, [1403.1862].
  60. A. Bhattacharya, R. Gandhi and A. Gupta, The Direct Detection of Boosted Dark Matter at High Energies and PeV events at IceCube, JCAP 1503 (2015) 027, [1407.3280].
  61. Y. Ema, R. Jinno and T. Moroi, Cosmological Implications of High-Energy Neutrino Emission from the Decay of Long-Lived Particle, JHEP 10 (2014) 150, [1408.1745].
  62. J. F. Cherry, A. Friedland and I. M. Shoemaker, Neutrino Portal Dark Matter: From Dwarf Galaxies to IceCube, 1411.1071.
  63. J. Kopp, J. Liu and X.-P. Wang, Boosted Dark Matter in IceCube and at the Galactic Center, JHEP 04 (2015) 105, [1503.02669].
  64. K. Murase, R. Laha, S. Ando and M. Ahlers, Testing the Dark Matter Scenario for PeV Neutrinos Observed in IceCube, Phys. Rev. Lett. 115 (2015) 071301, [1503.04663].
  65. A. Esmaili and P. D. Serpico, Gamma-ray bounds from EAS detectors and heavy decaying dark matter constraints, JCAP 1510 (2015) 014, [1505.06486].
  66. L. A. Anchordoqui, V. Barger, H. Goldberg, X. Huang, D. Marfatia, L. H. M. da Silva et al., IceCube neutrinos, decaying dark matter, and the Hubble constant, Phys. Rev. D92 (2015) 061301, [1506.08788].
  67. S. M. Boucenna, M. Chianese, G. Mangano, G. Miele, S. Morisi, O. Pisanti et al., Decaying Leptophilic Dark Matter at IceCube, JCAP 1512 (2015) 055, [1507.01000].
  68. P. Ko and Y. Tang, IceCube Events from Heavy DM decays through the Right-handed Neutrino Portal, Phys. Lett. B751 (2015) 81–88, [1508.02500].
  69. M. Chianese, G. Miele, S. Morisi and E. Vitagliano, Low energy IceCube data and a possible Dark Matter related excess, Phys. Lett. B757 (2016) 251–256, [1601.02934].
  70. P. S. B. Dev, D. Kazanas, R. N. Mohapatra, V. L. Teplitz and Y. Zhang, Heavy right-handed neutrino dark matter and PeV neutrinos at IceCube, JCAP 1608 (2016) 034, [1606.04517].
  71. M. Chianese and A. Merle, A Consistent Theory of Decaying Dark Matter Connecting IceCube to the Sesame Street, 1607.05283.
  72. V. Barger and W.-Y. Keung, Superheavy Particle Origin of IceCube PeV Neutrino Events, Phys. Lett. B727 (2013) 190–193, [1305.6907].
  73. A. N. Akay, U. Kaya and S. Sultansoy, Color octet neutrino as the source of the IceCube PeV energy neutrino events, 1402.1681.
  74. I. Alikhanov, The Glashow resonance in neutrino–photon scattering, Phys. Lett. B741 (2015) 295–300, [1402.6678].
  75. L. A. Anchordoqui, V. Barger, H. Goldberg, J. G. Learned, D. Marfatia, S. Pakvasa et al., End of the cosmic neutrino energy spectrum, Phys. Lett. B739 (2014) 99–101, [1404.0622].
  76. K. Ioka and K. Murase, IceCube PeV–EeV neutrinos and secret interactions of neutrinos, PTEP 2014 (2014) 061E01, [1404.2279].
  77. K. C. Y. Ng and J. F. Beacom, Cosmic neutrino cascades from secret neutrino interactions, Phys. Rev. D90 (2014) 065035, [1404.2288].
  78. J. Zavala, Galactic PeV neutrinos from dark matter annihilation, Phys. Rev. D89 (2014) 123516, [1404.2932].
  79. F. W. Stecker and S. T. Scully, Propagation of Superluminal PeV IceCube Neutrinos: A High Energy Spectral Cutoff or New Constraints on Lorentz Invariance Violation, Phys. Rev. D90 (2014) 043012, [1404.7025].
  80. M. Ibe and K. Kaneta, Cosmic neutrino background absorption line in the neutrino spectrum at IceCube, Phys. Rev. D90 (2014) 053011, [1407.2848].
  81. T. Araki, F. Kaneko, Y. Konishi, T. Ota, J. Sato and T. Shimomura, Cosmic neutrino spectrum and the muon anomalous magnetic moment in the gauged model, Phys. Rev. D91 (2015) 037301, [1409.4180].
  82. A. N. Akay, O. Cakir, Y. O. Gunaydin, U. Kaya, M. Sahin and S. Sultansoy, New IceCube data and color octet neutrino interpretation of the PeV energy events, Int. J. Mod. Phys. A30 (2015) 1550163, [1409.5896].
  83. E. Aeikens, H. Päs, S. Pakvasa and P. Sicking, Flavor ratios of extragalactic neutrinos and neutrino shortcuts in extra dimensions, JCAP 1510 (2015) 005, [1410.0408].
  84. J. I. Illana, M. Masip and D. Meloni, A new physics interpretation of the IceCube data, Astropart. Phys. 65 (2015) 64–68, [1410.3208].
  85. C. S. Fong, H. Minakata, B. Panes and R. Zukanovich Funchal, Possible Interpretations of IceCube High-Energy Neutrino Events, JHEP 02 (2015) 189, [1411.5318].
  86. F. W. Stecker, S. T. Scully, S. Liberati and D. Mattingly, Searching for Traces of Planck-Scale Physics with High Energy Neutrinos, Phys. Rev. D91 (2015) 045009, [1411.5889].
  87. A. DiFranzo and D. Hooper, Searching for MeV-Scale Gauge Bosons with IceCube, Phys. Rev. D92 (2015) 095007, [1507.03015].
  88. G. Tomar, S. Mohanty and S. Pakvasa, Lorentz Invariance Violation and IceCube Neutrino Events, JHEP 11 (2015) 022, [1507.03193].
  89. P. Di Bari, P. O. Ludl and S. Palomares-Ruiz, Unifying leptogenesis, dark matter and high-energy neutrinos with right-handed neutrino mixing via Higgs portal, 1606.06238.
  90. U. K. Dey, S. Mohanty and G. Tomar, Leptoquarks: 750 GeV Diphoton Resonance and IceCube Events, 1606.07903.
  91. P. S. B. Dev, D. K. Ghosh and W. Rodejohann, R-parity Violating Supersymmetry at IceCube, Phys. Lett. B762 (2016) 116–123, [1605.09743].
  92. E. Fermi, On the Origin of the Cosmic Radiation, Phys. Rev. 75 (1949) 1169–1174.
  93. T. Kashti and E. Waxman, Astrophysical Neutrinos: Flavor Ratios Depend on Energy, Physical Review Letters 95 (Oct., 2005) 181101, [astro-ph/0507599].
  94. J. G. Learned and S. Pakvasa, Detecting tau-neutrino oscillations at PeV energies, Astropart. Phys. 3 (1995) 267–274, [hep-ph/9405296].
  95. A. Bhattacharya, R. Gandhi, W. Rodejohann and A. Watanabe, The Glashow resonance at IceCube: signatures, event rates and vs. interactions, JCAP 1110 (2011) 017, [1108.3163].
  96. S. L. Glashow, Resonant scattering of antineutrinos, Phys. Rev. 118 (Apr, 1960) 316–317.
  97. V. S. Berezinsky and A. Z. Gazizov, Neutrino - electron scattering at energies above the W boson production threshold, Sov. J. Nucl. Phys. 33 (1981) 120–125.
  98. R. Gandhi, C. Quigg, M. H. Reno and I. Sarcevic, Ultrahigh-energy neutrino interactions, Astropart. Phys. 5 (1996) 81–110, [hep-ph/9512364].
  99. H. Athar, C. S. Kim and J. Lee, The Intrinsic and oscillated astrophysical neutrino flavor ratios, Mod. Phys. Lett. A21 (2006) 1049–1066, [hep-ph/0505017].
  100. J. F. Beacom and J. Candia, Shower power: Isolating the prompt atmospheric neutrino flux using electron neutrinos, JCAP 0411 (2004) 009, [hep-ph/0409046].
  101. S. L. Glashow, Resonant Scattering of Antineutrinos, Phys. Rev. 118 (1960) 316–317.
  102. V. Barger, L. Fu, J. G. Learned, D. Marfatia, S. Pakvasa and T. J. Weiler, Glashow resonance as a window into cosmic neutrino sources, Phys. Rev. D90 (2014) 121301, [1407.3255].
  103. F. Halzen, “Particle physics beyond laboratory energies.” Phenomenology 2016 Symposium, University of Pittsburgh, USA, 2016.
  104. S. Razzaque, The Galactic Center Origin of a Subset of IceCube Neutrino Events, Phys. Rev. D88 (2013) 081302, [1309.2756].
  105. M. Ahlers and K. Murase, Probing the Galactic Origin of the IceCube Excess with Gamma-Rays, Phys. Rev. D90 (2014) 023010, [1309.4077].
  106. Y. Bai, R. Lu and J. Salvado, Geometric Compatibility of IceCube TeV-PeV Neutrino Excess and its Galactic Dark Matter Origin, JHEP 01 (2016) 161, [1311.5864].
  107. C. Lunardini, S. Razzaque, K. T. Theodoseau and L. Yang, Neutrino Events at IceCube and the Fermi Bubbles, Phys. Rev. D90 (2014) 023016, [1311.7188].
  108. P. Padovani and E. Resconi, Are both BL Lacs and pulsar wind nebulae the astrophysical counterparts of IceCube neutrino events?, Mon. Not. Roy. Astron. Soc. 443 (2014) 474–484, [1406.0376].
  109. M. Ahlers and F. Halzen, Pinpointing Extragalactic Neutrino Sources in Light of Recent IceCube Observations, Phys. Rev. D90 (2014) 043005, [1406.2160].
  110. Y. Bai, A. J. Barger, V. Barger, R. Lu, A. D. Peterson and J. Salvado, Neutrino Lighthouse at Sagittarius A*, Phys. Rev. D90 (2014) 063012, [1407.2243].
  111. R. Moharana and S. Razzaque, Angular correlation of cosmic neutrinos with ultrahigh-energy cosmic rays and implications for their sources, JCAP 1508 (2015) 014, [1501.05158].
  112. K. Emig, C. Lunardini and R. Windhorst, Do high energy astrophysical neutrinos trace star formation?, JCAP 1512 (2015) 029, [1507.05711].
  113. IceCube, VERITAS collaboration, M. Santander, Searching for TeV gamma-ray emission associated with IceCube high-energy neutrinos using VERITAS, PoS ICRC2015 (2016) 785, [1509.00517].
  114. A. Neronov and D. V. Semikoz, Evidence the Galactic contribution to the IceCube astrophysical neutrino flux, Astropart. Phys. 75 (2016) 60–63, [1509.03522].
  115. L. S. Miranda, A. R. de León and S. Sahu, Blazar origin of some IceCube events, Eur. Phys. J. C76 (2016) 402, [1510.00048].
  116. IceCube, Pierre Auger, Telescope Array collaboration, M. G. Aartsen et al., Search for correlations between the arrival directions of IceCube neutrino events and ultrahigh-energy cosmic rays detected by the Pierre Auger Observatory and the Telescope Array, JCAP 1601 (2016) 037, [1511.09408].
  117. A. Neronov and D. Semikoz, Galactic and extragalactic contributions to the astrophysical muon neutrino signal, Phys. Rev. D93 (2016) 123002, [1603.06733].
  118. L. A. Anchordoqui, M. M. Block, L. Durand, P. Ha, J. F. Soriano and T. J. Weiler, Evidence for a break in the spectrum of astrophysical neutrinos, 1611.07905.
  119. IceCube collaboration, M. G. Aartsen et al., Atmospheric and astrophysical neutrinos above 1 TeV interacting in IceCube, Phys. Rev. D91 (2015) 022001, [1410.1749].
  120. M. Chianese, G. Miele and S. Morisi, Dark Matter interpretation of low energy IceCube MESE excess, 1610.04612.
  121. E. Waxman and J. N. Bahcall, High-energy neutrinos from astrophysical sources: An Upper bound, Phys. Rev. D59 (1999) 023002, [hep-ph/9807282].
  122. J. N. Bahcall and E. Waxman, High-energy astrophysical neutrinos: The Upper bound is robust, Phys. Rev. D64 (2001) 023002, [hep-ph/9902383].
  123. S. M. Boucenna, M. Chianese, G. Mangano, G. Miele, S. Morisi, O. Pisanti et al., Decaying Leptophilic Dark Matter at IceCube, JCAP 1512 (2015) 055, [1507.01000].
  124. K. Harigaya, M. Kawasaki, K. Mukaida and M. Yamada, Dark Matter Production in Late Time Reheating, Phys. Rev. D89 (2014) 083532, [1402.2846].
  125. Y. Kurata and N. Maekawa, Averaged Number of the Lightest Supersymmetric Particles in Decay of Superheavy Particle with Long Lifetime, Prog. Theor. Phys. 127 (2012) 657–664, [1201.3696].
  126. R. Allahverdi and M. Drees, Production of massive stable particles in inflaton decay, Phys. Rev. Lett. 89 (2002) 091302, [hep-ph/0203118].
  127. R. Allahverdi and M. Drees, Thermalization after inflation and production of massive stable particles, Phys. Rev. D66 (2002) 063513, [hep-ph/0205246].
  128. K. Griest and M. Kamionkowski, Unitarity Limits on the Mass and Radius of Dark Matter Particles, Phys. Rev. Lett. 64 (1990) 615.
  129. K. Agashe, Y. Cui, L. Necib and J. Thaler, (In)direct Detection of Boosted Dark Matter, JCAP 1410 (2014) 062, [1405.7370].
  130. J. Berger, Y. Cui and Y. Zhao, Detecting Boosted Dark Matter from the Sun with Large Volume Neutrino Detectors, JCAP 1502 (2015) 005, [1410.2246].
  131. K. Kong, G. Mohlabeng and J.-C. Park, Boosted dark matter signals uplifted with self-interaction, Phys. Lett. B743 (2015) 256–266, [1411.6632].
  132. H. Alhazmi, K. Kong, G. Mohlabeng and J.-C. Park, Boosted Dark Matter at the Deep Underground Neutrino Experiment, 1611.09866.
  133. IceCube collaboration, M. G. Aartsen et al., Observation and Characterization of a Cosmic Muon Neutrino Flux from the Northern Hemisphere using six years of IceCube data, 1607.08006.
  134. Fermi-LAT collaboration, M. Ajello et al., Fermi-LAT Observations of High-Energy -Ray Emission Toward the Galactic Center, Astrophys. J. 819 (2016) 44, [1511.02938].
  135. Fermi-LAT collaboration, M. Ackermann et al., The spectrum of isotropic diffuse gamma-ray emission between 100 MeV and 820 GeV, Astrophys. J. 799 (2015) 86, [1410.3696].
  136. KASCADE-Grande collaboration, Z. Feng, D. Kang and A. Haungs, Limits on the isotropic diffuse -rays at ultra high energies measured with KASCADE, PoS ICRC2015 (2016) 823.
  137. AGRAPES-3 collaboration, S. K. Gupta et al., The current status of the GRAPES-3 extensive air shower experiment, Nucl. Phys. Proc. Suppl. 196 (2009) 153–158.
  138. M. Ahlers and K. Murase, Probing the Galactic Origin of the IceCube Excess with Gamma-Rays, Phys. Rev. D90 (2014) 023010, [1309.4077].
  139. T. Cohen, K. Murase, N. L. Rodd, B. R. Safdi and Y. Soreq, Gamma-ray Constraints on Decaying Dark Matter and Implications for IceCube, 1612.05638.
  140. E. Izaguirre, G. Krnjaic and B. Shuve, The Galactic Center Excess from the Bottom Up, Phys. Rev. D90 (2014) 055002, [1404.2018].
  141. G. D’Ambrosio, G. F. Giudice, G. Isidori and A. Strumia, Minimal flavor violation: An Effective field theory approach, Nucl. Phys. B645 (2002) 155–187, [hep-ph/0207036].
  142. F. Kahlhoefer, K. Schmidt-Hoberg, T. Schwetz and S. Vogl, Implications of unitarity and gauge invariance for simplified dark matter models, JHEP 02 (2016) 016, [1510.02110].
  143. Planck collaboration, P. A. R. Ade et al., Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. 594 (2016) A13, [1502.01589].
  144. M. R. Buckley, D. Feld and D. Goncalves, Scalar Simplified Models for Dark Matter, Phys. Rev. D91 (2015) 015017, [1410.6497].
  145. M. Chala, F. Kahlhoefer, M. McCullough, G. Nardini and K. Schmidt-Hoberg, Constraining Dark Sectors with Monojets and Dijets, JHEP 07 (2015) 089, [1503.05916].
  146. ATLAS collaboration, T. A. collaboration, Search for New Phenomena in Dijet Events with the ATLAS Detector at =13 TeV with 2015 and 2016 data, .
  147. M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, P. Panci et al., PPPC 4 DM ID: A Poor Particle Physicist Cookbook for Dark Matter Indirect Detection, JCAP 1103 (2011) 051, [1012.4515].
  148. M. Cirelli, E. Moulin, P. Panci, P. D. Serpico and A. Viana, Gamma ray constraints on Decaying Dark Matter, Phys. Rev. D86 (2012) 083506, [1205.5283].
  149. J. F. Navarro, C. S. Frenk and S. D. M. White, The Structure of cold dark matter halos, Astrophys. J. 462 (1996) 563–575, [astro-ph/9508025].
  150. M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, P. Panci et al., PPPC 4 DM ID: A Poor Particle Physicist Cookbook for Dark Matter Indirect Detection, JCAP 1103 (2011) 051, [1012.4515].
  151. A. Belyaev, N. D. Christensen and A. Pukhov, CalcHEP 3.4 for collider physics within and beyond the Standard Model, Comput. Phys. Commun. 184 (2013) 1729–1769, [1207.6082].
  152. Planck collaboration, P. A. R. Ade et al., Planck 2013 results. XVI. Cosmological parameters, Astron. Astrophys. 571 (2014) A16, [1303.5076].
  153. O. Mena, S. Palomares-Ruiz and A. C. Vincent, Flavor Composition of the High-Energy Neutrino Events in IceCube, Phys. Rev. Lett. 113 (2014) 091103, [1404.0017].
  154. A. Esmaili, A. Ibarra and O. L. G. Peres, Probing the stability of superheavy dark matter particles with high-energy neutrinos, JCAP 1211 (2012) 034, [1205.5281].
  155. T. Sjostrand, S. Mrenna and P. Z. Skands, A Brief Introduction to PYTHIA 8.1, Comput. Phys. Commun. 178 (2008) 852–867, [0710.3820].
  156. H.-L. Lai, M. Guzzi, J. Huston, Z. Li, P. M. Nadolsky, J. Pumplin et al., New parton distributions for collider physics, Phys. Rev. D82 (2010) 074024, [1007.2241].
  157. R. Gandhi, C. Quigg, M. H. Reno and I. Sarcevic, Neutrino interactions at ultrahigh-energies, Phys. Rev. D58 (1998) 093009, [hep-ph/9807264].
  158. A. Cooper-Sarkar, P. Mertsch and S. Sarkar, The high energy neutrino cross-section in the Standard Model and its uncertainty, JHEP 08 (2011) 042, [1106.3723].
  159. M. D. Kistler and R. Laha, Multi-PeV Signals from a New Astrophysical Neutrino Flux Beyond the Glashow Resonance, 1605.08781.
  160. A. Ishihara and I. Collaboration, Extremely high energy neutrinos in six years of icecube data, Journal of Physics: Conference Series 718 (2016) 062027.