Implications of a Electroweak Triplet Scalar Leptoquark on the Ultra-High Energy Neutrino Events at IceCube

We study the production of scalar leptoquarks at IceCube, in particular, a particle transforming as a triplet under the weak interaction. The existence of electroweak-triplet scalars is highly motivated by models of grand unification and also within radiative seesaw models for neutrino mass generation. In our framework, we extend the Standard Model by a single colored electroweak-triplet scalar leptoquark and analyze its implications on the excess of ultra-high energy neutrino events observed by the IceCube collaboration. We consider only couplings between the leptoquark to first generation leptons and quarks and carry out a statistical analysis to determine the parameters that best describe the IceCube data as well as set $95\%$ CL upper bounds. We analyze whether this study is still consistent with most up-to-date LHC data and various low energy observables.


Introduction
In this work we study the implications of a colored electroweak-triplet scalar leptoquark (LQ) on the ultra high energy (UHE) neutrino spectrum observed at IceCube, focusing particularly on the range above PeV, where a bit higher than expected event rate has been reported [1]. The potential of the IceCube facility to probe LQ models has been exploited in many works. In ref. [2], for example, the inelasticity distribution of the events detected at IceCube are used to test LQ production; in refs. [3,4] electroweak-singlet scalar LQs, with different flavor structure for its couplings, are introduced to fit the neutrino flux at the PeV range. In this regard, besides the many explanations that incorporate new physics effects, other possibilities within the picture of the Standard Model (SM) have also been proposed [5,6].
Leptoquarks (LQs) are fields that arise naturally from the unification of quarks and leptons in extensions of the SM [7][8][9]. In particular, unification of quarks and leptons into simple groups of SU (5) requires the unification of LQs with the SM-like Higgs boson. However, one main obstacle that arises from the introduction of LQs is how they can mediate proton decay at tree level, specially in the case of LQs that violate lepton and baryon numbers, if those quantum numbers are indeed assigned. Unification schemes to accommodate very heavy LQs to avoid proton decay bounds have also been studied, in particular a scheme based on a flipped SU(5) framework where SM fields are embedded into representations of a SU(5)×U(1) gauge group has proved successful [10][11][12][13]. However, in view of the current experimental effort to produce particles beyond the SM, most studies have focused on two particular scalar LQ representations out of the six possible ones [14], where phenomenologically light LQs are natural. These fields transform as (3, 2, 1/6) and (3, 2, 7/6) under the SM SU(3) c ×SU(2) W ×U(1) Y gauge group and have been implemented to address several hints of new physics beyond the SM, in particular the excess reported by IceCube [15] and the anomalous LHC same-sign lepton events [16] such as [17,18]. These two weak doublets do not couple to baryon number violating operators at tree level. However, effects of higher dimensional operators can cause baryon number violation. In this regard, the authors in [19] discuss a framework where one can naturally suppress these operators. Despite the fact that the two representations mentioned above are the most frequently used, other LQ models with diquark operators have been also considered to address other reported anomalies. One recent work, for example, uses a LQ with the quantum numbers (3, 1, −1/3) to address deviations on R D * , R K , and the (g − 2) of the muon [20] and similarly using the electroweak doublets introduced above [21][22][23][24]. Another work uses this electroweak singlet LQ to explain the excess of high energy neutrino events [3]. Only two other scalar LQ representations can couple to SM neutrinos and quarks and are thus relevant in the explanation of the IceCube excess, these are the (3, 1, −1/3) mentioned earlier and a weak triplet (3, 3, −1/3). In contrast to the former, the latter has not been probed through the UHE neutrino spectrum observed at Icecube. Both LQs couple to diquarks and can induce proton decay at tree level. The authors in [25] discuss a scenario to suppress the diquark operators by embedding the weak triplet and singlet into a 45 H -dimensional Higgs representation of a SU(5) GUT model. It is therefore plausible to consider light weak triplet and singlet LQs, with masses accessible at colliders, as a possible source of the UHE neutrino events observed at IceCube. The study of LQs has become very active; with a focus also on R-parity violating scenarios of supersymmetry (SUSY) which yield couplings of scalar superpartners to quarks and leptons. As far as the UHE neutrino events observed at IceCube is concerned, these can be used to constrain R-parity violating supersymmetric models [26,27]. LQ have a rich phenomenology and for this reason we direct the reader to a recent review [28] for more an in depth discussion and references therein.
In this work we focus on the weak triplet since this class of particles has recently been used to mediate the generation of neutrinos masses radiatively and at three loops [29,30]. The model considered by these authors also includes a heavy Majorana neutrino dark matter candidate. The work focuses solely on the phenomenology of a LQ coupling right-handed up-type quarks to the Majorana neutrino, yielding a mechanism for its relic abundance. In addition, a monotop search strategy was introduced and limits were placed on the model using current LHC data. In this work, we wish to go one step further, that is, analyze the phenomenology of the weak triplet, originally with masses set at the TeV scale, and introduce a mechanism to produce high energy neutrino events in detectors such as IceCube. Our model is very attractive since the coupling of the LQ to up-and down-type quarks is the same. It also allows us to directly connect the observations by the IceCube collaboration to the mechanism of neutrino mass generation and specific GUT scenarios where tree-level baryon number violating operators are absent.
The remainder of this paper is organized as follows. In sec. 2 we review the model specifying those aspects related to the UHE neutrino events at IceCube. In sec. 3 we probe the proposed weak triplet with the IceCube detector. We review the SM neutrino-nucleon scattering cross section and compute the respective LQ contribution in sec. 3.1. We then obtain in sec. 3.2 the new physics contribution to the rate of events expected at IceCube and study its behaviour with respect to the LQ masses. In sec. 3.3 we perform an statistical analysis in order to determine the parameters that best fit the IceCube data as well as set upper limits as a function of the LQ mass. Secs. 3 and 4 are dedicated to the analisis of constraints arising from the LHC experiments with the 8 and 13 TeV data sets, lepton flavor violation (LFV), and low energy precision measurements such as atomic parity violation. Finally, in sec. 6, we compare the results obtained from the analysis of the IceCube data with the constraints derived in secs. 4 and 5 and provide some concluding remarks. The Appendix gives some further details about the attenuation effects on upward-going neutrinos resulting from their passage through the Earth.

