Phenomenology of colored radiative neutrino mass model and its implications on 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 Z2 symmetry that stabilizes the DM, and generates a tiny neutrino mass at the two-loop level with the color seesaw mechanism. After investigating the DM and flavor phenomenology of this model systematically, we further focus on its imprint on two 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 the relic density, direct detections, and dwarf spheroidal galaxies, in the Milky Way. Although the UHE neutrino events at the IceCube could be accounted for by the resonance production of a TeV-scale leptoquark, the relevant Yukawa couplings have been severely limited by the current low-energy flavor experiments. We subsequently derive the IceCube limits on the Yukawa couplings by employing its latest six-year data.


Introduction
Dark matter (DM) and tiny neutrino masses pose an outstanding challenge to both theoretical and experimental particle physicists. Although the current studies 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). Meanwhile, observations from high-energy cosmic rays (CR) may offer another angle to face the challenge. We herein 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 that generates the radiative neutrino mass, and has a cold DM particle built in. Nev-ertheless, we will first briefly review the current status of the two observations.
The GCE was first reported in Ref. [1] through analyzing the Fermi-LAT data, and the signal significance was confirmed by subsequent analyses [2][3][4][5][6][7][8]. While astrophysical interpretations such as millisecond pulsars or unresolved gamma-ray point sources [5,6,[9][10][11] are plausible, DM annihilation remains a popular interpretation because its thermally averaged cross section and morphology of density distribution match the standard WIMP scenario. In particular, Ref. [8] provides 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 pulsar-like sources located in the galactic bulge, which they referred to as the galactic bulge population, while the DM interpretation is disfavored, because its distribution is inconsistent 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 restrict their contribution only up to 4%−23% of the observed gamma-ray excess [14]. Moreover, the analyses of the spatial distribution and luminosity function of those sources were inconclusive regarding the presence of such galactic bulge population [15]. Therefore, the DM interpretation of GCE is still competitive. When using a model-independent fitting with DM directly annihilated into a pair of SM particles, the GCE spectrum is best fitted by the bb final state [7]. The other final states (τ + τ − , qq cc, gg, W + W − , ZZ, hh, and tt) with different DM masses and annihilation cross sections are also acceptable [16][17][18][19]. Additionally, when considering the 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, such as the 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 significant 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 the annihilation patterns: 1) DM annihilates directly into the SM final states, 2) DM annihilates into some intermediate particles, which subsequently cascade decay into SM particles.
While the first scenario typically suffers from stringent constraints from DM direct detections and collider searches, the second is advantageous in 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 semiannihilation 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 four-year data set released in 2015, 54 UHE neutrino events are collected (including 39 cascade events and 14 muon track events) with a 7σ excess over the expected atmospheric background [47]. Particularly, three events with an energy above PeV present excess on the SM prediction [48][49][50]. Very recently, the IceCube Collaboration has published the preliminary six-year result [51], with the total number of events increased to 82, with 28 of them observed in the recent two years. It is noteworthy that all of the new events exhibit 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 significant interest in both astrophysics and particle physics communities. While the astrophysics community focuses on various astrophysical sources [52][53][54], the particle physics community attempts 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] * , a DM particle of PeV mass is required to reproduce the desired UHE neutrino events. Such superheavy particles are extremely 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 typical belief that new physics should appear thereof. This latter scenario appears phenomenologically advantageous and could be examined using other methods, 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 renders it challenging to explain them in a single framework. We herein 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 on-shell 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 masses at two loops. In the next section, we describe the model and discuss the relevant experimental constraints on its parameter space. Sections 3 and 4 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 3, 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 subsequently presented incorporating all these limits. In section 4.1, we calculate the SM and LQ contributions to the neutrino-nucleon scattering cross section. Subsequently, in section 4.2, 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 present our conclusion in section 5.

