Neutrino Structure Functions from GeV to EeV Energies

The interpretation of present and future neutrino experiments requires accurate theoretical predictions for neutrino-nucleus scattering rates. Neutrino structure functions can be reliably evaluated in the deep-inelastic scattering regime within the perturbative QCD (pQCD) framework. At low momentum transfers ($Q^2 \le {\rm few}$ GeV$^2$), inelastic structure functions are however affected by large uncertainties which distort event rate predictions for neutrino energies $E_\nu$ up to the TeV scale. Here we present a determination of neutrino inelastic structure functions valid for the complete range of energies relevant for phenomenology, from the GeV region entering oscillation analyses to the multi-EeV region accessible at neutrino telescopes. Our NNSF$\nu$ approach combines a machine-learning parametrisation of experimental data with pQCD calculations based on state-of-the-art analyses of proton and nuclear parton distributions (PDFs). We compare our determination to other calculations, in particular to the popular Bodek-Yang model. We provide updated predictions for inclusive cross sections for a range of energies and target nuclei, including those relevant for LHC far-forward neutrino experiments such as FASER$\nu$, SND@LHC, and the Forward Physics Facility. The NNSF$\nu$ determination is made available as fast interpolation LHAPDF grids, and can be accessed both through an independent driver code and directly interfaced to neutrino event generators such as GENIE.

