IceCube Events from Heavy DM decays through the Right-handed Neutrino Portal

The recently observed IceCube PeV events could be due to heavy dark matter (DM) decay. In this paper, we propose a simple DM model with extra $U(1)_X$ gauge symmetry and bridge it with standard model particles through heavy right-handed neutrino. The Dirac fermion DM $\chi$ with mass ~5 PeV can dominantly decay into a dark Higgs ($\phi$), the SM Higgs ($h$) and a neutrino ($\nu$). If the lifetime of $\chi$ is ~O($10^{28}$) sec, the resulting neutrino flux can fit data consistently. The neutrino flux from $\chi \rightarrow \phi h \nu$ in our model is softer than the one predicted from $\chi \rightarrow \nu h$, for example. We also discuss a possible mechanism to produce DM with the right relic abundance.


I. INTRODUCTION
Recently, the IceCube Collaboration has reported the detection of 37 neutrino events with energy between 30 TeV − 2 PeV, of which three have energy above 1 PeV [1][2][3]. According to the recent analyses [4][5][6][7], the three-year data are consistent with equal fluxes of all three neutrino flavors and with isotropic arrival directions. However, the neutrino flux required to fit the data in 100TeV − PeV range is around 10 −8 GeV cm −2 s −1 sr −1 per flavor, and rejects a purely atmospheric explanation at 5.7σ. Therefore astrophysical and/or new physics explanations have been pursued for the origin of these high energy neutrinos.
Possible astrophysical sources are involved with supernova remnants (SNR) [8][9][10][11][12], active galactic nuclei (AGN) [13][14][15][16], and gamma-ray bursts (GRB) [17,18], all of which assume some specific emission spectra due to different production environments. In a model independent analysis, the best-fit power law spectrum from the IceCube analysis [1] is It is not very straightforward to fit such a spectrum by astrophysical sources. In all cases, extragalactic sources are needed due to the isotropic feature and galactic constraints [19][20][21][22]. Then identifying such astrophysical sources will be crucial for further understanding of IceCube events.
In this paper, we propose a simple DM model to explain the IceCube PeV events. A dark sector with new U (1) X gauge symmetry is introduced and can have connection with the standard model (SM) sector through neutrino-portal interactions as well as the Higgs portal interaction. Fermionic DM (χ) has ∼ PeV mass and mostly decays into three-body final state with dark Higgs (φ), SM Higgs (h) and active neutrino. The produced neutrinos from primary χ decay and the secondary h and φ decays can explain the observed PeV event spectra, while the atmospheric and astrophysical neutrinos are included for the low-energy part.
This paper is organized as follows. In Section. II we introduce our DM model with dark U (1) X gauge symmetry, heavy right-handed neutrino portal and Higgs portal interactions. In Section. III we outline the general formalism for calculating neutrino flux from galactic and extragalactic DM decay. In Section. IV we present both total and differential decay width for the relevant three-body decay in our model and compare the numerical results with IceCube data. In Section. V we discuss a possible mechanism to generate the correct relic density for DM and direct/indirect detection constraints. Finally, we give our conclusions.

