Seeking leptoquarks in IceCube

We investigate the sensitivity of IceCube(-Gen2) to a scalar leptoquark scenario with couplings only to heavy quark flavors which may be connected to solving discrepancies in B-meson semileptonic decays. We take into account, for the first time, the non-negligible neutrino-gluon cross section induced by leptoquarks, and we systematically account for indirect and direct constraints which have been overlooked in previous studies. We conclude that IceCube(-Gen2) can only probe the light LQ regime, already disfavored by the combination of flavor physics constraints, electroweak precision data and the direct searches at the LHC.


Introduction
It is remarkable that the two types of constituents in the flavor sector of the Standard Model (SM), quarks and leptons, although independent, share a number of common features. They come in three generations that conspire to exactly cancel the triangle anomalies of the gauge interactions, they exhibit mixing between mass and flavor eigenstates resulting in oscillations, and flavor changing neutral interactions do not occur at tree-level. These common features strongly suggest that quarks and leptons should be somehow interrelated. In fact, leptoquarks (LQs), bosons that can mediate quark-lepton transitions, appear in many extensions of the SM, in particular, in Grand Unified Theories [1,2], and composite Higgs models [3][4][5].

JHEP06(2018)032
From a phenomenological point of view, LQs are SU(3) c colored particles that can come as SU(2) L singlets, doublets or triplets [6,7]. Their masses, strength and flavor structure of their couplings to quarks and leptons are all undetermined when the ultraviolet completion is unknown. Since LQs with masses of the order of the electroweak scale are not disallowed by any fundamental reason, it is important to look for them in all possible ways. In fact, many collider experiments have been searching for them in pair as well as in single production modes. From HERA to the LHC, experiments have been setting stringent limits on their masses and couplings, cf. [8,9] and references therein. These limits depend, unfortunately, unavoidably on the assumed flavor structure, some of which are difficult to explore in accelerators. Other bounds deduced from the direct searches have been discussed in ref. [10], see also references therein.
Besides the purely theoretical reasons that compel experimentalists to look for these particles it seems there might be an experimental one too, coming from B-physics. The excellent agreement between the SM predictions and experimental data seems to be disrupted by flavor physics. Recent data from the LHC, Belle and BaBar suggest lepton flavor universality violation (LFUV) in B-meson semileptonic decays. The LHCb experiment, for instance, has reported about 2.6σ discrepancy in the ratio of the partial branching fractions of the loop induced transitions (denoted by B ), and the measured values [11,12] appear to be smaller than predicted in the SM [13,14]. Another intriguing indication of LFUV comes from the tree-level ratio [15] R where = e, µ, obtained by combining Babar and Belle results [16][17][18][19][20] and 2.2σ above R SM D = 0.300 (8), the SM prediction [15,21,22]. This feature was confirmed with R exp D * = 0.304 (15) [16,17,23], which also appears about 3σ larger than the SM prediction, R SM D * = 0.260 (8) [24,25]. Another confirmation of this phenomenon was reported by LHCb who measured R J/ψ = B(B c → J/ψτ ν)/B(B c → J/ψµν) = 0.71 ± 0.25 [26], also larger than the SM estimates although the SM value in this case is less reliable. Many proposed solutions to these anomalies discussed in the literature so far rely on existence of LQ states , notice that similar solutions can be implemented by considering supersymmetric models including R-Parity violating squarks [54]. The couplings of LQs to heavy quarks can be as peculiar as to make them difficult for direct observations in colliders and it is therefore useful to try and look for them in other ways.
Here we study a possibility of using the IceCube neutrino telescope in the South Pole as a complementary way to search for those states. Many authors have explored the potential of IceCube to test LQ models, in particular, after the observation of an excess of high energy stating events consistent with a flux of high-energy astrophysical neutrinos from outside our galaxy [55]. In ref. [56] electroweak singlet LQ that couple to either second or third generation of quarks and/or leptons were studied, emphasizing the use JHEP06(2018)032 of the inelasticity distribution of events to improve background rejection. Refs. [57,58] and [59] evoke LQ couplings to the first generation of quarks and leptons in order to try to accommodate IceCube high-energy neutrino data. Similarly, in [63] a supersymmetric model including R-parity violating squark resonances is considered to explain IceCube data. The limits that IceCube data can impose on a specific low-scale quark-lepton unification model with scalar LQs was studied in [60]. Scalar LQ doublets that couple the first and third family of quarks to leptons, with a flavor structure that can address the anomalies in R K ( * ) , pass other LHC constraints and give a sizable contribution to IceCube neutrino data were considered in [61] and by a more elaborated analysis in [62]. The LQ states in all of these works, with exception of ref. [56], couple to the first generation of quarks. Couplings to the first generation suffer, however, from severe experimental constraints as they have to comply with limits from direct searches at colliders, atomic parity violation experiments and flavor physics observables in the kaon/pion sectors. Moreover, all previous studies of the neutrino-nucleon cross-section at IceCube only consider the resonant s-and t-channel interactions of high energy neutrinos with quarks (valence and/or sea) via a LQ exchange. Nonetheless, we should not forget that in order to compute neutrino-nucleon cross sections we must also include the neutrino interaction with gluons mediated by a LQ, particularly important in accessing LQs that only couple to heavy quarks. Those are, in fact, the kinds of LQs that can take part in solving the discrepancies between data and heavy meson semileptonic decays. In this paper, we include for the first time the neutrino-gluon cross-section in order to consistently study the sensitivity of IceCube to LQ scenarios. Furthermore, we systematically compare the future IceCube future sensitivity with the indirect constraints coming from flavor physics and Z-pole observables, as well as the direct limits from the LHC, which are very efficient in constraining the LQ couplings.
Our paper is organized as follows. In section 2 we describe a LQ model with a particular simple flavor structure that can solve R D ( * ) and give rise to new neutrino-nucleon interactions which can be probed by IceCube. Furthermore, we perform a complete study of the indirect bounds on LQ couplings stemming from flavor physics and Z-pole observables. In section 3 we review the High Energy Starting Events sample collected by IceCube after six years and we discuss the current status of the SM fit. We explain the procedure to simulate both the SM and LQ events at IceCube in order to compare with experimental data. In section 4 we present results of our analysis and discuss what we can learn about our model from six years of IceCube data as well as the prospects for what can be expected in the future. We finally conclude in section 5. All LQ contributions to the neutrinonucleon scattering cross section, including the new ones, are collected in appendix A, while in appendix B we describe the transport equation for Earth's attenuation. The decay rate expressions for B → D ν are provided in appendix C.