Depending on the neutrino energies E ν involved, different regimes are relevant in the corresponding theoretical calculations. In the sub-GeV region, the dominant interaction process is charged current quasielastic scattering (e.g.ν µ p → µ + n), and then as E ν is increased resonance scattering processes (e.g. ν µ n → µ + Λ − → µ + nπ − ) become the leading contribution. At higher energies and above the resonance region, starting at E ν = O(10) GeV and final state invariant masses of W ∼ > 2 GeV, inelastic scattering dominates (e.g.ν µ p → µ + X, with X being the hadronic final state). Inelastic neutrino scattering is further divided into shallow inelastic scattering (SIS) and deep-inelastic scattering (DIS) involving momentum transfers Q 2 below and above the threshold Q 2 ≃ few GeV 2 which separates the non-perturbative and perturbative regions, respectively. Other interaction processes, typically subdominant but relevant in specific phase space regions, include (in)elastic scattering off the photon field of nucleons, coherent scattering off the photon field of nuclei, and scattering on atomic electrons via the Glashow resonance. Theoretical models of neutrino scattering are implemented in various neutrinos event generators [19], such as GENIE [20,21] and its high-energy module HEDIS [22], GiBUU [23], and NuWro [24]. Most of these generators are tailored to specific energies and setups and cannot be straightforwardly applied to the whole range of experiments listed above.
Differential cross sections in inelastic neutrino-nucleon scattering are decomposed [25,26] in terms of structure functions, F νA i (x, Q 2 ), with x being the Bjorken variable and Q 2 the momentum transfer squared between the neutrino and the target nucleon. In the DIS regime, these structure functions can be expressed in the framework of perturbative QCD as the factorised convolution of parton distribution functions (PDFs) [27][28][29] and hard-scattering partonic cross sections. Their state-of-the-art calculation is based on PDFs and hard-scattering coefficient functions evaluated at next-to-next-to-leading order (NNLO) in the strong coupling α s expansion, with partial and exact [30][31][32] results also available one perturbative order higher (N 3 LO) and used in [33] to extract the proton PDFs. Furthermore, heavy quark (charm, bottom, and top) mass effects can be accounted for by means of general-mass variable-flavour-number (GM-VFN) schemes [34][35][36][37]. In addition, the applicability of fixed-order perturbative QCD calculations can be extended to the large (small) x kinematic region by means of all-order threshold [38] (BFKL [39,40]) resummation.
While DIS dominates inclusive neutrino-nucleon event rates for energies E ν ∼ > few TeV, at lower energies these rates receive a significant contribution from the SIS region. For instance, at E ν ∼ 100 GeV up to 20% of the inclusive cross section can arise from the Q ∼ < 2 GeV region [41]. Theoretical predictions of neutrino structure functions in the SIS region are therefore affected by much larger uncertainties than their DIS counterparts, given that they are sensitive to low momentum transfers, Q 2 ∼ < few GeV 2 [42], where the QCD perturbative and twist expansions break down and the factorisation theorems stop being applicable.
In order to bypass the limitations of perturbative QCD in the SIS region, phenomenological models of low-Q 2 neutrino structure functions have been developed and implemented in various neutrino event generators. One of the most popular is the Bodek-Yang (BY) model [43][44][45][46][47][48], based on structure functions with effective leading-order (LO) PDFs from the GRV98 analysis [49] with modified scaling variables and K-factors to approximate mass and higher-order QCD corrections. Drawbacks of the Bodek-Yang approach include the reliance on an obsolete set of PDFs that neglects constraints on the proton and nuclear structure obtained in the last 25 years, ignoring available higher-order QCD calculations, and the lack of a systematic estimate of the uncertainties associated to their predictions. Another restriction of the BY structure functions is that they cannot be consistently matched to calculations of high-energy neutrino scattering based on modern PDFs and higher-order QCD calculations [41,[50][51][52][53], introducing an unnecessary separation between the modelling of neutrino interactions for experiments sensitive to different energy regions.
Here we present the first determination of inelastic neutrino-nucleon structure functions which is valid for the full range of momentum transfers Q 2 relevant for phenomenology, from oscillation measurements involving multi-GeV neutrinos to ultra-high energy scattering experiments at EeV energies. This determination is based on the NNSFν approach, which combines a data-driven machine-learning parametrisation of inelastic structure functions at low and intermediate Q 2 matched to perturbative QCD calculations at larger Q 2 values. Following the NNPDF methodology [54][55][56][57][58][59][60][61], already applied to neutral-current structure functions in [62,63], neural networks are adopted as universal unbiased interpolants and trained on available accelerator data on neutrino structure functions and differential cross sections, with the Monte Carlo replica method to estimate and propagate uncertainties. The perturbative QCD calculations are provided by YADISM, a new framework for the evaluation of DIS structure functions from the EKO [64] family. Their input is the nNNPDF3.0 determination of nuclear PDFs [65], which reduces in the A = 1 limit to the proton NNPDF3.1 fit [55], specifically to the variant with LHCb D-meson data included [41]. This YADISM perturbative calculation is applicable up to EeV neutrino energies and accounts for exact top mass effects in charged-current scattering.
The NNSFν determination hence provides a parametrisation of the structure functions F νA i (x, Q 2 , A) (i = 2, 3, L) reliable for arbitrary values of the three inputs x, Q 2 , and A together with a comprehensive uncertainty estimate. Our strategy provides a robust estimate of all relevant sources of uncertainty, including those related to experimental errors, functional forms, (nuclear) PDFs, and missing higher order (MHO) uncertainties in the perturbative QCD calculation. We demonstrate the stability of NNSFν with respect to variations of both the input data and methodological settings, and show how it correctly inter-and extrapolates for A values not directly constrained in the fit. Upon integration of the structure functions over the kinematically allowed ranges in x and Q 2 , we obtain predictions for inclusive inelastic cross sections, together with the associated uncertainties, without restrictions on the values of the neutrino energy E ν and mass number A of the target nuclei.
We compare the NNSFν structure functions and inclusive cross sections to related calculations available in the literature, in particular to the Bodek-Yang (as implemented in GENIE), BGR18 [41], and CSMS11 [52] predictions. We assess the dependence of our results with respect to the values of E ν and A and provide dedicated predictions for the energy ranges and target materials relevant for the LHC far-forward neutrino scattering experiments, namely FASERν, SND@LHC, and the FPF, quantifying in each case the role played by the various sources of uncertainty. The NNSFν determination is made available as stand-alone fast interpolation grids in the LHAPDF [66] format which can be accessed either through an independent driver code and or directly interfaced to neutrino event generators such as GENIE. We also provide look-up tables for the inclusive cross sections and their uncertainties as a function of E ν and A.
The outline of this paper is the following. In Sect. 2 we review the theoretical formalism underlying neutrino-nucleus inelastic scattering on proton and nuclear targets, assess the perturbative stability of the QCD calculation, and compare different existing predictions among them. Sect. 3 describes the NNSFν approach to construct a data-driven parametrisation of the neutrino structure functions matched to pQCD calculations. Sect. 4 presents the NNSFν determination, verifies that it reproduces the input data and theory calculations, studies its robustness and stability, compares it with existing analyses, and evaluates the Gross-Llewellyn Smith sum rule. Predictions for inclusive neutrino cross sections are provided in Sect. 5, in particular we study their E ν and A dependence, the agreement with experimental data, and the sensitivity to kinematic cuts, and we conclude by outlining some possible future developments in Sect. 6.
Technical details are collected in three appendices. App. A describes the software framework underlying the NNSFν determination and how to access and use its results for the inelastic structure functions and cross sections. App. B outlines the main features of the YADISM code used to calculate neutrino structure functions in perturbative QCD and reports the outcome of representative benchmark comparison with APFEL. App. C summarises the impact of nuclear effects in neutrino inelastic scattering at the level of the parton distributions provided by nNNPDF3.0.

Neutrino inelastic structure functions
In this section we summarise the theoretical formalism underpinning the evaluation of neutrino-nucleon inelastic scattering cross sections in terms of structure functions and the calculation of the latter in the framework of perturbative QCD. We then present results neutrino structure functions evaluated with YADISM, quantify their perturbative stability and their dependence of the input PDFs, and compare them with the Bodek-Yang and BGR18 predictions.

DIS structure functions in perturbative QCD
The double-differential cross section for neutrino-nucleus scattering can be decomposed in terms of three independent structure functions F νA i (x, Q 2 ) with i = 1, 2, 3. Focusing on the charged-current (CC) scattering case mediated by the exchange of a W + weak boson, this differential cross section reads where s = 2m n E ν is the neutrino-nucleon center of mass energy squared, m n the nucleon mass, A the atomic mass number of the target nucleus, E ν the incoming neutrino energy, and the inelasticity y is defined as An analogous expression holds for antineutrino scattering, mediated now by the exchange of a W − weak boson, with the only difference being a sign change in front of the parity-violating structure function xF 3 , While the differential cross sections are a function of three kinematic variables, (x, Q 2 , y), the structure functions themselves depend only on x and Q 2 . Furthermore, both the cross sections and the structure functions depend on the atomic mass number A of the target nucleus only through the nuclear modifications of the free-nucleon structure functions. Kinematic considerations indicate that inelastic structure functions vanish in the elastic limit x → 1, that is, Alternatively, Eq. (2.1) can be expressed in terms of the longitudinal structure function F νA L (x, Q 2 ) defined by F L = F 2 − 2xF 1 , leading to where Y ± = 1 ± (1 − y) 2 and with the counterpart expression for anti-neutrino scattering, Expressing the differential cross section as in Eqns. (2.5)-(2.6) is advantageous because in the parton model (and in perturbative QCD at leading order) the longitudinal structure function vanishes, and hence F νA L (x, Q 2 ) ̸ = 0 starting only at NLO. The combination of neutrino and antineutrino measurements makes it possible to disentangle the different structure functions, for example the cross-section difference is proportional to the parity-violating structure function xF 3 averaged over neutrinos and antineutrinos. As discussed in the introduction, depending on the values of the momentum transfer squared Q 2 and of the hadronic final-state invariant mass W , (2.8) different processes contribute to these neutrino structure functions. In this work we consider only inelastic scattering, defined by the condition that the hadronic state invariant mass satisfies W 2 ∼ > 3.5 GeV 2 to avoid the resonance region. In the DIS regime, where Q 2 ∼ > few GeV 2 , neutrino structure functions can be evaluated in perturbative QCD in terms of a factorised convolution of process-dependent partonic scattering cross sections and of process-independent parton distribution functions, x z , Q 2 , i = 2, 3, L , (2.9) where j is an index that runs over all possible partonic initial states, C νN i,j is the process-dependent (but target-independent) coefficient function, and f (A) j indicates the PDFs of the average nucleon bounded into a nucleus with mass number A.
DIS coefficient functions can be expressed as a series expansion in powers of the strong coupling α s (Q 2 ), (2.10) The leading-order (k = 0) term in the coefficient function expansion Eq. (2.10) is independent of α s for i = 2, 3 since the Born scattering is mediated by the weak interaction. For massless quarks, charged-current neutrino DIS coefficient functions have been evaluated up to N3LO (third-order, m = 3) in [30,31]. For massive quarks, the calculation of strange-to-charm transitions with charm mass effects has been performed at NNLO (second-order, m = 2) in [67]. Mass effects can be incorporated in the massless calculation by means of a general-mass variable-flavour-number scheme [34][35][36][37]. For neutrino structure functions the expansion Eq. (2.10) displays good perturbative converge unless either Q 2 approaches the boundary of the non-perturbative region, Q 2 ≃ 1 GeV 2 , or the Bjorken-x variable becomes small enough to be sensitive to BFKL corrections. As well know, both PDFs and coefficient function are scheme-dependent objects and only their combination Eq. (2.9) is scheme-independent (up to higher orders). Each of the neutrino and anti-neutrino structure functions in Eq. (2.9) depends on a different combination of quark and antiquark PDFs, bringing in unique sensitivity to quark flavour separation in nucleons and nuclei. To illustrate this, if we consider a LO calculation on a proton target with n f = 4 active quark flavours, neglect heavy quark mass effects, and assume a diagonal CKM matrix, one can express the F νp 2 and xF νp 3 structure functions as where f q indicates the proton PDFs. The corresponding expressions for a neutron target and or isoscalar target are obtained from isospin symmetry, for instance the neutrino-neutron structure functions are expressed in terms of the proton PDFs as Different combinations of the neutrino structure functions in Eqns. (2.11) and (2.12) are sensitive to different PDF combinations. For instance, assume an isoscalar target and neglect nuclear corrections. In such scenario one has that expressed in terms of the valence PDF combinations, e.g. f u V = f u − fū, and hence the cross-section difference Eq. (2.7) yields (2.14) The same result is obtained if the target is not isoscalar but rather a purely hydrogen target, indicating how the difference between neutrino and antineutrino parity-violating structure functions xF 3 in Eq. (2.7) is a sensitive probe of the valence quark content in protons and nuclei. The neutrino structure function xF νA 3 must also satisfy the Gross-Llewellyn Smith (GLS) sum rule [68] calculable in perturbative QCD. For an isoscalar target the GLS sum rule is given by where n f is the number of active flavours at the scale Q 2 and the coefficients c k have been computed. The same expression holds for the anti-neutrino counterpart. The leading-order contribution to Eq. (2.16) follows from the partonic decomposition of the isoscalar xF νA 3 in terms of the valence quark PDFs, Eq. (2.13). In this work we do not impose the GLS sum rule in the data-driven fit and instead verify a posteriori that it is satisfied within uncertainties in the region of applicability of perturbative QCD. We note that experimentally one cannot access the x → 0 region, and hence the evaluation of Eq. (2.16) depends on the modelling of the small-x extrapolation region for the neutrino structure functions.

PDF dependence and perturbative stability
Here we study the PDF dependence and perturbative stability of neutrino DIS structure functions. We focus on the x region relevant for scatterings involving neutrino energies of E ν ∼ < 1 TeV and momentum transfers of Q ∼ > 2 GeV. Recall that from DIS kinematics, for given values of Q 2 and E ν the Bjorken-x variable satisfies (2.17) and hence it suffices to consider x ∼ > 10 −3 . We compare the following structure function calculations: • YADISM. The YADISM package described in App. B evaluates DIS charged-lepton and neutrino inclusive and heavy quark structure functions up to NNLO, and up to N 3 LO whenever available. YADISM has been benchmarked with APFEL [69] finding good agreement. Heavy quark mass effects are implemented in various schemes, including the zero-mass variable-flavour-number (ZM-VFN) scheme, the fixed-flavour-number (FFN) scheme, and the FONLL general-mass variable-flavournumber scheme [34,70].
For the purpose of the benchmarking comparisons shown in this section, YADISM inclusive neutrino structure functions are evaluated using either LO, NLO, and NNLO coefficient functions and in all cases NNPDF4.0 NNLO as input PDF set, with FONLL to account for heavy quark mass effects at NLO accuracy. We denote these calculations as YADISM-LO, YADISM-NLO, and YADISM-NNLO respectively in the following.
As will be discussed in Sect. 3, for the NNSFν determination of neutrino structure functions the perturbative QCD baseline calculation from YADISM will instead use NLO coefficient functions, nNNPDF3.0 NLO as input PDF sets for all targets including hydrogen, and a 5FNS where top quark mass effects are accounted for exactly while charm and bottom mass effects are neglected.
• BGR18. This is the calculation of neutrino DIS structure functions first presented in [41] in the context of predictions for UHE neutrino-nucleus cross sections and then updated in [22] when evaluating attenuation rates for UHE neutrinos propagating within Earth matter. BGR18 is based on APFEL with either NNPDF3.1 [55] or NNPDF3.1+LHCb [71] as input PDF sets. Nuclear corrections from the nNNPDF2.0 determination [72] were used in [22] to account for deviations with respect to the free-nucleon calculation. Several variants of the BGR18 calculation are available, both at fixed-order QCD (NLO and NNLO) and with BFKL resummation (NLO+NLLx and NNLO+NLLx).
In this work we consider the variant of BGR18 based on NLO coefficient functions and NNPDF3.1 NLO as input PDF, as implemented in the HEDIS module of GENIE. Since NNPDF3.1 and NNPDF4.0 are consistent within PDF uncertainties, one expects agreement between BGR18 and YADISM-NLO.
• Bodek-Yang. The BY calculation [43][44][45][46][47][48] is a phenomenological model for inelastic neutrino-and electron-nucleon scattering cross sections based on effective leading order PDFs which can be applied from intermediate-Q 2 in the DIS region down to the photo-production region Q 2 ≃ 0. The starting point is the GRV98 LO PDF set [49] evaluated at a modified scaling variable ξ w replacing the standard Bjorken-x, which allows BY to be extended to the low-Q 2 non-perturbative region. Several phenomenological corrections are applied to approximate NLO QCD and nuclear effects. In the following, we display the BY structure functions as implemented in the GENIE event generator.
• LO-SF. This calculation is defined by the LO expression of neutrino DIS structure functions on a proton target, Eq. (2.11), with PDFs accessed directly from the LHAPDF interface [66]. We consider two variants. First, LO-SF-NNPDF4.0, which uses NNPDF4.0 NNLO as input and that should coincide with the YADISM-LO calculation. Second, LO-SF-GRV98, which adopts GRV98LO as input PDF set and that should reproduce the Bodek-Yang calculation in the large Q 2 limit where its scaling variable ξ w reduces to x. The LO-SF-NNPDF4.0 and LO-SF-GRV98 calculations are only meant for benchmarking purposes and will not be used beyond this section.
The settings of the inelastic structure functions calculations that we just described are summarised in  functions, the software tool used for its evaluation, the treatment of heavy quark mass effects and its region of applicability. Mass effects and higher-order QCD corrections are approximated in the Bodek-Yang calculation by means of phenomenological model parameters.
Perturbative stability. Figs. 2.1 and 2.2 display the YADISM-LO, YADISM-NLO, and YADISM-NNLO calculations of F νp 2 , xF νp 3 , and F νp L and their antineutrino counterparts as a function of x for Q = 2 GeV and 10 GeV respectively. As indicated in Table 2.1, in all cases the common PDF set NNPDF4.0 NNLO is used. We display both the absolute structure functions and their ratios to the NLO calculation, and focus on the x region relevant for DIS structure functions with E ν ∼ < 1 TeV. It is found that for Q = 2 GeV higher-order QCD corrections are in general significant and exhibit a similar pattern both for neutrinos and for antineutrinos. For the dominant F 2 structure function, the LO calculation underestimates at large-x the NNLO result for up to 25%, while for x ∼ < 0.2 it overestimates it by more than 20%. NLO corrections reduce these differences at medium and large-x, but for x ≃ 10 −3 the NLO structure functions still overshoot the NNLO result by 20%. Since PDF uncertainties in F 2 are at the few percent level at most (see below), at Q = 2 GeV missing higher order uncertainties (MHOU) are the main source of theory errors. Concerning the parity-violating structure function xF 3 , for neutrino beams the LO calculation overestimates the NNLO result by 50% at x = 0.1 and by a factor 2 for x = 10 −2 . A similar pattern is observed for antineutrinos, with now LO becoming more negative at small-x while NLO is similar to NNLO. The longitudinal structure function F L vanishes at LO and displays large NNLO corrections, up to +40% as compared to the NLO result.
Once we increase the scale to Q = 10 GeV, the perturbative expansion exhibits an improved convergence, and in particular differences between NLO and NNLO structure functions (for a fixed common PDF set) are moderate in all cases. The only exception is F L at large-x, a region anyway not relevant for phenomenology since F 2 is much larger there. Nevertheless, neglecting NLO and NNLO coefficients functions still leads to sizable differences, specially for F 2 , of up to 10% at small-x and −20% at large-x.
Benchmarking. Fig. 2.3 compares the neutrino structure functions on a proton target in the YADISM-LO, Bodek-Yang, LO-SF-NNPDF4.0, and LO-SF-GRV98 calculations as a function of x for Q = 2 GeV and Q = 10 GeV and then as a function of Q (in the perturbative region) for x = 0.0126 and x = 0.25. In the case of LO-SF-NNPDF4.0, we show the 68% CL uncertainties, evaluated over N rep = 100 Monte Carlo replicas.
As expected, YADISM-LO coincides with the central value of LO-SF-NNPDF4.0 for all values of Q. Residual differences are found only at large-x and small-Q and are explained in terms of the target mass corrections (TMCs) accounted for in the YADISM calculation. Likewise, Bodek-Yang reduces to LO-SF-GRV98 at large-Q, and in particular at Q ≃ 10 GeV the two calculations are almost identical. This agreement   indicates that the phenomenological corrections to the GRV98 LO PDFs in the Bodek-Yang model have a negligible effect for Q ∼ > 10 GeV, while they become instead important at low-Q. For instance, at Q = 2 GeV, the differences between Bodek-Yang and LO-SF-GRV98 range between 10% and 30% depending on the structure function and the value of x.
The comparisons of Fig. 2.3 also highlight the typical behaviours of the neutrino structure functions in different regions of the (x, Q 2 ) plane. Concerning the x dependence, that of F νp 2 and Fν p 2 is similar and displays the valence peak at x ≃ 0.3 followed by the rise at small-x driven by DGLAP evolution, which is more marked the higher the value of Q. In the case of xF 3 , being a non-singlet structure function, the valence peak is followed by a small-x behaviour that depends on the value of Q and of whether the beam is composed by neutrinos or anti-neutrinos, since in each case the quark flavour combinations, Eq. (2.11), are different. In terms of the Q dependence, whether higher Q values lead to an increase or a decrease of the structure function depends on the value of x. For F νp 2 and Fν p 2 the structure functions grow (decrease) with Q for x = 0.0126 (x = 0.25). For xF 3 also for x = 0.25 there is a similar decrease with Q, while for x = 0.0126 the Q dependence varies strongly with the choice of input PDF set.
Given that in Fig. 2  differences between Bodek-Yang and YADISM-LO as Q is increased can only be attributed to those at the input PDFs level, GRV98LO and NNPDF4.0 NNLO respectively. For instance, for F νp 2 the GRV98LO calculation undershoots the YADISM-LO one based on NNPDF4.0 by 40% at x ≃ 10 −2 and Q = 10 GeV and by 10% at x = 0.25 and Q = 6 GeV. In the case of xF 3 , the small-x behaviour is qualitatively different between GRV98 and NNPDF4.0, for example at Q = 10 GeV for xF νp 3 the former predicts a steep rise while a flat extrapolation is preferred by the latter. This indicates that predictions based on the Bodek-Yang model, and thus on the obsolete GRV98LO PDF set, will in general disagree with those based on modern PDF determinations.
PDF dependence. We compare in Fig. 2.4 the YADISM-NLO predictions with those from the Bodek-Yang and BGR18 calculations. The YADISM-NLO and BGR18 predictions are very similar, consistent with the agreement within uncertainties of the underlying NNPDF4.0 and NNPDF3.1 PDF fits respectively. Differences between YADISM-NNLO and Bodek-Yang are significant, specially for F 2 in the x ∼ < 0.1 region and for xF 3 at small-x. These differences are explained by the reliance of BY on the obsolete GRV98LO PDF set (cfr Fig. 2.3) and due to neglecting higher-order QCD corrections (cfr Figs. 2.1 and 2.2).
The longitudinal structure function. In Fig. 2.3 we display the F 2 and xF 3 structure functions which provide the dominant contribution to the double-differential cross section Eq. (2.5). While the longitudinal    structure function F L vanishes at LO, it becomes non-zero at NLO and in specific kinematic regions can lead to non-negligible contributions to the scattering cross sections. To illustrate this hierarchy in the relative magnitude of the different structure functions, Fig. 2.5 displays their ratio to F 2 for neutrinos and antineutrinos for Q = 2 GeV and Q = 10 GeV in the YADISM-NNLO calculation.
From Fig. 2.5 we observe how in the large-x valence region the F 2 and xF 3 structure functions are of comparable magnitude, with F L being much smaller. Since xF 3 is a valence structure function, it is suppressed as x decreases and indeed for x = 10 −2 it becomes at most 20% of the value of the dominant F 2 . The relative contribution from F L is similar for neutrinos and neutrinos, since it is dominated by the gluon contribution, and becomes more important as both x and Q 2 decrease. In particular, for x ∼ < 10 −2 the magnitude of F L becomes larger than that of xF 3 . At Q = 2 GeV, F L can be up to 30% the value of F 2 , indicating a contribution to the double differential cross section larger than the typical experimental uncertainties and that hence must be accounted for.
Whenever possible, we fit data for the double-differential cross section rather than for the individual F 2 and xF 3 structure functions separately, since the former provides also sensitivity to F L .   Fν p 2 (x, Q) Bodek-Yang BGR18 YADISM-NLO      neural network parametrisation of neutrino structure functions, how it is trained on both the data and the QCD predictions, and the uncertainty estimate based on the Monte Carlo replica method.

General strategy
An schematic representation of the NNSFν strategy to determine neutrino structure functions is displayed in Fig. 3.1 with the Bodek-Yang predictions at x = 0.0126 for illustration. The (x, Q) plane is divided into three disjoint regions, with complementary methods to evaluate the structure functions in each of them: • Region I. At low momentum transfers Q ∼ < Q dat , with Q dat ≃ 5 GeV, the perturbative calculation of neutrino structure functions in Eq. (2.9) is either invalid or affected by significant theory uncertainties related to higher twists, missing higher perturbative orders, and large-x resummation effects.
In this region we parametrise the structure functions in terms of the information provided by the available experimental data on neutrino-nucleus inelastic scattering summarised in Sect. 3.2. Following the NNPDF fitting methodology, this parametrisation combines neural networks as universal unbiased interpolants with the Monte Carlo replica method for the uncertainty estimate.
• Region II. The region of intermediate momentum transfers, Q dat ∼ < Q ∼ < Q thr with Q thr ≃ 25 GeV, is well described by the perturbative QCD formalism. DIS structure functions are computed at NLO by YADISM with nNNPDF3.0 as input for all targets. In this region the neural network parametrisation is fitted to these QCD predictions rather than to the data as in Region I. The figure of merit is given in terms of the theory covariance matrix with PDF and the MHO uncertainties, see Sect. 3.4.
The nNNPDF3.0 determination already includes information from neutrino measurements, in particular from CHORUS (inclusive) and NuTeV (charm) structure functions, and therefore in Region II no neutrino data needs to be explicitely used to further constrain the parametrisation.
• Region III. For large momentum transfers, Q ∼ > Q thr , the neural network predictions are replaced by the direct outcome of the same YADISM calculation used to constrain the fit in Region II. Hence, in Region III the central prediction and uncertainties of NNSFν coincide with the YADISM ones, which extend up to Q = 10 TeV and down to x = 10 −9 to cover the entire kinematic region relevant for neutrino phenomenology including the UHE scattering.
Furthermore, Region III is extended in the small-x region with x ≤ 10 −5 to cover also momentum transfers of Q min ≤ Q ≤ Q thr , with Q min = 2 GeV. The reason for this choice is that for x ≤ 10 −5   Fig. 3.1. We indicate the publication reference, the range of x and Q 2 , the observables included, the scattering target, the final state measured, and the number of data points before (after) applying kinematic cuts. The coverage of these datasets in the (x, Q 2 ) plane is displayed in Fig. 3.2. Some of experiments provide additional observables, which are excluded to prevent double counting.
the neural network extrapolation trained in Regions I and II exhibits large uncertainties and using the QCD calculation is preferred on theoretical grounds [73].
As we will show in the subsequent sections, such a strategy allows us to consistently extend the state-ofthe-art perturbative QCD computations into the non-perturbative region and provide predictions that are valid across a wide range of energy relevant for neutrino phenomenology.
We have verified that the NNSFν determination is stable with respect to moderate variations of the values of the Q dat and Q thr hyperparameters. App. A provides additional details on the implementation of the NNSFν procedure including the prescriptions to evaluate and match inelastic structure functions in the various regions of the (x, Q 2 ) kinematic plane.

Experimental data
The parametrisation of neutrino structure functions applicable in Regions I and II defined in Fig. 3.1 requires two different inputs: experimental data in Region I and the corresponding QCD calculations for Region II. For the former, we consider all available data on inelastic neutrino structure functions and double differential cross sections. We restrict our analysis to those measurements where the incoming neutrino energy E ν is sufficiently large to ensure that the contribution from the inelastic region dominates. For this reason, we do not consider neutrino measurements from experiments such as ArgoNeuT [74], MicroBooNE [75], T2K [76], or MINERνA [77], where E ν is too low to cleanly access inelastic scattering.
Two kinematic cuts are applied to the data used as input to the NNSFν fit. First, a cut in the invariant mass of the final hadronic state W 2 ≥ 3.5 GeV 2 filters away points in the quasi-elastic and resonant scattering regions. Second, data points with Q ≥ Q dat are excluded according to the definition of Region I in Fig. 3.1. This cut does not result on a net information loss in the fit, since as mentioned above nNNPDF3.0 already includes the constraints from neutrino DIS data present in the Q ≥ Q dat region. Table 3.1 lists the datasets used to constrain neutrino structure functions in Region I. For each dataset, we indicate the publication reference, the range of x and Q 2 covered, the observables included, the scattering target, the final state measured, and the number of available data points n dat before and after applying kinematic cuts. In total we have 6224 (4184) data points in the fit before (after) cuts. The corresponding kinematic coverage of these datasets in the (x, Q 2 ) plane is displayed in Fig. 3.2, which also indicates the regions excluded by cuts. Some of the datasets from Table 3.1 provide additional observables on top of those indicated there, which however need to be excluded from the fit to prevent double counting. In particular, the same measurement is often presented in terms of both the differential cross section d 2 σ/dxdQ 2 and of the individual structure functions F 2 and xF 3 . We always select observables that are closer to the actual measurements, in this case the double differential cross-sections as they also constrain the longitudinal structure function F L .   Table 3.1. The region covered in grey is excluded from the fit from W 2 ≥ 3.5 GeV 2 cut required to isolate inelastic scattering, while the one in light red is excluded from the Q ≤ Q dat condition that defines Region I in Fig. 3.1. Table 3.1 and Fig. 3.2 indicate that the NNSFν fit is sensitive to inelastic neutrino structure functions for momentum transfers down to Q ≃ 400 MeV, well in the non-perturbative region. The range of x covered reaches x min = 0.015, and as a consequence of the DIS kinematics the values of Q being probed increase with x. Measurements in the non-perturbative region with Q ∼ < 1 GeV are provided by several experiments and cover momentum fractions up to x ∼ < 0.3. The datasets with the largest number of points (CHORUS, NuTeV, and CDHSW) present their measurements in terms of the double differential cross section d 2 σ/dxdy, while BEBCWA59, CCFR, and CHARM only provide data for separate structure functions F 2 and xF 3 .
Concerning nuclear effects, Table 3.1 shows that available data have sensitivity to neutrino scattering on Ne (A = 20), Fe (A = 56), and Pb (A = 208) targets. For CaCO3, the target used in the CHARM experiment, we assume A = 20 as the average atomic mass number of the nuclei that form this compound. Neutrino structure function measurements are not available on hydrogen or deuteron targets, and hence in Region I the low-A behaviour is extrapolated from the measurements with A ≥ 20. The inter-and extrapolation to A values not included in the fit is provided by the smoothness of the neural network output, as we validate in Sect. 4.
For momentum transfers Q > Q dat (Regions II and III) the dependence with the atomic mass number A of NNSFν follows that provided by nNNPDF3.0, and is hence constrained by other types of processes beyond neutrino DIS, such as charged-lepton fixed-target nuclear DIS and weak boson, dijet, and D-meson production cross sections in proton-lead collisions at the LHC.

Structure function parametrisation
The NNSFν parametrisation of neutrino structure functions in Regions I and II is obtained by training a machine learning model to experimental data and to the QCD predictions, respectively. It follows the NNPDF fitting methodology based on the combination of neural networks as universal unbiased interpolator with the Monte Carlo replica method for error estimate and propagation. This methodology was originally developed for DIS neutral-current structure functions [62,63] and subsequently extended to proton PDFs [54,56,[59][60][61]84], helicity PDFs [85,86], nuclear PDFs [65,72,87], and fragmentation functions [88,89].
Here we apply for the first time the NNPDF approach to i) the determination of neutrino (charged- Flowchart summarizing the NNSFν fitting framework. The three inputs are the momentum fraction x, the momentum transfer squared Q 2 , and the atomic mass number A, suitably preprocessed to lie in a common range. The network output is supplemented by a power-like term to facilitate the learning of the small-x region and the endpoint behaviour at x = 1 is subtracted to reproduce the elastic limit. The architecture of this neural network is 3-70-55-40-20-20-6 and hence composed by five hidden layers. The parameters of the neural network are optimised to reproduce the experimental data in Region I and the QCD predictions in Region II. current) structure functions and to ii) the parametrisation of a three-dimensional function, with neural networks receiving (x, Q 2 , A) as inputs. This is achieved by means of a stand-alone open-source NNSFν framework described in App. A. This framework shares many similarities with the NNPDF4.0 codebase, in particular it is also built upon TensorFlow [90] and uses adaptative stochastic gradient descent (SGD) methods for the minimisation such as Adam [91]. The double-differential neutrino-nucleus cross sections, Eqns. (2.1) and (2.3), are expressed in terms of three independent structure functions and hence one needs to parametrise six independent quantities which we choose to be each of them function of the three inputs (x, Q 2 , A). In Eq. (3.1) the limit A = 1 is to be understood as that of the structure functions on an isoscalar, free-nucleon target, rather than those for a proton target. That is, F ν 2 (x, Q 2 , A = 1) coincides with F νp 2 (x, Q 2 ) + F νn 2 (x, Q 2 ) /2, and likewise for the other structure functions.
The mapping between the inputs (x, Q 2 , A) and the outputs F ν i , Fν i with i = 2, 3, L in Eq. (3.1) is provided by a deep neural network as illustrated in Fig. 3.3. The free parameters of this neural network, its weights and thresholds, are determined by training the parametrisation to the experimental data (in Region I) and to the QCD predictions (in Region II) for neutrino structure functions and differential cross sections. Hyperparameters like the network architecture are determined by means of a dedicated optimisation procedure. The best choice for the architecture of this network is found to be 3-70-55-40-20-20-6, hence composed by five hidden layers with 70, 55, 40, 20, and 20 neurons in each of them.
As indicated by Fig. 3.3, two corrections are applied to the network output before it can be identified with the neutrino structure functions. First, we supplement the network output with a preprocessing factor x 1−α i that facilitates the learning and extrapolation of structure functions in the small-x region [73]. Second, we subtract the endpoint behaviour at x = 1 to reproduce the elastic limit where structure functions vanish due to kinematic constraints. With these considerations, the relation between the network output and the structure functions is given by