Model and relevant constraints
2.1 The model The particle contents and their charge assignments are shown in Table 1. In addition to the LQ ψ and diquark ω, we further introduce two singlet scalars, ϕ with Table 1. Particle contents and their charge assignments. The double vertical line separates the SM particles from the new ones. lepton number L=2 and S with L= 1 2 . Herein, ϕ is used to break the U (1) B−L gauge symmetry spontaneously; thus, generating the L-breaking trilinear term ψ * ψ * ω required for the radiative neutrino masses. Notably, owing 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 an ad hoc discrete symmetry [84][85][86]. Such that U (1) B−L is 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 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. Further, 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. It is noteworthy that owing 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 assumed to be positive, and the trace is over the color indices. Hence, 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. Owing to the B−L charge assignment of S, one can still obtain S = 0 after a spontaneous symmetry breaking, such that a residual Z 2 symmetry remains under which only S is odd. This blocks all potential decays of S, rendering 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 where g BL is the gauge coupling of U (1) B−L . The LEP bound requires that [95] yielding a lower limit on v ϕ 3.5 TeV. Meanwhile, the direct searches for the Z ′ -boson at the 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, such that v ϕ =20 TeV in our following discussion. The masses of the DM S, LQ ψ, and diquark ω can be obtained 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 regions, respectively. The constraints from the relic density and direct detections will be discussed in section 3. Assuming the LQ ψ decaying exclusively into eq, µq, and τ q, CMS (ATLAS) has excluded M ψ < 1010, 1165, 850 GeV [102][103][104] (M ψ < 1100, 1050, 534 GeV [105][106][107][108]). However, both ψ → ℓq and ψ → ν ℓ q ′ exist in our model. The maximum exclusion limits by CMS (ATLAS) for the first and second generation LQs are 850, 960 GeV [102,103] (900, 830 GeV [105]) when assuming BR(ψ → ℓq) = 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 primarily consider M ψ 1 TeV and M ω =7 TeV to adhere to 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,H 0 ) 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 the 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 the 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 θ, and the VEVs v φ,ϕ :

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 downtype quarks are much lighter than the colored scalars, it can be simplified for the order of magnitude estimate to where Typically, a neutrino mass m ν ∼ 0.01 eV can be realized 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 the trilinear coupling λv ϕ ψ * ψ * ω, and the choice of λ ∼ 0.1; further, v ϕ = 20 TeV satisfies the perturbativity requirement λv ϕ 5min(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 the proper parameterization [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, such 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 will be discussed below, we follow the typical phenomenological practice to use Yukawa components y ij L as the input parameters whose values will be constrained by the IceCube data and lowenergy experiments.

Flavor constraints
The LQ ψ can induce various flavor-violating processes at the tree level. To minimize such processes, one typically assumes y R =0 [94,127], because the y R term is less important to the neutrino masses. This also fits our interest in the IceCube UHE neutrino events that may be induced by y L but not y R couplings. As the 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: The relevant upper limits on ǫ ijkn in the colored Zee-Babu model are summarized in Table 3 of Ref. [94]. In particular, two ǫ ijkn 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] One method to satisfy these bounds is to assume, e.g., y ℓu L 0.001 and y ℓc L 0.1 at M ψ ∼ 1 TeV. The constraints on other components of ǫ ijkn are loose, and can be readily avoided by, e.g., y ℓq L O(0.1) for a TeV scale M ψ [130].
A by-product of lepton radiative decays is the LQ contribution to the anomalous magnetic moment of the charged lepton ℓ [134,135] Under the constraints from LFV, the predicted values are ∆a e = −2×10 −19 , ∆a µ = −1×10 −14 , and ∆a τ = −2×10 −12 for the 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, because 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 ecm, and 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] 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.