II. MODEL
We consider a dark sector with a dark Higgs field Φ and a Dirac fermion DM χ associated U (1) X gauge symmetry. Their U (1) X charges are assigned as follows 2 : We begin with the following renormalizable and gauge invariant Lagrangian including just one singlet right-handed (RH) neutrino N and one lepton flavor (more N s and/or flavors can be easily generalized): and is the kinetic mixing parameter. Two types of Yukawa couplings, y and f , can be taken as real parameters, ignoring CP violation for simplicity. We define covariant derivative as D µ = ∂ µ − ig X X µ . Since we are interested in explaining the IceCube PeV events in terms of DM χ decay, we shall take m χ ∼ PeV. Other parameters in our model are free variables.
The scalar potential V of this model is given by Both electroweak and dark gauge symmetries are spontaneously broken by the nonzero vacuum expectations values of H and Φ: Here v H 246GeV is the same as SM value but v φ might be taken as a free parameter. In the unitarity gauge, we can replace the scalar fields with Note that h and φ shall mix with each other thanks to the Higgs-portal operator, (the λ φH term) 3 . Through this mixing, φ can decay into SM particles. Another important mixing happens among three neutral gauge bosons, photon A µ , Z µ and X µ . Such a mixture would enable an extra mass eigenstate Z µ (mostly X µ ) to decay SM fermion pairs. Then DM χ scattering off nucleus is possible by the Z exchange, and the cross section essentially depends on , v φ , m Z . It is easy to choose small , or heavy masses to evade the constraints from DM direct detection [54]. When the right-handed neutrino N is much heavier than χ, we can integrate it out and obtain an effective operator, yf m Nχ ΦH † L + h.c., (2.4) 2 A similar setup with different dark charge assignments has been considered for the AMS02 positron excess [48]. One may also use discrete symmetries, see Ref. [49] for example. 3 The λ hφ term can also help to stabilize the electroweak vacuum [50][51][52][53]. which would make χ decay possible but long lived. After spontaneous gauge symmetry breaking, we have several higher dimensional effective operators from the aforementioned operator Eq. (2.4) as follows: with the common factor yf 2 for all these operators. If kinematically allowed, all the above operators induce χ decays into different channels with fixed relative branching ratios. Within the heavy χ limit, m χ m φ , m Z , m h , m Z , m W , the mass operatorχν in Eq. (2.5) would induce a tiny mixing between χ and ν with the mixing angle β approximately given by Then the gauge interactions for χ and ν will generate the decay channels, with their branching ratios being proportional to ∼ v 2 H : v 2 φ : 2v 2 φ . Two dim-4 operators, χhν andχφν, would lead χ to the following decays, with their branching ratios being proportional to ∼ v 2 φ : v 2 H . It is also straightforward to get the following relation for the branching ratios, Therefore, all the decay branching ratios are basically calculable and completely fixed in this model 4 . Note that the decay modes with Z or φ are unique features of DM models with dark gauge symmetries 5 . Another interesting phenomenon in this model is that three body decay channel χ → φhν is dominant over all other channels when m χ v φ : since we actually have an enhancement from heavy m χ even though there is a phase space suppression from three-body final states. There are another three-body decay channels that are equally important: with branching ratios 1 : 1 : 2 due to the Goldstone boson equivalence theorem. In the following, if not otherwise stated explicitly, we use χ → φhν to represent all these channels and in numerical calculations we take all of them into account.
To give a rough impression for the relevant parameter ranges, we can perform an orderof-magnitude estimation (complete formulas and details of calculation are given in the Appendix): which shows the required value for the combination yf /m N . If one additionally assumes that N should be responsible for neutrino mass through the usual Type-I see-saw mechanism, then we have y ∼ 10 −5 m N /PeV. Therefore we would have y ∼ 1 and f ∼ 10 −22 for m N ∼ 10 14 GeV and, y ∼ 10 −5 and f ∼ 10 −25 for m N ∼ PeV. In any case, f is very tiny and seems unnaturally small. However, it is still natural a la 't Hooft [55] since taking f = 0 enhances symmetries of the theory, namely DM number conservation 6 . In later discussion, we shall take y, f, m N as free parameters, unless specified.