2)
dat . The structure functions parametrised by neural networks according to Eq. (3.2) and Fig. 3.3 are combined to construct the observables that input (O Inp ) the fit (and composed by either experimental data or QCD calculations) e.g. by means of Eqns. (2.5) and (2.6). The predicted observables are compared to the corresponding data points to evaluate the χ 2 entering the optimisation process and the crossvalidation stopping. χ 2 tr and χ 2 vl represent the training and validation losses, respectively.
for i = 1, 2, 3 corresponding to F 2 , xF 3 , F L , and where NN j indicates the activation state of the j-th neuron in the output layer of the network. By construction, structure functions parametrised this way vanish in the elastic limit x = 1 for all values of Q 2 and A without restricting the behaviour in the x < 1 region. The small-x preprocessing exponents (α i , α i ) in Eq. (3.2) are constrained from the data as part of the training procedure at the same time as the neural network parameters. Furthermore, the neural network inputs are rescaled to a common range, logarithmically in x and linearly in Q 2 and A to ensure that no specific kinematic region is arbitrarily privileged by the training. We do not enforce the positivity of the F 2 and F L structure functions, since it is found that the data together with the QCD constraints included are sufficient to avoid the unphysical negative region.

Fitting and error propagation
The neural network parametrisation of neutrino structure functions from Fig. 3.3 is constrained by the experimental data from Table 3.1 in Region I and by the QCD predictions in Region II. Here we discuss the figures of merit used in the minimisation, the error propagation strategy based on the Monte Carlo replica method, and the settings and performance of the training procedure. Fig. 3.4 provides a diagrammatic representation illustrating the evaluation of the figure of merit in the fit, χ 2 , as a function of the kinematic inputs z i = x i , Q 2 i , A i with i labelling the fitted data points. The structure functions parametrised by neural networks according to Eq. (3.2) and Fig. 3.3 are combined to construct the observables that input the fit e.g. by means of Eqns. (2.5) and (2.6). These fit inputs are classified into experimental data and QCD calculations depending on the range of Q 2 being considered. The predicted observables from the neural network parametrisation are then compared to the corresponding input data points to evaluate the χ 2 entering both the optimisation process and the cross-validation stopping.
Constraints from experimental data. In the Monte Carlo replica method one first generates a sample of N rep artificial replicas of the experimental data as follows. Given n dat experimental measurements of neutrino structure functions or differential cross sections characterised by central value F i,α , and n norm normalisation uncertainties σ (norm) i,n , artificial replicas of these measurements are generated by i,α , and r (k) i,n indicate univariate Gaussian random numbers generated such that experimental correlations between systematic and normalisation errors are accounted for. This procedure is formally equivalent to generating data replicas according to the statistical model provided by the experimental covariance matrix, which is reproduced by averaging over replicas, in the large replica limit, N rep → ∞. For the neutrino structure function data considered in NNSFν, cov exp contains only additive systematic uncertainties and hence it is not affected by the D'Agostini bias which would require introducing a t 0 covariance matrix [92] based on a previous iteration of the fit. Subsequently, for each of the N rep replicas generated according to Eq. (3.3), a separate neural network with the structure of Fig. 3.3 is trained by minimising the error function defined as in terms of the experimental covariance matrix. For each replica, the training is stopped once the crossvalidation stopping criterion described below is satisfied. The overall goodness-of-fit between the model predictions and the experimental data is then quantified by the χ 2 defined in a similar manner as Eq. (3.5) now in terms of the average neural network prediction, with averages over replica sample are evaluated as The N rep trained neutral network parametrisations provide a representation of the probability density in the space of neutrino structure functions, from which expectation values and other statistical estimators can be computed. For instance, the 1σ uncertainty δF in a structure function at generic (x, Q 2 , A) values can be computed by evaluating the standard deviation over the replica ensemble, where F represents any of the three structure functions F 2 , xF 3 , and F L . Similar considerations apply to other statistical estimators such as correlations coefficients, higher moments, and confidence level intervals. Both the experimental covariance matrix and the associated correlation matrix given by are displayed in the left panels of Fig. 3.5 for the data points listed in Table 3.1 after cuts. The matrices are block-diagonal since the different experiments are uncorrelated among them. In most cases, neutrino inelastic scattering experiments are limited by the correlated systematic uncertainties rather than by the statistical errors.
QCD constraints. In Region II, the same neural network parameterisation is trained now on QCD predictions based on YADISM and nNNPDF3.0 rather than to the experimental data. Taking the central nNNPDF3.0 set as baseline, we generate n th = 100 × n sf structure function QCD predictions distributed in x, Q 2 , and A covering Region II, with n sf = 6 denoting the number of independent structure functions being parametrised. As indicated in Eq. (3.10) we include a QCD boundary condition for A = 1, understood as an isoscalar free-nucleon target. By analogy with Eq. (3.3), also in Region II one generates N rep (independent) Monte Carlo replicas, this time starting from the central QCD predictions rather than from the experimental data. If we denote these central QCD structure functions by F where cov ⋆ th indicates the transpose of the Cholesky decomposition of the theory covariance matrix and r (th) j , r (th)(k) j represent stochastic noise sampled from a standard normal distribution. In order to treat the QCD data on the same footing as the experimental data, it is necessary to add two levels of stochastic noise to the central predictions F (qcd)(0) i : first to account the statistical fluctuations around the true underlying law, and second to generate the Monte Carlo replicas themselves. In the language of the closure test formalism [59], the first and second lines of Eq. (3.11) correspond to level-1 and level-2 pseudodata generation, while F (qcd)(0) i is the level-0 underlying law.
The theory covariance matrix entering Eq. (3.11) is constructed [93,94] as the sum in quadrature of the contributions from the MHO and PDF uncertainties, where the MHOU contribution is evaluated using the NLO scheme-B with the 9-point scale variation prescription as implemented in YADISM, and the PDF contribution is evaluated from the N rep = 200 replicas of the nNNPDF3.0 determination. In constructing the theory covariance matrix Eq. (3.12), the correlations between different nuclear targets are neglected.
In Region II the figure of merit used for the neural network training should be, instead of Eq. (3.13), with goodness-of-fit after the training of all replicas quantified now by the counterpart of Eq. (3.6), (3.14) Analogously to the experimental case, we can define the theoretical correlation matrix as which is displayed, together with the corresponding theory covariance matrix, in the right panels of Fig. 3.5, where QCD data points are sorted in increasing order of x, Q 2 , and A. Recall that there is no kinematic overlap between the entries of the experimental (Region I) and the theoretical (Region II) correlation matrices. One difference between ρ (th) and ρ (exp) is that in the former case QCD induces correlations between all data points (for a given value of A), arising because the same underlying nPDF determination is used as well as due to the correlation of factorisation scale variations entering DGLAP evolution.
Minimisation. Adding up the contributions from the experimental data in Region I and from the QCD predictions in Region II, the total figure of merit used for the training of the neural network parametrisation of neutrino structure functions is given by In minimising Eq. (3.16), it is paramount to achieve a balanced description of the two regions, that is, the average error function in both regions should be similar. Furthermore, one expects ⟨E exp ⟩ rep ∼ ⟨E th ⟩ rep ∼ 2 in the absence of tensions or inconsistencies in the data [63]. A consequence of this requirement is that the fit quality to the data should not be distorted by the inclusion of the QCD constraints, and hence in a fit variant using E exp as figure or merit the description of the experimental data should be comparable to that in the fits based on Eq. (3.16). We will demonstrate this stability of the NNSFν procedure in Sect. 4.
For the overall goodness-of-fit the total χ 2 is evaluated as which again for a balanced fit should satisfy χ 2 exp ∼ χ 2 th ∼ 1. In the following we will only quote the values of the experimental contribution χ 2 exp , since we verify that a comparable fit quality is obtained for the QCD component of the error function.
Eqns. (3.16) and (3.17) can also be understood as adding to the experimental χ 2 an extra contribution in the form of a Lagrange multiplier that enforces an external theory constraint, in this case that the neural network extrapolation in Q 2 reproduces the QCD prediction. Such Lagrange multiplier method is commonly used in the NNPDF framework to account for theory constraints, such as positivity and integrability in NNPDF4.0 and the A = 1 free-nucleon boundary condition in nNNPDF3.0. One benefit of the approach adopted here is that the theory covariance matrix Eq. (3.12) provides automatically the appropriate normalisation for the Lagrange multiplier contribution.
The minimisation of Eq. (3.16) is carried out by means of the adaptive SGD methods available in Tensor-Flow, specifically by Adam. The choice of optimisation algorithm, as well as that of other hyperparameters defining the methodology, has been determined by inspection of the fit results and performance. Table 3.2 lists the values of the hyperparameter configuration used in the baseline NNSFν determination, from the network architecture to the learning rate and the settings of the cross-validation stopping described below. A gradient clipping procedure is used to normalize gradient tensors such that their L 2 -norm is less than or equal to the clipnorm value. Training fraction 0.75 Once all N rep neural network replicas have been trained, a post-fit selection procedure is carried out to filter out eventual outliers associated to e.g. minimisation inefficiencies. Specifically, we define (and remove) an outlier replica as that exhibiting E val values 4σ away from the mean of the associated distributions. We also filter out replicas for which E (k) val never reaches below a threshold of 3.
Stopping criterion. To avoid overfitting, the same early-stopping cross-validation algorithm as used in NNDPF4.0 determination is adopted. For each dataset, 75% of the datapoints are randomly sampled to generate a training dataset while the other 25% of datapoints constitute the validation set. The SGD minimisation algorithm is trained of E tr while simultaneously monitoring E val . The optimal state of the neural network corresponds to the training iteration at which the value of E val has the lowest value, and once training has ended the state of the model is reverted to this point before producing the final outputs. Training can end when one of two conditions is met, whichever comes earliest; either the validation χ 2 has not improved for a given number of steps (known as stopping patience), or a threshold number of training epochs is reached. See Table 3.2 for the values of the associated hyperparameters.
Performance. For the baseline hyperparameter configuration summarised in Table 3.2, the training of a NNSFν replica takes on average around 10 hours. While the number of fitted data points n dat is of the same order to that of the NNPDF4.0 analysis, the fit takes a factor 25 more to converge, see Table 3.4 of [60], despite the network output being directly compared to the data without the intermediate requirement of the FK-table convolution present in PDF fits. This behaviour has a two-fold explanation. First of all, one is now exploring a three-dimensional parameter space in (x, Q 2 , A), as compared to the one-dimensional space relevant for a PDF determination where only the x dependence is constrained by the data. Second, the fit needs to reproduce not only the experimental data but also the QCD constraints which impose boundary conditions also in a three-dimensional (x, Q 2 , A) space.

