Neutrino Telescopes as QCD Microscopes

We present state-of-the-art predictions for the ultra-high energy (UHE) neutrino-nucleus cross-sections in charged- and neutral-current scattering. The calculation is performed in the framework of collinear factorisation at NNLO, extended to include the resummation of small-$x$ BFKL effects. Further improvements are made by accounting for the free-nucleon PDF constraints provided by $D$-meson data from LHCb and assessing the impact of nuclear corrections and heavy-quark mass effects. The calculations presented here should play an important role in the interpretation of future data from neutrino telescopes such as IceCube and KM3NET, and highlight the opportunities that astroparticle experiments offer to study the strong interactions.


Introduction
Ultra-high energy (UHE) neutrinos represent a unique messenger for a wide variety of astrophysical and cosmological phenomena [1][2][3]. Being charge-neutral, neutrinos are not deflected by galactic magnetic fields, while being weakly interacting implies that they are not attenuated e.g. by dust. Therefore, neutrinos point directly to their original source of production, offering crucial information which would be not available from other messengers such as electromagnetic waves or cosmic rays. This potential is illustrated by the recent report of the detection of UHE neutrinos from the same direction as a known Blazar [4].
UHE neutrinos also represent a very promising probe of physics beyond the Standard Model [5], for instance testing scenarios where PeV neutrinos arise from heavy dark matter decays [6] or those where neutrinos have non-standard interactions (NSI) [7]. Moreover, neutrino-nucleon scattering at very high energies also provides unique opportunities to test Quantum Chromodynamics (QCD) in an extreme kinematic regime far from those accessible at colliders [8].
At very high energies, the astrophysical neutrino flux Φ ν decreases steeply with the neutrino energy E ν with a power-law of approximately E −2 ν [9]. The fact that the incoming UHE neutrino flux is so small, combined with the very feeble interactions of neutrinos with matter, makes the detection of UHE neutrinos extremely challenging. To bypass these limitations, neutrino telescopes such as IceCube [10], KM3NET [11], and Baikal [12] instrument large volumes of ice or water, which then act as detectors with an effective volume of up to 1 km 3 . This way, by accumulating data taken during several years, neutrino telescopes are expected to be able to detect neutrinos with energies as high as E ν 10 12 GeV. In addition, UHE neutrinos can also be studied by means of an array of radio antennas [13] or with balloon experiments such as ANITA [14].
Irrespective of the specific detection method adopted, the physical interpretation of the UHE neutrino event rates relies on an accurate theoretical prediction for the relevant signal process. In the case of the neutrino telescopes described above, the signal arises from the interaction [8] of a highly energetic neutrino with a nucleon from the target material, typically an H 2 O molecule. The prediction for this scattering process is sensitive to both perturbative and non-perturbative QCD effects including the quark and gluon content of the target nucleon [15,16], the treatment of heavy-quark mass effects [17], and the stability of the perturbative expansion at small values of the partonic momentum fraction x [18]. The sensitivity of the prediction to QCD dynamics represents both a challenge and an opportunity, given that the strong interactions must be probed in an unexplored kinematic regime. It is the main goal of this work to present a state-of-the-art calculation, taking into account recent developments in our understanding of perturbative and non-perturbative QCD, for the UHE neutrino-nucleon cross-section relevant for signal detection at neutrino telescopes. While a number of calculations have been provided by different groups in the past, both in the framework of collinear DGLAP factorisation [8,[19][20][21][22][23][24] and beyond it [25][26][27][28][29][30], there are now several significant motivations to revisit this calculation. The improvements made in this work, and their motivation, are detailed below.
Firstly, we account for the effects of small-x (BFKL) resummation up to the nextto-leading logarithm (NLLx) accuracy in the calculation of deep-inelastic scattering (DIS) structure functions and in the collinear evolution of parton distribution functions (PDFs). Specifically, the fixed-order calculations at next-to-next-to-leading order (NNLO) are consistently matched to the corresponding resummed results [31,32]. Small-x corrections to the DIS structure functions are provided in the framework of the FONLL generalmass scheme [17,33]. Recently, evidence for the onset of small-x resummation has been found [34,35] in the inclusive and charm HERA data following the approach proposed in Refs. [36,37]. It is thus of the utmost importance to include these effects also in the computation of the UHE cross-sections, which probe much smaller values of x than at HERA, down to x 10 −8 .
Secondly, we include the constraints on PDFs arising from the LHCb D-meson production measurements in pp collisions at 5, 7, and 13 TeV [38][39][40] following the approach of Ref. [41]. As demonstrated in this work, and in the previous studies [42][43][44], the LHCb data provides important constraints at small-x which are in turn relevant for the UHE crosssection predictions. Here we revisit the study of Ref. [41] to account for the impact that this data has on the PDFs extracted from a global analysis where small-x resummation effects have been included [34]. As will be shown, including the LHCb D-meson measurements leads to a significant reduction of the uncertainties on the UHE cross-sections associated to PDFs. It is worth mentioning that the same PDF constraints are also an important ingredient in the prediction of the so-called prompt neutrino flux [45], induced by the decays of D mesons produced in cosmic ray collisions in the atmosphere. The neutrinos produced in this process are a dominant background source in the search for UHE neutrinos of astrophysical origin. See Refs. [43,46,47] for more details as well as a number of related studies [48][49][50][51][52][53][54].
Finally, as the signal process of interest is neutrino scattering on an H 2 O target (as opposed to a free nucleon), we account for nuclear modifications in the UHE cross-section predictions. There have been a number of studies quantifying the impact that nuclear modifications have on the distribution of quarks and gluons within bound nuclei in Refs. [55][56][57][58], which can be used to study the impact that nuclear corrections have on neutrino-nucleon scattering. These modifications (and the associated uncertainty) are absent from many of benchmark UHE cross-section calculations.
Our results represent a significant improvement over previous calculations of the UHE cross-sections both in terms of the perturbative (NNLO accuracy, small-x resummation) and non-perturbative (PDF constraints from D-meson data, nuclear corrections) content. This allows us to provide a more reliable comparison to the very recent direct measurements of the neutrino-nucleus interaction cross-sections at high energies [59,60] extracted from the IceCube data, and to extrapolate the calculation to as of yet unexplored energies. Our calculations are made publicly available at the level of the total cross-section for a large range of (anti)neutrino values up to 10 12 GeV. In addition, we also provide predictions for the DIS structure functions in the form numerical grids. These grids (with the relevant interpolation routines) may be used to construct a double-differential cross-section and can also be integrated numerically to obtain the total cross-section.
The outline of this paper is as follows. In Sect. 2 we present the theoretical formalism adopted for the calculation of the cross-sections. In Sect. 3 we revisit the impact of the LHCb D-meson data on the small-x PDFs once the effects of small-x resummation are accounted for. The main results of this work are presented in Sect. 4, where we also discuss the various sources of theoretical uncertainties associated to our calculation. We conclude in Sect. 5 and outline possible future developments. Appendix A contains a tabulated version of our calculation for the UHE cross-sections with the corresponding uncertainties, and also includes a description of the structure function grids and how they may be used to construct UHE cross-section predictions at both differential and integrated levels.

Theoretical formalism
In this section we present an overview of the theoretical settings used throughout this work to compute the UHE neutrino-nucleus cross-sections. We first review the basic ingredients required to evaluate predictions for neutrino-induced DIS cross-sections, focussing on the kinematic region relevant to UHE neutrino-nucleus scattering. We then provide details on the theoretical accuracy and implementation of our calculation. Finally, we discuss the role played by heavy-quark mass effects and the treatment of nuclear corrections.

Neutrino-nucleon deep-inelastic scattering
As illustrated in Fig. 2.1, a high-energy neutrino interacts with a nucleon from the target through the exchange of an electroweak gauge-boson, either neutral (Z) or charged (W ± ). The energy of the incoming neutrino can then be reconstructed by the energy deposited in the detector by either the outgoing lepton, in charged-current (CC) scattering, or the recoiling hadronic system in the case of neutral-current (NC) scattering.
The formalism for describing both NC and CC processes is equivalent, and in the following we focus on the CC case, ν(k) + N (p) → ± (k ) + X(W ). Details on the implementation of the NC case are discussed in Appendix A. The differential cross-section for this process can be written in terms of DIS structure functions as follows: where the sign of the xF 3 term is positive (negative) for (anti)neutrino scattering, and N represents the struck nucleon. 1 In Eq. (2.1) we have defined Y ± = 1 ± (1 − y) 2 and the DIS kinematic variables are given by where some of the Lorentz-invariant quantities are also given in the laboratory centre-ofmass frame where the target nucleon is at rest. In these expressions, m N is the nucleon mass, E ν and E are the incoming-and outgoing-lepton energies, and k, q, and p are the four-momenta of the neutrino, the gauge boson, and the nucleon (see Fig. 2.1). The main ingredients for the theoretical prediction of the double-differential crosssection defined in Eq. (2.1) are the structure functions F ν(ν)N i=2,3,L (x, Q 2 ), which describe the underlying QCD dynamics of the scattering process. Structure functions factorise as follows: corresponding to the convolution of the universal PDFs (f a ) with the process-dependent coefficient functions (C i,a ). Typically, coefficient functions are computed in perturbation theory as a truncated expansion in powers of the strong coupling α s . Before discussing the theoretical accuracy of the structure function predictions, it is useful to illustrate the kinematic coverage of the UHE neutrino-nucleon cross-section for some representative values of the incoming neutrino energy E ν . To do so, we compute the total CC neutrino-nucleon cross-section integrating Eq. (2.1) over the relevant kinematic region. The computation is performed at NLO using the NLO fixed-order NNPDF3.1sx central PDF set [34] assuming an isoscalar target. 2 Fig. 2.2 shows the integration regions on the (log 10 x vs. log 10 Q 2 ) plane that contribute most (20%, 40%, 60%, and 80%) to the total cross-section for E ν = 5 × 10 8 GeV (left) and E ν = 5 × 10 10 GeV (right). The dot-dashed diagonal line indicates the upper bound of the integration of Q 2 as a function of x, that is given by Q 2 max (x) = 2m N E ν x. For E ν = 5 × 10 8 (5 × 10 10 ) GeV, the main contribution to the cross-section arises from the region with Q 2 M 2 W and x 2 × 10 −5 (2 × 10 −7 ), with sensitivity to the small-x region down to x 10 −6 (10 −8 ). This pattern is the consequence of two separate effects. For Q 2 M 2 W the contribution to the total cross-section is power-suppressed by the W -boson propagator, see Eq. (2.1). For Q 2 M 2 W , instead, the suppression of the cross-section is driven by the decrease of PDFs as Q 2 decreases.
In the kinematic regions highlighted in Fig. 2.2, fixed-order calculations are affected by the presence of large logaritms of x that spoil the perturbative convergence, see Refs. [31,32,61] and references therein. However, it has been shown that fixed-order calculations   The percentage of the total UHE neutrino-isoscalar cross-section corresponding to the contribution of different kinematic regions in the (x, Q 2 ) plane in CC scattering. Predictions are computed in the FONLL scheme at NLO using the NNPDF3.1sx NLO PDF set. Results are shown for E ν = 5 × 10 8 GeV (left) and E ν = 5 × 10 10 GeV (right). The dot-dashed diagonal lines correspond to the kinematic limit Q 2 can be complemented by small-x resummation effects to extend their validity down to very small values of x. Remarkably, this leads to a marked improvement in the description of the precise HERA 1+2 combined data [34,35]. It is therefore important to account for these effects when computing predictions for the UHE neutrino-nucleus cross-section. As discussed below, predictions for the DIS structure functions in this work are accurate to NNLO+NLLx.

Numerical implementation
The DIS structure functions appearing in Eq. (2.3) are computed at fixed-order with the APFEL program [62]. APFEL has been interfaced to the HELL code [31,32,61] that provides small-x resummation corrections to the DIS coefficient functions and the DGLAP splitting functions. The computation of the structure functions is performed with the NNPDF3.1sx sets obtained in the corresponding global PDF analysis [34]. Heavy-quark mass effects are included in the calculation of the structure functions using the FONLL general-mass variable flavour number scheme [17,33]. In the NC case, the FONLL-B and -C variants have been used for the NLO and NNLO computations, respectively. In the case of CC interactions, NNLO corrections to the massive coefficient functions have recently been computed [63] (see also Ref. [64]). However, the results of this calculation have not been made publicly available yet in a format suitable for the present calculation. Mass effects to the CC structure functions are therefore taken into account only up to NLO accuracy [65]. For a number of phenomenological applications, such as the calculation of event rates at neutrino telescopes using Monte Carlo simulations, the double-differential distribution in Eq. (2.1) is the most relevant physical quantity. To this end, predictions for both CC and NC structure functions differential in x and Q 2 are made publicly available in the form of interpolation grids (see Appendix A).
However, the total neutrino-nucleon cross-section is also a highly relevant quantity for a number of analyses at neutrino telescopes. To construct the theoretical predictions for this quantity, it is necessary to integrate Eq. (2.1) over the allowed kinematic range in x and Q 2 . The total cross-section is therefore computed as This integral is evaluated numerically using Gauss quadrature as implemented in GSL with relative precision set to 10 −4 . Results were cross-checked by comparing to those obtained by performing a two-dimensional Monte Carlo integration with the VEGAS algorithm (as implemented in the CUBA library [66]) and requiring a relative precision of 5 × 10 −3 . The integration limits of the integrals in Eq. (2.4) are obtained by inspection of the DIS kinematics in Eq. (2.2). The inelasticity y is bounded within the range of y ∈ [0, 1], and consequently the maximum value that Q 2 can take for a given neutrino projectile energy is given by Q 2 max = 2m N E ν . While the integration in Q 2 in Eq. (2.2) should extend all the way down to zero, it is in practice necessary to impose a cut so that the structure functions are evaluated in the region of validity of perturbation theory. We choose Q 2 min = Q 2 0 , where Q 0 is the scale at which the input PDFs are parameterised. In the case of NNPDF3.1sx, Q 0 = 1.64 GeV. As far as the integration over x is concerned, the upper integration bound is naturally set to one, while the lower bound for each particular value of Q 2 should be set to x 0 (Q 2 ) = Q 2 /(2m N E ν ). However, when E ν becomes very large and for values of Q 2 close to Q 2 min , x 0 can potentially become very small. Since we cannot access PDFs down to indefinitely small values of x, we need to impose a cut on the minimum value that x 0 can take. To this end, we replace the lower bound of the integration in x in Eq. (2.4) with x 0 (Q 2 ) = max Q 2 /(2m N E ν ), x min . For our nominal results we choose x min = 10 −9 .
In order to investigate the dependence of the results on the integration limits Q 2 min and x min , we have computed the total CC neutrino cross-section on an isoscalar target as a function of incoming neutrino energy E ν for different values of Q 2 min and x min . The results are shown in Fig. 2.3 where the cross-section is computed at NLO for different values of Q 2 min (left) and x min (right), presented as a ratio to that obtained with the nominal values of Q min = 1.64 GeV and x min = 10 −9 . Regarding the dependence on the value of x min , we find that the total cross-section becomes sensitive to the precise value of x min only for extremely large energies (E ν 10 10 GeV). For instance, at E ν = 10 12 GeV the total cross-section is reduced by around 15% if the integration is restricted to x min = 10 −7 as compared to the nominal value x min = 10 −9 . With respect to the choice of Q min , we find that the total cross-section receives non-negligible contributions from the region Q ∈ [1.64, 2.2] GeV for neutrino energies below E ν 5 × 10 3 GeV. For instance, raising -7 - Q min from the nominal value of 1.64 GeV up to 2.2 GeV leads to a suppression of the total cross-section by around 15% at E ν = 100 GeV. The predictions provided in this work can thus be considered reliable for neutrino energies E ν 5 × 10 3 GeV. Accurate predictions in the region E ν 5 × 10 3 GeV could be achieved by matching the perturbative QCD calculation of the structure functions to a parameterisation of low-Q 2 neutrino structure functions [67], following for instance the approach developed in Ref. [68].

Heavy-quark mass effects
Previous calculations of the UHE neutrino-nucleon cross-sections have been performed in the zero-mass variable-flavour-number scheme (ZM-VFNS), where finite mass effects in the computation of the partonic cross-sections are neglected. Nonetheless, it is now well known that heavy-quark mass effects are required in order to provide an adequate description of the HERA data in the low-Q 2 region. Since the computation of the total neutrino crosssections requires integrating over a wide range in Q 2 that includes the regions close to the charm-and bottom-quark masses, heavy-quark mass corrections can be relevant. As mentioned above, our predictions for the DIS structure functions are computed using the FONLL scheme that accounts for heavy-quark mass effects. In the following, we assess the importance of including these effects in the UHE cross-section predictions. We also discuss the relevance of top-quark mass effects.
We first consider the impact of charm-and bottom-quark mass effects on both the CC and NC calculations at the level of the single-differential cross-section dσ/dQ 2 . This is obtained by integrating Eq. (2.1) over the allowed kinematic range of x, namely Comparison between the FONLL and ZM-VFN heavy-quark mass schemes at NLO for the calculation of the single-differential cross-section dσ νI /dQ 2 , Eq. (2.5), for CC (left) and NC (right) scattering at E ν = 10 7 GeV. In both cases the limit m t → ∞ is taken to decouple the effects of the top quark.
display the absolute distributions while the lower panels show the ratio to the ZM-VFNS. In these plots the limit m t → ∞ is taken to decouple the effect of the top quark. 3 In the upper panel of these plots, the total cross-section obtained by integrating over the allowed range in Q 2 is also reported.
In the case of CC scattering, the inclusion of charm-mass effects introduces a positive shift of up to 5% in the region of Q 2 10 GeV 2 but has a negligible impact on the total cross-section. The impact of bottom-mass effects is entirely negligible. This is due to the strong suppression for the subprocesses involving the CKM matrix elements V ub and V cb . Heavy-quark mass corrections in NC case display a similar behaviour at low energies but introduce a negative correction in the intermediate energy range of Q 2 ∈ [20, 500] GeV 2 . The net effect is a reduction of the total cross-section by (1-2)%. While not shown here, similar results are also obtained for neutrino energies in the range E ν ∈ [10 3 , 10 12 ] GeV. From this comparison, we conclude that charm-and bottom-quark mass effects are negligibly small at the level of the total cross-section, but may still be relevant for the calculation of the double-differential distributions in Eq. (2.1). These effects are included in our predictions unless specifically stated.
In addition to charm-and bottom-quark mass effects, it is also important to consider the mass effects related to the top quark (see also Ref. [69]). The Born processes for massive top-quark production in neutrino-nucleon scattering are where in the CC process q d denotes a down-type quark. Kinematically, the CC (NC) process may only proceed if the physical threshold for the production of a top quark (pair) is reached, namely where the target nucleon mass is neglected. An accurate prediction for these processes relies on the inclusion of mass effects in the partonic cross-section. We account for these effects by adopting the FONLL scheme with a maximum of n f = 6 active flavours (n f = 6 scheme, in short). In this scheme the top quark is treated as a massive parton for energies close to the top-quark mass, and becomes effectively massless at much larger energies, Q 2 m 2 t . This has the consequence of introducing top-quark PDFs whose evolution resums potentially large logarithms of the form ln(Q 2 /m 2 t ). As discussed in Sect. 2.2, heavy-quark mass effects are accounted for at O α 2 s for NC scattering and O(α s ) for CC scattering. We note that the scheme choice in this work contrasts with the approach commonly used in the extraction of PDFs from experimental data, such as the NNPDF3.1sx sets used in this work, in which the n f = 5 scheme is preferred. In order to employ the n f = 6 scheme, we have generated variants of the NNPDF3.1sx sets in which the top-quark PDFs are dynamically generated at Q 2 = m 2 t [70]. These sets coincide with the original ones for Q 2 < m 2 t and are made publicly available (see Appendix A). In Fig. 2.5 we assess the impact of top-quark production on the total cross-section for both CC (left) and NC (right) scattering. To do so, cross-sections are computed at NLO in the m t → ∞ limit and compared to the default n f = 6 FONLL calculation with m t = 172.5 GeV. Taking the limit of m t → ∞ effectively increases the production threshold for top production such that it becomes kinematically inaccessible. Comparing these predictions is therefore useful to single out the relative contribution for the top-quark production to the total cross-section.
For CC scattering (left plot of Fig. 2.5), it is found that the top-quark contribution becomes relevant for neutrino energies of E ν 10 7 GeV, while for NC scattering (right plot) the contribution is negligible. Inspection of Eq. (2.7) reveals that the production of a single top-quark in CC scattering is kinematically possible for Q 2 /x m 2 t . Therefore, for neutrino energies of E ν 10 7 GeV top-quark production becomes accessible at relatively small values of Q 2 in the low-x region. This process is dominated by the partonic subprocess W + b → t as a consequence of the steep growth of the b-quark PDF (generated at Q = m b ) at low-x, and of the fact that |V tb | ≈ 1. The combination of these effects results in a nonnegligible contribution to the total CC cross-section which represent up to 5% of the total cross-section at E ν ∼ 10 12 GeV.

Impact of nuclear modifications
At neutrino telescopes, the incoming UHE neutrinos scatter upon nucleons which may be bound within a nucleus of the target material. Therefore, in principle it is necessary to account for the presence of nuclear-binding effects. The most relevant of these effects for the calculation of the UHE cross-sections is that of shadowing [71], namely the depletion of nuclear structure functions as compared to their free-nucleon counterparts. In experiments such as IceCube, KM3NET, and Baikal the nuclear target is a molecule of H 2 O, and the scattering is dominated by bound nucleons inside the oxygen nuclei with a small admixture of the free protons from the hydrogen. The mass-averaged structure function for a water/ice target may thus be written as where F p is the structure function of a free proton, and F p(n),A represents that of a proton (neutron) which is bound within a nucleus with mass number A containing Z protons and N neutrons, with A = Z + N . In this case, the bound nucleus is that of oxygen which is an isoscalar target (N = Z) with A = 16. Predictions for these structure functions can be computed according to Eq. (2.3), where the PDFs correspond to those of the (bound) nucleon. This implies that one should account for the effects that quark and gluon PDFs of a nucleon experience inside a heavy nucleus. These corrections have been quantified in a number of analyses of nuclear parton distribution functions (nPDFs) [55][56][57][58].
To assess the impact of nuclear corrections on the UHE neutrino-nucleon cross-sections, we have computed the mass-averaged total cross-section for an H 2 O molecule using the nPDFs from the EPPS16 global analysis [56]. These predictions are obtained by first computing the cross-section for both a free proton and an oxygen target, and then performing a combination according to The EPPS16 fit is constructed taking the CT14 NLO free-nucleon PDFs [72] as a baseline. Therefore, we use the central value of this PDF set for the free-nucleon predictions. 4 Results for both CC and NC scattering cross-sections (for the sum of neutrino-and antineutrinoinduced processes) are shown in Fig. 2.6. Distributions are presented as normalised to those of the corresponding free-nucleon predictions. The quoted uncertainty bands represent the 1σ uncertainty of the EPPS16 set (excluding free-nucleon uncertainties) evaluated using the asymmetric Hessian prescription.
The CC (left) and NC (right) neutrino-nucleon cross-sections (adding the contributions from neutrinos and antineutrinos) for an H 2 O molecule computed with the EPPS16 nPDF set and presented as a function of the neutrino projectile energy E ν(ν) . The quoted 1σ confidence-level (CL) uncertainty bands include only the uncertainty from the nuclear PDF fit, and have been evaluated using the asymmetric Hessian prescription. Each distribution has been normalised with respect to the baseline free-nucleon prediction.
We find a suppression of the cross-section for (anti)neutrino energies E ν(ν) 10 6 GeV due to the shadowing effect present in the nPDFs at small-x. The central value is reduced by 3% for E ν(ν) = 10 6 GeV, and by as much as 10% for E ν(ν) = 10 10 GeV. However, it should be noted that while a suppression of the cross-section is preferred, the uncertainty of the nuclear corrections are almost as large as the shift of the central value. Therefore, using EPPS16, the significance of nuclear modifications is mostly within the 1σ level. At lower energies, instead, the impact of the nuclear corrections becomes less important. The results of Fig. 2.6 indicate that nuclear corrections represent a large source of theoretical uncertainty in the predictions of the UHE neutrino-nucleon cross-section. Therefore, it is necessary to account for such effects to provide reliable predictions.
In Sect. 4, where predictions are provided for the total UHE cross-sections, the impact of nuclear corrections is accounted for in a factorised form. To illustrate this procedure, we consider here the construction of the cross-section for the neutrino-induced scattering on an oxygen nucleus. First, the nuclear modification factor R νO (E ν ) is computed with the EPPS16 nPDFs as follows: where σ free νI (E ν ) is the cross-section for an isoscalar target computed with the central CT14 NLO free-nucleon PDFs, and the normalisation is such that R νO (E ν ) → A = 16 in the limit of vanishing nuclear effects. Note that the flavour symmetry of PDFs at small-x implies that Eq. (2.10) gives essentially the same results for both NC and CC scattering, as also seen from Fig. 2

.6.
This modification factor is then applied to the cross-section for an isoscalar target computed with a different set of free-nucleon PDFs according to for either NC or CC scattering. The mass-averaged cross-section for neutrino scattering on a molecule of H 2 O is obtained according to Eq. (2.9). In fact, for the highest neutrino energies, any departures from non-isoscalarity can be ignored and the following approximation can be made The total uncertainty in Eq. (2.12) can then be computed by adding in quadrature the free-nucleon PDF uncertainties (arising from σ νI ) with those of the nPDFs (from R νO ). The factorised expression in Eq. (2.12) makes it straightforward to improve the prediction of the total cross-section when a more precise determination of R νO (E ν ) becomes available. 5 The large uncertainties associated to R νO (E ν ) in the current calculation can be related to the lack of experimental data used to determine nPDFs (and in particular those sensitive to the gluon) in the region of x 10 −2 . Indeed, to compute R νO (E ν ) at large E ν values it is necessary to extrapolate nPDFs to small values of x by several orders of magnitude. In the extrapolation region results are driven to a large extent by the particular methodological choices made to extract nPDFs. Examples are the parameterisation (i.e. the functional form assumed to parameterise the x and A dependence) or χ 2 tolerances that define the 1σ PDF uncertainties. 6 As will become apparent in Sect. 4, the uncertainty associated to the nuclear corrections is a limiting factor in the calculation of the UHE neutrino-nucleon cross-section. This provides a strong motivation to improve global fits of nPDFs by extending the kinematic coverage of the input data set. One possibility is to include LHC data collected in p+P b collisions which are sensitive to the small-x region. The nPDFs are typically parameterised as continuous functions of the nucleus mass number A. Therefore, constraints obtained for nucleons bound within a P b nucleus (A = 208) are relevant also for lighter nuclei such as oxygen. Progress in this direction may be possible by studying forward D-meson production in p+P b collisions [73], but these data have not yet been included in any nPDF fits (see Ref. [74] for initial work in this direction). In the longer term, stringent constraints would also be provided by possible future lepton-ion colliders such as the EIC [75] and the LHeC [76].

Constraining small-x resummed PDFs with D-meson data
The analysis of Ref. [41] quantified the impact of the LHCb D-meson cross-section measurements on the NNPDF3.0 NLO PDFs at small x. Here we revisit this analysis, applying it to the NNPDF3.1sx sets [34,77] which have been extracted either with or without including the effects of small-x resummation. In this section, we review the fit settings and discuss the experimental inputs along with the corresponding theoretical calculations used to include the LHCb data into the fit. We then present the fit results and describe the tests performed to assess their robustness.

Fit settings, experimental data, and theory calculations
The kinematic coverage of the UHE neutrino-nucleon cross-section as shown in Fig. 2.2 illustrates the sensitivity of this observable to PDFs in a region of extremely small x (x 10 −8 ). This region is outside the coverage of the input datasets used in the current PDF fits. In particular, it is beyond the coverage of the HERA data [78] that only reach values of x as small as x 2 × 10 −5 for Q 2 1 GeV 2 . However, recent work has demonstrated that smaller values of x can be probed using D-meson production measurements as provided by the LHCb experiment. Several groups have shown [42][43][44] that forward D-meson crosssection measurements can provide important information on both the normalisation and the shape of the small-x PDFs, especially that of the gluon. In particular, in Ref. [41] it has been shown that the combination of the LHCb measurements at √ s =5, 7, and 13 TeV allows to substantially reduce the uncertainties of the gluon PDF at small x.
In order to illustrate the interplay between the small-x gluon PDF and the UHE neutrino-nucleon cross-sections, in Fig. 3.1 we show the correlation coefficient ρ [79] evaluated between the gluon PDF at Q 2 = 4 GeV 2 and the CC neutrino-isoscalar cross-section σ CC νI (E ν ) for fixed values of E ν . The correlation is shown as a function of the momentum fraction x carried by the gluon. Note that although there is no direct coupling between the electroweak bosons and the gluon, this correlation mostly arises from the impact of the latter on the sea quarks by means of the DGLAP evolution. From this comparison, we find a large correlation (ρ 0.6) in the small-x region (x 10 −4 ) for the cross-section evaluated at large values of the neutrino energy E ν 5 × 10 8 GeV. Therefore, a better understanding of the small-x gluon PDF is necessary to obtain reliable predictions for the UHE neutrino-nucleon cross-sections in this energy range.
In this work we aim at providing predictions for the UHE neutrino-nucleon crosssections using a state-of-the-art calculation based on structure functions accurate at NNLO+NLLx. This requires using PDFs exatracted with the same accuracy. The main limitation in achieving this goal is that theoretical predictions for D-meson production are not available at this accuracy. Firstly, because NNLO differential cross-sections for heavy-quark pair production are only available for top quarks [80]. Secondly, because small-x resummed hard cross-sections [81] are not yet available in a form suitable for phenomenology.
Following Ref. [41], the impact of the LHCb D-meson production data on PDFs at small-x can be quantified by means of the following observables: where X = 5 or 7, p D T and y D are the transverse momentum and rapidity of the D mesons, and y D ref denotes a reference rapidity bin. Eq. (3.1) implies that we use D-meson differential distributions always in a normalised form, taking as a reference either a given rapidity bin or data taken at a different centre-of-mass energy. An important implication of this approach is that missing higher-order corrections to partonic cross-sections partially cancel in the ratios, reducing the sensitivity of these observables to such corrections. Crucially, numerators and denominators in Eq. (3.1) probe different regions in x, and therefore these observables still provide constraints on the PDFs at small x.
We apply Bayesian reweighting [82,83] to prior PDF sets which have previously been extracted with different theory settings. Specifically, we use the NNPDF3.1sx sets based on: NLO, small-x resummation matched to NLO (i.e. NLO+NLLx), and small-x resummation matched to NNLO (i.e. NNLO+NLLx). In all cases, the partonic cross-sections for Dmeson production are computed at NLO. We discuss below the stability of our results with respect to this choice.
For the PDF reweighting, we exploit the LHCb measurements of D-meson cross-section measurements in pp collisions at √ s = 5, 7, and 13 TeV [38][39][40]. These data are provided double differentially with respect to p D T and y D for different types of D meson. The kinematic coverage of these cross-sections is approximately p D T ∈ [0, 8] GeV and y D ∈ [2, 4.5], with small differences depending on the specific D-meson species and the hadronic centreof-mass energy. From these measurements we construct the two normalised observables defined in Eq. In what follows we adopt as a baseline the results of the fit with the normalised cross-sections N ij X , denoted by N 5+7+13 . The LHCb data is restricted to the region p D T ∈ [1,8] GeV. The reason for this cut is that in our calculation the factorisation and renormalisation scales are set equal to the transverse mass of the outgoing heavy quark. Since the NNPDF3.1sx fits are determined at the input scale of Q 0 = 1.64 GeV, restricting to the region p D T ∈ [1,8] GeV avoids sampling the PDFs outside their validity range. Practically, we set µ min = 1.64 GeV in the calculation which is relevant for a very small fraction of events with p c T 0 GeV and p D T > 1 GeV. As summarised in Table 3.1, with this kinematic cut, a total of 78, 72, and 119 data points at 5, 7, and 13 TeV, respectively, are included in the fit.
Applying the reweighting procedure requires computing the theory predictions for the observables in Eq. (3.1) using the N rep = 100 replicas of the NNPDF3.1sx sets. These predictions are obtained at NLO+PS accuracy using POWHEG [84][85][86][87] to match the fixed-order calculation [88] to the Pythia8 shower [89,90] with the default Monash 2013 tune [91]. The calculation is performed with an input value for the charm-quark mass of m c = 1.5 GeV.

Results and validation
Following the procedure outlined above, we have produced three variants of the NNPDF3.1sx+LHCb fits, based on the NLO, NLO+NLLx, and NNLO+NLLx theory settings. In all cases, the matrix-elements for the partonic process are evaluated at NLO [88]. In Fig. 3.2 we compare the gluon PDF from the three prior NNPDF3.1sx sets at Q 2 = 4 GeV 2 with the corresponding results after the LHCb D-meson cross-sections have been included in the fit. For completeness, we also compare the NLO results which have been obtained with the NNPDF3.0+LHCb set [41]. PDF uncertainties are computed as 1σ intervals.
From the comparisons in Fig. 3.2 we find a significant reduction of the PDF uncertainties due to the inclusion of the LHCb D-meson cross-section data. The magnitude of this reduction turns out to be similar for all the three theory settings considered. We also observe that at NLO, consistent results are obtained for the two different priors, NNPDF3.0 and NNPDF3.1sx. This stability is reassuring taking into account the differences between the two fits in terms of input dataset, treatment of the charm quark PDF, and the values of the heavy-quark masses.
In the lower-right plot of Fig. 3.2 we also display the comparison of the three NNPDF3.1sx+LHCb sets based on NLO, NLO+NLLx, and NNLO+NLLx theory. The stability of the perturbative expansion once small-x resummation effects are accounted for is evident. Indeed, the differences between the central values of the NLO+NLLx and NNLO+NLLx gluon PDFs are always smaller than the corresponding uncertainties in the entire range of x considered. We also observe that the central value of the small-x gluon PDF is larger in the NNLx case as compared to fixed-order and that the effective behaviour at small x is a moderate rise rather than a constant behaviour as at NLO.
In Table 3.1 we report the values of χ 2 /N dat for each of the LHCb D-meson datasets considered in this analysis. For each of the three theory settings considered, we show the results both before (χ 2 orig ) and after (χ 2 new ) adding the LHCb data into the fit. In the first column of this table, we also indicate the values of N dat for each dataset. The results for the N 5+7+13 combination, corresponding to N dat = 269 data points, represent the baseline fit of this work. For completeness, we also provide the χ 2 values for the ratio between the cross-sections at 13 and 5 TeV, R 13/5 . From Table 3.1 one finds that an excellent description is obtained for the normalised D-meson cross-sections, with similar values of the χ 2 /N dat for the three different theory settings, and with the NLO+NLLx fit leading to the smallest χ 2 new . As discussed in Ref. [41],  in the calculation of these χ 2 values the experimental bin-by-bin correlation matrices are included for R 13/5 , while for the normalised cross-section data these correlations (which are only available for N 5 and N 13 ) are not included. In this analysis the partonic cross-sections for charm-quark production, convoluted with NLO and (N)NLO+NLLx accurate PDFs, are accurate to NLO in all cases. In order to assess the robustness of these results with respect to this approximation, the analysis of the N 5+7+13 normalised cross-section dataset has been repeated in the following way. A modified χ 2 is introduced to quantify the agreement between theory and data, defined as where O i corresponds to the i-th bin value of the observable O and δO i to its uncertainty. Note that in our baseline analysis the theoretical uncertainty δO th i is not accounted for. This estimator is then used to repeat the NLO+NLLx analysis, and in this case we define  the theory error δO th i to be the following shift This additional source of theoretical uncertainty δO th i is introduced in the χ 2 in Eq. (3.2) in order to estimate the possible impact of the missing contributions in the evaluation of the hard cross-sections for charm production. It effectively reduces the contribution to the χ 2 due to those bins that are most sensitive to the difference between (N)NLO+NLLx and (N)NLO PDF evolution. Consequently, the weight of these bins is reduced in the reweighting procedure.
In Fig. 3.3 we show the same comparison as in the upper plots of Fig. 3.2, now normalised to the central value of the NNPDF3.1sx+LHCb baseline result and adding the results (indicated by a *) of the fits obtained using the modified definition for the χ 2 in Eq. (3.2). One finds that adding the additional theory uncertainty in the χ 2 leads to slightly larger PDF uncertainties together with a small positive shift of the central value, both at NLO+NLLx and at NNLO+NLLx. These results suggest that missing NLLx corrections in the prediction of the D-meson production cross-sections can be neglected as compared to the gluon PDF uncertainties, thus justifying the approximations used in this work.   The results of the reweighting analysis are affected by further theoretical uncertainties, such as the choice of renormalisation and factorisation scales (set equal in this analysis) and of the value of the charm-quark mass used in the calculation of Ref. [41]. For illustration purposes, we show here the impact of performing scale and charm-quark mass variations in the determination of the gluon PDF at small x. As a baseline we take the NNPDF3.0+LHCb set of Ref. [41] restricting the input data to within p T ∈ [1.0, 8.0] GeV to match the current analysis. This comparison is shown in Fig. 3 From the comparison in Fig. 3.4 one observes that in all cases the results of the fits varying either the scale or the value of m c are consistent with the baseline within uncertainties. In particular, the effects of charm-quark mass variations are much smaller than the PDF uncertainties, since they partially cancel out in the ratios in Eq.

The neutrino-nucleon cross-section at ultra-high energies
In this section we present the main results of this work, namely the predictions for the UHE (anti)neutrino-nucleon cross-section for CC and NC scattering. Firstly, the perturbative stability of our calculation is assessed by studying the convergence of results obtained with different theoretical accuracies. We then provide a comparison of our baseline predictions to previous calculations and discuss the origin of differences and similarities. Finally, we compare our predictions to recent measurements from IceCube [59] and assess the impact of nuclear corrections following the strategy outlined in Sect. 2.4.

Impact of theory settings and perturbative stability
In order to assess the perturbative stability of our calculation, we compute the total crosssection in the n f = 6 FONLL scheme at NLO, NLO+NLLx, and NNLO+NLLx accuracy using the corresponding NNPDF3.1sx+LHCb PDF sets presented in Sect. 3. The results are shown in Fig. 4.1, where the total cross-section for both the CC (left) and NC (right) scattering is shown for the sum of neutrino-and antineutrino-induced processes, presented as a function of the (anti)neutrino energy E ν(ν) . Predictions are normalised to the central value of those obtained with the NNPDF3.1sx NLO set (i.e. without the LHCb D-meson data) and the quoted uncertainties indicate the 1σ PDF uncertainty. For the NLO+NLLx predictions, we only show the central value, as the results are almost identical to those obtained at NNLO+NLLx, both in terms of central value and uncertainty. Here we assume an isoscalar target without nuclear modifications. Fig. 4.1 reveals that the impact of the LHCb D-meson data is significant in the region E ν(ν) 10 8 GeV. For instance, at E ν(ν) 10 12 GeV the PDF uncertainties decrease from about 30% to below 10%. The same qualitative behaviour is observed for both CC and NC processes. A further interesting observation is that the inclusion of the LHCb D-meson data leads to a suppression of the central value of the total UHE neutrino cross-section by around 10% at the highest energies. Although this is shown here only for the NLO case, the LHCb data have a similar impact when different theory settings, such as (N)NLO+NLLx, are adopted. This is a consequence of the impact of the LHCb D-meson data on the behaviour of the gluon PDF at small x. As shown in Fig. 3.2, the LHCb data leads to relatively lower values of the gluon PDF in the small-x region (irrespective of the theory settings).
Concerning the perturbative stability, Fig. 4.1 shows that the NNPDF3.1sx+LHCb predictions are consistent within the 1σ PDF uncertainties at all considered perturbative accuracies. The central values of the NLO+NLLx and NNLO+NLLx calculations are remarkably consistent, in agreement within 1% for both CC and NC scattering across the entire range of E ν(ν) values. This difference is negligible as compared to PDFs uncertainties and nuclear corrections (see Sect. 2.4).
The perturbative stability of the UHE cross-sections upon inclusion of small-x resummation effects is a direct consequence of the stability of the structure functions. To illustrate this point, in Fig. 4.2 we show the predictions for the F ν+ν 2 (x, Q 2 ) structure function for both CC (left) and NC (right) scattering on an isoscalar target at Q = 100 GeV, computed with the NNPDF3.1sx+LHCb sets for the three different theoretical settings: NLO, NLO+NLLx, and NNLO+NLLx. The central values of the NLO+NLLx and NNLO+NLLx structure functions differ typically by around (1-2)%, and 4% at most for NC scattering. Differences between the fixed-order and the resummed calculations are instead more pronounced and can be as large as 15% at x 10 −8 . We also note that differences between the corresponding input PDF sets for each of these three theoretical settings are typically larger (see Fig. 3.2). Nonetheless, predictions for physical observables are in better agreement due to the partial compensation between PDF evolution and DIS coefficient functions.
The agreement between the NLO+NLLx and NNLO+NLLx predictions shown in Figs. 4.1 and 4.2 demonstrates the excellent convergence of the perturbative expansion at small x after resummation effects are included. This indicates that missing higher-order (MHO) corrections are likely to be small as compared to other sources of theoretical un-certainty. Altogether, we find that the combination of the constraints from the LHCb D-meson data and the inclusion of NLLx small-x resummation leads to robust predictions for CC and NC neutrino-nucleon cross-section predictions up to the highest energies, with PDF uncertainties below 10% and negligible uncertainties due to MHO corrections. At this level of precision, other sources of theoretical uncertainty, such as nuclear corrections, cannot be neglected.

Comparison to previous calculations
As discussed in Sect. 1, predictions for the UHE neutrino-nucleus cross-sections based on a variety of the different theoretical setups have been provided in the past by a number of groups [8,[19][20][21][22][23][24][25][26][27][28][29][30]. In the following, we provide a comparison of our results (labelled as BGR18) to a number of calculations of the UHE cross-sections present in the literature. The BGR18 predictions shown here correspond to the NNLO+NLLx calculation in the n f = 6 FONLL scheme with the corresponding NNPDF3.1sx+LHCb set as an input PDF set. As nuclear corrections are absent from the selected benchmark calculations, for the purpose of comparison we do not include them in our predictions.
We begin by comparing our results to the calculations from Gandhi et al. (GQRS98) [8], Connolly, Thorne, and Waters (CTW11) [22], and Cooper-Sarkar, Mertsch, and Sarkar (CMS11) [23], all of which have been obtained in the framework of collinear factorisation. The GQRS98 and CTW11 predictions are obtained with LO-accurate coefficient functions with CTEQ4M [92] and MSTW08 [93] PDFs respectively, while the CMS11 calculation is performed at NLO with the HERAPDF1.5 PDFs [94]. The comparison for the sum of neutrino and antineutrino cross-sections is presented in Fig. 4.3 for both the CC (left) and NC (right) processes. For GQRS98 and CTW11 we only show the central values, while for BGR18 and CMS11 we also include the PDF uncertainties. In the case of the CMS11 predictions we have added the central values of the neutrino and antineutrino induced crosssections, and show the relative uncertainty of the neutrino induced process. The lighter and darker uncertainty bands of this prediction correspond to different prescriptions for estimating the PDF uncertainties. The darker band is obtained by excluding one particular member of the HERAPDF1.5 PDF set. This member dominates in the computation of the PDF uncertainty at small x and thus removing it results in a much smaller uncertainty band, see Ref. [23] for more details.
In general, there are marked differences between the BGR18 and the other calculations. For the case of CC scattering, we find that the BGR18 and GQRS98 predictions are in agreement in the region 10 5 GeV E ν(ν) 10 8 GeV, but differ significantly outside this range. At E ν(ν) 10 12 GeV the GQRS98 calculation is larger by around 30%, which corresponds to a deviation of more than 3σ in units of the PDF uncertainty of the BGR18 calculation. The origin of this difference can be understood by considering that the CTEQ4M PDF set used in the GQRS98 calculation must be extrapolated beyond its region of validity, x ∈ [10 −5 , 1]. Applying the extrapolation adopted in Ref. [8], we find the CTEQ4M gluon and quark PDFs evaluated at x 10 −7 and Q ∼ 100 GeV overshoot those of NNPDF3.1sx+LHCb by around 25%. While the CMS11 and CTW11 calculations are broadly consistent with one another, we find that at intermediate energies these predictions are approximately (8-10)% larger than BGR18. In the case of the CTW11 prediction, this can be partly attributed to the absence of the O(α s ) corrections to the coefficient functions, which are negative and amount to (4-5)% for E ν(ν) ∈ [10 4 , 10 8 ] GeV. The origin of the difference with respect to the CMS11 prediction (which includes these corrections) is instead likely to originate from the treatment of top-quark production. Inspection of Fig. 11 (left) of Ref. [23] suggests that the contribution to the total cross-section from b-quark initiated diagrams (top-quark production) amount to +10% at E ν = 10 6 GeV. This calculation is performed in the ZM-VFNS where heavy-quark mass effects are absent. As discussed in Sect. 2.3, our calculation of top-quark production for CC scattering includes heavy-quark mass correction to NLO. We find that at this energy the relative contribution of top-quark production is below 1% (see Fig. 2.5). At high (anti)neutrino energies, E ν(ν) 10 11 GeV, the three calculations are instead consistent within the BGR18 uncertainties, although the CTW11 central value is suppressed by around 10% as compared to BGR18.
In the NC scattering case, the GQRS98 calculation agrees with BGR18 at intermediate values of E ν(ν) but overshoots it at higher and lower energies. The same behaviour was observed for CC scattering and these differences can again be primarily attributed to the behaviour of the input PDFs (see the discussion above). We find reasonable agreement between BGR18 and the CMS11 and CTW11 results for NC scattering. For the highest values of neutrino energies, both of these predictions tend to undershoot the BGR18 predictions. For instance, the CTW11 predictions are suppressed by around a factor 20% at E ν 10 12 GeV as compared to the BGR18 calculation. This can be partly traced back to differences at the level of input PDFs. However, when the PDF uncertainties of the CTW11 are accounted for (see Ref. [22]), this behaviour is not significant.
The calculations displayed in Fig. 4.3 are all based on collinear factorisation. In order to assess the sensitivity of the UHE cross-sections to other QCD theoretical frameworks, in   Fig. (4.3), now for the sum of NC and CC neutrino-isoscalar cross-sections. In addition, the AIS15 calculation based upon the non-linear QCD (saturation) framework is included.  [26]. This calculation is based on the Balitsky-Kovchegov equation with running coupling, which incorporates non-linear effects from gluon recombination (saturation).
In the AIS15 framework, the steep growth of the PDFs at small x is tamed by nonlinear effects, leading to a suppression of the UHE neutrino-nucleon cross-section. Indeed, Fig. 4.4 shows that at the highest energies the AIS15 calculation is in general suppressed by about 25% as compared to the BGR18 result. However, for large neutrino energies the AIS15 prediction is affected by a large theoretical uncertainty arising from the limited information on some of the input parameters that enter the calculation. Given the small theory uncertainties of the BGR18 calculation, a possible suppression of the measured UHE cross-sections at high energies as compared to our predictions may indicate the onset of saturation effects.

Comparison to IceCube data and impact of nuclear corrections
To conclude the discussion of our results, in Fig. 4.5 we show predictions for the sum of neutrino and antineutrino scattering on different targets. We consider the case of a free isoscalar target and that of H 2 O molecule, where nuclear modifications have been evaluated using the EPPS16 nPDF set (see Sect. 2.4). Results are also shown for an "isoscalar" H 2 O molecule (N ≈ H 2 O), according to Eq. (2.12). These predictions are obtained from the baseline NNLO+NLLx accurate results in the n f = 6 FONLL scheme and are presented as a function of E ν for both CC (left) and NC (right) processes. The upper panels of Fig. 4.5 display the absolute cross-section, while in the lower panels the predictions are shown normalised to the central value of the BGR18 calculation on H 2 O. In the CC case, we also show the recent IceCube measurements based on the 6-year HESE (high-energy showers) dataset [59]. This dataset is based on high-energy starting events (or containedvertex events). The experimental uncertainties are the sum in quadrature of the statistical and systematic errors. Note that for the rightmost data point only a lower limit on the magnitude of the cross-section can be derived.  The results in Fig. 4.5 highlight that, given the current precision of free-nucleon calculations, effects due to the nuclear modifications of the PDFs of nucleons bound inside a H 2 O molecule are significant and cannot be neglected. As already discussed in Sect. 2.4, nuclear modifications induce a suppression of the UHE cross-section by up to 10% as compared to the free-nucleon case. However, uncertainties associated to these nuclear modifications are large and currently dominate over PDF and other theoretical uncertainties. The fact that nuclear modifications are now the dominant source of theoretical uncertainty on the UHE neutrino cross-section predictions provides a strong motivation to improve on the knowledge of nPDFs. This goal could be achieved by exploiting the LHC p+P b collision data as well as, in the near future, data collected at the EIC. It is worth noting that the prediction for N ≈ H 2 O provides an excellent approximation for E ν values. The difference between the full and approximate results at the lowest E ν values is ≈ 1%.
Focussing now on the CC scattering cross-section of Fig. 4.5, we observe that the IceCube measurements extend up to neutrino energies of around E ν 10 6 GeV. At the current level of precision, these measurements cannot discriminate between the different theoretical predictions. Nonetheless, future data based on a much larger sample from neutrino telescopes, such as IceCube and KM3NET, and from other experiments sensitive to very high-energy neutrinos should improve the precision of the current measurements at intermediate energies and extend the measurement to larger E ν values. As discussed in Sect. 4.2, the different theoretical predictions give rise to small differences at intermediate energies. The origin of these differences are well understood, and are related to the perturbative accuracy of the calculation. A cross-section measurement at higher energies would instead probe QCD in the very small-x region providing a test of the assumptions of non-perturbative information related to the distribution of quarks and gluons within bound nuclei, and may also test for the presence of saturation effects.
- 25 -In this work, we have presented state-of-the-art predictions for the cross-sections of highenergy neutrino scattering on nucleons, with particular attention to a target material composed of H 2 O molecules. With respect to previous calculations, we have made a number of major improvements.
Firstly, we have extended the calculation of deep-inelastic structure functions to NNLO matched with small-x resummation corrections at next-to-leading logarithmic (NLLx) accuracy. We find that small-x resummation corrections stabilise the perturbative convergence of the UHE neutrino cross-sections. This feature is highlighted by the similarity of the NLO+NLLx and NNLO+NLLx calculations.
The second main improvement is the inclusion of the D-meson production data from LHCb in the determination of the PDF sets used in our calculations. This dataset imposes a stringent constraint on the PDFs at small x, with a consequent significant impact on the UHE neutrino cross-section. This was achieved by producing dedicated PDF sets based on the NNPDF3.1sx sets and including the LHCb data by means of Bayesian reweighting. Due to the unavailability of small-x resummation corrections to the partonic cross-section for D-meson production, predictions for this process have instead been obtained at NLO+PS accuracy convoluted with PDFs determined including small-x resummation effects. We have demonstrated that the LHCb measurements lead to a reduction of the uncertainties on the UHE neutrino cross-sections related to PDFs by up to a factor three.
As compared to previous studies, we have introduced further improvements in our calculations. One of these is the inclusion of charm-, bottom-, and top-quark mass effects using the FONLL general-mass scheme for both CC and NC structure functions. We have also provided an estimate of the impact on the UHE cross-sections related to nuclear effects. We find that these corrections can be as large as 10% at high neutrino energies but are also affected by large uncertainties, which in turn impact the quoted precision of the calculation.
We found notable differences between our predictions and previous benchmark calculations, which can be traced back to differences in both the perturbative and non-perturbative inputs to the calculation. For example, in the energy range of E ν(ν) ∈ [6 × 10 3 , 10 5 ] GeV we find a prediction which is ≈ 5% lower than the quoted SM prediction in the IceCube measurement [59]. When more precise cross-section measurements are obtained at neutrino telescopes, these effects, as well as a reliable estimate of the nuclear corrections, will eventually become relevant to interpret the data. For the moment, the recent data from IceCube are still affected by large uncertainties and do not extend to large enough energies to discriminate between the different predictions.
We foresee that our calculations could be improved in at least two ways. Firstly, we found that uncertainties attributed to nuclear effects represent one of the dominant sources of theoretical uncertainty. This uncertainty could be reduced by including in the nPDF fits measurements from the LHC in p+P b collisions to constrain the distributions in smallx region. Secondly, small-x resummation effects and NNLO corrections, once available, should also be included in the partonic cross-section for D-meson production. Our study indicates that these corrections are likely to be small at the level of a normalised crosssection, but at the level of precision that the calculation has achieved it might be necessary to account for them.
To summarise, our analysis demonstrates how measurements of the neutrino-nucleus cross-section will represent a unique probe to test the strong interaction in an extreme regime where new dynamics are expected to arise, such as BFKL or non-linear effects. In this context, our calculations provide a robust building block for the data analysis and interpretation of neutrino telescopes in the coming years. Our results should also be relevant for other phenomenological applications, for instance to compute the attenuation of the high-energy neutrino flux as they pass through the Earth [95].

Acknowledgments
We are grateful to Aart Heijboer and Alfonso García for discussions about the KM3NET experiment, to Marco Bonvini for discussions regarding HELL, to Subir Sarkar for information about the CMS11 calculation, and to Javier Albacete and Alba Soto for providing us with the AIS15 results. We thank Hannu Paukkunen and Aleksander Kusina for providing us with the EPPS16 and nCTEQ15 Oxygen nPDF grids. We thank Rabah Abdul Khalek and Luca Rottoli for providing cross-checks relevant for some of the results presented in this paper. V. B. and J. R. are supported by an European Research Council Starting Grant "PDF4BSM". J. R. is also partially supported by the Netherlands Organization for Scientific Research (NWO). The research of R. G. is funded by the ERC Advanced Grant "MC@NNLO" (340983).

A The BGR18 UHE neutrino-nucleus cross-sections
In this appendix, we discuss the delivery of the results presented in this paper. We provide a tabulation of the total UHE cross-sections for a range of E ν values together with the corresponding theoretical uncertainties. Results are obtained with our baseline theory settings, i.e. the FONLL scheme at NNLO with a maximum of n f = 6 active flavours augmented by small-x resummation corrections to NLLx accuracy, and based on the NNPDF3.1sx+LHCb NNLO+NLLx PDF set.
Structure functions. In order to allow the users to reproduce our results for the doubledifferential cross-sections, we provide predictions for the neutrino structure functions F 2 , xF 3 , and F L both for NC and CC scattering. Structure functions are made available in the form of interpolation grids in the LHAPDF6 format [96]. Grids for NC and CC, neutrino and antineutrino structure functions are provided for a free isoscalar target. The grids contain the N rep = 40 Monte Carlo replicas resulting from the PDF reweighting analysis. 7 Mean and standard deviation of a structure function F are obtained according to Note that some of the replicas are equivalent, as a consequence of the unweighting procedure [83].
where F (k) is the value of the structure function computed with the k-th Monte Carlo replica.
Total cross-sections. The structure function grids discussed above allow constructing the double-differential UHE neutrino cross-section with the associated PDF uncertainties. Predictions for the total cross-section can then be obtained by integrating the doubledifferential cross-section by means of Eq. (2.4). While in Eq. (2.1) we presented the explicit formula for the CC cross-section, here we discuss in more detail the NC cross-section. Throughout this paper we have implicitly assumed that in our computations higher-order electroweak (EW) corrections can be neglected. However, for the computation of the NC cross-section we find advantageous to employ an "improved" scheme that includes part of the higher-order EW corrections. A pure LO treatment of the EW effects entails relations between the relevant parameters, such as G F , M W , M Z and sin 2 θ W , that lead some of them to be in strong disagreement with the measured values. In order to overcome this limitation, we employ the prescription of Ref. [97] to include in our computation the leading universal corrections, so that: where ρ = 1 + ∆ρ, with ∆ρ given in Eq. (8.22) of Ref. [97]. Note that, in contrast with the CC case, the structure functions in Eq. (A.2) for neutrino and antineutrino scattering are the same. Therefore, the only difference between neutrino and antineutrino cross-sections is the sign of the term proportional to xF 3 . In addition to the free-nucleon cross-sections, we also provide the values of the nuclear correction factor R νO /A as defined in Eq. (2.10). This allows one to evaluate the central value and the associated uncertainty to the total cross-sections including the effects of nuclear modifications. Presenting the results in this format has the advantage that predictions for scattering off a molecule of H 2 O can be easily updated once improved predictions for R νO become available, see Sect. 2.4. Since proton PDF uncertainties and nPDF uncertainties are uncorrelated, they can be combined by adding them in quadrature.
The BGR18 CC and NC (anti)neutrino total cross-sections, σ CC ν(ν)I and σ NC ν(ν)I , as functions of the energy E ν are tabulated in Tables A.1 and A.2, respectively. Results for the scattering off an isoscalar target without nuclear effects are shown along with the corresponding PDF uncertainty δσ CC ν(ν)I . The values for the nuclear correction factor R ν(ν)O /A, Eq. (2.10), and the corresponding nPDF uncertainty δR CC ν(ν)O are also provided. The LHAPDF structure function grids and an example code that computes the total cross-sections tabulated in Tables A.1