Phenomenology of Colored Radiative Neutrino Mass Model and Its Implications on the Cosmic-ray Observations

We extend the colored Zee-Babu model with a gauged $U(1)_{B-L}$ symmetry and a scalar singlet dark matter (DM) candidate $S$. The spontaneous breaking of $U(1)_{B-L}$ leaves a residual $Z_2$ symmetry that stabilizes the DM and generates tiny neutrino mass at the two-loop level with the color seesaw mechanism. After investigating dark matter and flavor phenomenology of this model systematically, we further focus on its imprint on two of cosmic-ray anomalies: the Fermi-LAT gamma-ray excess at the Galactic Center (GCE) and the PeV ultra-high energy (UHE) neutrino events at the IceCube. We found that the Fermi-LAT GCE spectrum can be well fitted by DM annihilation into a pair of on-shell singlet Higgs mediators while being compatible with the constraints from relic density, direct detections as well as dwarf spheroidal galaxies in the Milky Way. Although the UHE neutrino events at the IceCube could be accounted for by resonance production of a TeV-scale leptoquark, the relevant Yukawa couplings have been severely limited by current low energy flavor experiments. We then derive the IceCube limits on the Yukawa couplings by employing its latest 6-year data.


I. INTRODUCTION
The existence of dark matter (DM) and tiny neutrino mass poses an outstanding challenge to both theoretical and experimental particle physics. Although current searches coming from the Large Hadron Collider (LHC) and DM direct detections have imposed stringent limits, their null results have not yet provided powerful guidance to physics beyond the standard model (SM). On the other hand, observations from high energy cosmic rays (CR) may offer another angle to face the challenge. In this paper, we will focus on two of them, i.e., the Fermi-LAT gamma-ray excess at the Galactic Center (GCE) and the PeV ultra-high energy (UHE) neutrino events at the IceCube. We will attempt to interpret the two observations in a colored seesaw extension of the SM which generates radiative neutrino mass and has a cold DM particle built in. But before we embark on that, let us briefly review the current status of the two observations.
The GCE was first reported in Ref. [1] through analysing the Fermi-LAT data, and the signal significance was confirmed by subsequent analyses [2][3][4][5][6][7][8]. While astrophysical interpretations like millisecond pulsars or unresolved gamma-ray point sources [5,6,[9][10][11] are plausible, DM annihilation remains one of popular interpretations because its thermally averaged cross section and morphology of density distribution match the standard WIMP scenario. In particular, Ref. [8] gives a comprehensive and systematic analysis with multiple Galactic gamma ray diffuse emission (GDE) models. Very recently, the Fermi-LAT Collaboration has released their updated analysis [12,13] and concluded that GCE can be caused by an unresolved pulsarlike sources located in the Galactic bulge which they referred to as Galactic bulge population, while the dark matter interpretation is disfavored since its distribution is not consistent with the morphology detected in their analysis. However, a large population of pulsars should be accompanied with a large population of low-mass X-ray binaries in the same region, which turns out to restrict their contribution only up to 4 − 23% of the observed gamma-ray excess [14]. Moreover, analyses of spatial distribution and luminosity function of those sources were inconclusive about the presence of such Galactic bulge population [15]. Therefore, dark matter interpretation of GCE is still competitive.
When using model independent fitting with DM directly annihilated into a pair of SM particles, the GCE spectrum is best fit by the bb final state [7]. The other final states (τ + τ − , qq cc, gg, W + W − , ZZ, hh and tt) with different DM mass and annihilation cross section are also acceptable [16][17][18][19]. Additionally, when taking into account uncertainties in DM halo profiles and propagation models, the annihilation cross section required by GCE is compatible with the limits from other indirect DM searches like dwarf spheroidal galaxies (dSphs) of the Milky Way and the antiproton and CMB observations [19][20][21][22][23]. The DM annihilation explanation of the GCE has attracted great interest in the past few years and has been extensively explored in various new physics models . These models can be classified into two scenarios from annihilation patterns: • DM annihilates directly into SM final states, • DM annihilates into some intermediate particles, which subsequently cascade decay into SM particles.
While the first scenario usually suffers from stringent constraints from DM direct detections and collider searches, the second has the advantage that cascade decays can soften and broaden the resulting photon spectrum, thus considerably enlarging the parameter space and relaxing the experimental constraints. More interestingly, GCE can also be interpreted in DM models with a global or local Z 3 symmetry by invoking semi-annihilation channels [33,42,44].
The IceCube observatory is a neutrino telescope located at the South Pole, and holds the unique window to cosmic UHE neutrinos. In the 4-year data set released in year 2015, a total of 54 UHE neutrino events are collected (including 39 cascade events and 14 muon track events) with 7σ excess over the expected atmospheric background [47]. Particularly, three events with an energy above PeV present a bit of excess on the SM prediction [48][49][50]. Very recently, the IceCube Collaboration has published the preliminary 6-year result [51], with the total number of events increased to 82 with 28 of them being observed in the recent two years. Note that all of new events have energies below 200 TeV, and the excess in the PeV range still exists. The origin of these PeV UHE neutrino events remains mysterious and immediately causes great interest in both astrophysics and particle physics communities. While the astrophysics community focuses on various astrophysical sources [52][53][54], the particle physics community tries to relate them to new physics phenomena. For instance, in the models of decaying superheavy DM [55][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70][71] 1 , a DM particle of PeV mass is required in order to reproduce the desired UHE neutrino events. Such superheavy particles are very difficult to probe in other experiments and thus phenomenologically less interesting. Another possible explanation invokes a new particle resonance in the TeV region [75][76][77][78][79][80][81][82], in accord with the common belief that new physics should appear there. This latter scenario appears phenomenologically advantageous and could be examined with other means, in particular by direct searches at the LHC.
The six orders of magnitude difference in the energy scale between the GCE (GeV) and IceCube (PeV) events makes it challenging to explain them in a single framework. Here we present a novel example for this issue. We extend the colored Zee-Babu model [83] with a U (1) B−L gauge symmetry and a singlet scalar DM candidate. Another singlet Higgs scalar associated with the U (1) B−L symmetry serves as an onshell mediator for DM annihilation resulting in the GCE spectrum, while the leptoquark (LQ) is responsible for the resonance production of extra UHE neutrino events. The same singlet Higgs scalar and leptoquark generates tiny neutrino mass at two loops. In the next section we describe the model and discuss relevant experimental constraints on its parameter space. Sections III and IV include the core contents of this work, in which the DM properties, GCE spectrum and UHE neutrino event rate at IceCube are systematically investigated. In section III, we explore the vast parameter space that satisfies the constraints from relic abundance and direct detections, and discuss the dominant annihilation channels. A comprehensive fit to the GCE spectrum is then presented incorporating all these limits. In section IV A, we calculate the SM and LQ contributions to the neutrino-nucleon scattering cross section. Then in section IV B, we estimate the LQ contribution to the UHE neutrino event rate at IceCube and perform a likelihood analysis to determine the parameter space. Finally, we draw our conclusion in section V.

II. MODEL AND RELEVANT CONSTRAINTS
A. The Model The particle contents and their charge assignments are shown in Table. I. In addition to the LQ ψ and diquark ω, we further introduce two singlet scalars, ϕ with lepton number L = 2 and S with L = 1 2 . Here, ϕ is used to break the U (1) B−L gauge symmetry spontaneously, thus generating the L-breaking trilinear term ψ * ψ * ω required for radiative neutrino masses. Notably, due to the proper charge assignment of S, the U (1) B−L symmetry forbids any gauge invariant terms that would allow S to decay, promoting S a DM candidate without imposing ad hoc discrete symmetry [84][85][86]. In order to make U (1) B−L anomaly free, some fermions neutral under the SM gauge group but with exotic B − L charges other than −1 could be employed [87][88][89][90][91][92].  The relevant Yukawa interactions involving the LQ ψ and the diquark ω are given by (1) where σ 2 is the second Pauli matrix, ij refers to the SM generations, and the color indices are suppressed.
Here, y ω is a symmetric matrix, while y L,R and y ψ are general complex matrices. The neutrinos interact with the LQ only through the y L term, which induces neutrino masses at the two-loop level as shown in Fig. 1. Compared to the original Zee-Babu model, no antisymmetric Yukawa couplings are involved in neutrino mass generation so that all neutrino masses can be non-zero in this colored Zee-babu model. And the y ψ together with the y L,R terms can lead to the tree-level proton decay [93]. In principle, this y ψ term can be forbidden by some discrete symmetry [94]. For simplicity, we will assume y ψ = 0 in the following discussion. Note that due to the charge assignments the two scalar singlets ϕ and S do not couple to fermions at the Lagrangian level.
The gauge invariant scalar potential is described by where µ 2 X (X = Φ, ϕ, S, ψ, ω) are all taken to be positive, and the trace is over the color indices. In this way, the SU (2) L ×U (1) Y and U (1) B−L gauge symmetries are spontaneously broken by the vacuum expectation values of Φ and ϕ, respectively. Due to the B − L charge assignment of S, one can still have S = 0 after spontaneous symmetry breaking, so that a residual Z 2 symmetry remains under which only S is odd. This blocks all potential decays of S, making it a viable DM candidate [84][85][86].
In unitary gauge the scalar fields Φ and ϕ are denoted as Here v φ = 246 GeV is the electroweak scale, and the vacuum expectation value (VEV) v ϕ generates the mass for the new gauge boson Z of U (1) B−L , where g BL is the gauge coupling of U (1) B−L . The LEP bound requires that [95] M Z /g BL = 2v ϕ 7 TeV, yielding a lower limit on v ϕ 3.5 TeV. On the other hand, the direct searches for the Z -boson at LHC in the dilepton channel have excluded M Z 4 TeV [96][97][98], and recasting these searches in the gauged U (1) B−L model has been performed in Refs. [99][100][101] to acquire the exclusion region in the M Z − g BL plane. Considering these bounds, we choose to work with M Z = 4 TeV and g BL = 0.1, so that v ϕ = 20 TeV in our following discussion. The masses of the DM S, LQ ψ and diquark ω can be figured out from the scalar potential in Eq. (2): In this work, we will consider M S in the interval [5,150] GeV and [500, 1500] GeV for the low and high mass region, respectively. The constraints from relic density and direct detections will be discussed in BR(ψ → ν q ) = 0.5 with = e or µ, respectively. ATLAS has also excluded M ψ < 625 GeV when BR(ψ → ν τ b) = 1 for the third generation LQ [109]. As for the scalar diquark ω, CMS has excluded M ω 7 TeV [110,111]. In the following, we will mainly consider M ψ 1 TeV and M ω = 7 TeV to respect these collider limits.
The λ Φϕ term induces mixing between φ 0 and ϕ 0 , with the squared mass matrix given by which is diagonalized to the mass eigenstates (h, by an angle θ determined by with −π/4 < θ < π/4. The masses of h, H 0 are Here h is regarded as the Higgs boson with M h = 125 GeV discovered at LHC [112][113][114]. According to previous studies on scalar singlets, in the high mass region M H 0 > 500 GeV [115][116][117][118], a small mixing angle | sin θ| 0.2 is allowed by various experimental bounds. In light of the recent Fermi-LAT GCE, we will also consider the low mass region M H 0 ∈ [5,150] GeV. In this region, the LHC SM Higgs signal rate measurement has excluded | sin θ| 0.36 [118][119][120], and the LEP search for ZH 0 associated production has excluded | sin θ| 0.2 when H 0 → bb dominates [121]. Thus, it is safe to consider | sin θ| 0.1 in the following discussion. For convenience, we express the Lagrangian parameters λ Φ,ϕ,Φϕ and µ Φ,ϕ in terms of the physical scalar masses M h,H 0 , mixing angle θ as well as the VEVs v φ,ϕ :

B. Neutrino Mass
As shown in FIG. 1, the neutrino masses are induced at two loops [94]: where the full analytical form for the loop function I mn can be found in Ref. [122]. Considering that the down-type quarks are much lighter than the colored scalars, it can be simplified for order of magnitude estimate to where Typically, a neutrino mass m ν ∼ 0.01 eV can be realised with λ ∼ 0.1, y L ∼ y ω ∼ 0.01 when v ϕ = 20 TeV, M b = 4.7 GeV, M ψ = 1 TeV, and M ω = 7 TeV. The radiative correction to the masses M ψ and M ω involves also the trilinear coupling λv ϕ ψ * ψ * ω, the choice of λ ∼ 0.1 and v ϕ = 20 TeV also satisfies the perturbativity requirement λv ϕ 5 min(M ψ , M ω ) for M ψ ∼ 1 TeV and M ω ∼ 7 TeV [123,124]. The neutrino mass in Eq. (20) can be written in a compact form where In principle, by adopting a proper parametrization [125,126], the Yukawa coupling y L can be solved in terms of the neutrino masses, mixing angles and a generalized orthogonal matrix with three free parameters, so that the neutrino oscillation data can be automatically incorporated. Following this approach, a benchmark point has been suggested in Ref. [127]; see Ref. [94] for more details. As to be discussed below, in this work we follow the usual phenomenological practice to take Yukawa components y ij L as input parameters whose values will be constrained by IceCube data and low-energy experiments.
Two-loop generation of neutrino mass.

C. Flavor Constraints
The LQ ψ can induce various flavor violating processes at the tree level. To minimize such processes, one usually assumes y R = 0 [94,127], since the y R term is less important to neutrino masses as well. This also fits our interest in the IceCube UHE neutrino events which may be induced by y L but not y R couplings.
Since LQ is heavy, its effects can be incorporated into effective four-fermion operators of the SM leptons and quarks. The constraints on these operators have been studied in Ref. [128] for the normalized Wilson coefficients: particular, there are two ijkn that are strongly constrained: one is eµuu < 8.5 × 10 −7 from µ-e conversion in nuclei, and the other is uc < 9.4 × 10 −6 from the K-meson decay. This indicates that [129] y eu L y µu One way to satisfy these bounds is to assume, e.g., y u for a TeV scale M ψ [130].
The Yukawa coupling y ij L (L Li ) C iσ 2 Q Lj ψ * is also responsible for lepton flavor violation (LFV) processes at one loop. According to Ref. [94], the constraints from the radiative decay → γ are usually more stringent than other LFV processes, and the branching ratio is calculated as [94,131] where N C = 3 and the LQ-quark loop yields Here r q = M 2 q /M 2 ψ , A L = A R | y L ↔y R , and the loop functions are [94] In the limit x → 0, the loop functions behave as F 1 (x) → 1/12 and F 2 (x) → (7 + 4 ln x)/6 < 0. If y q L ∼ y q R , the second term in Eq. (28) is expected to be dominant, since |M F 1 (r q )| |M q F 2 (r q )|. Hence we assume y R = 0 in numerical analysis partly for minimizing the LQ contribution to lepton radiative decays. With this assumption, A R dominates over A L considering M M , and Eq. (27) simplifies to Currently, the most stringent limits on lepton radiative decays are BR(µ → eγ) < 4.2 × 10 −13 [132], [133], and BR(τ → eγ) < 3.3×10 −8 [133]. They translate into the constraints on the Yukawa couplings q=u,c,t y eq * L y µq q=u,c,t q=u,c,t For a flavor universal structure, the above requires |y q L | 0.02 at M ψ ∼ 1 TeV. On the other hand, a hierarchal structure |y eq L | |y µq L | ∼ |y τ q L | ∼ O(0.1) is still allowed at M ψ ∼ 1 TeV, because radiative τ decays are less stringently constrained [130].
A by-product of lepton radiative decays is the LQ contribution to the anomalous magnetic moment of the charged lepton [134,135] Under constraints from LFV, the predicted values are ∆a e = −2 × 10 −19 , ∆a µ = −1 × 10 −14 , and ∆a τ = −2 × 10 −12 for universal Yukawa couplings |y q L | ∼ 0.01 at M ψ ∼ 1 TeV and assuming y R = 0, which are far below the current experimental limits [136,137]. It is also clear that with the assumption of y R = 0 the observed discrepancy ∆a µ = (27.8 ± 8.8) × 10 −10 [137] cannot be explained, since the contribution of the |y q L | 2 term is negative. To resolve the discrepancy, a nonzero y R is necessary, e.g., with y µc R ∼ y µt R ∼ 0.01, y µc L ∼ 2.4, y µt L ∼ 0.5, and M ψ ∼ 1 TeV [138]. If Im(y q * L y q R ) is nonzero, the LQ also contributes to the electric dipole moment (EDM) of the charged lepton at one loop [135] Typically for |y eq L | ∼ |y eq R | ∼ 0.01, M ψ ∼ 1 TeV, and an order one CP phase, the top quark would dominate and contribute to the electron EDM |d e | ∼ 10 −24 e-cm, which has already been excluded by the current limit |d e | < 8.7 × 10 −29 e-cm [139]. If we still assume y R = 0, the EDM will arise at three loops, whose order of magnitude is [94] d For |y eq L | ∼ 0.01, M ψ ∼ 1 TeV, and an order one combined CP phase, one has |d e | ∼ 10 −37 e-cm, which is much smaller than the current limit.
As for the diquark ω, the Yukawa couplings y ij ω are tightly constrained by neutral mesons mixings [140].

III. DM PHENOMENOLOGY AND GCE SPECTRUM FITTING
For the purpose of illustrating the effects of various annihilation processes on relic abundance and direct detection, we list all important annihilation channels in Fig. 2. A DM pair can annihilate into (1) a b quark  Table III. Several features learned from these results are summarized as follows. 2 For the low mass DM scenario: • There are much less survived samples than for the high mass scenario. This is due to the fact that the coupling between the DM and SM Higgs, λ Sh , is tightly constrained by relic abundance and direct detections. As a consequence, only the channels mediated by the singlet scalar H 0 and t-channel   • Both W + W − and ψψ * channels could be dominant when M S 1.3 TeV, while only the ψψ * channel dominates for M S 1.3 TeV. The reason is that we have chosen the corresponding couplings λ Sh , λ Sψ < 0.5 in our scan.
while the ψψ * channel in the same plane is much scattered.
We now turn to GCE spectrum fitting in our model. The hard photons due to DM annihilation arise mainly from subsequent decays of SM particles, since their direct production is typically loop-suppressed.
The continuous gamma-ray spectrum results from light mesons produced through hadronization and decay of SM fermions. The gamma-ray flux due to DM annihilation in the Galaxy can be expressed as where f sums over all quark and lepton annihilation channels. σv f halo is the thermally averaged annihi-  lation cross section in the Galactic halo, and dN γ f /dE γ the prompt photon spectrum per annihilation for a given final state f . The astrophysical factorJ is expressed as where r(s, ψ) = r 2 + s 2 − 2r s cos ψ. Here r = 8.5 kpc is the Sun-Galactic Center distance, s is the line of sight (l.o.s) distance, and ψ is the angle between the observation direction and the Galactic Center.
In terms of the Galactic latitude and longitude coordinate (b, l), one has cos ψ = cos b cos l.
For a DM interpretation of the GCE, the angular region of interest for the Fermi-LAT is, ∆Ω: 2 • ≤ |b| ≤ 20 • and |l| ≤ 20 • . In our calculation, we take the generalized Navarro-Frenk-White (gNFW) profile for the DM halo distribution [149] ρ(r) = ρ r r where the scale radius r s = 20 kpc. Based on the analyses of Refs. [16,22,29,30], the local DM density ρ and index γ are estimated to be ρ = (0.4 ± 0.2) GeV/cm 3 and γ = 1.2 ± 0.1. We thus choose their central values (ρ , γ) = (0.4 GeV/cm 3 , 1.2) for the benchmark halo profile, which yields the valuē J ben forJ. The uncertainties of (ρ , γ) then translate intoJ ≡ JJ ben , where the factor J ∈ [0. 14, 4.4] parameterizes the allowed range for DM distribution. We will do the GCE scan for J in the above range and J = 1 for the benchmark profile.
To fit the GCE we use the results in Ref. [8], which explored in detail multiple galactic diffuse emission (GDE) models. We employ micrOMEGAs and PPPC4DMID [150] to generate the photon spectrum and perform global fitting by using where dΦ th,obs i /dE γ are respectively the theoretical and observed gamma-ray flux in the i-th energy bin.
Σ ij is the covariance matrix provided by Ref. [8] which includes both statistical and correlated systematic errors. Here we focus on the on-shell mediator scenario, in which DM annihilates into a pair of on-shell singlet scalars H 0 , which in turn decay to the SM quarks and leptons. The decay branching ratios of H 0 are presented in Fig. 6 versus its mass, which have a similar pattern to those of the SM Higgs due to the φ 0 − ϕ 0 mixing. We vary M S , M H 0 in the GCE scan while fixing other parameters as shown in Table IV.
In addition to relic abundance and direct detections, one must take into account the constraint from dwarf spheroidal galaxies (dSphs) in the Milky Way. The lack of gamma-ray excess from dSphs imposes a tight bound on the DM annihilation cross section in the galactic halo, and also gives a stringent constraint on the DM interpretation of GCE for various annihilation channels. Here we adopt dSphs limits provided in Ref. [151], which performed a model-independent and comprehensive analysis on various two-body and The cyan region corresponds to the 2σ ranges allowed by GCE fitting, i.e., for J ∈ [0.14, 4.4], and the green region is for the benchmark halo profile, i.e., J = 1. Scan samples that satisfy the R+D constraints cover the blue region, and those passing all of the R+D+dSph constraints are highlighted in red. Moreover, we show three benchmarks for GCE spectrum fitting in Fig. 8 and in Table V      • The "track-like" events, which are induced by muons produced in charged-current (CC) interactions of ν µ .
• The "shower-like" events, which are induced by neutral-current (NC) interactions of all neutrino flavors, and by CC interactions of ν e in all energy ranges and ν τ with E ντ ≤ 100 TeV.
• The "double-bang" events, which are generated by high energy ν τ . In this case its displaced vertices between the hadronic shower at the τ generation and the shower produced at the τ decay can reach tens of meters.
For the Yukawa structure in Eq. (55) that we will employ for illustration, only "track-like" CC and "showerlike" NC events have to be taken into account in our calculation.
In the SM, the neutrino-nucleon (νN ) interactions are mediated by the W, Z bosons: where = e, µ, τ denotes the SU (2) L lepton flavor, N = (n + p)/2 is an isoscalar nucleon, and X is the corresponding hadronic final state. At leading order (LO), the differential cross sections are [159,160] momentum transfer squared, and G F is the Fermi constant. The Bjorken variables x and y are defined as, where E ν (E ) is the energy of the incoming neutrino (outgoing lepton). The quark and anti-quark parton distribution functions (PDFs) f q , fq (f q 0 , fq0) are summed over all flavors of valence and sea quarks which are involved in CC (NC) interactions [159,160]: where ing angle. The cross sections for antineutrino-nucleon interactions (νN ) are obtained by the following replacements, The neutrino-electron interactions (in the target material) can generally be neglected compared to the neutrino-nucleon interactions due to the fact that m e M N [160]. The only important exception arises when the incoming neutrino has an energy of E ν ∼ 4 − 10 PeV. In this case, the resonance production of the W boson [161] enhances theν e e cross section significantly with the peak at E ν = M 2 W /2m e = 6.3 PeV. Since this energy is higher than most of the shower events observed at IceCube, we do not include neutrinoelectron interactions in our analysis; for a detailed discussion on this issue, see Ref. [160].
With differential cross sections in Eqs. (50) and (53), the total cross section is obtained by In Fig. 9, we present the total SM cross section as a function of the incoming neutrino energy E ν for both νN andν N interactions using the NNPDF2.3 PDF sets [162] at LO, NLO, and NNLO respectively. Due to the large uncertainty in small x grids, we have set the lower limit of x to be 10 −6 in numerical integration to reach a reliable result, which is in good agreement with Ref. [160].
Now we compute the cross section due to LQ interactions. The neutrino-nucleon CC and NC processes are mediated by an sand u-channel exchange of the LQ through Yukawa couplings in Eq. (1), and in addition there is interference between the LQ and SM amplitudes. Nevertheless we have numerically verified that both u-channel exchange and interference are negligible compared with the resonant s-channel LQ exchange. It is therefore sufficiently accurate to calculate the LQ contribution in the narrow width approximation (NWA) which only takes into account the s-channel resonance process. In order to keep at least two massive neutrinos as required by oscillation experiment, we assume a simple Yukawa structure in which only the first two generations of quarks and leptons are involved: In the NWA, the differential cross section for the NC or CC process can be written as where NC (CC) means L j = ν j ( j ), i, j, k, k = 1, 2 refer to the first two generations of quarks and leptons, and s = 2M N E ν . Neglecting the final state fermion masses, the total decay width of the LQ ψ The Bjorken scaling variable x has been integrated out in the NWA, so that the distribution functions are evaluated at x = M 2 ψ /s and Q 2 = xys = M 2 ψ y. The expressions forνN scattering can be obtained from Eq. (56) by f q ↔ fq.
For the purpose of illustration we plot in Fig. 10 the total νN cross section due to the LQ resonance for typical values of M ψ . We have assumed y 11 L , y 21 L = 1 and others vanishing, and included both NC and CC contributions. Comparing with the relatively smooth variation of the SM cross sections in Fig. 9, one finds that the LQ resonance contribution is triggered and rises rapidly once the incoming neutrino energy goes above the threshold E th ν = M 2 ψ /(2M N ). Since E th ν is in the multi TeV to PeV range in the current IceCube data, one expects that it is sensitive to the LQ in the mass range of M ψ ∼ 100 GeV − 2 TeV. With the above preparation, we move on to evaluate the event rate at the IceCube which includes the LQ contribution and perform a statistical analysis to constrain model parameters.

B. Event rate at IceCube and constraint on the model parameters
The distribution of number of events with respect to the incoming neutrino energy and the inelasticity parameter is estimated as where T is the exposure time, Ω(E ν ) is the effective solid angle of coverage, N eff (E ν ) = N A V eff (E ν ) with N A = 6.022 × 10 23 /cm 3 the water equivalent Avogadro number and V eff (E ν ) the effective target volume   of the detector, dΦ ν /dE ν the incoming neutrino flux, and dσ/dy the differential νN cross section shown In the above equation, E dep is always smaller than E ν and their relation depends on the interaction channel. In this paper, we follow the method in Ref. [160]. For NC events, the neutrino final state leads to missing energy, and the hadronic final state carries energy E X = yE ν . Thus the total EM equivalent deposited energy for ν e,µ is given by where the factor F X is the ratio of the number of photoelectrons originated from the hadronic shower to that from the equivalent-energy electromagnetic shower, which is a function of E X and parameterized as [163] where the parameters E 0 , m, f 0 are extracted from the simulations of a hadronic vertex cascade with the best-fit values E 0 = 0.399 GeV, m = 0.130, and f 0 = 0.467 [164]. On the other hand, for CC events, the leptonic final states e, µ entirely deposit their energy E e,µ = (1 − y)E ν into the EM shower. Together with the accompanying hadronic shower, the total EM equivalent deposited energy yields The remaining parameters in Eq. (58) are determined as follows: • Exposure time T = 2078 days, corresponding to the IceCube data-taking period from year 2010 to 2016 [47].
• The effective target volume V eff (E ν ) = M eff /ρ ice , where ρ ice = 0.9167 g/cm 3 is the density of ice, and M eff is the effective target mass. M eff depends on the incoming neutrino energy and reaches the maximum value 400 Mton above 100 TeV for ν e CC events (corresponding to V max eff 0.44 km 3 water equivalent), and above 1 PeV for ν µ,τ CC and NC events [49]. Here we choose V eff = 0.44 km 3 water equivalent in the calculation.
• The solid angle of coverage Ω is different for neutrino events coming from the southern hemisphere (downgoing events) and northern hemisphere (upgoing events). While for isotropic downgoing events Ω is essentially equal to 2π sr, for isotropic upgoing events Ω is generally smaller by a shadow factor S(E ν ) due to the Earth attenuation effects [159,165]. The total solid angle of coverage is then given by Ω tot (E ν ) = 2π(1 + S(E ν )) sr. In the extreme case of a fully neutrinoopaque (neutrino-transparent) Earth, one has Ω tot = 2π sr (4π sr), and for the realistic Earth one has Ω tot ∈ [2π, 4π] sr. The LQ could have a potential impact on the shadow factor through modification of the interaction length, but it has been shown in Ref. [81] that this effect is small enough to be negligible. For simplicity, we will work with the above limiting values of Ω tot in our numerical analysis, and this will yield the two edges of the upper limit band on the Yukawa couplings y ij L for a given LQ mass.
• The incoming neutrino flux dΦ ν /dE ν is assumed to be an isotropic, single power-law spectrum for each neutrino flavor i: where Φ 0 is the flux normalization at E ν = 10 5 GeV for all neutrino flavors, f i is the fraction for the ith flavor at the Earth, and γ the spectral index. Typical astrophysical processes yield source neutrinos with a flavor ratio of ν e : ν µ : ν τ = 1 : 2 : 0 when they are produced by the decay of pions. Since the distance to the source is much larger than the neutrino oscillation length, one actually observes at the Earth an oscillation-averaged flavor composition, which tends to be in a ratio of 1 : 1 : 1 [166]. We will thus use f i = 1/3 for i = e, µ, τ . For flux normalization Φ 0 and spectral index γ, we assume the best-fit values in Ref. [167]: which were obtained by performing maximum likelihood combination of different IceCube results.
In order to investigate the number of events coming from the LQ contribution and its effect at the IceCube, we use Eq. (58) to calculate all of the 14 deposited energy bins in the IceCube data points. In the left panels of Fig. 11, we present the numbers of NC and CC events due to LQ as a function of the deposited energy. The plots are done for various Yukawa components in Eq. (55) and typical LQ mass M ψ = 500, 1000 GeV, respectively. Here we simply assume a universal Yukawa coupling |y L | for the nonzero components and the legends in the figure are understood as follows: for instance, (y 11 L , y 21 L ) indicates y 11 L = y 21 L = |y L | while others vanishing. It is straightforward to extend our analysis to non-universal cases by assuming specific relations for the Yukawa components in Eq. (55). For comparison, the corresponding total numbers of events (SM + background + LQ) for the same Yukawa components and 6-year IceCube data points are presented in the right panels, where both IceCube data and SM + background fit are taken from Ref. [51]. Some important information can be observed from Fig. 11: • The resonance peak broadens and shifts according to the threshold incoming neutrino energy E th ν = M 2 ψ /(2M N ) for both NC and CC events.
• The CC events are distributed only in the deposited energy bins above the threshold energy, while the NC events are spread in all of bins. This arises from the fact that NC and CC processes deposit different amounts of energy according to Eqs. (59,61), respectively.
• The numbers of events obey the sequence N bin (y 11 L , y 21 L ) > N bin (y 11 L , y 22 L ) > N bin (y 12 L , y 22 L ), which clearly reflects the effects of PDF dependence. Since the u and d quarks are the dominant constituents of the nucleon, Yukawa components involving only the first generation of quarks give the most significant contribution while that involving the second generation of quarks is suppressed.
The interpretation of the IceCube excess in the energy interval 1 − 3 PeV generically demands a LQ mass above TeV, where the production cross section and the neutrino flux are significantly suppressed.
This may require a large Yukawa coupling beyond perturbation theory, for instance, |y L | = 3 for M ψ = 1 TeV as shown in the lower panels of Fig. 11. Nevertheless, one expects that a small fraction of the LQ contribution with a perturbative Yukawa coupling could relax the tension between the IceCube data and the SM prediction thus marginally improving the SM + background fit, which is also a part of motivation for this paper. Alternatively, one can also treat the current IceCube result as a complementary constraint, which allows us to put an upper bound on the Yukawa coupling for a given LQ mass. Along this way, we preform a binned statistical analysis with the Poisson likelihood function [82,168], where n obs, th i are respectively the observed and theory counts in the i-th bin. We then use the test statistics −2∆ ln L = −2(ln L − ln L max ) , to derive upper limits on y ij L at 90% C.L. (corresponding to −2∆ ln L = 2.71) in the LQ mass region M ψ ∈ [100,2000] GeV. Here L max is the likelihood value assuming y ij L = 0. Our results are presented in Fig. 12 for the same Yukawa structure discussed above. As expected, the most stringent bound is set on the (y 11 L , y 21 L ) components, while that on (y 12 L , y 22 L ) is relatively weak due to subdominant PDFs of the second generation of quarks in the proton.
There also exist stringent limits on y ij L from flavor physics and on M ψ from LHC direct searches. For the former, according to our discussion in section II C, the components (y 11 L , y 21 L ) and (y 11 L , y 22 L ) components are most sensitive to the K-meson decay K + → π +ν ν, while (y 12 L , y 22 L ) are sensitive to the LFV decay µ → eγ. As an illustration of the collider constraints, we use the ATLAS limits on the LQ mass at 13 TeV [105]. These limits are also shown in Fig. 12 for comparison. In all the cases, the limits derived from K + → π +ν ν and µ → eγ decays are much stronger than that from the IceCube in the entire mass range considered. This severely restricts the LQ interpretation of the IceCube excess in the 6-year data.
However, it is worthwhile to treat the excess as a supplementary constraint although it is highly limited by current statistics. With the increase of exposure time and data collection, one expects that the IceCube limit will improve and that the distribution of data in the bins may even change remarkably. In that case better agreement or more severe discrepancy with the SM prediction will serve as a complementary limit or hint of new physics. acts as an on-shell mediator. We found that the GCE spectrum is well fitted when the H 0 mass is close to the DM mass which is consistent with constraints coming from relic abundance, direct detections as well as dSphs in the Milky Way. We studied the feasibility that the resonance LQ production is responsible for the extra UHE neutrino events at the IceCube. Using the 6-year dataset in the multi TeV to PeV energy range, we derived upper limits on the LQ Yukawa couplings as a function of its mass. Although the fraction of the LQ contribution to the IceCube excess is tightly limited by flavor physics constraints at low energies, we expect that better limits will be possible with more statistics in the near future. Together with the limits from LHC direct searches, the parameter space will be explored complementarily by multi-experiments.