A simple leptoquark scenario to explain R D ( * )
In this paper we focus on the so called R 2 -model which involves an electroweak doublet of scalar LQs with hypercharge Y = 7/6, described by the Yukawa Lagrangian, where in the first line the Lagrangian is given in the flavor basis and the superscripts in the LQ fields refer to the hypercharge (Y ), while in the other two lines the Lagrangian is given in the mass eigenstate basis in which the superscripts of LQs refer to the electric charge Q = Y + T 3 , where Y is the hypercharge and T 3 the third component of weak isospin. Note that ∆ = iτ 2 ∆ * is the conjugate of the doublet of mass degenerate LQs, while V and U stand for the Cabibbo-Kobayashi-Maskawa (CKM) and the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrices, respectively. Furthermore, P L,R = (1 ∓ γ 5 )/2 are the chiral projectors, T denote quark and lepton SU(2) L doublets, whereas u L , d L , L and ν L are the fermion mass eigenstates. We opt for the minimalistic structure of the Yukawa coupling matrices g L,R and assume that allows ∆ (2/3) to mediate a tree-level contribution to R D ( * ) and generates a contribution to the neutrino-nucleon scattering in IceCube.

R 2 explanation of R D ( * )
In order to confront the LQ contributions with the experimental data involving the B → Dτ ν decay we consider the effective Hamiltonian where g S,T are the Wilson coefficients induced by the LQ state mediating the semileptonic decay via a tree-level contribution shown in figure 1. After integrating out ∆ (2/3) , the expression for g S,T , at the matching scale µ = m ∆ , reads: clearly enhanced by V −1 cb . By virtue of the leading order QCD running from µ = m ∆ ≈ 1 TeV down to µ = m b , the above relation between g S and g T becomes (2.5) Furthermore, we included the small electroweak corrections discussed in ref. [64] which can induce mixing between the scalar and tensor contributions.  We focus on B → D semileptonic decays for which the hadronic form factors have been computed by means of numerical simulations of QCD on the lattice in refs. [21,22]. 1 To fill the discrepancy between R SM D and R exp D one then allows for g S = 0. In this way we find the region of allowed values of Re[g S (m b )] and Im[g S (m b )] which is depicted by the blue regions in figure 2 to 1, 2 and 3σ accuracy. We decided not to include R D * in the same fit because the lattice QCD determination of the form factors at non-zero recoil has not been made, in addition to the fact that the relevant pseudoscalar form factor has never been computed on the lattice. To verify that our solution to R D also gives an improvement of R D * , despite the uncertainties mentioned above which need to be clarified by lattice simulations, we show in the same plot the 2σ region consistent with R exp D * . The latter region was obtained by using the B → D * ν form factors extracted from experimental results [15], combined with the ratio A 0 (q 2 )/A 1 (q 2 ) and T 1−3 (q 2 )/A 1 (q 2 ) computed in ref. [65]. By using this approach, we conclude that, in the R 2 -model, only the complex g S solutions to R D can provide a reasonable agreement with R exp D ( * ) , a fact that agrees with findings of ref. [66]. Therefore, we will focus our analysis on the solutions to R D with complex Wilson coefficients obtained to 1σ accuracy. This relation implies that the magnitude of LQ couplings should be such that where we used |V cb | = 0.0417(20) [67]. We see that the couplings needed to explain R D for m ∆ 1 TeV should not be too large, mostly due to the CKM enhancement in eq. (2.4).
In the following we show that the explanation of R exp D ( * ) > R SM D ( * ) described above is consistent with other limits from flavor physics and the direct searches for pair-produced LQs at the LHC. To that purpose, we will assume that the coupling g bτ R is purely imaginary.
JHEP06(2018)032 Figure 2. Regions of allowed values for g S (µ = m b ) on the complex plane, compatible with experimentally measured R D to 1−, 2− and 3σ (from darker to lighter blue). We assumed g S (µ = m ∆ ) = 4×g T (µ = m ∆ ), as predicted by the R 2 LQ scenario, and considered the QCD running from µ ≈ 1 TeV to m b , cf. eq. (2.5). We also show in black the region compatible to 2σ accuracy with R exp D * , which indicates that the effective couplings must be complex in this scenario. See discussion in the text.