Model
The authors in [29,30] investigated a model that incorporates LQs, one triplet under the SU(2) W gauge interaction and a singlet. In addition, the model contains a single Majorana right-handed neutrino used to both explain the nature of dark matter and the mechanism for neutrino mass generation. Within this framework Majorana masses for the active neutrinos were made possible via a radiative process which involves a three loop diagram. In order for the mechanism to work, two representations of weak triplets were implemented: a lepton and baryon number violating LQ transforming as a (3, 3, −1/3) under the SM gauge group and a weak triplet transforming as a (3, 3, 2/3) with no treelevel coupling to fermions. In this work we are primarily interested in the former because it couples quarks to leptons and therefore affects the neutrino-nucleon cross section, which may lead to new features in the spectrum of UHE neutrinos observed by the IceCube collaboration. In the following we will refer to the field transforming as a (3, 3, −1/3) under the SU(2) W × U(1) Y SM interactions as χ. Given its quantum numbers, one may choose to write χ as a 2 × 2 matrix with the following transformation property where U = exp(iω j τ j /2) and τ j is the j-th Pauli matrix. We then represent the weak triplet χ with the following matrix: and parametrize its interactions with left-handed quarks and leptons with the following Lagrangian: where P R = (1 + γ 5 )/2, λ i j represents the coupling between the i-th generation of quarks and the j-th generation of leptons with i, j = 1, 2, 3, and ψ c denotes the conjugate field of ψ.
The terms in the above Lagrangian are not the only ones allowed by gauge invariance. One can incorporate a quark bilinear operator coupling to the weak triplet given by where the indices i and j run over the three quark generations and Q L denotes the quark weak doublet. The interaction in eq. (2.4) induces rapid proton decay and a symmetry needs to be imposed in order to suppress the strength of these interactions. However, as shown in [31], the above operator also induces a Planck scale suppressed dimension five operator that gives the decay modes p → π + ν, p → K + ν and p → K + π + l − . In order to generate two-body nucleon decay partial rates near the present limits one would require m χ ∼ (3k) 1/4 (y · Y 5 ) 1/2 10 7 GeV, (2.5) where Y 5 denotes the coefficient of the dimension five operator and k is in the range 0.17 ≤ k ≤ 6.7. With this in mind, one can obtain LQ masses within the reach of particle colliders with couplings of order 10 −5 to 10 −3 . Of course allowing for the above diquark operator will make χ not a genuine LQ in the sense that not only the operators in eq. (2.3) are present. However, the diquark operators can be suppressed or neglected by imposing a symmetry, in particular a GUT symmetry in a supersymmetric framework. This case has been discussed in [25] where one embeds χ in a 45 H -dimensional Higgs representation and has different contractions leading to the quark-lepton interaction and the diquark interaction. Allowing only for the lepton-quark contraction will also lead to the absence of any mixing induced proton decay. In what follows, we will assume that χ is a genuine LQ in the sense that either the diquark operator is suppressed or it is altogether absent. In this case, one can assign a lepton and baryon number to χ such that the accidental lepton and baryon number symmetries of the SM are conserved.
Another aspect relevant in the study of the impact of our model in the UHE neutrino spectrum observed in the IceCube detector is the flavor structure of the interactions in eq. (2.3), and its consistency with measurements looking for deviations from the minimal flavor structure of the SM, specially the 2.6σ deviation from lepton universality presented by the LHCb collaboration on the measurement of R K [32], the ratio between the branching fractions of B → Kµµ and B → Kee. In addition, there are hints of LFV reported by the CMS collaboration on the decay h → µe [33]. The authors in [34] have used the weak triplet introduced in this work to explain these two measurements by adapting frameworks with non-abelian flavor symmetries that predict the leptonic mixing matrices. Even though the simplest scenarios are those for which the LQ couples to a single generation of leptons, the authors use a data-driven approach to constrain the LQ Yukawa couplings in a generalized scenario using a hierarchical pattern consistent with the observed quark mass pattern. They then analyze various flavor models that lead to different textures of the LQ Yukawa matrix consistent with LFV decays, rare meson decays and lepton universality.
In what follows, we will assume that the LQ couples primarily the first family of quarks to the electron and the muon and the correponding neutrinos. We will also consider the couplings of the LQ to the second and third families of quarks to be suppressed in order to make the collider phenomenology more tractable and, at the same time, to simplify the computation of the rate of UHE neutrino events arising from the LQ component. The study of the viable parameter space consistent with collider constraints and low energy measurements such as LFV decays and atomic parity violation is postponed to secs. 4 and 5.