NNSFν structure functions
Here we present the main results of this work, the NNSFν determination of inelastic neutrino structure functions valid for momentum transfers Q 2 in Regions I and II as defined in Fig. 3.1. The corresponding implications for inclusive neutrino scattering cross-sections are then presented in Sect. 5, while the matching procedure of the NNSFν outcome with the YADISM QCD calculations appropriate for Region III is described in App. A.
First of all we assess the quality of the fit to the experimental data on neutrino structure functions and compare NNSFν with representative measurements. Then we study the dependence with x and Q 2 of the NNSFν determination, and in particular demonstrate that in Region II it correctly reproduces the QCD predictions defining the theoretical boundary condition of the fit. We also compare the NNSFν determination with the Bodek-Yang and BGR18 calculations. We present a number of alternative NNSFν fits with dataset or methodology variations in order to assess the stability of our results. Finally, we study the implications of our analysis for the Gross-Llewellyn Smith sum rule and verify the agreement with the perturbative QCD expectations.

Fit quality and performance
Here we assess the fit quality to the experimental data and the QCD boundary conditions, quantify the fit performance including the small-x preprocessing, and provide representative comparisons between the NNSFν predictions and some of the fitted observables.
Fit quality. Table 4.1 we display the values of the experimental χ 2 exp per data point, defined by Eq. (3.6), for the individual datasets entering the NNSFν determination as well as for the total dataset. For each dataset we indicate the nuclear target, the fitted observables, the number of data points before and after kinematic cuts, and the values of the χ 2 exp . The latter are also provided by a fit variant (labelled as "wo QCD") where only the experimental data (in Region I), but not the QCD predictions (in Region II), are included in the fitted error function. The datasets in brackets are not included in the baseline fit to avoid double counting. The quoted χ 2 exp values are instead evaluated a posteriori using the outcome of the baseline NNSFν fit.
From Table 4.1 one finds that in the baseline fit a good description of the experimental data is obtained, with a total of χ 2 exp = 1.287 per data point for the 4089 points of considered. The fit quality is in general similar among the input datasets, without major outliers. This feature holds specially true for the three datasets that contribute the most in terms of statistical weight in the fit, namely CDHSW, CHORUS, and NuTeV, which are also the ones provided in terms of the cleaner double-differential cross-sections. In addition, a balanced descriptions of neutrino and antineutrino data is obtained whenever measurements for the two initial states are provided. A somewhat worse fit quality is obtained for the F 2 data from BEBCWA59 and CCFR, which in any case carry less weight in the fit and are potentially affected by the model-dependent separation from the measured cross-section. For the CHORUS observables, we find χ 2 exp = 1.185 and 0.797 per data point for the neutrino and antineutrino data respectively in the baseline fit. These values can be compared with the corresponding ones of 1.12 and 1.06 obtained in the NNPDF4.0 NNLO fit, where more stringent kinematic cuts are applied resulting in n dat = 832 data points as compared to 1580 in the NNSFν analysis.
As discussed in Sect. 3, in the context of a matched analysis such as NNSFν it is crucial to achieve a balanced description of the experimental data (in Region I) and the QCD predictions (Region II) during the fit. In this respect, we have verified that in the baseline fit the contribution to the total χ 2 tot arising from the QCD constraints, χ 2 th in Eq. (3.17), is of similar size as that for the experimental component χ 2 exp as expected for a balanced training.
Furthermore, by comparing the last two columns of Table 4.1, one can assess the impact of accounting for the QCD structure function constraints in the fit quality to the experimental data. As expected, the total fit quality is improved, χ 2 exp = 1.187 down from χ 2 exp = 1.287 in the baseline, in the fit where only  , for the individual datasets entering the NNSFν determination in Region I as well as for the total dataset. The datasets in brackets are not included in the baseline to avoid double counting, and their χ 2 values are evaluated a posteriori using the outcome of the NNSFν fit. For each dataset we indicate the nuclear target, type of observable, number of data points before and after kinematic cuts, and the resulting values of χ 2 exp . The latter are also provided by the fit variant where only the experimental data (Region I), but not the QCD predictions (Region II), are included in the fit and labelled as "wo QCD". the experimental data enters the figure of merit, E tot = E th instead of Eq. (3.16). Indeed, a better (or comparable) description of the data is generically expected once theoretical constraints are removed from the figure of merit, given that the functional form space available for F i (x, Q 2 , A) becomes less restricted. Nevertheless, this improvement remains moderate confirming that the addition of the QCD constraints in Region II does not significantly distort the description of the experimental data in Region I. Furthermore, this improvement in the χ 2 is homogeneously spread among the input datasets, rather than being associated to specific ones. We hereby conclude that the NNSFν fit is dominated in Region I by the experimental data constraints, with QCD boundary conditions providing a smooth transition to Region II. Table 4.1 also indicates the values of the χ 2 exp corresponding to datasets not included in the baseline fit, but rather evaluated a posteriori using the NNSFν predictions. Specifically, we list the values of the χ 2 exp for the individual structure functions from CDHSW, CHORUS, NuTeV experiments. These structure function datasets are not considered in the baseline fit since they would overlap with the corresponding reduced crosssections. The NNSFν predictions for the separate F 2 and xF 3 data excluded from the fit lead to a poor χ 2 , indicating a potential internal inconsistency between the reduced cross-section and separate structure function data. In Sect. 4.3 we investigate this issue by assessing the stability of the NNSFν fit results when the differential cross-section data is actually replaced by these separate structure function measurements.    Table 3.2), demonstrating that in all cases convergence is reached in a way that satisfies the cross-validation stopping criterion. Furthermore, the distribution is approximately Gaussian, with a mean of N ep ∼ 10 5 epochs, and does not present any significant long tails. These two considerations point to a stable training which is evenly distributed among the replica distribution. Fig. 4.2 displays the posterior probability distributions associated to the preprocessing exponents α i ,ᾱ i as defined in Eq. (3.2) for each of the fitted neutrino and antineutrino structure functions. Recall that these exponents are part of the fitted parameters and their values are restricted to the interval α i ∈ [0, 2]. As in Fig. 4.1, the distribution is sampled from the N rep = 200 replicas entering NNSFν. For each structure function, the distributions for the small-x neutrino and antineutrino exponents, α i andᾱ i respectively, turn out to be similar, indicating that the small-x behaviour of NNSFν depends only mildly on the neutrino flavour. The fact that despite α i andᾱ i being fitted separately, the same distributions are obtained, is another indication of the fit stability, given that QCD predicts that asymmetries between neutrino and antineutrino structure functions are washed out in the small-x region. The three structure functions also prefer similar small-x preprocessing exponents, with median values of α 1 ∼ 1.25, α 2 ∼ 1.1, and α 3 ∼ 1.05 for F 2 , xF 3 , and F L respectively, and in agreement within uncertainties. The distributions of α i ,ᾱ i are non Gaussian and exhibit a skewed tail towards smaller values of the exponents.

Comparison with data and previous calculations
Here we compare the NNSFν results with other calculations of neutrino structure functions, in particular with BGR18 and Bodek-Yang, as well as with the YADISM predictions based on nNNPDF3.0 which enter the fit as the QCD boundary condition. We then show the agreement between NNSFν and representative datasets used in the fit. Subsequently, we study the NNSFν uncertainties and their dependence with respect to variations in the x, Q 2 and A inputs of the parametrisation.
Comparison with Bodek-Yang and BGR18. First, we compare the NNSFν predictions with those from the Bodek-Yang and BGR18 calculations described in Sect. 2.2. In both cases, we access their predictions by means of their implementation in the GENIE event generator. For BGR18, we consider the variant based on NLO coefficient functions and NNPDF3.1 NLO as input PDF. We also display the YADISM predictions based on nNNPDF3.0 which enter the fit as the QCD boundary condition in Region II. The corresponding comparisons at the inclusive neutrino cross-section level will be presented in Sect. 5. We consider predictions for an isoscalar free-nucleon target, and also in Sect. 5 we will compare results for other nuclear targets for inclusive cross-sections. Figs. 4.3 and 4.4 display the NNSFν predictions as function of x for Q 2 = 4 GeV 2 (Region I) and Q 2 = 100 GeV 2 (Region II) and then as a function of Q for x = 0.0126 and x = 0.25, respectively. We display the F 2 and xF 3 structure functions for neutrinos, antineutrinos, and for their sum for an isoscalar free-nucleon target ( 2 H). The error band on the NNSFν predictions indicates the 68% confidence level intervals evaluated over the N rep = 200 Monte Carlo replicas. We also consider the central values of the Bodek-Yang and BGR18 calculations, as well as the YADISM prediction including PDF uncertainties. For the BGR18 and YADISM calculations we only display results corresponding to the perturbative region with Q 2 > 3.5 GeV 2 . In Fig. 4.4, the area covered in light gray indicates the Q 2 coverage of Region II, where the NNSFν parametrisation is constrained to reproduce the YADISM QCD boundary condition.
From the comparisons in Figs. 4.3 and 4.4 one can observe how the NNSFν predictions reproduce the YADISM boundary conditions at Q 2 = 100 GeV 2 in the relevant region of x. We verify that within the whole Region II there is agreement within uncertainties between NNSFν and YADISM, demonstrating that as required the QCD boundary condition is being reproduced by the structure function parametrisation. At medium and small-x, there is a good agreement between the BGR18 calculation and the NNSFν predictions in the region of validity of the former (for Q 2 ≥ 4 GeV 2 ), with some differences in the large-x region. It is interesting to note that the agreement found between NNSFν and YADISM in Region II is not automatic: for instance at Q = 2 GeV (Region I), we find that for x ∈ [0.05, 0.3] the two results disagree within uncertainties, showing that the experimental neutrino data (rather than the QCD boundary condition) is driving the fit results there. Furthermore, one observes from these comparisons how the Bodek-Yang calculation falls outside the 1σ error band of the NNSFν prediction for a significant region of the relevant x and Q 2 values, in particular for F 2 at small-and large-x. For instance, for x ≃ 0.25 in the low-Q region, the Bodek-Yang prediction is around 25% smaller than the NNSFν one. The agreement between the NNSFν and Bodek-Yang structure functions improves for heavier nuclei, as we will demonstrate in Sect. 5 when evaluating the inclusive neutrino cross-sections on iron, tungsten, and lead targets.
Another interesting feature of Figs. 4.3 and 4.4 is the behaviour of NNSFν in the extrapolation regions. Concerning the x dependence, the NNSFν uncertainties increase at small-x specially at low-Q due to the lack of direct experimental data, while in the same x region at higher values of Q these uncertainties are reduced due to the information provided by the QCD boundary condition. Note that in the case of the xF 3 structure function, the small-x behaviour is fixed in the case of the ν +ν combination, rather than for the individual ν andν structure functions. Concerning the extrapolation in Q 2 , the NNSFν uncertainties increase as Q 2 decreases due to the lack of data, and decrease as Q 2 increases as a consequence of the constraints from the QCD boundary condition.
Comparisons with experimental data. The values of the χ 2 exp reported in Table 4.1 indicate good agreement between the experimental data and the NNSFν parametrisation. This agreement can be further illustrated by comparing the NNSFν predictions with representative datasets entering the fit in Region I in selected kinematic regions as a function of Q 2 , as done in Fig. 4.5. For each dataset we also indicate the values of x and A for the bin shown. For the experimental data points, the error band corresponds to the diagonal entry of the associated covariance matrix. Specifically, we show the F 2 and xF 3 structure functions (averaged over ν andν) for the BEBCWA59, CHARM, and CCFR experiments. A similar level of agreement is obtained for the various regions of x, Q 2 , or A considered in this analysis, again indicating a well-balanced fit where the different kinematic regions are satisfactorily described.  Uncertainty estimate and kinematic dependence. Concerning the uncertainty estimate of the NNSFν determination, in general its 1σ errors are found to be rather smaller as compared those of the corresponding experimental measurements, see also Fig. 4.5. This behaviour is expected, since effectively the neural network parametrisation is averaging over the input data [63] which displays partly overlapping kinematic coverage. In addition, one has to account for the effects of the QCD boundary conditions, and indeed one can verify that as the lower boundary of Region II (Q dat ) is approached, the NNSFν uncertainties decrease as a consequence of these constraints. These effect are illustrated in Figs. 4.6 and 4.7, which display the absolute 68% CL relative uncertainties in the NNSFν structure functions as Q 2 is varied for x = 0.25 and as x is varied for Q 2 = 2 GeV 2 , respectively. We compare the uncertainties for the isoscalar free nucleon 2 H target with those for Ne, Fe, and Pb targets, for the three structure functions and different initial states. The uncertainties of the NNSFν determination stabilize in Region II, where they approach those of the QCD boundary condition based on YADISM and nNNPDF3.0. In Region II, the Q 2 dependence of the absolute uncertainties is moderate and arising from scaling violations. In the low-Q 2 extrapolation region, the NNSFν uncertainties increase as a consequence of the limited experimental constraints, as also mentioned above. Uncertainties also increase as x decreases, both for F 2 (which rises at small-x) as for xF 3 (which being a non-singlet does not). Lastly, the NNSFν uncertainties become small in the large-x region where structure functions vanish due to the elastic limit.
Another noticeable feature from Figs. 4.6 and 4.7 is the dependence of the NNSFν uncertainties with respect to the atomic mass number A. At low Q values (Region I), uncertainties are the largest for 2 H, consistent with the fact that there is no experimental information in this region. As Q is increased, the NNSFν uncertainties for all nuclei become similar, specially for the F 2 and F L structure functions. In general, the most precise NNSFν prediction is obtained for an iron target, at least in the region of x and Q 2 being shown, which is consistent with the fact that Fe is the most abundant target in the input dataset. However, we point out that this is not the case in the small-x region with x ∼ < 10 −5 , since there are no constraints from D-meson production on an Fe target, and hence for the UHE neutrino cross-sections the uncertainties on a Fe target are higher than those of either a 2 H or a Pb target as shown in Sect. 5.2. The 68% CL (absolute) uncertainties in the NNSFν baseline fit as a function of Q 2 in the region between 1 GeV 2 and 100 GeV 2 for x = 0.25. We compare the results of an isoscalar nucleus 2 H with those corresponding to various nuclear targets entering the fit, namely 20 Ne, 56 Fe, and 208 Pb. We display results for the three structure functions F 2 , xF 3 , and F L , separately for the ν,ν and ν +ν initial states.