General constraints
The choice of g cτ L = 0 and g bτ R = 0 in eq. (2.2) allows us to avoid many flavor physics constraints, making this explanation of R exp D ( * ) > R SM D ( * ) particularly simple. Processes mediated by the flavor-changing neutral ∆F = 1 and ∆F = 2 currents are not altered in this scenario. Only the charged-current transitions b → cτν are modified, as desired. The most significant constraints to this scenario actually come from the Z boson decay widths.

Constraints from Z decays
Our LQ model induces modifications to Z-boson decays through loop contributions, as depicted in figure 3. Those effects have to be consistent with the LEP measurements of B(Z → + − ) [68]. The Lagrangian describing the Zff coupling can be written as  where g is the SU(2) L gauge coupling, and The corresponding branching ratio is then given by where m f is the fermion mass and β f = 1 − 4m 2 f /m 2 Z . The LQ couplings lead to up-and down-type quark loop contributions to δg τ L(R) , which are given by and where x t = m 2 t /m 2 ∆ and x z = m 2 Z /m 2 ∆ , and N c = 3 is the number of colors. 2 The dominant term in δg τ R comes from the top quark since it enhances the loop function. The effective 2 One can easily identify the top and bottom quark contributions in the first two lines of eq. (2.11), while the ones coming from the up and charm quarks have been recast in the same equation by using the unitarity of VCKM, i.e. |V ub g bτ

JHEP06(2018)032
coupling δg τ L does not exhibit a similar mass enhancement, and it is for this reason less relevant for phenomenology. Similarly, the LQ interaction with the charm quarks also induces an effective coupling to τ neutrinos, which reads The effective couplings given above are constrained by the LEP measurements of the Z decay widths and other electroweak observables [68]. The measurements of the effective Z couplings give [69], These results translate into useful bounds on the LQ contributions to δg τ L,R . Furthermore, the LEP bound on the effective number of neutrinos [69] N ν = 2.9840 (82) , will constrain the LQ couplings via the expression, (2.16) As we shall see below, the coupling g bτ R is tightly limited by B(Z → τ τ ) due to the enhancement of the loop function by the top quark mass. On the other hand, the constraints on g cτ L derived from N ν would be relevant only for very light LQ states, already excluded by the constraints stemming from g τ V,A . Finally, we have also checked that the bound on the LQ couplings arising from B(Z → bb) is not significant with current experimental precision.