IceCube and PeV Neutrinos
In this section, we study the impact of the model of LQs proposed in sec. 2 in the spectrum of PeV neutrinos measured at IceCube. In the first place, we revisit the computation of the neutrino-nucleon scattering cross section within the SM, and then we derive the corresponding LQ contribution. The addition of this new physics component leads to particular features in the spectrum, which is studied by computing the expected rate of events. Finally, by adding the rate of events expected from the LQ contribution on top of the SM we determine the masses and couplings that best accommodate the observed spectrum.

Neutrino−nucleon scattering cross section
At the IceCube detector, the ultra-high energy (UHE) neutrinos coming from outside the atmosphere are detected by observing the Cherenkov light emitted by the secondary particles produced in their interactions with the nucleons present in the ice. In the standard model (SM), there are charged current (CC) as well as neutral current (NC) neutrino-nucleon interactions, which are mediated by a W or a Z boson, respectively. The topology of the events observed at IceCube depends on the interaction channel as well as on the flavor of the incoming neutrino. The track-like events are induced by CC ν µ interactions, while the shower-like events are induced by CC ν e and ν τ interactions and NC interactions of neutrinos of all flavors.
The SM differential cross section for the generic CC interaction ν N → X, with = e, µ, τ , N the target nucleon and X the hadronic final state, can be written as, where M W and M N are the masses of the W and the nucleon respectively, −Q 2 is the invariant momentum transferred by the intermediate boson to the hadronic system, and G F is the Fermi constant. The Bjorken scaling variable x and the inelasticity y used in eq. (3.1) are defined as where E ν and E are the energies carried by the incoming neutrino and by the outgoing lepton in the laboratory frame respectively. Finally, in the case of an isoscalar nucleon N ≡ (n + p)/2 1 , the quark distribution functions in the differential cross section are given by [35], where u, d, s, c, t, b denote the distributions corresponding to the various quark flavors in a proton, and the subscripts v and s indicate the valence and sea contributions.
Similarly to the CC case, we can write the differential cross section corresponding to the NC process ν + N → ν + X in terms of the variables x and y, with the following quark distribution functions, with the quiral couplings given by . Likewise, for the NCνN differential cross section, the corresponding expression is obtained from eq. (3.5) with the replacement q 0 ↔q 0 . In addition to their interactions with nucleons, the UHE neutrinos can also interact with electrons in the detection volume. These interactions are proportional to the electron mass and then can be generally neglected compared to the neutrino-nucleon interactions. The only exception is the resonant production of W − inν e e interactions, which occurs at 6.3 PeV. Since this energy is high compared to the most energetic showers observed at IceCube, we will not enter in details regarding the neutrinoelectron interactions. The expressions for the differential cross section for these interactions can be found for example in [35].
For the numerical computations performed in this paper, we have used the NNPDF2.3 PDF sets [36]. In particular, we use the central values of the PDF sets with α s (M Z ) = 0.118 at NLO. The NNPDF2.3 sets provide a grid division that can go up to Q 2 max = 10 8 GeV 2 in the Q 2 axis, and down to x min = 10 −9 in the x axis. However, given the large uncertainties in the grids for low x, we have taken in most of the computations 10 −6 as the lower limit for the x-integration. For illustration, we show in figure 1 the total νN andνN cross sections in terms of the incoming neutrino energy E ν for the SM contributions.
In order to study the impact of the proposed LQ in the energy distribution of the events expected at IceCube, we must compute the LQ contribution to the neutrino-nucleon cross section. From eq. (2.3), we see that only χ 1 and χ 2 give contributions to this cross section; the analogue to the SM's NC processes are provided by both χ 1 and χ 2 , whereas the final states corresponding to the SM's CC processes are produced only through χ 2 . The corresponding Feynman diagrams are depicted in figure 2, where U and D denote up-and down-type quarks, and the indices i, i and j, j indicate the number of family for quarks and leptons respectively.  Figure 2: Feynman diagrams contributing to the νN interaction. The two first rows correspond to NC processes and the last one to CC processes.
Since the cross sections corresponding to the s-channel processes shown in the first column of figure 2 are resonance enhanced, we will assume that they dominate the neutrino-nucleon interaction and use the narrow width approximation for both χ 1 and χ 2 2 . Also, we consider a scenario in which the LQ triplet couples only to the first generation of quarks and, for the sake of simplicity, to the first and second generations of leptons; hence, we have λ i j = 0 for i = 1 and/or j = 3. Within this scenario, the differential cross section for the NC and CC processes can be written as follows where j, j = 1, 2, Γ χ 1 and Γ χ 2 are the total widths of χ 1,2 , s = 2M N E ν is the center-of-mass energy squared, and f u,d are the distribution functions of the up and down quarks in the nucleon, respectively. In the case of an isoscalar nucleon, these functions turns out to be equal and given by, In eqs. (3.8) and (3.9), the fractional momentum x has been integrated out by using the narrow width approximation for both χ 1 and χ 2 . As a consequence of this, the distribution functions are evaluated at x = m 2 χ 1,2 /s and Q 2 = xys = m 2 χ 1,2 y. In order to compute the widths of χ 1 and χ 2 we assume that Γ χ 1 is saturated by the decay χ 1 → ν e,µ u while Γ χ 2 is saturated by χ 2 → ν e,µ d and χ 2 → (= e, µ) u. Thus, these widths are written in terms of the couplings as follows By combining eqs. (3.11) and (3.12) with eqs. (3.8) and (3.9), the differential cross sections for NC and CC νN scattering are expressed in terms of the couplings λ 1 1 and λ 1 2 and the LQ masses m χ 1 and m χ 2 . In the case of antineutrino-nucleon scattering the expressions for the NC and CC processes are the same but with the quark distributions replaced by the respective antiquark distributions, i.e., f u → fū and f d → fd. In figure 3 we show the LQ contribution to the total νN cross section for the NC and CC processes. For concreteness, we have considered the case in which χ 1 and χ 2 are degenerate in mass with m χ 1 = m χ 2 = 800 GeV and the LQ couplings are such that |λ 1 1 | = |λ 1 2 | = 1. From figure 3 and comparing with the SM cross sections in figure 1, we see that the LQ contribution turns on when the incoming neutrino energy is enough to produce the resonances χ 1,2 . This occurs when the center-of-mass energy is such that √ s m χ 1,2 or, equivalently, when E ν m 2 [GeV] In contrast to the SM cross section, the total cross section induced by the resonant LQ production is higher for the NC reactions, which involve both χ 1 and χ 2 , than for the CC reactions, which proceed only via χ 2 . In the case of non-degenerate masses, we note that the splitting in mass cannot be greater than ∼ 50 GeV due to the constraints arising from the oblique parameters [28,37]. For such small difference in mass between χ 1 and χ 2 , the behaviour of the cross section with respect to the incoming neutrino energy is quite similar, with the actual values being slightly smaller or higher depending on wether the mass of χ 1 or χ 2 is increased or decreased with respect to the degenerate case. In what follows, we will focus then in the case m χ 1 = m χ 2 .