Stability and validation
We now study the stability of the NNSFν determination by comparing the baseline fit with variants where either the input dataset or some aspect of the fitting methodology are modified. In particular, we assess the dependence of the fit on the value of the threshold scale Q 2 dat separating Regions I and II; quantify its stability when double-differential cross-section data is replaced with their structure functions (F 2 , xF 3 ) counterparts; and assess the quality of the NNSFν interpolation to values of the atomic mass number A not considered in the fit.
Dependence on the matching scale. Fig. 4.8 compares the NNSFν baseline results with a fit variant in which the matching scale between Regions I and II has been increased from Q 2 dat = 30 GeV 2 (just above the bottom mass) to Q 2 dat = 50 GeV 2 . Results are shown at Q 2 = 40 GeV 2 as a function of x for both Fe and Pb targets, and are normalised to the central NNSFν baseline. We display results for the three structure functions separately for neutrinos, antineutrinos, and for their sum. One finds that, in all cases considered, this NNSFν variant is in agreement with the baseline fit within uncertainties, demonstrating the fit stability with respect to variations of the Q 2 dat threshold value. This stability is particularly visible for F 2 and xF 3 , while somewhat larger effects are observed for F L in the large-x region. The likely explanation of this effect is that experimental constraints on F L are limited and hence this structure function is more sensitive to the settings of the matching to the QCD predictions. Some differences are also observed for the lead structure function F 2 for x ∼ < 0.01. In the case of this target the direct experimental constraints end at x ≈ 0.025, so again the matching scale has some influence on the results. Nevertheless, agreement within uncertainties is preserved in all cases, demonstrating that the NNSFν analysis is robust with respect to moderate variations of the hyperparameter Q 2 dat . This said, as  Fig. 4.6 now as a function x for Q 2 = 2 GeV 2 discussed in Sect. 3.1, Q 2 dat can neither be arbitrarily reduced, which would cut away most of the neutrino data, nor increased, since a gap would arise between the data region and the region where QCD boundary conditions are imposed.
Reduced cross-sections vs structure functions. Fig. 4.9 presents the same comparison as that in Fig. 4.8 now for a NNSFν fit variant in which the data on the double-differential cross-sections for the CDHSW, CHORUS, and NuTeV experiments has been replaced by the corresponding measurements at the level of the individual structure function F 2 and xF 3 , see also the discussion in Sect. 3.2. We note that since in these experiments the separate structure functions are extracted from the differential cross-sections by means of a theory-assisted averaging procedure, this replacement entails removing 3805 cross-section data points and replace them by a much smaller number, 490, of structure functions data points.
In terms of the fit quality, from the outcome of this fit variant we find a marked deterioration of the χ 2 exp , which increases to 1.450 as compared to the value of 1.287 for the baseline (see Table 4.1). This fit quality worsening can be traced back to the effect of the CHORUS and NuTeV experiments in particular, where for instance the F 2 structure function data has χ 2 exp = 2.073 and 2.445 per data point respectively. We note that the poor fit quality to the NuTeV F 2 and xF 3 structure function data has already been reported and studied in the literature (see [95][96][97] and references therein), and concerns about the internal consistency of this dataset have been raised. While our choice of reduced cross-sections d 2 σ/dxdQ 2 for the baseline dataset is motivated by a priori considerations based on it being a more robust, less theory-dependent observable, the poor description of the F 2 data from CHORUS and NuTeV provides a further argument in favour of the choice adopted here.
Concerning the impact that this dataset variation has at the level of the NNSFν output, from Fig. 4.9 we find that in general results are compatible at the one-sigma level, specially when considering structure functions averaged over neutrinos and antineutrinos. Nevertheless, non-negligible differences, not covered by the respective uncertainties, are observed for instance for Fν 2 at large-x both for Fe and Pb, for xF 3   for A = 56 (the NuTeV target) when separated into neutrino and antineutrino predictions, and for F L for intermediate x and also for A = 56. We note that in the latter case, this fit variant does not include direct experimental constraints on F L and hence the only information provided by the fit comes from the QCD boundary conditions in Region II. Taking into account both the deterioration at the fitted χ 2 exp level, the lack of agreement in the fitted structure functions for specific regions of (x, Q 2 , A), and their poor description when the baseline NNSFν predictions are used, one concludes the separate F 2 and xF 3 structure functions are not equivalent, and may be inconsistent, as compared with their differential cross-sections counterparts which is our default choice. Even so, we note that at the level of F 2 and xF 3 averaged over neutrinos and antineutrinos, except for F 2 at x ∼ > 0.2, the baseline and the variant fits are in agreement at the 68% CL and hence predictions obtained from them, for instance for inclusive cross-sections, are also likely to be in agreement, the only possible exception being low E ν values where the large-x region dominates.
Interpolation in atomic number A. One of the key advantages of the NNSFν strategy is the ability to interpolate the predictions for neutrino structure functions to other targets, with different atomic mass numbers, beyond those directly considered in the fit and that may be relevant for neutrino phenomenology, such as oxygen (A = 16), argon (A = 39), calcium (A = 40), and tungsten (A = 184) among others.
Here we validate the NNSFν interpolation in A by comparing the baseline fit with two variants. First, we compare the baseline fit results with those of a variant in which the datasets with atomic mass number A ≃ 20, namely BEBCWA59 (Ne) and CHARM (CaCO 3 ), are excluded. In such case, the predictions for A = 20 of the fit variant in Region I will be obtained from a extrapolation of the constraints provided by the fitted data with A ≥ 56 for iron and lead targets. Second, we compare the NNSFν predictions from the baseline fit in Region II with the YADISM+nNNPDF3.0 calculations for calcium (A = 40), a nuclear species which is not contained in the baseline analysis. These two tests make possible validating the interpolation (and extrapolation) of the NNSFν predictions to new values of A not used in the fit neither in Region I nor in Region II.
First, Fig. 4.10 displays the results comparing the baseline fit with the variant in which the A ≃ 20 datasets, namely BEBCWA59 (Ne) and CHARM (CaCO 3 ), are excluded. We display the F 2 and xF 3 structure functions for A = 20 for Q 2 = 4 GeV 2 as a function of x, separately for the neutrino, antineutrino, and their sum initial states. In this kinematic region, the predictions for A = 20 of the NNSFν variant arise entirely from extrapolating the constraints provided by the fitted data with A ≥ 56 for iron and lead targets. As one can observe from the comparison, the two fits are compatible within the respective uncertainties, with the possible exception of F νA 2 for x ≃ 0.2 where the baseline and its variant overlap at the 2σ level. The good agreement between the baseline fit and its variant, specially for F 2 for x ∼ < 0.1 and for xF 3 , confirms that the methodology extrapolates in the atomic mass number A in a way that is consistent with the available experimental constraints in Region I, and that the predicted structure functions are stable upon removal of a subset of the fitted datasets. Fig. 4.11 then compares the NNSFν predictions for A = 40 with the corresponding YADISM calculations, with nNNPDF3.0 as input, for a calcium 40 Ca target. As in Fig. 4.10, we compare the F 2 and xF 3 neutrino structure functions but this time at Q 2 = 100 GeV 2 , which corresponds to Region II. Since neither the fitted data (in Region I) nor the QCD boundary conditions (in Region II) include information on A = 40, this comparison assesses the interpolation capabilities of NNSFν to atomic mass numbers not considered in the fit. The error band in the YADISM predictions of Fig. 4.11 contains the PDF uncertainties but not the MHOU ones.
One finds that the interpolated predictions of NNSFν for A = 40 are in very good agreement with the YADISM calculation, which do not enter the fit. While excellent agreement can be seen at the averaged level, some differences are noticed when looking at large-x of the separate ν-andν-results for xF 3 . These discrepancies could however be attributed to the small constraints on the sea-quark distributions that cancel out when taking the average. This result implies that the NNSFν parametrisation interpolates in A in a manner consistent with the underlying behaviour of the QCD prediction in Region II, and furthermore that the fit results would not be affected in the case that QCD data with A = 40 is added to the NNSFν fit. These results demonstrate that NNSFν is able to reliably interpolate between nuclear targets and therefore capable of providing predictions for A values for which no direct experimental measurements are available. We use this feature to provide LHAPDF sets for all values of A relevant for neutrino phenomenology as listed in App. A.

The Gross-Llewellyn Smith sum rule
Finally, we use the NNSFν determination to evaluate the Gross-Llewellyn-Smith sum rule of neutrino structure functions, Eq. (2.16), and assess its agreement with respect to the perturbative QCD calculation.   The latter prediction indicates that for an isoscalar target A one expects to find with n f being the number of active flavours at the scale Q 2 and the coefficients c k known up to third order in perturbation theory. Evaluating Eq. (4.1) requires extrapolating down to the x → 0 region which is not accessible experimentally, and hence instead we compute a truncated variant of the GLS sum rule, for different values of the lower integration limit x min and the atomic mass number A. The truncated sum rule Eq. (4.2) should converge to the QCD prediction in Eq. (4.1) in the x min → 0 limit. Furthermore, the GLS sum rule should hold true irrespective of the value of the mass number A entering the calculation. Fig. 4.12 displays the outcome of the calculation of the (truncated) Gross-Llewellyn-Smith sum rule, Eq. (4.2), evaluated with the NNSFν baseline fit as a function of Q 2 . We display results corresponding to lower integration limits of x min = 10 −3 and of 10 −4 for A = 1, A = 40, and A = 56 nuclei. We compare the NNSFν baseline truncated results with the corresponding QCD predictions for the exact GLS sum rule. Note that the latter is the same in all panels, since it is independent of both x min and A, and that its Q 2 dependence is entirely dictated by that of the running of the strong coupling α s (Q 2 ). As in the rest of this section, the uncertainty in the NNSFν results is computed as the 68% CL intervals from the N rep = 200 replicas that constitute the representation of its probability density.
From this comparison, one can conclude that there is agreement within uncertainties between the (truncated) calculation of the GLS sum rule from the phenomenological NNSFν determination of the xF 3 neutrino structure function and the corresponding perturbative QCD prediction. The agreement in central values, especially in the large-Q 2 region, slightly deteriorates once x min is reduced from 10 −3 to 10 −4 . But at the same time, the uncertainties in the truncated sum rule calculation increase as x min decreases, as expected since NNSFν does not have direct experimental constraints below x ≃ 0.01, see also Table 3.1. Even more remarkably, not only the value but also the slope with Q 2 , which in the QCD calculation is dictated by the α s (Q) running, is correctly reproduced by the data-driven NNSFν analysis. The analysis of Fig. 4.12 hence indicates that the NNSFν data-driven determination of xF 3 is consistent with the QCD expectations, and in particular that extrapolates to the small-x, low-Q 2 region in a manner which improves the agreement with the N 3 LO result given by Eq. (4.1).
The agreement between the truncated GLS sum rule computed with NNSFν and the QCD calculation holds true for the three values of A considered, namely A = 1, 40, and 56. The main effect of the A dependence appears to be the increase in the uncertainties in the A = 1 case as compared to the heavier nuclei, due to the lack of direct experimental constraints on light nuclear targets. This insensitivity of the results of the truncated GLS sum rule with respect to the value of A is consistent with the expectation that this sum rule is related to the nucleon valence sum rules, which are satisfied and take the same values irrespective of the value of A entering the calculation of the neutrino structure functions.

Inclusive neutrino cross sections
Here we deploy the NNSFν determination of neutrino structure functions to evaluate inclusive neutrino scattering cross sections as a function of E ν for different projectiles and targets. Specifically, we provide predictions for the complete range of energies relevant for neutrino phenomenology, from the GeV region entering accelerator experiments up to the multi-EeV energies of cosmic neutrinos. Our predictions come accompanied with an estimate of the dominant uncertainties relevant in each energy region.
First, we compare the NNSFν determination with previous results in the literature, in particular with the Bodek-Yang, BGR18, and CSMS11 calculations. We evaluate the mean inelasticity ⟨y⟩ for different neutrino energies, compare with experimental measurements of the inclusive cross-section from NuTeV, assess the sensitivity of our calculation to the low-Q region, and quantify the impact of nuclear corrections for the most relevant target materials. Furthermore, we provide dedicated predictions for the energy range and target materials required for the interpretation of far-forward neutrino scattering experiments at the LHC.