Limits from direct detection
Direct searches at the LHC for pair-produced LQs provide useful limits on the couplings g cτ L and g bτ R as a function of the LQ mass m ∆ . For the flavor ansatz given in eq. (2.2) the allowed LQ decay modes are The CMS collaboration recently improved bounds of LQs decaying to bτ , obtaining a lower bound on m ∆ of 900 GeV, with the assumption B(∆ (2/3) → bτ ) = 1 [70,71]. Similarly, CMS excludes LQs decaying into tτ with masses m ∆ 850 GeV if B(∆ (5/3) → tτ ) = 1 [72,73]. Furthermore, a search for pair-produced LQs decaying into a light quark (q = u, d, s, c) and a neutrino has been reported by CMS, allowing us to set the limit m ∆ 980 GeV if B(∆ (2/3) → cν) = 1 [74].
In any realistic scenario, the limits mentioned above have to be reinterpreted to account for values of the branching ratios smaller than one. The reinterpretation of these exclusion and where, for simplicity, we assumed |V tb | ≈ 1 and neglected the light fermion masses (m c , m b and m τ ). From the above formulas we conclude that scenarios with m ∆ 1 TeV can comply with the limits of figure 4 if the branching ratios are diluted by a large coupling to the charm quark. In other words, the ratio |g bτ R |/|g cτ L | should be of O(1), which also agrees with the constraints on g bτ R derived from Z → τ + τ − , cf. eq. (2.11). We checked that m ∆ can be as low as ≈ 650 GeV for some combinations of Yukawa couplings.
In figure 5, we confront the allowed region by R D with the indirect and direct constraints discussed above by assuming g bτ R to be purely imaginary, as discussed in section 2.1. For instance, we see in these plots that both |g bτ R | and g cτ L are bounded from above by the Z-pole observables. Furthermore, the couplings needed to explain R D for m ∆ 1 TeV are not very large, mostly due to the CKM enhancement in eq. (2.4), and they can therefore be perfectly consistent with phenomenological constraints discussed in this section.
In what follows, we will discuss the IceCube phenomenology of the LQ scenario considered in this paper. In practice, we will focus on the analysis of the parameter space defined by m ∆ and g cτ L since they are the relevant parameters to determine IceCube signals. In this plane the region allowed by the Z-pole observables is approximately limited by a straight line, such that the allowed values are g cτ L 1.6 + 1.3(m ∆ /300 GeV − 1), which is also illustrated in this paper in figures 12 and 14.

High energy neutrino events in IceCube
The cubic kilometer IceCube Neutrino Observatory in the South Pole observed between 2010 and 2012 the first evidence of a diffuse flux of astrophysical neutrinos above 100 TeV [55], which was further confirmed in the updated analyses [75][76][77]. In this paper we consider the 6-year sample [77] which contains 80 High Energy Starting Events (HESE). In particular, we focus on the 28 events with deposited energies above 100 TeV, which are divided in 23 showers and 5 tracks. This flux, of yet unknown astrophysical origin, seems to be isotropic and consistent with a single power law behavior followed by neutrinos and antineutrinos of all flavors [78]. In other words, the neutrino energy spectrum can be parameterized by a spectral index, γ, and a normalization constant C 0 as so that at the Earth the flavor fluxes are given by where f ⊕ α , α = {ν e , ν µ , ν τ ,ν e ,ν µ ,ν τ }, are the normalized neutrino and antineutrino flavor fractions at the Earth.
Assuming the flavor composition to be the standard one expected from neutrino oscillations, i.e., (f ⊕ νe : f ⊕ νµ : f ⊕ ντ ) = (f ⊕ νe : f ⊕ νµ : f ⊕ ντ ) = (1 : 1 : 1), IceCube best fit of 6-year sample yields γ = 2.92 +0. 33 −0.29 and C 0 = 2.46±0.8 GeV cm −2 str −1 s −1 [77]. To test our fitting JHEP06(2018)032 Figure 5. The allowed region on the plane g cτ L × |g bτ R | by current experimental constraints is shown in yellow for m ∆ ∈ {0.65, 0.8, 0.9} TeV. We assume g bτ R to be a purely imaginary coupling, as described in the text. The separate allowed regions both from the direct searches at the LHC, the Z-pole observables and R D are shown from lighter to darker blue, respectively. Thus, the yellow region is given by their intersection. See text for details.