Event rate at IceCube and LQ contribution
The LQ contribution to the total number of events induced by neutrinos of a certain flavor can be written as follows, where T is the exposure time, Ω the solid angle of coverage, N eff the effective number of target nucleons, dφ/dE ν is the flux of the incoming neutrinos and dσ/dy is the differential neutrino-nucleon cross section corresponding to the sum of eqs. (3.8) and (3.9). From eq. (3.13), the distribution of the number of events with respect to the incoming neutrino energy and the inelasticity is (3.14) We note that, in order to compare with the rate of events observed at IceCube, we must use the distribution of the number of events with respect to the deposited energy E, which is always smaller than the incoming neutrino energy E ν . The predicted number of events due to the LQ contribution in the deposited energy interval ∆ ≡ (E i , E f ) is given by (3.15) By changing variables from the deposited energy E to the incoming neutrino energy E ν in eq. (3.15), we obtain The relation between the deposited energy and the incoming neutrino energy depends on the interaction channel. In this study we follow the approach used in ref. [5], which we summarize in the following. For NC events, the outgoing hadrons carry an energy E X = yE ν , and the corresponding deposited energy is given by E had = F X yE ν , where F X is the ratio of the number of photo-electrons yielded by the hadronic shower to that produced by an equivalent-energy electromagnetic shower. This quantity is parameterized as [5,38] where E 0 = 0.399 GeV, m = 0.130 and f 0 = 0.467 are the best-fit values obtained from simulations in ref. [38]. The energy carried by the final state neutrino is missed and thus the total deposited energy for NC ν e -and ν µ -events is E NC = F X yE ν . In the case of CC events, in contrast, the energy of the final state lepton, E e,µ = (1 − y)E ν , is completely deposited giving rise to a total deposited energy given by E CC = E e,µ + E had . The remaining ingredients appearing in eq. (3.16) are set in the following manner: • For the time of exposure we take T = 1347 days, corresponding to four years of IceCube data between 2010 and 2014 [1].
• The solid angle of coverage is Ω = 2π sr for events coming from the southern hemisphere (downward-going neutrino events). Due to attenuation effects in the Earth, the effective solid angle for northern events turns out to be smaller by a shadow factor that depends on E ν . This factor can be written as [35] S(E ν ) = 1 2π where the function z(θ) gives the thickness of the Earth as a function of the angle of incidence of the incoming neutrinos and L int (E ν ) is the interaction length, which depends on the flavor of the incoming neutrino. Thus, for an isotropic neutrino flux, the total solid angle of coverage is given by Ω tot = 2π(1+S(E ν )) sr, from which we see that for a fully opaque Earth Ω tot = 2π sr, while for a transparent Earth, Ω tot = 4π sr. The LQ contribution modify in principle the interaction length and therefore the shadow factor. However, the deviation from the SM expectation for the total solid angle turns out to be small, so that the LQ contribution can be neglected in the computation of the shadow factor. On the other hand, S(E ν ) is a monotonically decreasing function of the incoming energy E ν that, in the range of energies relevant at Icecube (10 TeV − 10 4 TeV), varies between 1 and ∼ 0.15. For simplicity, we will cosider the total solid angle of coverage as a constant and present in the following the results for both limiting cases mentioned above (for further details on attenuation effects in the Earth see App. A).
• The effective number of target nucleons depends on the energy of the incident neutrinos, N eff = N A V eff (E ν ), where N A = 6.022 × 10 23 cm −3 water equivalent (we) is Avogadro's number. The effective target volume can be written as V eff (E ν ) = M eff /ρ ice with M eff the effective target mass and ρ ice the density of ice. The effective target mass increases with E ν and reaches a maximum value of 400 Mton above 100 TeV, in the case of ν e CC events, and above 1 PeV for NC events and CC events induced by ν µ and ν τ [39]. For the computation of the LQ contribution to the event rate observed at IceCube, we use the maximum value for M eff which corresponds to V eff = 0.44 km 3 we.
• For each neutrino flavor i, we assume an isotropic, single power-law flux that is parameterized as follows, where f i is the fraction of neutrinos of the i-th flavor at Earth, γ is the power law spectral index and φ 0 is the all-flavor neutrino flux at 100 TeV. We use the most commonly considered scenario in which the flux is dominated by the decay of pions and their daugther muons giving rise to a flavor ratio of (1/3, 2/3, 0) at source. This ratio tends to equalize at Earth due to neutrino oscillations averaged over astronomical distances. Hence, in eq. (3.19), we set f i = 1/3 for i = e, µ, τ . Also, an equal ν andν flux is used [40]. Regarding the spectral parameters φ 0 and γ, we take the best-fit values obtained in ref. [41] by performing a maximum-likelihood combination of the results from six different IceCube searches. The spectral parameters resulting from this analysis in the case of the single-power law are given by In order to illustrate the LQ component expected for the number of events, we have applied eq. (3.16) to 15 bins of deposited energy in the range [10 TeV, 10 PeV]. In figure 4 we show the LQ component to shower-and track-like events along with their sum for |λ 1 1 | = |λ 1 2 | = 1 and for different values of the LQ masses ranging between 500 GeV and 1 TeV. We note that, in our scenario, the LQ contributes to shower-like events via NC processes initiated by ν e,µ (andν e,µ ) or CC ν e N (andν e N ) interactions; in the case of track-like events, the LQ contribution arises only from the CC process ν µ N → µ − X (andν µ N → µ + X). An important feature of the distributions in figure 4 is that the regions of deposited energy at which they peak increase with the LQ mass. This general behaviour is inherited from the distribution of the number of events with respect to E ν . On the other hand, due to the fact that NC and CC processes deposit different ammounts of energy, the distributions of track-like events exhibit the threshold at m 2 χ /2M N for resonant production of χ 2 (middle panels of figure 4), while those corresponding to shower-like events keep different from zero for all the considered bins (upper panels of figure 4).   (= m χ 2 ). The upper and middle panels correspond to shower-and track-like events, respectively, and the lower panels display the total number of events.
In the following subsection, we add the LQ contribution in top of the spectrum expected from the SM + background hypothesis and study its implications on the spectrum actually observed by IceCube.