The NNSFν cross sections
Starting from the double-differential neutrino-nucleus scattering cross sections of Eqns. (2.5) and (2.6), the inclusive cross section is obtained [41] by integrating over the kinematically allowed range in x and Q 2 , where the integration limits are given by and with the inelasticity related to the neutrino energy by y = Q 2 /(2m N xE ν ). In an inclusive calculation, the lower integration limit Q 2 min should go all the way down to Q 2 min = 0. While most previous studies impose a cut in Q 2 min to restrict the integration in Eq. (5.1) to the perturbative region, this is not required within the NNSFν approach since our structure function predictions are valid for all Q 2 values.
For the calculations presented in this section we take Q min = 0.03 GeV and study the dependence of the predictions with respect to variations of this choice. A kinematic cut in the final-state invariant mass of W ≥ 2 GeV is applied to restrict the calculation to the inelastic scattering region. See App. A for technical details about the implementation of Eq. (5.1) and the estimate of the associated uncertainties.
Comparison with Bodek-Yang, BGR18, and CSMS11. σ νN /E ν is provided in units of 10 −38 cm 2 /GeV per nucleon, hence assuming the proton and neutron content of the average nucleon in the target nuclei. We display neutrino energies from E ν ≃ 10 GeV up to E ν ≃ 10 11 GeV = 100 EeV. The NNSFν determination is the only available theory prediction applicable in the complete range of E ν relevant for neutrino phenomenology. The region of formal applicability of the BY calculation is E ν ≤ 10 5 GeV, while CMS11 and BGR18 are restricted to E ν > 100 GeV. Furthermore, BGR18 was optimised to high-energy scattering and does not provide a good description of the energy region sensitive to low-Q values. In the comparisons of Fig. 5.1, the CMS11 and BGR18 predictions are provided for a free isoscalar target and hence do not account for nuclear modification effects. The latter could be included by means of an external nuclear PDF analysis within a factorised approach, but only in the perturbative QCD region. Nuclear corrections are data-driven and model-independent in NNSFν, while they are model-dependent and neglect the constraints from proton-ion collisions at the LHC in the Bodek-Yang case. These nuclear effects are specially significant for heavy nuclei such as lead, as we discuss below.
The uncertainty band of the NNSFν prediction varies from a few percent up to a maximum of 15%, depending on E ν and the nuclear target. At low energies, it is the largest for 2 H, given the lack of direct experimental constraints on low-Q neutrino-hydrogen structure functions. At very high energies, it is the largest in iron since both for a free nucleon and for a lead target the nNNPDF3.0-based calculation accounts for the constraints on small-x PDFs provided by charm production at LHCb. On the intermediate energy region with 100 GeV ∼ < E ν ∼ < 100 TeV, uncertainties in the NNSFν calculation of the inclusive cross-sections are at the few percent level at most.
Concerning the comparison between NNSFν and the theory predictions from other groups, starting with the case of lead nuclei one observes agreement within uncertainties with the Bodek-Yang calculation for E ν ∼ < 10 TeV, except for antineutrinos with very low energies. For higher energies, NNSFν agrees with BGR18 for energies above 1 TeV up to 10 6 TeV, with the latter overshooting the former for E ν ≥ 10 6 TeV due to the missing nuclear corrections related to the strong small-x quark and gluon shadowing in lead nuclei found in nNNPDF3.0. The differences with CSMS11 at high energies are explained from the absence of nuclear effects, the choice of proton PDFs and the corresponding small-x behaviour, and the treatment of top quark mass effects, see the discussion of [22,41].
In the case of an iron target, there is also good agreement between NNSFν and the Bodek-Yang predic- tions for neutrinos, while for antineutrinos the latter is suppressed by a factor around 10%. As demonstrated below in Fig. 5.2, the NNSFν antineutrino predictions on Fe are preferred by the NuTev measurements of the inclusive cross-sections. In the case of the ν +ν sum, the differences between NNSFν and BY are relatively moderate. The NNSFν predictions are consistent with BGR18 for E ν ∼ > 100 TeV, as expected from the use of common settings for the QCD structure functions and the similarities in the input (n)PDFs, given that no strong nuclear small-x shadowing effects have been identified in the case of lead. The CSMS11 calculation agrees with NNSFν until around 10 TeV, and then overshoots it for the reasons discussed in the case of lead [22,41].
Finally, concerning the isoscalar free-nucleon target 2 H, NNSFν predicts a larger cross-section as compared to Bodek-Yang by a factor between 10% and 20%, depending on the E ν value. For this target, there exist no direct experimental constraints on neutrino-nucleon scattering, and hence in NNSFν the A = 1, low-Q behaviour arises from the extrapolation from the data with A ∼ > 20 included in the fit and from the QCD constraints in Region II for Q ≥ Q thr with A = 1. Despite the absence of direct constraints, thanks to the NNSFν matching procedure for energies within E ν ∼ > 1 TeV a reliable prediction can be obtained by means of the perturbative QCD calculation using the reasonably accurate proton PDFs. In this high-energy region, the NNSFν prediction is bracketed by the CSMS11 from above and BGR18 from below. The differences with BGR18 are explained by a combination of the input PDFs and the treatment of top quark mass effects [22], since NNSFν uses FFNS5 and BGR18 is based on FONLL instead, see also Sect. 2.
Nevertheless, neutrino cross-sections on a free-nucleon target are of limited phenomenological interest due to the absence of data for this target material. This said, the large spread in the theory predictions for a free-nucleon target motivates future experimental analyses of neutrino scattering in the region E ν ∼ < 1 TeV with data taken on hydrogen, deuterium, or other light nuclei as targets.
Comparison with experimental data. Several experiments have measured the inclusive neutrinonucleus interaction cross section for different nuclear targets and energy ranges. To further validate the NNSFν predictions on neutrino cross-section measurements, Fig. 5.2 compares NNSFν with the NuTeV [83] experimental data on inclusive neutrino cross sections an iron target, separately for neutrinos and antineutrinos. We focus on the energy region (between 20 GeV and 380 GeV) relevant for the interpretation of the NuTeV data. In addition to NNSFν, we also include the predictions from the Bodek-Yang, CSMS11 and BGR18 calculations. For the energy region in which inelastic scattering dominates, E ν ∼ > 100 GeV, we find excellent agreement between the NuTeV data and the NNSFν predictions, separately for neutrinos and antineutrinos. The Bodek-Yang calculation undershoots the NuTeV measurements on antineutrinos by up to 10%. CSMS11 is in good agreement with NNSFν and the NuTeV data, while the BGR18 calculation, optimized for highenergy scattering, undershoots it for both neutrinos and antineutrinos. We note that the invariant mass cut W ≥ 2 GeV is applied only to the theoretical calculations but not to the NuTeV measurement. This difference explains the disagreement for E ν ∼ < 100 GeV where other contributions in addition to inelastic scattering become significant, in particular resonant scattering. In this low-energy region, NNSFν and Bodek-Yang exhibit the same qualitative behaviour.
Inelasticity. In neutrino-hadron scattering, the inelasticity variable y = Q 2 /(2m N xE ν ) represents the fraction of the incoming neutrino energy which is transferred to the hadronic final state and hence to the associated hadronic shower. The higher the value of y, the more energy that is transferred to the hadronic shower, facilitating the experimental measurement and characterisation of the latter. At intermediate energies, with E ν ∼ < 100 TeV, neutrino interactions are expected to produce on average more energetic hadronic showers than their antineutrino counterparts. Therefore, a measurement of the event-by-event inelasticity provides a useful handle in order to statistically separate ν-fromν-bar initiated interactions, a feature exploited in the by IceCube analysis of [98]. The expected mean value of the inelasticity as a function of the neutrino energy can be computed as using as input a given prediction for the double-differential scattering cross-section. Fig. 5.3 displays the mean value of the inelasticity, Eq. (5.3), as a function of the neutrino energy for both iron and lead targets. The NNSFν prediction for ⟨y⟩, together with the associated 68% CL uncertainty band, is compared with the central values of the Bodek-Yang, BGR18, and CSMS11 calculations. The bottom panels display the ratio of the different calculations with respect to the NNSFν central value. As mentioned above, for E ν ∼ < 100 TeV the average inelasticity is markedly larger for neutrino-initiated scattering, while at higher energies the prediction for ⟨y⟩ becomes projectile-independent. An overall good agreement is  Fig. 5.5 for an iron target and restricted to the energy region with E ν ≤ 10 7 GeV. We assess the dependence of the inclusive cross-section with respect to variations in the lower integration limit Q min in Eq. (5.1). The default value of Q min = 0.03 GeV is compared to calculations based on Q min = 1, 1.41, and 2 GeV. observed between the different predictions, with differences up to the 10% level depending on the target and E ν ranges. Specifically, at high neutrino energies, the NNSFν prediction for ⟨y⟩ in a lead target is around 10% smaller as compared to CSMS11, partially explained by the nuclear corrections accounted for in the former but not in the latter. Good agreement between the Bodek-Yang and the NNSFν calculations of ⟨y⟩ is obtained in the region of applicability of the former, with residual differences at the few percent level.
Sensitivity to the low-Q region. Fig. 5.4 displays the NNSFν inclusive neutrino scattering cross-sections on an iron target for energies between 10 GeV and 10 7 GeV, for different values of he lower integration limit Q min in Eq. (5.1). Specifically, the default value of Q min = 0.03 GeV is compared to calculations based on Q min = 1, 1.41, and 2 GeV respectively. This comparison allows one to determine the relative contribution of the low-Q region to the inclusive cross-sections as a function of E ν , an effect which is similar for neutrinos and antineutrinos.
From the analysis of Fig. 5.4 one determines that neutrinos with energies around E ν ≃ 100 GeV receive a contribution of up to 10% from the region with Q ≤ 1.41 GeV and up to 25% from the region with Q ≤ 2 GeV. As the neutrino energy is increased, the contribution from the low-Q region decreases, but even for E ν ≈ 1 TeV neutrinos, up to 10% of the inclusive cross-section can arise from momentum transfers of Q ≤ 2 GeV. As discussed in Sect. 2, this low-Q region is affected by sizable theory uncertainties (MHOUs, higher twists, mass effects, factorisation breakdown, ....), and hence  cross-sections in the region E ν ≤ few TeV. It is interesting to highlight that for neutrino energies between 300 GeV and a few TeV, relevant for the interpretation of LHC far-forward neutrino scattering experiments, neither a data-driven, theory-agnostic calculation nor a purely QCD calculation can accurately predict the inclusive cross-sections. For higher neutrino energies, E ν ∼ > 10 TeV, the NNSFν prediction becomes independent of structure functions in the Q ≤ 2 GeV region and hence a perturbative QCD calculation of the cross-section such as that provided by YADISM can be reliably deployed.

Impact of nuclear effects
As illustrated by Fig. 5.1, the NNSFν scattering cross section (per nucleon) on targets such as iron or lead is in general different to that evaluated on a isoscalar free-nucleon target. To further scrutinise the role that nuclear effects have on the inclusive neutrino cross-sections, Fig. 5.5 displays the NNSFν predictions for the inelastic cross-section (divided by E ν ) comparing the results of neutrino and antineutrino projectiles on deuterium (isoscalar free-nucleon), iron, and lead targets. The top panel displays the absolute cross-sections while the bottom ones show the ratio to the iron target baseline, sequentially for neutrinos, antineutrinos, and their sum. As in previous comparisons, the band in the NNSFν predictions indicates the associated 68% CL uncertainties.
To facilitate the identification of the different energy regions of interest for neutrino phenomenology, Fig. 5.5 also displays grey bands indicating the E ν regions relevant for different experiments and detection techniques. Specifically, from low to high energies, we indicate the coverage of: • 10 GeV ∼ < E ν ∼ < 300 GeV: KM3NET-ORCA and DeepCore/IceCube-upgrade.
This list is not exhaustive, and has been added to the plot for illustrative purposes only. As further discussed in App. C, two different types of nuclear effects are responsible for differences between scattering rates on a heavy nucleus and on a free proton (hydrogen target). The first is related to the different content in protons and neutrons, with iron (lead) containing Z = 26 (Z = 82) protons and A − Z = 30 (126) neutrons. Therefore, as compared to a hydrogen target, the heavy nuclear targets display an enhanced (suppressed) content of down (up) valence quarks which leads in turn to an enhancement (suppression) of the associated structure functions for neutrino (antineutrino) scattering. This effect is most relevant in the valence structure function region, while as E ν increases the cross section becomes dominated by small-x scattering involving isospin-symmetric sea quarks and gluons. In the comparisons in Fig. 5.5 this effect is partially factored out since we normalise to an isoscalar 2 H target, though the non-isoscalarity of heavy nuclei, specially in lead, still plays a role.
The second type of nuclear effect is the modifications of the structure of bound nucleons as compared to their free-nucleon counterparts. In the perturbative QCD region, these modifications are encoded by the nuclear PDFs, here taken from the nNNPDF3.0 determination. The corrections, quantified in App. C, are similar for neutrinos and antineutrinos and become most important in the shadowing region at medium and small-x, where a strong suppression in heavy nuclei is preferred. Nuclear structure modification effects are also present for E ν ∼ < 1 TeV, a region where they are partially described by the DIS calculations but which also receives contributions from the SIS region which cannot be expressed in terms of factorised nuclear PDFs. Within the NNSFν framework, one estimates the nuclear modifications in Region I by allowing a free dependence of F i (x, Q 2 , A) to be directly constrained from the data, and then matched to the QCD calculation in Region II.
From the comparisons in Fig. 5.5, one observes that differences between inclusive cross-sections on Fe and Pb targets and those on 2 H targets are moderate (a few percent at most) for energies in the intermediate region between 10 TeV and a few PeV. In the high energy region, for E ν ∼ > 10 TeV, the most marked effect is the strong suppression of the cross-sections in lead as compared to those in a free-nucleon target, reaching up to 20% at 100 EeV. This suppression is a consequence of the small-x shadowing of quark and gluons in nNNPDF3.0. Given the large nPDF uncertainties associated to an iron target, high-energy cross-sections for Fe are compatible within errors with those for both 2 H and Pb.
In the region where the non-DIS contribution is sizable, E ν ∼ < 1 TeV, differences between Fe and Pb targets for the ν +ν cross-section remain at the few percent level, with again the latter being suppressed in comparison to the former. The cross-sections on a 2 H target are enhanced at low energies, with the difference being up to 20% in the E ν = 10 GeV where the inelastic component is very small. In the same energy region, differences between the three nuclear targets appear to be less (more) significant once we consider separately the predictions for neutrinos (antineutrinos).
While in this section we discuss explicitely nuclear effects only for Pb and Fe targets, NNSFν structure functions and inclusive cross-section predictions are provided for all targets relevant for neutrino phenomenology and summarised in Table A.1. In addition, below we provide dedicated predictions for the energy ranges and target materials relevant for far-forward neutrino experiments at the LHC.