JHEP06(2018)032
procedure we used the publicly available IceCube data, fit them to the SM and obtained 3 γ = 2.98 +0. 42 −0.38 , with a p-value of 0.25, in full agreement with the above-mentioned IceCube results [77]. The best fit curves for shower and track spectra are shown in figure 6. 4 For simplicity, in the fit we have used a fixed atmospheric background spectrum which is derived from the central values for the total amount of atmospheric neutrinos estimated in [77], which are N µ = 25.2 ± 7.3 and N atm = 15.6 +11.4 −03.9 . To reduce the effects of uncertainties we focus on the energy range [100 TeV, 10 PeV], in which the background contribution is sub-leading. Indeed, the computation of the best fit considering the minimum (maximum) expectations for the atmospheric background yields γ = 3.0 (2.9) and C 0 = 2.7 (2.4) with a p-value of 0.24 (0.23), which lies perfectly inside the allowed 1σ region. We will apply the same strategy to analyze LQ scenarios. Before passing onto that discussion we examine some aspects of the simulation of IceCube events, with special attention to the new contributions introduced by LQ interactions.

Calculation of the event rate in IceCube
Icecube detects individual neutrino events through the Cherenkov radiation emitted by charged secondary particles produced in the neutrino interaction with the Antarctic ice.
where c = sh, tr indicates the type of topology, h = S, N the events entering the detector from the south (S) or north (N) hemisphere, and X = NC , CC , e-scattering stands for the different type of processes we consider. Here, T is the exposure time in seconds, 3 The details about our statistical analysis considering the SM and several LQ scenarios are explained in section 4. We advance the results from the best fit of the SM in this section both to visualize graphically the IceCube 6-year dataset, figure 6, and also to establish some general but important remarks about the analysis. 4 Notice that the contribution of the atmospheric background to the neutrino spectrum, which we obtain by using the modelling of ref. [80], is consistent with the results shown in ref. [77] concerning absolute numbers and spectrum shape. However, the validity of the model is challenged at low energies, for instance by the discrepancy in the number of tracks around 30 TeV, which can be treated in a way discussed in refs. [62,90]. Note that in this work we focus on the events above 100 TeV and rely on the validity of the atmospheric background model at high energies.  N A = 6.022 × 10 23 g −1 is the Avogadro's number times the number of moles per each gram of protons/neutrons, and A h α accounts for the effect of Earth's attenuation, cf. eq.(3.3). The effective mass of the detector in grams is denoted by M eff (E t ), which corresponds to the mass of the target material convoluted with the efficiency of converting an event with true deposited energy E t into the observed signal. The resolution function, R(E t , E dep , σ Et ), is taken to be a Gaussian distribution with mean equal to the true deposited energy E t and standard deviation σ Et , both clearly dependent on the particular process taking place. Finally, dσ X α /dy is the differential cross section for process X in the inelasticity interval [y, y + dy]. Our calculations were performed in the same way as described in refs. [79] and [80], and we refer to these papers for further details.

Leptoquark contribution to the event rate
In our model, only the coupling g cτ L affects the neutrino-nucleon production cross section at IceCube, cf. eq. (2.1). The other coupling, g bτ R , only contributes to the total decay width of the LQ bosons. In this case, since the new interactions only involves quarks, neutrinos and τ -leptons, this scenario will only produce nonstandard shower-like events in IceCube. The event rate with ν τ in eq. (3.4) will then receive the New Physics contribution where the differential cross-section expressions dσ ∆ ντ /dy can be found in appendix A, with y being the inelasticiy variable. The same expression applies mutatis mutandis to theν τ event rate. Figure 7. Feynman diagrams contributing to ν + q u → ν + q u and ν +q u → ν +q u via a LQ exchange. Here ν , ν are generic neutrino flavors and q u , q u are generic up-type quarks. Figure 8. Feynman diagrams contributing to the process g + ν → q u + ∆, where the LQ particle decays intoq u + ν .

JHEP06(2018)032
In previous LQ studies [56][57][58][59][60][61][62], only the parton level processes ν + q were considered in the cross-section, cf. figure 7. In addition to the interaction with quarks, high-energy neutrinos can also interact with gluons via the transition ν + g → q u + ∆, with the LQ decaying inside the detector, as shown in figure 8. This contribution to the cross-section can be of the same order of the ν + q contribution when dealing with heavy quarks, q ∈ {c, b}, and thus cannot be neglected, as we illustrate in figure 9. This can be understood from the suppression of the s-and t-channel interactions with heavy quarks by the small values of the parton distribution functions (PDFs). Furthermore, the gluon PDF shows a fast growth when the Bjorken scale variable x approaches zero, giving a non-negligible contribution to heavy quark cross-sections.
To estimate the impact of LQ interactions to the IceCube data we compare in figure 10 the distribution of shower and track events generated by the SM and two representative LQ scenarios. These results were obtained under the restriction of a fixed total number of events. We see that the nonstandard effects appear only at high neutrino energies (E ν 500 TeV), anticipating a modest sensitivity to LQs with current IceCube data. This instigates us to investigate the situation with maximal future sensitivity of IceCube and IceCube-Gen2 [89] which we postpone to section 4.