Statistical Analysis and Results
We consider the number of events in the i-th bin of deposited energy, n i , as a poisson variable and parameterize the respective expected number of events as where µ ≡ |λ 1 1 | 2 + |λ 1 2 | 2 , µ y s i is the number of events arising from the LQ contribution and y b i the number of events expected from the SM + background model. The values y s i were computed by using eq. (3.16), whereas the y b i 's were taken from ref. [1]. In order to estimate the µ parameter, we minimize the following statistic test, where L(µ) is the likelihood function. We note that the minimization of χ 2 (µ) is equivalent to the maximization of the likelihood function and that the last term of eq. (3.23) can be dropped during the minimization. The results obtained for LQ masses between 500 GeV and 1.5 TeV are shown in table 1, whereμ denotes the value of µ that minimizes χ 2 (µ), and where the cases Ω = 4π and 2π sr have been considered.
For the two lowest masses considered here, namely 500 and 600 GeV, we obtain non-physical values forμ. This result can be understood from the fact that the event distributions corresponding to these masses peak in the region between 100 and 1000 TeV, where the SM + background prediction is already above the observed spectrum. For LQ masses between 700-1200 GeV, we obtain increasing positive values ofμ, which is in fact expected since the LQ contribution deacreases with increasing m χ (see figure 4). Finally, for m χ > 1200 GeV, the maximum of the corresponding event distributions lies at the right end or even beyond the range of deposited energy considered at IceCube. Accordingly, LQs with these masses contribute mainly to the most energetic bins in which no event have been observed, forcing theμ's to decrease and even to become negative for m χ = 1500 GeV.
We have also used the estimates in table 1 to obtain upper limits for the parameter µ. In this case we use the following statistic for testing values of µ such that µ ≥μ [42,43], where λ(µ) is the profile likelihood ratio, andν i is obtained from eq. (3.22) with the replacement µ →μ. Then, the 95% CL upper limit is defined as the maximum value of µ for which p µ ≥ 0.05, with the p µ -value computed as follows where q µ,obs is the value of q µ obtained from the data and f (q µ |µ) is the probability density function (pdf) of q µ assuming the data correspond to the value µ. In the |λ 1 1 |-|λ 1 2 | plane the 95% CL contour is simply a circle of radius √ µ. For this reason, we list in table 2 the 95% CL upper limits on the quantity √ µ rather than µ.
In order to further specify the improvement in the fit obtained by adding the LQ contribution, we quantify the level of disagreement between the data and the hypothesis µ = 0. For this purpose, we use the statistic test q 0 = −2ln(λ(0)) forμ ≥ 0, and compute the respective p-value as where f (q 0 |0) is the pdf of q 0 assuming the SM + background hypothesis (µ = 0). We note that the data is consider to show lack of agreement with the hypothesis µ = 0 only ifμ > 0. Thus, we apply this test only to the LQ masses for which a physical value ofμ was obtained. In figure 5 we show the p 0 -value as a function of the LQ mass in the range [700 GeV, 1.4 TeV]. We see that the hypothesis µ = 0 cannot be rejected conclusively in any of the considered cases. For m χ = 700 GeV, the level of disagreement between data and the SM+backgound hypothesis is such that the latter could be rejected at a confidence level of 56%. This confidence level increases with the LQ mass and attains its maximum (minimum p 0 ), given by ∼ 69.5%, at m χ 1025 GeV. In the lower panel of figure 6, we show the total number of events observed at Icecube along with the predictions from the SM + background component and when the LQ contribution corresponding to m χ = 1025 GeV is added on top of it. For masses deacreasing from 1025 GeV to 700 GeV, the p 0 -value increases, indicating that the fit worsen (see the upper left panel of figure 6). Since the LQ contribution of the smaller masses affects mainly the bins of deposited energy where the majority of the events appear, the correspondingμ is forced to small values which leads to a negligible impact on the two bins that exhibit a weaker agreement with the SM + background explanation. As long as the LQ mass increases, their contributions become maximum at the region around these two bins, improving the fit with higher values ofμ (see table 1). The p 0 -value also increases for masses higher than the best fit value because the maximum of the respective LQ contributions moves away from the two bins between 2-3 PeV and start to affect the most energetic bins in which no events have been observed (see, for example, the upper right panel of figure 6). This leads to worse fits and for m χ > 1200 GeV pushes theμ again towards smaller values (see table 1).