DM phenomenology and GCE spectrum fitting
To investigate the DM phenomenology, we use FeynRules [141] to generate the CalcHEP [142] model file and implement it into the micrOMEGAs4.3.2 package [143] to calculate the DM relic abundance and DMnucleon scattering cross section. We performed random scans for the parameter space in both the low and high mass DM scenarios (with 3×10 5 samples for each), with the input parameters shown in Table 2. The constraints from the DM relic abundance and direct detection experiments were imposed on each sample. For the DM relic abundance, we adopted the combined Planck+WP+highL+BAO result in the 2σ range, 0.1153 < Ω DM h 2 < 0.1221 [144]. For direct detections, we used the latest spin-independent limits obtained by LUX [145], XENON1T [146] and PandaXII [147] Collaborations.
To illustrate the effects of various annihilation processes on relic abundance and direct detection, we list all important annihilation channels in  Table 3. Several features learned from these results are summarized as follows 1) .
For the low-mass DM scenario: 1) There are much less survived samples than for the high mass scenario. This is because the coupling between the DM and SM Higgs, λ Sh , is tightly constrained     by relic abundance and direct detections. Consequently, only the channels mediated by the singlet scalar H 0 and t-channel DM can survive. On the contrary, the annihilation channels SS * → W + W − /H 0 H 0 are available in a wide parameter region. Specifically, regions of λ Sh 0.03 and λ SH 0 0.01 are favored for the low-mass DM scenario.
2) The bb and W + W − channels are respectively dominant when M S 75 GeV and 75 GeV. This is the typical behavior of the Higgs (h/H 0 ) portal DM [148]. In addition, although the H 0 H 0 channel could satisfy the relic abundance requirement in broad DM mass regions, only samples with M S >100 GeV could escape the direct detection bounds.
For the high-mass DM scenario: 1) 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. This is because we have chosen the corresponding couplings λ Sh ,λ Sψ <0.5 in our scan.
2) λ Sh (λ Sψ ) 0.2 is required when the W + W − (ψψ * ) channel dominates. Moreover, the W + W − channel fills a narrow band in the λ Sh −M S plane where λ Sh increases with the increase in M S , while the ψψ * channel in the same plane is much scattered.
We now discuss the GCE spectrum fitting in our model. The hard photons due to DM annihilation arise primarily from the subsequent decays of the SM particles, because their direct production is typically loop suppressed. The continuous gamma-ray spectrum results from the light mesons produced through the 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 annihilation cross section in the galactic halo, and dN γ f /dE γ is 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 ⊙ scosψ. 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ψ=cosbcosl.
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 use the generalized Navarro-Frenk-White (gNFW) profile for the DM halo distribution [149] 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 that yields the valueJ ben forJ. The uncertainties of (ρ ⊙ , γ) subsequently translate intoJ ≡ JJ ben , where the factor J ∈[0.14,4.4] parameterizes the allowed range for the DM distribution. We performed the GCE scan for J in the range above, and J = 1 for the benchmark profile.
To fit the GCE, we used the results in Ref. [8], which explored in detail the multiple GDE models. We employ micrOMEGAs and PPPC4DMID [150] to generate the photon spectrum and perform global fitting using where dΦ th,obs i /dE γ are the theoretical and observed gamma-ray flux in the i-th energy bin, respectively. Σ ij is the covariance matrix provided by Ref. [8] that includes both statistical and correlated systematic errors. Herein, we focus on the on-shell mediator scenario, in which the DM annihilates into a pair of on-shell singlet scalars H 0 , which in turn decays into the SM quarks and leptons. The decay branching ratios of H 0 are presented in Fig. 6 versus its mass; they exhibit a similar pattern to those of the SM Higgs owing to the φ 0 −ϕ 0 mixing. We vary M S , M H 0 in the GCE scan while setting other parameters as shown in Table 4. In addition to relic abundance and direct detections, one must consider 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 imposes a stringent constraint on the DM interpretation of GCE for various annihilation channels. Herein, we adopt the dSphs limits provided in Ref. [151], which performed a modelindependent and comprehensive analysis on various twobody and four-body annihilation channels based on the Planck [21] (CMB), Fermi-LAT [152][153][154][155][156] (dSphs) and AMS-02 [157] (antiproton) results. For our model, The most relevant are the 4b, 4τ , and 2b2τ channels. During the scan, we translated the corresponding limits into each M H 0 sample weighted by Br(H 0 → bb/τ + τ − ) and subsequently extracted the strictest one. We present our results in Fig. 7, where the allowed parameter regions for fitting the GCE spectrum and fulfilling various constraints are displayed in the M S −M H 0 (left panel) and M S − σv halo (right) plane. 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 5. Among them, benchmark1 (benchmark2) is the best fit point of the GCE spectrum for J ∈ [0.14, 4.4] (J = 1) in the total samples, while benchmark3 is the best fit point in the R+D+dSph samples. Except for benchmark1, the other two nearly degenerated H 0 and S with M H 0 ≈M S ∈ [40,50] GeV. This feature can be understood by a simple analysis of kinematics. For nearly degenerated H 0 and S, the H 0 pair is produced almost at rest, and each decay final state of H 0 carries an energy M H 0 /2 ≈ M S /2, thereby producing a spectrum similar to the two-body annihilation process with a doubled Table 4. The ranges or values of the input parameters used in GCE scan. All masses are in units of GeV.  Table 5. Three benchmarks for GCE spectrum fitting. Here the benchmark1 (benchmark2) is the best fit point of GCE spectrum in the total samples for factor J ∈ [0.14, 4.4] (J = 1), and benchmark3 is the best fit point in the R+D+dSph samples.      Table 5. The GCE data with statistical and systematic errors (cyan ) in Ref. [8].
number of injection fermions that reproduces the best fit result as the two-body bb final state. Finally, the exception of benchmark1 can be understood because it only occasionally produces a minimal χ 2 by using a marginal value of J , and will yield an unacceptably large χ 2 when J =J ben .