Attenuation and regeneration effects
The neutrinos that arrive in IceCube from the south hemisphere do not suffer Earth effects, having then an attenuation factor of A S α (E ν ) = 1 2 for an isotropic flux in eq. (3.4). The neutrinos entering from the north hemisphere have to go through the Earth before reaching JHEP06(2018)032

JHEP06(2018)032
the IceCube detector with the attenuation factor of the resulting flux given by where θ is the nadir angle, and is simply a fraction of the initial flux after propagating though a column depth X(θ), where ρ(s) is the matter density at the point s, L(θ) = 2R ⊕ cos θ, and R ⊕ is the Earth's radius. In the above formulas, the function F α is defined by To calculate A N α (E ν ) we need to solve the coupled integro-differential transport equations given in appendix B which include absorption and regeneration effects via particle decays. We have followed the prescription described in ref. [81] to numerically solve these equations. Furthermore, we assumed that the absorption of τ and ∆ (2/3) is only given by their decays into lighter particles, neglecting the subdominant scattering contributions.
The average attenuation factor A N α (E ν ) for each (anti-)neutrino flavor is shown in figure 11 for the SM (left panel), and for a LQ scenario with m ∆ = 400 GeV and g cτ L = 2 (righ panel). As expected, the attenuation factors A N ντ (A N ντ ) are highly affected by LQ interactions for E ν 500 TeV. Note, however, that since the incoming flux falls rapidly with E ν , this effect turns out to be smaller than suggested by figure 11. Furthermore, we included the regeneration effects in our simulation which can lower the attenuation factors by at most 10%.

Statistical data analysis at IceCube, present and future
As stated in previous sections, we use the 6-year HESE sample in the energy range [100 TeV, 10 PeV] [77]. To quantify the goodness-of-fit of each model we follow the maximum likelihood method described in ref. [82]. The test statistic is defined by  Figure 11. Average attenuation factor A N α (E ν ) for each (anti-)neutrino flavor for the SM (left panel) and for a LQ model with m ∆ = 400 GeV and g cτ L = 2 (right panel). These results has been obtained by fixing γ = 2.9, but they only show a mild dependence on this parameter.
The total number of bins, N bins = 20, is comparable to the 28 events detected above 100 TeV. We consider a track to shower misidentification factor (missID) of 20%, which is reasonable considering that IceCube has measured it to be approximately 30% for low energy tracks [78] and since it is expected to decrease with energy. Nevertheless, we have verified that the final results are unaffected by using a 10% or a 30% missID factor, similarly to what has been obtained in ref. [62].
We reiterate that in our statistical analysis we consider 5 parameters: which we vary as m ∆ ∈ [300 GeV, 700 GeV], |g cτ L | ∈ [1,4], and we set g bτ R ≈ 0 since it does not directly contribute to the cross-sections relevant to this study.

Current data analysis
To evaluate the quality of the fit of a hypothesis H with respect to current data, we compute the minimum of the test statistic (4.1) and its corresponding p-value, obtained by simulating pseudo-data with the best fit values of γ and C 0 . As already mentioned in the text, we find that the SM best fit point yields γ = 2.98 and C 0 = 2.63 GeV cm −2 str −1 s −1 with a test statistic value of 17.05, and a p-value of 0.25. This indicates that the SM compatibility with the current data is quite acceptable.
The results obtained in the LQ scenario are shown in figure 12 for different values of LQ mass m ∆ and coupling |g cτ L |. We find that the best fit point is given by (m ∆ , |g cτ L |) = (500 GeV, 2.5) with a minimum test statistic of 16.7, and therefore the event distribution of the SM and the best LQ model are indistinguishable from the IceCube data, also shown JHEP06(2018)032 in left panel of figure 13. From the right panel of figure 13 we learn that only LQ models with very large g cτ L ≈ 3.5 and small masses m ∆ ≈ 300 GeV can be excluded by current IceCube data to 95% CL (p ≤ 0.05). These parameters, however, are already excluded by the limits arising from the Z-pole observables, as discussed in section 2.