LHC Constraints
In this section we discuss the most up-to-date LHC constraints on our colored electroweak-triplet scalar. We begin our discussion by filtering the model through the latest 8 TeV data. Our framework leads to five distinct final state topologies that we classify as follows We simulate the (a) and (b) topologies separately using MadGraph 5 [44]. We implement PYTHIA [45] for the parton shower and hadronization and the detector simulation is carried out using Delphes 3 [46]. We simulate the two topologies in a separate manner since the lepton misidentification rate is very small and events with final states containing only jets and missing energy will not significantly contribute to final state topologies containing leptons. In fact, the electron fake rate can be anywhere between 10 −4 and 10 −5 [47] while a recent study finds a muon fake rate of 2×10 −5 [48]. The events are generated for masses in the range 600 < m χ < 1200 GeV, for different combinations of the couplings assuming that λ i j = 0 if i = 1 and/or j = 3. In order to set bounds on the parameter space of the model, we use the latest CheckMATE validated analyses [49]: • ATLAS search for squarks and gluinos with jets and missing momentum [50], • ATLAS search for third generation squarks via charm quarks or compressed supersymmetric scenarios [51], The IceCube data as well as the SM + background fit were taken from ref. [1].
for the (a) topologies, and • ATLAS search for direct top-squark pair production in final states with two leptons [53], • ATLAS search for top squark pair production with one isolated lepton and missing transverse momentum [54], • ATLAS search for supersymmetry in events containing a same-sign dilepton pair, jets and large missing transverse momentum [55], • ATLAS search for direct slepton and chargino production in final states with two opposite-sign leptons, missing energy and no jets [56] for the (b) topologies. In figure 7 we show results for LQ masses in the range 600-1200 GeV in the λ 1 1 -λ 1 2 plane after applying all of the 8 TeV LHC results listed above. We compare our results to the 95% upper confidence limits on the number of signal events using the variable r defined in [49] given by where the numerator parametrizes the 95% lower limit on the number of signal events determined by CheckMATE and the denominator the 95% experimental limit on the number of signal events.
Regions of parameter space are excluded if r ≥ 1. In figure 7 we depict, for all LQ masses, the r = 1 contour with a black solid line. We do not show results for masses below 600 GeV, since for couplings λ 1 j > 0.1, which is the case for the simulations performed in this work, these masses are not allowed by current experimental constraints.  We also note that ATLAS and CMS have dedicated searches for first and second generation LQs with the 8 TeV [57,58] and 13 TeV [59,60] data sets, with a slight improvement on the limits with the latter. The searches target LQ pair production. The CMS collaboration focuses primarily on the second generation and place limits of 1165 and 960 GeV for LQ branching fractions of 0.5 and 1 respectively using 2.7 fb −1 of data, while the ATLAS collaboration places limits of 1100 and 1050 GeV for first and second generation LQs respectively when the branching ratio is 100% to a lepton and a quark. In addition, the ATLAS collaboration obtains limits varying the branching ratio into electrons and muons which are shown in figure 7 of [59].
In order to apply these LQ dedicated searches to our model we follow a conservative approach since the only component of χ that decays purely to a charged lepton and a quark is χ 3 ; hence limits on the mass of χ 1 and χ 2 will turn out to be much weaker. However, since we are assuming mass degeneracy to avoid tensions with electroweak precision data (EWPD) [28,37], the limits apply across the components of χ. In addition, since we are assuming that χ primarily couples the first family of quarks to electrons and muons, the decay width must be saturated with these two decay modes. As a consequence, the constraints given in figure 7 of [59] only imply that our LQ must lie above 900 GeV. Therefore, our model is basically unconstrained by the LQ searches and the limits derived from the more general searches described above dominate.

Low energy physics observables
The renormalizable interactions introduced in eq. (2.3) can lead to rare flavor changing and CP violating processes both at tree-level and at the one-loop level. Our working assumption is that χ couples primarily the first family of quarks to the electron and the muon and helps us to avoid the most stringent bounds arising from tree level semi-leptonic and leptonic meson decays as well as semileptonic τ decays. However, our LQ can yield new contributions to muon rare decays such as µ → eγ, the magnetic dipole moment of the muon, and atomic parity violation measurements. We discuss these constraints below.