III. NEUTRINO FLUX FROM DM DECAY
The neutrino flux from dark matter decay is composed of galactic and extragalactic contributions which are equally important as we shall see below. Galactic neutrino flux at kinetic energy E from DM decay in our Milky Way dark halo is given by where Γ i is partial width for decay channel i, dN i ν /dE ν is the neutrino spectrum at production, r = r 2 + r 2 − 2r r cos θ, r is the distance to earth from the DM decay point, r 8.5kpc for the solar system and θ is the observation angle between the line-of-sight and the center of the Milky Way. For the galactic DM density distribution, we use the following standard NFW profile [56], with parameters r c 20 kpc and ρ 0.4 GeV/cm 3 . For decaying dark matter the flux is not very sensitive to DM density profile, so our discussions and results will still apply if another different profile is used.
We can also get the extragalactic or cosmic contribution from a formula similar to the above one, by taking cosmic expansion into account, namely the red-shift effect [24]: where E is red-shifted to E as E = (1 + z)E, the red-shift z is defined as 1 + z = a 0 /a with present scale factor a 0 being normalized to 1, the critical energy density ρ c = 5.5 ×  10 −6 GeV/cm 3 and Ω χ 0.27 is DM χ's fraction. The Hubble parameter H is related to its present value through Ω λ , Ω m and Ω r are energy fractions of dark energy, all matter, and radiations, respectively. We shall use the latest results from Planck [57] for numerical evaluation.

IV. NUMERICAL RESULTS
To compute the neutrino flux from DM decay, we first need to calculate the total and differential three-body decay width for χ → φ+h+ν. In the heavy χ limit, we have obtained the total width and normalized differential decay widths The details of the calculation are given in the Appendix where complete formulas with nonzero mass parameters are also presented. The above differential widths are essential ingredients to get the final neutrino flux. For example, νs from different decay final states are given by where x = ν, h, W, Z, Z , φ. Note that dN ν (E x ) /dE in the integrand can be calculated with Pythia [58] or PPPC4DMID [59].
In Fig. 1 we show the neutrino spectra. Just for illustration, we choose the mass of DM χ around 5 PeV and its lifetime 2 × 10 28 s. The neutrino spectra are multiplied by E 2 ν in order to account for the energy dependence of the neutrino-nucleus cross section so that it might be compared with IceCube data more easily. The blue dot-dashed curve indicates the spectrum from final ν in χ's decay and red dashed one marks h/φ's contribution. For simplicity, the mass of φ has been chosen to be just as the SM Higgs mass. However, other choice does not affect our result much, since in the left panel of Fig. 1 we see that in the high energy part h/φ's contributions are basically negligible and it is the high energy part that explains the IceCube PeV events. In the right panel, we include the red-shifted effects and show both galactic (blue dashed) and extragalactic (red dot-dashed) contributions.
Next, we compare our model predictions with IceCube three-year data [1]. To parameterize the possible astrophysical neutrino fluxes at low energy, we consider either a broken power law (BPL) or unbroken power law (UPL) [32], where the first one has an exponential cut-off at energy scale E 0 which is chosen to be 125 TeV in agreement with the SNR results [8].
We illustrate in Fig. 2 with two different astrophysical flux choices. The low energy data are best fitted by varying J 0 and the spectral index γ, but the high energy PeV data points are fit by DM decay. In the left panel, we have used J BPL 0 = 4.1 × 10 −8 GeV/cm 2 /s/sr and γ 1 = 0(red dot-dashed curve). In the right panel, we have used J UPL 0 = 1.3 × 10 −8 GeV/cm 2 /s/sr and γ 2 = 0.7(red dot-dashed curve). DM's contributions and total flux are labeled with purple dashed and blue solid curves, respectively. As we can see in the figure, our model can agree with the PeV data. Also a gap could appear around 400 TeV although the current data can not tell existence of the gap is statistically significant. The feature for the DM decay spectrum is that there would be a sudden drop around E ν m χ /2. As more data are accumulated, it should be possible to test our model in the future.
We should note that in our discussion J 0 and γ are just adjusted visually to be consistent with the low-energy data points. Dedicated investigation would require global fitting, which is beyond our scope here. Just for comparison, IceCube [60] gives the best-fit parameters of a single unbroken power law for neutrinos energies between 25TeV and 2.8PeV without DM contribution, J 0 = 6.7 +1.1 −1.2 × 10 −8 GeV/cm 2 /s/sr and γ = 0.5 ± 0.09. In Fig. 3, we also compare our model with the preliminary updated results [61] based on IceCube 4-year data which has already filled the gap a bit. Here we have only shifted to J BPL 0 = 5.6 × 10 −8 GeV/cm 2 /s/sr and τ χ ∼ 1.5 × 10 28 s (left panel), and J UPL 0 = 2.1 × 10 −8 GeV/cm 2 /s/sr and τ χ ∼ 2 × 10 28 s (right panel). In Fig. 4, we illustrate two cases, one with m χ = 8PeV, τ χ ∼ 1.7 × 10 28 s, and the other m χ = 10PeV, τ χ ∼ 1.5 × 10 28 s. There are some crucial differences between our model and some others in the literature. For example, the authors in Ref. [23,29] considered the effective operator, yL Hχ with y ∼ 10 −30 , which induces mainly two-body decay of DM χ, χ → νh, νZ, l ± W µ .
In this scenario, the neutrino spectrum shows that there should be no gap between 400 TeV ∼ 1 PeV [26]. Our model predicts that the dominant decay mode are . which is a consequence of U (1) X dark gauge symmetry and the dark charge assignments of the dark Higgs and dark matter fermion χ. The neutrino spectra from primary χ decay and the secondary decays of h and φ have different shapes and could account for the possible gap. However, we should note that the current data can not favor one over another yet due to its low statistics. Also the neutrino flux in our model is softer than the one predicted in Ref. [23,29], for example. In Ref. [32], leptophilic three-body decay induced by dimension-sixL α l βLγ χ was considered with global U (1) or A 4 flavor symmetries. Besides the neutrino spectrum difference, our model involves an additional gauge boson which mediates the DM-nucleon scattering, and could be tested by DM direct searches.
Our scenario is also different from those in which DM decay is also responsible for the low-energy flux [24]. The DM lifetime in Ref. [24] should be around 2 × 10 27 s, as mainly determined by the low energy part of events. This is partly due to the reason that the branching ratio into neutrinos and bb there should be about 10% and 90%, respectively, to account for the possible gap. On the other hand, in our scenario 1/2 of the decay channels have prompt neutrinos. Another main difference is that three-body-decay usually gives broader spectra at PeV range than two-body-decay considered in Ref. [24], but more data is required in order to discriminate this difference.
Assuming the dark photon X µ is much heavier than Z, the DM-nucleon scattering cross section can be roughly estimated as 10 −39 cm 2 is the typical cross section value for SM Z-mediating DM-nucleon process. Comparing it with the direct detection bound for 100GeV DM, we should have for heavy m χ ∼ 5 PeV. This can be easily satisfied, for example, with m X ∼ TeV and sin 0.1.