Future IceCube(-Gen2) projections
To investigate the projected sensitivity of IceCube to our LQ scenario, we compute the projected exclusion regions in the plane (m ∆ , |g cτ L |) by assuming the SM is the null hypothesis. We considered the test statistic defined by the likelihood ratio between the LQ and the SM hypotheses, given by q = −2 log L(LQ)/L(SM) with L(LQ) and L(SM) defined analogously to eq. (4.1). The value of q is obtained after the minimization of each hypothesis with respect to (γ, C 0 ) by using projected data, which is generated from Monte Carlo simulations of the SM best fit to current data. More specifically, we evaluate the luminosity L, at which the projected data is simulated, by using multiples of current 6-year exposure. To establish a criteria for exclusion we consider the p-value associated to each tested value of q. We simply assume that q follows a χ 2 -distribution with k degrees of freedom, where k is the difference in the number of parameters between the LQ and SM hypotheses, which in our case is k = 2. Thus, for a given value of q, the p-value is given by p = exp(−q/2).
In figure 14 we plot our results as the projected 95% CL exclusion regions in the (m ∆ , |g cτ L |) plane, along with the limits arising from the Z-pole observables, as well as the region eliminated by the reinterpretation of LHC searches for pair-produced LQs decaying  Figure 13. Comparison between the SM and LQ best fit curves to the IceCube 6-year dataset. In the left panel we show the best fit curve for (m ∆ , |g cτ L |) = (500 GeV, 2.5), which is one representative case of scenarios with −2 ln L = 16.7. In the right panel, we show the best fit curve for (m ∆ , |g cτ L |) = (300 GeV, 4), which has a p-value around 0.05 and therefore it is in tension with current data. In the plot legend we quote within parentheses the best fit values of γ and C 0 for each model. into bτ [70,71], cν [74] and tτ [72,73], cf. discussion in section 2. We see that the parameter region accessible to IceCube increases at a modest pace with the luminosity L, reaching the corners (m ∆ = 300, |g cτ L | = 1) and (m ∆ = 600, |g cτ L | = 4) only with about 20 and 40 times more events than the current data. An exposure equivalent to 10 (40) times the current one might be within reach of IceCube-Gen2 [89] after about 5 (20) years of data taking. From that plot we also see that the region of parameters accessible to IceCube(-Gen2) is already disfavored by a combination of constraints arising from the direct (LHC) and indirect searches (flavor physics and LEP).

Conclusions
Several authors have inquired about the potential of IceCube detector in the South Pole to search for LQs [56][57][58][59][60][61][62]. In most of these works LQs couple to the first generation of quarks. On the one hand this allows LQ production to profit from the valence quark contribution to the neutrino-nucleon cross section. On the other hand limits from direct searches, atomic parity violation and flavor physics are at odds with these possibilities.
Here, instead, we focus on LQs that couple to heavy quark flavors, a scenario that recently became of particular interest in connection with B-meson semileptonic decays. The hope that IceCube could contribute to the search of this type of LQ models springs from a few considerations. First, these models cannot be dismissed by current direct or indirect detection searches. Second, their production cross section has been underestimated in the literature as the neutrino-gluon contribution has been neglected up to now. Third, attenuation effects when crossing the Earth are modified by the new LQ interactions which was JHEP06(2018)032 Figure 14. Estimation of the 95% CL exclusion regions in the parameter space (m ∆ , |g cτ L |) considering several projections of IceCube data at exposures measured in units of the current 6year data sample. We also show the constraints on large LQ couplings derived from the Z-pole observables (red line) and the lower limits on the LQ mass stemming from the reinterpretation of LHC limits on pair-produced LQs (dashed line), see section 2.2 for details. not considered previously. Finally, the proposal of IceCube-Gen2 [89] allows for exploring a much higher statistics of HESE than at the present time.
In this paper we focused on a simple LQ scenario that can explain R exp and involves only three new parameters: the LQ mass (m ∆ ) and two couplings (g cτ L and g bτ R ). Parameter space is constrained in such a way as to solve R D ( * ) and satisfy experimental constraints coming from the Z-decay modes. With parameters constrained in this way we estimate the present (future) sensitivity of IceCube (IceCube-Gen2) to this model. Our results are summarized in figure 14 from which we see that the parameter space accessible to IceCube/IceCube-Gen2 is already excluded by the direct searches (see figure 4) and the low-energy flavor physics observables. That conclusion is based on our current understanding of the HESE data and the value of the spectral index γ in particular. If that value turns out to be lower, then also figure 14 would change because there would be more data corresponding to heavier LQs but with the mass still lower than about 650 GeV, to comply with the results of the direct searches at LHC.
We should also stress that in our calculation we consistently used a democratic flavor distribution of the incoming flux, which is plausible and compatible with the current data.