µ → eγ and (g − 2) µ
Our LQ, a colored electroweak-triplet scalar, can give rise to lepton flavor violating decays such as µ → eγ as well as a contribution to the muon anomalous magnetic moment, a µ . Both contributions come in at the 1-loop level. The Feynman diagrams for the µ → eγ decay process are depicted in figure 8.
To place constraints in our model we follow the conventions used in [28] where the relevant parts of eq. (2.3) contributing to the µ → eγ decay and the muon's (g − 2) can be expressed as where V is the Cabibbo-Kobayashi-Maskawa mixing matrix. The above expression was obtained by starting from a mass-ordered mass eigenstate basis for the down type quarks and charged leptons and applying the following transformations: u iL → (V † ) ik u kL , d iL → d iL , and e jL → e jL . Since V 12 is roughly 20% of V 11 , we will assume that the coupling of both the muon and the electron to downand up-type quarks is the same. With this working assumption, and using the following effective Lagrangian for the µ → eγ decay In the same way, loops of LQs can modify the anomalous magnetic moments of leptons, a l . The effective Lagrangian parameterizing modifications to a l can be written as L a l = e ·l a l 4m l σ µν F µν l.
The contribution to a µ in the m q /m χ → 0 limit from χ 2 and χ 3 is given by [28] a µ ≈ 9 32π 2 m 2 µ m 2 The most precise experimental result on (g − 2) µ was obtained by the E821 experiment carried out at BNL [62,63]. The deviation from the SM value is given by δa µ = (2.8 ± 0.9) × 10 −9 where the SM value is given by a SM µ = 1.16591803(70) × 10 −3 [64]. Using this result we can directly constrain the value of λ 1 2 : From the two constraints discussed above, one can see that one scenario of interest could lead to a very suppressed value of λ 1 1 compared to λ 1 2 . In particular, for LQ masses in the TeV range, one needs λ 1 1 ∼ 10 −3 for O (1) λ 1 2 couplings. These scenarios are not unnatural if one takes into account specific flavor models where quarks transform as different non-trivial singlets of A 4 [34]. Below we will discuss how this specific scenario is also consistent with low energy precision measurements such as atomic parity violation.

Atomic Parity Violation
Below the electroweak scale the parity violation such as in the Cesium 133 atom can be studied with the following effective Lagrangian The SM maximally violates parity and one can calculate very precisely the values of C 1u and C 1d with C SM 1u = (−1/2 + 4/3 sin 2 θ W ) and C SM 1d = (1/2 + 2/3 sin 2 θ W ), where θ W denotes the Weinberg angle of the SM. Using these values one can define a nuclear weak charge by (5.11) where Z and N are the number of protons and neutrons respectively. For cesium, the experimentally measured value of Q W is −73.20(35) [65]. Using this measurement one can extract strong constraints on the LQ couplings to the first quark family. In particular, one can parametrize the contributions arising from χ by δC 1u and δC 1d . Given that in the SM Q W = −73.15 (35) [66], using eq. (5.11) with C 1u = C SM 1u + δC 1u and C 1d = C SM 1d + δC 1d and assuming that the coupling of the electron to the upand down-type quarks is the same, as discussed in the previous section, one can extract the following matching contribution to δC 1u = δC 1d = δC 1 [28]: where G F denotes the Fermi constant. With the above result and the experimentally measured value of Q W one has the following bound on λ 1 1 : which is roughly four times stronger than the bound on λ 1 2 derived from the measurement of the muon anomalous magnetic moment. In light of this result and the bound arising from the µ → eγ rare decay, our framework leans towards values of λ 1 1 which are suppressed relative to λ 1 2 .