V. RELIC ABUNDANCE AND CONSTRAINTS
In our above investigation, we have not discussed the relic abundance for DM χ yet. Since χ is very heavy, the unitary bound on its annihilation makes the DM χ impossible to be thermally produced and a non-thermal process is needed (see Ref. [62] for a recent review on such topics). Here, we discuss one possible non-thermal production mechanism for DM χ in our model. We assume that the dark Higgs φ 7 once shared a common temperature with SM particles and had a thermal distribution when its temperature T was larger than m χ . Here we do not specify the mechanism how φ reached such a temperature; it could be due to reheating after inflation or some heavy particle decays.
In the thermal bath, χ could be produced through φ + φ → χ +χ, whose thermal cross section is given by Here we considered the m N T case only. We can calculate the χ's yield, Y χ ≡ n χ+χ /s, where s ∼ g * s (T )T 3 is the total entropy density in the Universe ( g * s (T ) ∼ 100), In the above derivation we have used the Hubble parameter H g * (T )T 2 /M pl , g * (T ) ∼ 100 is the total effective number of degree of freedom when φ's temperature is T and Planck mass M pl = 3/8πG 4.2 × 10 18 GeV. Since the yield has a positive power dependence on temperature, χ is mostly produced at φ's highest temperature T φ max , .