JHEP06(2018)032
Other possibilities are, however, possible too which is one reason more why settling that issue is of major importance for searching the effects of physics beyond the SM [77].

Acknowledgments
This work was partially supported by Fundação de Amparoà Pesquisa do Estado de São Paulo (FAPESP) and Conselho Nacional de Ciência e Tecnologia (CNPq). R.Z.F. is grateful for the hospitality of the Laboratoire de Physique Théorique of Université Paris-Sud during the final part of this work and acknowledges the CNRS for support during this time. D.B. acknowledges the INP-CNRS support via"Inphyniti". This project has also received support from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement N • 690575 and 674896.

A Leptoquark contributions to the neutrino-nucleon scattering cross section
For completeness we give here the expressions for the neutrino-nucleon scattering cross section at the parton level. The total differential cross section induced by LQs can be written as for an incident neutrino ν τ . Here g(x,ŝ) (q(x,ŝ)) is the gluon (quark q) PDF of the nucleon. In our computaions, we considered the NNPDF2.3 PDF sets [84], which are implemented in LHAPDF6 [85].

A.1 ν -gluon contribution
We give here the expressions for the gluon-neutrino differential cross sections induced by LQs. The kinematics of these processes can be described by the following variables with the inelasticity y and the Bjorken variable x defined in the intervals 5 , and x ∈ (m q + m ∆ ) 2 2m N E ν , 1 ,

JHEP06(2018)032
where E ν and E ∆ are the (anti-)neutrino and LQ energies in the laboratory frame, Q 2 = −t is the invariant squared momentum transfered in the t-channel, m q is the mass of the quark in the process, m ∆ the LQ mass, and m N = (m p +m n )/2 is the average mass of the nucleon. We have used the fact that water is a isoscalar target, so that N = Z. We also define the function λ(a, b, c) ≡ (a 2 − (b − c) 2 )(a 2 − (b + c) 2 ).
The two diagrams that contribute to the amplitudes that enter these cross sections can be seen in figure 8. The cross section for ν g → q u ∆ (−2/3) is given by 6 ∂ 2 σ gν ∂x∂y = (2m p E ν x) α s |g qu L | 2 8ŝ where N c is the number of colors q u and are a generic up-type quark and lepton, and α s = α s (ŝ) is the strong coupling constant at the relevant scale of the process.

A.2 ν -quark contribution
The LQ state contributes to the neutrino-quark cross-section via an exchange of a LQ in the t-channel as depicted in figure 7. The SM cross-section ammended with the LQ contribution reads where q u , q u are generic up-type quarks and , are generic leptons. We define the weak couplings g u V = T 3 u −2Q u sin 2 θ W = 1/2−4/3 sin 2 θ W and g u A = T 3 u = 1/2. The Mandelstam variables in our setup are taken to bê s = (p ν + p qu ) 2 ,t = (p ν − p ν ) 2 andû = (p ν − p q u ) 2 . (A.4)

JHEP06(2018)032
Similarly, the νq cross-section receives a s-channel contribution from the LQ state, as shown in figure 7. The corresponding expression is given by where Γ ∆ is the total LQ decay width.

B Transport equations including the leptoquark constribution
The set of transport equations that have to be solved in order to compute the attenuation factor for the neutrino flux is the following where E y = E/(1 − y), m τ (m ∆ ) is the τ (∆) mass, τ τ (τ ∆ ) is the τ (∆) lifetime, σ ∆ ν is the total LQ contribution to the neutrino-nucleon cross section, dσ ∆−q ν /dy (dσ ∆−g ν /dy) is the neutrino-quark (neutrino-gluon) part of the differential cross section and σ CC,NC ν N are the SM CC and NC neutrino-nucleon cross sections. The distribution of the τ partial decay widths dΓ τ →ν X /dy can be found in refs. [81,86] while the distribution of the ∆ partial decay width dΓ ∆→ν X /dy is given by (B.4) where γ = E ∆ /m ∆ , Θ(x) is the Heavyside function.