Far-forward neutrino scattering at the LHC
Precise predictions for neutrino scattering rates are a key ingredient for the interpretation of data from experiments aiming to detect and study the far-forward neutrinos produced in LHC collisions. These include the current SND@LHC and FASERν experiments as well as the proposed AdvSND@LHC, FLArE, and FASERν2, which would be installed in the Forward Physics Facility operating concurrently with the HL-LHC. The dominant component of these far-forward LHC neutrino fluxes lies in the region between a few hundreds of GeV and a few TeV [99,100], with the high-energy component most sensitive to neutrinos produced from charm meson decays. Here we provide dedicated predictions for the inclusive neutrino scattering cross-sections in the FPF energy range, specifically between E ν = 100 GeV and 10 TeV, assuming a tungsten (W) nuclear target with A = 184, which is the current target of FASERν and the intended target for the FASERν2 experiment. For completeness, we also provide predictions for neutrino scattering on an oxygen (O) target, which are relevant for the intepretation of ongoing and future neutrino oscillation measurements taking place at the DeepCore/IceCube and KM3NET-ORCA experiments in a partially overlapping energy range as the FPF.
With this motivation, Fig. 5.6 (left) presents the same comparison as in Fig. 5.1 now on a tungsten target and restricted to the energy region relevant for the FPF experiments. In this energy region and for this specific nuclear target, one observes excellent agreement between the NNSFν determination and the Bodek-Yang model predictions, except perhaps for antineutrinos with E ν of several TeV. Given that the two calculations are completely different in terms of both data and QCD input as well as in methodological assumptions, such agreement at the 1% level is remarkable. The CSMS11 (BGR18) predictions tend to overestimate (underestimate) the FPF neutrino scattering rates in comparison to the NNSFν baseline, specially for BGR18 at low energies where the calculation falls outside its regime of applicability.
All in all, Fig. 5.6 demonstrates that the uncertainties on the inclusive neutrino cross-section predictions for a W target at FPF energies are estimated to be at most at the 5% level, and even less for the combined ν +ν cross-section. We conclude that state-of-the-art calculations of neutrino scattering in the SM at the FPF energies can achieve a precision of a few percent, thus providing an excellent starting point for further investigations of neutrino interactions in this unexplored energy range, for instance in terms of anomalous couplings or of EFT effects.
The right panel of Fig. 5.6 displays a similar comparison as that of Fig. 5.5, now for the NNSFν predictions for tungsten and oxygen targets. While neutrino cross-sections display almost no differences between the two targets, for the antineutrino case we observe sizable nuclear effects, of up to 20% and specially for energies below 1 TeV. As the energy is increased, the differences related to nuclear corrections wash away but even at 10 TeV they can represent a 10% effect. The predictions for the combined ν +ν crosssection, averaged over the two projectiles, lead to more moderate differences but still up to the 10% level at energies of few hundreds of GeV. Fig. 5.6 indicates that properly accounting for nuclear effects associated to the target material is required in order to translate constraints on neutrino interactions obtained from the FPF experiments to those from KM3NET and IceCube analyses and vice-versa, specially in those cases where one is sensitive to the separation between neutrinos and antineutrinos in the initial state.

Summary and outlook
We have presented a novel approach to the determination of neutrino inelastic structure functions based on the combination of a data-driven parametrisation at low and moderates values of Q 2 matched to perturbative QCD calculations at high Q 2 . The resulting structure functions, dubbed NNSFν, enable the evaluation of inclusive neutrino scattering cross sections over 12 orders of magnitude in energy E ν from a few GeV up to the multi-EeV region. In particular, it makes possible the accurate evaluation of scattering event rates for neutrinos with energies between 100 GeV and a few TeV, relevant for both far-forward neutrino detection at the LHC and for atmospheric neutrino oscillation experiments, without the need to impose acceptance cuts in Q 2 . The NNSFν determination comes accompanied with a faithful estimate of the associated uncertainties and accounts for the constraints provided by state-of-the-art determinations of nucleon and nuclear structure.
We have compared the NNSFν determination with existing calculations in the literature, such as the popular Bodek-Yang model, both at the level of structure functions and of inclusive cross sections. We have demonstrated the agreement between our parametrisation and the experimental data, highlighted the smooth matching between the data-driven extraction and the QCD boundary conditions, and assessed the behaviour of NNSFν in the extrapolation regions. We have quantified the impact of nuclear PDF modifications, finding that these are most significant for heavy nuclei such as tungsten and lead. We have also studied the implications of the NNSFν analysis for key observables in neutrino DIS, such as the Gross-Llewellyn-Smith sum rule.
The NNSFν determination is made available in terms of LHAPDF interpolation tables [66] for neutrino structure functions corresponding to the relevant nuclear targets, as described in App. A. Two different sets of grids are provided, one associated to the neural network determination and the other to the YADISM calculation, and a Python script showcasing how to use and combine these two grids is provided. We also provide a script that takes as input these LHAPDF structure function grids and integrates them to evaluate the inclusive cross section and the corresponding uncertainties. For completeness, we also make available look-up tables with the NNSFν predictions for inclusive cross sections σ νA (E ν ) in the relevant range of E ν and target nuclei. The LHAPDF grids, Python driver scripts, and look-up tables are available from the project GitHub repository https://nnpdf.github.io/nnusf together with the data, code, and theory tables that enter the structure function determination. The NNSFν structure functions have also been implemented in the GENIE [20] event generator, where it can be accessed in a similar manner as the BGR18 and Bodek-Yang calculations. It is hence readily available in a software framework already integrated in the analysis pipeline of many neutrino experiments.
The results presented in this work could be expanded in different directions. First, far-forward neutrino measurements taking place at LHC experiments such as FASERν and SND@LHC, and in the longer term the Forward Physics Facility, will provide further constraints on the modelling of neutrino scattering at low and moderate Q 2 . Second, the NNSFν calculation could be integrated into other neutrino event generators to streamline its usage, also for neutrino experiments where inelastic scattering, without dominating, still represents a sizable contribution to the total event rates. Third, by incorporating updated determinations of proton and nuclear PDFs benefitting from new experimental constraints and higher-order theory calculations, such as the information provided by charm production at LHCb combined with small-x BFKL resummation to improve predictions of neutrino structure functions at UHE energies [40,41,71]. Fourth, by considering the information provided by lattice QCD calculations [101] to constrain structure functions in the non-perturbative low Q 2 region.
All in all, the NNSFν analysis bridges a significant gap in current neutrino phenomenology by providing a consistent determination of structure functions enabling the calculation of cross sections in the full range of E ν relevant for inelastic scattering. For the first time, a determination of structure functions can be used for predictions involving accelerator neutrinos with E ν ∼ few GeV, atmospheric and collider neutrinos with energies between ∼ 100 GeV and few TeV, and ultra-high energy neutrino scattering at EeV energies. The availability of NNSFν allows bypassing a major limitation of current analyses where different cross-section calculations, inconsistent among them, have to be adopted depending on the range of energies involved.