Discussion and Concluding Remarks
Taking into account the analysis performed in sec. 3, we conclude that, in order to improve the explanation of the spectrum of UHE neutrinos observed at IceCube through the addition of a LQ triplet, higher values for the mass m χ 1 = m χ 2 = m χ are preferred. Additionally, since the rate of events expected from the LQ component decreases with the LQ mass, large values of µ = |λ 1 1 | 2 + |λ 1 2 | 2 are also required (see table 1). Specifically, under the hypotheses used in eq.(3.21) and described in sec. 3.2, we have found that the best fit of the four year IceCube data is achieved when the LQ mass is approximately 1025 GeV and the couplings are such that µ = 2.189 (see figure 6). We note that this mass is allowed by the dedicated searches of LQs in the LHC at 8 TeV and 13 TeV.
Regarding the 95% CL limits derived from the IceCube data in table 2, we see that these are considerably weaker than the constraints placed by the general searches at the LHC at 8 TeV listed in sec. 4 (see figure 7). This is mainly due to the lack of statistics in the most energetic bins of the IceCube spectrum, where the data is not sufficiently explained by the SM expectation and the LQ contribution may become more relevant.  Figure 9: Contours corresponding to the fit of the IceCube spectrum along with the 95% CL LHC upper limits for m χ = 800 GeV, 1000 GeV and 1200 GeV.
The estimates of the parameter µ for different LQ masses shown in table 1 can be confronted with the constraints arising from the 8 TeV LHC data. In figure 9, we display theμ values obtained from the fit to the IceCube data along with the 95% CL LHC limits for the masses 800 GeV, 1000 GeV and 1200 GeV. We see that in the case of a fully opaque Earth, Ω = 2π sr, the contours preferred by the IceCube data are excluded by the LHC upper limits at 95% CL. This is also the case when no attenuation effects are considered, Ω = 4π sr, for 1000 GeV and 1200 GeV, while for 800 GeV the ranges 0.20 < λ 1 1 < 0.40 and 0.55 < λ 1 1 < 0.67 are not ruled out. After analyzing the low energy constraints on our framework, we are led to conclude that a scenario of interest will include a very suppressed value of λ 1 1 compared to λ 1 2 . This resulted from a combination of the rare µ → eγ decay and atomic parity violation measurements and our working assumption that χ coupled primarily the first family of quarks to the electron and the muon. Therefore, by taking |λ 1 1 | sufficiently small and |λ 1 2 | ∼ 1.47, the parameters that give the best fit of the IceCube data, m χ = 1025 GeV andμ = 2.4, are compatible with the low energy constraints and also, as said above, with the dedicated searches of LQs at the LHC. However, as shown in figure 9, this scenario for the LQ triplet is clearly in tension with the 8 TeV LHC constraints. On the other hand, even though in the idealized case of Ω = 4π sr a LQ with mass around 800 GeV is not ruled out by these constraints, its contribution to the spectrum above PeV is not significant (see figure 6). Furthermore, such a value for the LQ mass is in conflict with the LHC 13 TeV dedicated searches if one requires the LQ to decay only to electrons and muons. Loosing this requirement with an additional decay mode will necessitate a more dedicated recast of LHC searches.
Acknowledgments This work has been partially supported by ANPCyT under grants No. PICT 2013-0433 and No. PICT 2013-2266, and by CONICET (NM, AS). ADP was supported in part by the Natural Sciences and Engineering Research Council of Canada (NSERC). In addition, the authors have benefited from conversations with Estefania Coluccio Leskow.

A Attenuation effects in the Earth
The rate of upward-going neutrinos is reduced due to the interactions of the incoming neutrinos with the nucleons in the Earth. Although the interactions with electrons can be important at E ν 6.3 PeV, where the resonant production of the W boson takes place, we will focus in this appendix on the neutrino-nucleon interactions. The water equivalent interaction length due to neutrino-nucleon interactions is given by where N A = 6.022 × 10 23 cm −3 (water equivalent) is Avogadro's number. We note that every neutrino (antineutrino) flavor has a different interaction length according to the specific cross section describing its interactions with nucleons. In order to study the attenuation effects, the interaction length should be compared with the ammount of material encountered by an upward-going neutrino, which is konwn as the column depth and depends on the angle of incidence of the incoming neutrinos. The thickness of the Earth as a function of the cosine of the angle of incidence is shown in figure 10, where the density profile of the Earth given in ref. [35] have been used. The maximum column depth is 11 kilotonnes/cm 2 , and correponds to a neutrino emerging from the nadir. By plugging the function z(θ) and the interaction length given in eq. (A.1) into the eq. (3.18), we obtain the shadow factor S(E ν ) (see ref. [35]).  Figure 10: Column depth as a function of the angle of incidence of the incoming neutrinos.
As mentioned in sec. 3.2, the shadow factor depends on the neutrino-nucleon cross section via the interaction length (see eq. (3.18)) and therefore the addition of the LQ contribution could have in principle an impact on the reduction of the rate of northern events. In order to study this possibility, let us define first the LQ contribution to the total neutrino-nucleon cross section for a flavor as follows where λ 1 = λ 1 1,2 for = e, µ respectively, andσ LQ νN is the same for the two flavors. A similar relation can be written as well for antineutrinos. Also, we denote the SM contribution as σ SM νN and define With these definitions the number of southern events for a given flavor can be written as and adding all flavors of neutrinos and antineutrinos we obtain  In figure (11), we show the ν e shadow factor for the SM hypotheses along with the deviations produced by adding different LQ contributions. These contributions correspond to the mass that gives the best fit to the Icecube data, m χ = 1025 GeV (see sec. 3.3), and |λ 1 1 | 2 = 1-6. Moreover, we display the relative difference between the effective solid angle for the SM hypothesis, Ω SM ≡ 2π(1 + S SM νe ), and for the SM+LQ hypothesis, Ω tot ≡ 2π(1 + S total νe ). From the left plot, we see that the shadow factor is a decreasing function of the incoming neutrino energy that, as expected, begins to deviate from the SM behaviour above the energy threshold associated to the specific LQ contribution, namely m 2 χ /2M N . However, the deviation due to the addition of the LQ contribution is not meaningful as can be concluded from the right plot in figure 11. Indeed, the relative difference between the effective solid angles is less than 9%, even for a squared coupling as large as 6. We note that the case of the muonic neutrino is entirely analogous (with the replacement λ 1 1 → λ 1 2 ). angle is very small, with the maximum value of the ratio being obtained for α = π/4. Thus, we have set |λ 1 1 | = |λ 1 2 | and computed the ratio R for µ in the range 0-10. From figure 12, we see that the ratio decreases as µ increases, but the deviation of the SM expectation is at most 17% for µ as large as 10. This result is consistent with the conclusions derived from figure 11; since the interaction length is dominated by the SM contribution, the impact of the LQ contribution in the disbalance between events coming from the two hemispheres is not significant and this is reflected in the plot of R as a function of µ.