Neutrino-nucleon scattering in SM and LQ contribution
The IceCube neutrino observatory is located at the South Pole. The overwhelming majority events recorded by IceCube are muons from CR air showers, and only approximately one in a million events results from neutrino interactions. In the latter case, the UHE neutrinos in CR penetrate the ice and scatter with nucleons through neutrino-nucleon deep inelastic scattering (DIS) inter-actions. The Cherenkov light emitted by the secondary particles produced during scattering is observed by the IceCube detector. Depending on the interaction channel and the incoming neutrino flavor, three types of signatures can be distinguished for the neutrino events [158]: 1) The "track-like" events that are induced by muons produced in the charged-current (CC) interactions of ν µ .
2) The "shower-like" events that 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.
3) The "double-bang" events that 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 the "track-like" CC and "shower-like" NC events must be considered 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] In the equations above, M N and M W, Z are the nucleon and W, Z boson masses, respectively; −Q 2 is the 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 antiquark parton distribution functions (PDFs) f q ,fq (f q 0 ,fq0) are added over all flavors of valence and sea quarks that are involved in CC (NC) interactions [159,160]: and L u = 1+R u with θ W the weak mixing 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 because 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. Because this energy is higher than most of the shower events observed at IceCube, we do not include neutrino-electron interactions in our analysis; for a detailed discussion on this issue, see Ref. [160].
With the 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. Owing to the large uncertainty in small x grids, we have set the lower limit of x to be 10 −6 in the numerical integration to yield a reliable result, which is in good agreement with Ref. [160]. We subsequently compute the cross section owing to LQ interactions. The neutrino-nucleon CC and NC processes are mediated by an s-and u-channel exchange of the LQ through Yukawa couplings in Eq. (1). In addition, interference occurs between the LQ and SM amplitudes. Nevertheless, we have numerically verified that both the 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) that only considers the s-channel resonance pro-   cess. To retain at least two massive neutrinos as required by the 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 refers 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 ψ is Γ ψ ≃M ψ /(8π) ij |y ij L | 2 . The Bjorken scaling variable x has been integrated out in the NWA, such 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.
To illustrate, 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 is above the threshold E th ν =M 2 ψ /(2M N ). Because 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 preparation above, we proceed to evaluate the event rate at the IceCube that includes the LQ contribution, and perform statistical analysis to constrain the model parameters.