A Usage and delivery of the NNSFν determination
The software framework used to produce the NNSFν determination of neutrino inelastic structure functions, together with the results obtained from it, can be obtained from the project website: https://nnpdf.github.io/nnusf Specifically, there one can find the following: • Installation instructions, code documentation, and user-friendly examples of the NNSFν software framework, including the data and theory runcards that allow reproducing the results of this paper.
• Fast interpolation grids in x and Q 2 for the inelastic structure functions F 2 , xF 3 , and F L for neutrinos, antineutrinos, and their sum. These structure function grids can be accessed by means of the LHAPDF interface and are provided for all nuclear targets of phenomenological relevance for neutrino scattering experiments.
• A driver code that evaluates, taking these structure functions grids as inputs, the corresponding central values, uncertainties, and correlations. As explained below, two structure function parametrisations with a different coverage in the (x, Q 2 ) plane are provided and the driver code takes care of their combination.
• A second Python driver code that evaluates the inclusive neutrino cross section Eq. (5.1) as a function of E ν from the LHAPDF grids and provides the corresponding uncertainty estimate.
• Look-up tables compiling the NNSFν predictions for the inclusive cross sections σ νA and σν A for phenomenologically relevant nuclei as a function of E ν , together with their uncertainty estimates.
In this appendix we provide details on these deliveries associated to the NNSFν determination, from the code and documentation to the inclusive cross-section look-up tables. In particular, we describe the prescriptions required to evaluate central values and uncertainties for the structure functions and inclusive cross sections from the LHAPDF grids provided.
The NNSFν framework. The framework developed in this work provides a stand-alone code to parametrise structure functions from experimental data in the presence of general theory constraints. It provides an independent implementation of the NNPDF fitting methodology, and shares the main methodological aspects with the NNPDF4.0 and nNNPDF3.0 proton and nuclear PDF determination such as the use of Stochastic Gradient Descent for the minimisation. For the first time in this context, neural networks with three independent inputs (x, Q 2 , A) are used to parametrise unknown functions from experimental data.
low-Q grid high-Q grid (1,2) [D] NNSFnu D lowQ NNSFnu D highQ (2,4) [He] NNSFnu He lowQ NNSFnu He highQ (3,6) [Li] NNSFnu Li lowQ NNSFnu Li highQ (4,9) [Be] NNSFnu Be lowQ NNSFnu Be highQ (6,12) [C] NNSFnu C lowQ NNSFnu C highQ (7,14)   for each nuclear target and the names of the corresponding low-Q and high-Q grids. Each grid provides, as a function of x and Q 2 , as outputs the neutrino and antineutrino structure functions F 2 , xF 3 , and F L as well as their sum. We note that both the number of replicas and the statistical interpretation thereof is different in the low-and high-Q grids, see text for more details.
Structure functions. We indicate the values of (Z, A) for each nuclear target and the names of the corresponding low-Q and high-Q grids. Each grid provides, as a function of x and Q 2 , as outputs the neutrino and antineutrino structure functions as well as their sum, namely with the PDG ID codes being respectively [1001,1002,1003,2001,2002,2003,3001,3002,3003]. We emphasize that the number of replicas, coverage on x and Q 2 , and statistical interpretation of the lowand high-Q grids is different, as discussed below. The nuclear targets listed in Table A.1 are selected for their relevance in the interpretation of past and present experimental data as well for theoretical predictions for upcoming experiments. For instance, tungsten (W) is relevant for the event rate predictions at FASERν and its eventual successor FASERν2, while oxygen (O) enters scattering rates for neutrinos on air, water, or ice targets. We also provide a grid for (A, Z) = (15,31), which correspond to the average nuclear numbers for Earth matter relevant for calculations of UHE neutrino attenuation when crossing the Earth [22]. The low-Q and high-Q NNSFν structure function grids are defined as follows: • low-Q grid. This grid encapsulates the direct output of the NNSFν neural networks trained on the neutrino structure function data and supplemented by the theoretical constraints from the QCD calculation based on nNNPDF3.0.
The LHAPDF grid is constituted by a central set and N rep = 200 Monte Carlo replicas, which represent the probability density in the space of fitted structure functions and from which statistical estimators such CL intervals can are evaluated with the usual NNPDF prescription. The central set is defined as the average over the N rep replicas. The uncertainty band associated with these replicas receives two contributions: the experimental uncertainties associated to the input fitted neutrino data and those associated to the methodology such as the functional uncertainty, in both cases subject to the constraints provided by the QCD boundary conditions imposed during the fit.
This low-Q grid can be used in the kinematic region 0.01 GeV ≤ Q ≤ 22 GeV and 10 −5 ≤ x ≤ 1 .
We note that extending this coverage could be achieved by enlarging the (x, Q 2 ) coverage of the QCD constraints added to the fit, but this may distort the fit quality by decreasing the weight given to the experimental data. If values of (x, Q 2 ) outside the region defined by Eq. (A.2) are requested, the output will be determined by the LHAPDF extrapolation algorithm.
• high-Q grid. This grid tabulates the direct output of the YADISM calculation of neutrino structure functions at NLO in QCD with nNNPDF3.0 together with the corresponding uncertainties. This calculation is independent of the output of the neural network parametrisation used for the low-Q grid, though the respective central values are matched due to their use of the same QCD calculation.
The LHAPDF grid is constituted by a central set, N rep = 200 replicas that correspond to the nNNPDF3.0 replicas and structure functions obtained with central factorisation and renormalisation scales, µ F = µ R = 1, and 9 additional "replicas" computed with the central nNNPDF3.0 set and with nine µ R and µ F scale variations following the procedure outlined in Sect. 3. The central set in this grid is evaluated as the average of the N rep = 200 replicas. Therefore this grid is constituted by 210 members (including the central predictions). PDF uncertainties and MHOUs can be obtained and combined by means of the procedure from [93], where the latter are obtained with the 9-point prescription and then added in quadrature with the PDF uncertainties. The scripts delivered with the grids illustrate how the uncertainty calculation is carried out.
This high-Q grid can be used in the kinematic region 2 GeV ≤ Q ≤ 10 TeV and 10 −9 ≤ x ≤ 1 , (A. 3) and again outside this region the grid output is obtained from LHAPDF extrapolation. Restricting the YADISM calculation to this region ensures the validity of the pQCD calculation.
We note that Eqns. (A.2) and (A.3) overlap in the region of (x, Q 2 ) where both the data-driven and pQCD approaches can be reliably applied. In this overlap region, we recommend using the outcome of the low-Q grid as baseline.
• Grid matching. Within their respective regions of applicability, one can evaluate the neutrino structure functions and the associated uncertainties using the appropriate prescription. If we denote by F νA i,lQ (x, Q 2 ) and F νA i,hQ (x, Q 2 ) the outcome of the F i structure function obtained from the low-Q and high-Q grids respectively, to evaluate the neutrino structure function one should use the prescription: with Q thr = 22 GeV, and the same holds for the corresponding uncertainties. We emphasize that the replicas of the low-Q grids can be treated as correlated among them, and the same holds for those of the low-Q grids, but that replicas of the low-Q grid are uncorrelated with those of the high-Q grid. The user can vary the matching parameters in Eq. (A.4), for instance by choosing a lower Q 2 threshold Q thr to switch to the high-Q grid calculation.
Since the same pQCD calculation that is tabulated in F νA i,hQ enters the data-driven fit of F νA l,hQ as theoretical constraint, central values obtained from the prescription of Eq. (A.4) should match within uncertainties as one crosses the threshold value Q thr . This does not necessarily hold for the structure function uncertainties, since these are typically larger in the data-driven fit than in the pQCD calculation. If required, the user may implement the matching of the structure function uncertainty across the threshold Q thr by means of a sliding window prescription.
As mentioned above, in the project website we provide a Python code that evaluates neutrino structure functions and their uncertainties for any value of x and Q 2 following the prescription outlined above.
Neutrino inclusive cross sections. Inclusive neutrino structure functions are evaluated by integrating the output of the structure function grids described above using Eq. (5.1). As highlighted by the analysis of Sect. 5.1, for neutrino energies of E ν ∼ > 10 TeV the contribution from the kinematic region involving Q ∼ < 2 GeV is negligible, and in such case one can evaluate the cross section in terms of only the high-Q structure function grid. Conversely, for lower neutrino energies, say E ν ∼ < 100 GeV, the integral in Eq. (5.1) will be dominated by momentum transfers with Q ∼ < 20 GeV and hence it can be evaluated in terms of entirely the output of the low-Q structure function grid. For intermediate values of the neutrino energy, one can combines the low-and high-Q structure function grid.
To further illustrate the kinematic coverage of the integrated neutrino cross-section calculation, Fig. A.1 displays the coverage in the (x, Q 2 ) plane of the inclusive cross section Eq. (5.1) in muon-neutrino inelastic scattering. We consider three values of the neutrino energy: from left to right, E ν = 11 GeV, 90 TeV, and 1 EeV. The blue bins indicate the regions in (x, Q 2 ) that contribute to the inclusive cross-section (normalised to the maximum value of the integrand), with darker bins dominating the integral. The region covered by red (green) indicates that the output of low-Q (high-Q) structure function grid is being used to evaluate the inclusive cross section. For E ν = 11 GeV, the cross section is determined entirely by the output of the low-Q grid. For E ν = 90 TeV, the output of the two grids are required, carrying a similar weight in the calculation. For E ν = 1 EeV, the cross section depends only on perturbative structure functions and the bulk of the contribution is associated to the region around x ≃ 10 −5 and Q ≃ 100 GeV. From this last case, one can see how at high neutrino energies it is important to also account for the region with x ∼ < 10 −5 and Q ∼ < Q thr , which is covered by the high-Q grid but not by the low-Q one.
Taking into account the contribution from the relevant kinematics, the integration is Eq. (5.1) is separated into three disjoint regions separated by the threshold Q 2 thr , in terms of structure functions computed with low-and high-Q grids, in analogy with Eq. (A.4), and with Q 2 min = 1.65 GeV to ensure that the high-Q grid is not extrapolated to the non-perturbative region. Depending on the E ν value, one we can simplify Eq. (A.5) to since at low and high energies only the low-Q and high-Q structure function grids are required respectively, with the other term leading to a negligible contribution. We note that, as also indicated by Fig. A.1, at high neutrino energies the high-Q grid is also used for Q ∼ < Q thr to determine the contribution for x ≤ 10 −5 which is not covered by the low-Q grid.
Given the different statistical interpretation of the uncertainties associated to the low-and high-Q structure function grids, the total uncertainty in Eq. (A.5) is evaluated by adding in quadrature the uncertainties associated to the low-and high-Q contributions, namely with δσ νN lQ and δσ νN hQ evaluated using the uncertainty prescription for the low-and high-Q structure function grids respectively. As opposed to the structure function case, the prescriptions Eq. (A.5), for central values, and (A.7), for the total uncertainties, lead to smooth predictions for neutrino energies E ν sensitive both to low-and high-Q structure functions. Furthermore, the outcome of the cross-section calculation is stable upon moderate variations of the threshold value Q 2 thr . The matching procedure described in this appendix to evaluate inclusive neutrino cross sections in terms of the NNSFν structure functions by means of the prescriptions of Eq. (A.5), for central values, and Eq. (A.7), for the total uncertainties, is illustrated in Fig. A.2. The NNSFν predictions of the inclusive cross sections for neutrino and antineutrino scattering off an iron target are compared to the corresponding individual contributions from the low-Q and high-Q structure function grids, denoted by σ νN lQ (E ν ) and σ νN hQ (E ν ) respectively in Eq. (A.5). For low (E ν ∼ < 1 TeV) and high (E ν ∼ > 100 PeV) values of the neutrino energy, the calculation is fully dominated by the contributions from the low-and high-Q structure function grids respectively, as indicated by Eq. (A.6), while for intermediate E ν values between 1 TeV and 100 PeV one must account for both contributions. It is worthwhile noting that the contribution from the low-Q grid remains relevant up to rather high neutrino energies. The smooth dependence in E ν both for central values and uncertainties further validates the NNSFν matching which connects the various energy regions.
B YADISM: DIS structure functions made easy Deep-inelastic structure functions can be evaluated with several public codes such as APFEL [69] and QCDNUM [102]. Available DIS codes differ in the accuracy with which structure functions can be computed, whether they are based on the x-space or the N -space formalism, the treatment of heavy quark mass effects and of target mass corrections, the availability of polarised and time-like coefficient functions, and the possible inclusion of QED corrections among other considerations.
This appendix introduces YADISM, a new framework for the evaluation of DIS structure functions from the same family as the EKO   One of the main advantages of YADISM is that it is integrated with the fast interpolation grid toolbox Pine-APPL [103], and hence DIS structure functions can be treated on the same footing as hadronic observables from the point of view of PDF fitting and related applications. PineAPPL provides a unique grid format, with application programming interfaces (APIs) for different programming languages and a command-line interface to manage the grid files. Furthermore, YADISM implements the available N 3 LO DIS coefficient functions, which combined with (approximate) N 3 LO evolution and heavy quark matching conditions available in EKO provide theoretical calculations required to carry out a N 3 LO PDF determination. YADISM will be described in an upcoming publication [104], here we summarise its main features, in particular those relevant to the present study, and highlight benchmarking studies carried out.
Grid formalism. As indicated by Eq. (2.9), in the perturbative regime DIS structure functions are given by the factorised convolution of process-dependent partonic scattering cross sections and of processindependent parton distribution functions, where j is an index that runs over all possible partonic initial states and C i,j is the process-dependent, but target-independent, coefficient function, given by an expansion in the QCD coupling α s (Q 2 ). In the third term of Eq. (B.1) and in the following, sum over repeated indices is implicit.
As standard for fast interpolation techniques developed in the context of PDF fits [103,[105][106][107], the PDFs can be expanded over an interpolation basis Already available as K-factors [67], now being integrated in the grid format. † Full calculation not available but an approximated expression can be constructed from partial results [108,109]. ‡ Calculation available, to be implemented. For each perturbative order (NLO, NLO, and N3LO) we indicate the light-to-light ("light"), light-to-heavy ("heavy"), heavy-to-light and heavy-to-heavy ("intrinsic") and "asymptotic" (Q 2 ≫ m 2 h limit) coefficients functions which have been implemented and benchmarked. The NNLO heavy quark coefficient functions for CC scattering are available in K-factor format and are being implemented into the YADISM grid formalism.
with p α (x) some suitable polynomial basis. This way the convolution in Eq. (B.1) can be replaced by a simple contraction F i = C j;i ⊗ f j = C jα;i · f α , C jα;i = C j;i ⊗ p α , (B.3) in terms of PDFs evaluated at fixed grid points ξ α and precomputed coefficients C jα;i . In YADISM the polynomial interpolation basis is provided by the EKO modules. The same grid structure can be generalised to accommodate extensions of the basic structure function calculation in Eq. (B.1) such as heavy quark mass effects, renormalisation and factorisation scale variations [93,94], and target mass corrections, among other effects. Isospin modifications, required to evaluate the neutron, deuteron, or heavy nuclear structure functions, can be accounted for either at the coefficient function level or at the input PDF level. The grid formalism summarised schematically in Eq. (B.3) requires as input the corresponding DIS coefficient functions. Table B.1 provides an overview of the different types and accuracy of the DIS coefficient functions currently implemented in YADISM. For each perturbative order (NLO, NLO, and N3LO) we indicate the neutral-current and charged-current light-to-light ("light"), light-to-heavy ("heavy"), heavy-tolight and heavy-to-heavy ("intrinsic") and "asymptotic" (Q 2 ≫ m 2 h limit) coefficients functions which have been implemented and benchmarked. The NNLO heavy quark coefficient functions for CC scattering are currently available in a K-factor format, and their implementation into the YADISM grid formalism is work in progress. We note that the full calculation of the N 3 LO NC heavy coefficient functions is not available, but that an approximated expression can be constructed from known partial results [108,109]. Heavy quark structure functions can be evaluated in the FONLL general-mass variable flavour number scheme (GM-VFN) [34], as well as in the fixed-flavour number (FFN) and zero-mass variable-flavour number (ZM-VFN) schemes. We point out that the list in Table B.1 is going to be updated as new features are added, and therefore the interested user is encouraged to consult the online documentation for an up-to-date states of available coefficient functions.
Scale variations. As done by other public DIS tools, YADISM also provides the option of varying the renormalisation and factorisation scales in the calculation. The code follows the definitions of scale variations from [110,111], which are consistent with the broader picture of scale variations relevant for PDF fits from [93] where they also affect the DGLAP evolution. There are two kinds of scale variations: renormalization scale Q R dependence, related to the ultraviolet renormalization scheme, and factorization scale Q F dependence, related to the subtraction of collinear logarithms in the adopted factorization scheme. The factorization scale Q F sets the boundary between the coefficient functions and the DGLAP-evolved PDFs. Scale variations at a given perturbative order can be constructed from combining ingredients already present at the previous perturbative order, and hence for this reason they represent a suitable predictor of potentially unknown missing higher orders. Within YADISM, the scale variation contributions to the DIS structure functions are stored in separate grids such that the values of the scale ratios µ 2 F = Q 2 F /Q 2 and µ 2 R = Q 2 R /Q 2 can be evaluated a posteriori.
As described in Sect. 3, in the NNSFν analysis the YADISM structure functions enter the fit to constrain the neural network parametrisation in Region II, with an error function defined in terms of a theory covariance matrix. This covariance matrix accounts both for the PDF and MHO theory uncertainties, the latter evaluated from scale variations using the 9-point prescription. The calculation of scale variations provided by YADISM and the subsequent determination of the MHOU theory covariance matrix has been benchmarked with the results of [93].
Benchmarking with APFEL. The calculations of DIS structure functions and reduced cross-sections provided by YADISM have been thoroughly benchmarked with those provided by APFEL and QCDNUM. Specifically, we have verified that YADISM reproduces the APFEL predictions for those of the DIS coefficient functions listed in Table B.1 which are also available in the latter. Excellent agreement is found in all cases considered, with some residual differences understood as will be discussed in more detail in [104].
To illustrate this agreement, Fig. B.1 displays the ratio between the YADISM and APFEL calculations of the neutrino-initiated structure function F νA 2 (x, Q 2 ) and of the corresponding double-differential crosssections at NNLO on a proton target for the same choice of input PDFs and theory settings. Specifically, we use in both cases the central replica of NNPDF4.0 NNLO with α s (m Z ) = 0.118 and FONLL-C for heavy quark mass effects. The benchmark comparison is presented both for fixed Q = 100 GeV and for fixed x = 0.01. The benchmark shows agreement at the ‰ level between the two calculations for the whole kinematic region in (x, Q 2 ) relevant for this study. A similar level of agreement is obtained for other DIS observables.

C Nuclear effects in neutrino scattering
As compared to baseline predictions on a proton target, there are two effects that modify the neutrino scattering rates on nuclear targets. The first is related to the different content of up and down quarks, since for an isoscalar nuclear target the average nucleon has on 25% less (more) up (down) valence quarks than in a proton target, an effect which is more marked for non-isoscalar nuclei where (A/2) > Z. The second is associated to the modifications of the bound nucleon structure that take place in heavy nuclei when compared with their free-nucleon counterparts, and in the DIS region is quantified by the nuclear PDFs. In this appendix we illustrate the impact of these two effects on neutrino inelastic structure functions.
For a general nuclear target, the F 2 structure functions per nucleon at LO in terms of nPDFs are The ratio between the YADISM and APFEL calculations of the neutrino-initiated structure function F νA 2 (x, Q 2 ) (top) and of the corresponding double-differential cross-sections (bottom panels) at NNLO for the same choice of input PDFs and theory settings. The benchmark comparison is presented both for fixed Q = 100 GeV (left) and for fixed x = 0.01 (right panels). A similar level of agreement is obtained for other DIS observables as well as for other regions of (x, Q 2 ) relevant for phenomenology. taking the sum of neutrino and antineutrino structure functions which coincides with the result for scattering on a free proton target. Once isospin effects are corrected for, neutrino structure functions on nuclear targets can still differ from those associated to free nucleons due to genuine nuclear modifications. These are typically quantified by nuclear modification ratios of the form which are equal to unity if the bound-nucleon structure functions equal the free-nucleon ones. In the DIS region we can express these nuclear modification ratios in terms of the nPDFs neutrons. Alternatively one can define , (C.6) which differ from unity due to both bound-nucleon modifications and due to the different valence quark content between the numerator and the denominator. for F 2 and xF 3 at Q = 10 GeV, where we show separately the ratios for the neutrino and antineutrino scattering as well as for their sum. Results are shown for two recent global nPDF determinations, nNNPDF3.0 (baseline in this work) and EPPS21 [112], which are in good agreement within uncertainties. In the case of nNNPDF3.0, we also display the nuclear modification ratio defined as in Eq. (C.4) where isospin effects are subtracted. The bands indicate the 90% CL intervals evaluated using the corresponding prescription for each nPDF set.
As compared to a proton target, in the valence region (x ∼ > 0.01) one finds that the F 2 and xF 3 structure functions in lead are enhanced (suppressed) for neutrino (antineutrino scattering), mostly as a consequence of the increased content of down quarks in lead as compared to the proton. These differences cancel out in the sum of neutrino and antineutrino structure functions due to Eq. (C.3), and there the ratio Eq. (C.6) reduces to Eq. (C.4), namely the genuine nuclear modifications of bound nucleons and compared to free nucleons. Bound-nucleon modification ratios, Eq. (C.4), are small in the valence peak region, except at rather large-x where structure functions are suppressed, and become more significant for F 2 in the sea quark region, x ∼ < 0.01, leading to a suppression of up to 25% as compared to the free-nucleon baseline. Therefore, from Fig. C.1 one concludes that in the valence region the dominant nuclear effects are the isospin-related ones, while in the sea region instead bound-nucleon modifications are important and need to be accounted for. These nuclear effects propagate to the inclusive cross section via Eq. (5.1), e.g. for neutrino scattering on lead one has an enhancement as compared to proton targets at low and intermediate energies and a suppression at high and ultra-high energies.
For completeness, we display in Fig. C.2 the corresponding nuclear modification ratios at the nPDF level in the case of the nNNPDF3.0 determination for both iron (Fe) and lead (Pb) targets. We compare the two definitions, Eqns. (C.5) and (C.6), differing in that isospin effects are accounted for in the former but not in the latter. As in Fig. C.1, the bands indicate the 90% CL intervals. In the valence region, one observes the expected suppression (enhancement) for up (down) quark nPDFs when using the ratio defined in Eq. (C.6). For isoscalar nPDF combinations such as the total quark singlet Σ and the gluon the two nuclear ratio definitions coincide. From Fig. C.2 we observe the quark shadowing in the sea region also reported by the structure function ratios in the case of lead target, while for iron targets in this region the nuclear ratios are consistent with unity within uncertainties.