(5.2)
Requiring χ give the correct relic density, we have a relation which puts a constraint on f and m N , When m N > T φ max > m χ and m χ ∼ PeV, we are able to give a lower bound for f |f | 10 −6 .
One more thing we can infer from Eq. (2.11) and (5.4) is that in this production mechanism y would be too small such that this right-hand neutrino N can not be fully responsible for active neutrino mass and mixing angles. This is because, in order to explain the active neutrino mass, we should have which can not be satisfied simultaneously with Eq. (2.11) and (5.4). This is not a problem for this model since we can expect there are additional right-handed neutrinos N i and they can couple toL H with large y i but toχΦ with tiny f i , so that they are just responsible for active neutrino mass and mixing angle but not for DM χ's production. DM direct detection can constrain DM-nucleon scattering cross section, whose value in our model is determined by the kinetic mixing parameter , gauge coupling g X and the mass of Z . For GeV-TeV DM, there is already plenty of viable parameter space to evade such a constraint, see Ref. [63] for example. For PeV DM, the constraint is even relaxed due to the low number density, see Eq. 4.9. Indirect detection from positron, anti-proton and γ rays also constrain DM χ's decay lifetime, see Refs. [64][65][66] and references therein for example. We have checked that τ χ ∼ 10 28 sec is still allowed by all such constraints.
As an illustration, in Fig. 5 we show the expected gamma-ray flux from DM decay with m χ ∼ 5PeV and lifetime τ χ ∼ 2 × 10 28 s. We have included the prompt gamma-rays and those from inverse compton scattering (ICS) of charged particles on CMB, starlight and dust-rescattered light. For extragalatic contribution, we have taken absorption factor into account. As we can easily see, the flux is well below the current constraints from Fermi-LAT [67] and KASCADE [68]. Our results are also consistent with the recent investigations about gamma-ray constraints on the lifetime of PeV DM in different scenarios [36,69], τ χ 3 × 10 27 .

VI. CONCLUSION
In this paper, we have proposed a dark matter (DM) model that can explain the IceCube PeV events in terms of DM decay. The model is based on an extra U (1) X dark gauge symmetry which is spontaneously broken by a dark Higgs field, Φ. One crucial bridge between DM χ and standard model (SM) particles is established by heavy right-handed (RH) neutrino portal interactions. This heavy neutrino can induce DM decays into SM particles, including the light neutrinos.
The dominant decay channel of DM is the three-body final state with SM Higgs, dark Higgs, and neutrino (χ → φ + h + ν) (and other channels due to the Goldstone boson equivalence theorem), and not the usual two-body decays such as χ → Zν, W ± l ∓ , etc.. This is a unique feature of the present model based on U (1) X dark gauge symmetry and the RH neutrino portal interactions. We have calculated both total and differential decay width to evaluate the galactic and extragalactic neutrino fluxes. We have found that neutrino flux from these decay products can agree well with the IceCube spectrum. Together with an astrophysical flux for lower energy events, we are able to fit IceCube data around 100 TeV ∼ 2 PeV if we assume DM mass is about m χ ∼ 5PeV and its lifetime is τ χ ∼ 2 × 10 28 sec.

APPENDIX
Here we show the complete differential decay width for χ → h + φ + ν. Throughout the calculation, we work in the rest frame of χ, so χ's momentum is (m χ , 0, 0, 0). For unpolarized χ, we have where we have used the averaged, squared matrix element, Then we get (with E ν < m χ /2) . and some other kinematic variables, where E * φ and E * ν are the energies of φ and ν in the m hφ rest frame, respectively. From the above differential decay width, we can easily get the total width Γ 1 24 (6.4) The differential decay width as function of E h or E φ can also be calculated similarly. For example, with similar definitions and kinematic bounds, φ and E * h are the energies of φ and h in the m φν rest frame, respectively.