Event rate at IceCube and constraint on the model parameters
The distribution of the number of events with respect to the incoming neutrino energy, and the inelasticity parameter is estimated as dN 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, 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 Eq. (56) In the equation above, E dep is always smaller than E ν and their relation depends on the interaction channel. Herein, we follow the method in Ref. [160]. For the 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 originating 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]. Meanwhile, for the 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: 1) Exposure time T = 2078 days, corresponding to the IceCube data acquisition period from year 2010 to 2016 [47].
2) 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]. Herein, we choose V eff =0.44 km 3 water equivalent in the calculation.
3) 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 ν ) owing to the Earth attenuation effects [159,165]. The total solid angle of coverage is subsequently given by Ω tot (E ν ) = 2π(1+S(E ν )) sr. In the extreme case of a fully neutrino-opaque (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 the modification of the interaction length; however, it has been shown in Ref. [81] that this effect is sufficiently small to be negligible. For simplicity, we will use the limiting values of Ω tot above 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. 4) 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. Because the distance to the source is much larger than the neutrino oscillation length, an oscillation-averaged flavor composition, which tends to be in a ratio of 1 : 1 : 1 [166]is observed on the Earth. Thus, we use f i = 1/3 for i = e,µ,τ . For flux normalization Φ 0 and spectral index γ, we assume the best-fit values as in Ref. [167]: These were obtained by performing maximum likelihood combination of different IceCube results.
To investigate the number of events from the LQ contribution and its effect at the IceCube, we used 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 were generated for various Yukawa components in Eq. (55), and typical LQ mass M ψ =500, 1000 GeV, separately. 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 ) indi-  cates y 11 L = y 21 L = |y L | while others are 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 six-year IceCube data points are presented in the right panels, where both the IceCube data and SM + background fit are extracted from Ref. [51]. Some important information can be observed from Fig. 11: 1) 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.
2) The CC events are distributed only in the deposited energy bins above the threshold energy, while the NC events are spread in all bins. This arises from the fact that NC and CC processes deposit different amounts of energy according to Eqs. (59,61), respectively.
3) 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. As the u and d quarks are the dominant constituents of the nucleon, the Yukawa components involving only the first generation of quarks present 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 an LQ mass above TeV, where the production cross section and the neutrino flux are significantly suppressed. This may require a large Yukawa coupling beyond the perturbation theory, e.g., |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 loosen the tension between the IceCube data and the SM prediction, thus marginally improving the SM + background fit, which is part of the motivation for this paper. Alternatively, one can also treat the current IceCube result as a complementary constraint that allows an upper bound to be placed on the Yukawa coupling for a given LQ mass. Further, we performed a binned statistical analysis with the Poisson likelihood function [82,168], where n obs, th i are the observed and theoretical counts in the i-th bin, respectively. We subsequently use the test statistics −2∆lnL=−2(lnL−lnL max ), to derive the upper limits on y ij L at 90% C.L. (corresponding to −2∆lnL = 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 owing to the subdominant PDFs of the second generation of quarks in the proton. Stringent limits also exist on y ij L from flavor physics, and on M ψ from the LHC direct searches. For the former, according to our discussion in section 2.3, the components (y 11 L , y 21 L ) and (y 11 L , y 22 L ) are the 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 cases, the limits derived from K + → π +ν ν and µ → eγ decays are much stronger than those from the IceCube in the entire mass range considered. This severely restricts the LQ interpretation of the IceCube excess in the six-year data. However, it is worthwhile to treat the excess as a supplementary constraint although it is highly limited by the current statistics. With the increase in 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 significantly. In that case, better agreement or more severe discrepancy with the SM prediction will serve as a complementary limit or hint of new physics.

Conclusion
We have investigated the phenomenology of the colored Zee-Babu model augmented with a U (1) B−L gauge symmetry and a singlet scalar DM S. The tiny neutrino masses were still generated via a two-loop radiative seesaw involving the SM quarks, a diquark, and an LQ; however, we have related to two high-energy CR observations: the Fermi-LAT GCE and the PeV UHE neutrino events at the IceCube. For the Fermi-LAT GCE, we focused on the annihilation channel in which the singlet (-dominating) Higgs H 0 acted as an on-shell mediator. We found that the GCE spectrum is well fitted when the H 0 mass was close to the DM mass, which is consistent with the constraints coming from relic abundance, direct detections, and dSphs in the Milky Way. We studied the feasibility of the resonance LQ production being responsible for the extra UHE neutrino events at the IceCube. Using the six-year dataset in the multi TeV to PeV energy range, we derived the 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 the 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 the LHC direct searches, the parameter space will be explored complementarily by multiple experiments.