Parton distribution function uncertainties in theoretical predictions for far-forward tau neutrinos at the Large Hadron Collider

New experiments to measure neutrinos in the far-forward region at the Large Hadron Collider (LHC) are under design or already in preparation. Two of them, FASER$\nu$ and SND@LHC, are expected to be active during Run 3 and have the potential to detect neutrinos that come from high-energy collisions in one of the LHC interaction points, extracted along the direction tangent to the beam line. Tau neutrinos and antineutrinos come predominantly from $D_s^\pm$ production in $pp$ collisions, followed by the leptonic decay of these mesons. Neutrino pseudorapidities in the range of $\eta>6.9$ and $\eta>8.9$ are relevant to these future experiments. At such pseudorapidities at high energies, theoretical predictions for the flux of tau neutrinos rely on parton distribution functions (PDFs) in a combination of very small and large parton$-x$ values. We evaluate PDF uncertainties in a next-to-leading order (NLO) QCD calculation of the flux of $\nu_\tau$ + $\bar{\nu}_\tau$ produced by $D_s^\pm$ decay in the far forward region at the LHC. The theoretical uncertainty associated with the 40 PDF sets of the PROSA19 group amounts to $\pm (20-30)$\% for the ($\nu_\tau$ + $\bar{\nu}_\tau$) number of charged-current (CC) events. Scale uncertainties are much larger, resulting in a range of CC event predictions from $\sim 70\%$ lower to $\sim 90\%$ higher than the central prediction. A comparison of the predictions with those obtained using as input the central PDFs from the 3-flavour NLO PDF sets of the CT14, ABMP16 and NNPDF3.1 collaborations show that far-forward neutrino energy distributions vary by as much as a factor of $\sim 2-4$ relative to the PROSA19 predictions at TeV neutrino energies. The Forward Physics Facility in the high luminosity LHC era will provide data capable of constraining NLO QCD evaluations with these PDF sets.

Abstract: New experiments dealing with neutrinos in the far-forward region at the Large Hadron Collider (LHC) are under design or already in preparation. Two of them, FASERν and SND@LHC, are expected to be active during Run 3 and have the potential to detect the interactions of ν andν that come from high-energy collisions in one of the LHC interaction points, extracted along the direction tangent to the beam line. Tau neutrinos and antineutrinos come predominantly from D ± s production in pp collisions, followed by the leptonic decay of these mesons. Neutrino pseudorapidities in the range of η > 6.9 and η > 8.9 are relevant to these future experiments. At such pseudorapidities at high energies, QCD theoretical predictions for the flux of ν τ plusν τ rely on parton distribution functions (PDFs) in a combination of very small and large parton−x values. We evaluate PDF uncertainties affecting the flux of ν τ +ν τ produced by D ± s decay in the far forward region at the LHC. Next-to-leading order (NLO) QCD uncertainties are included in the calculation of D ± s production and NLO PDF sets are used for consistency. The theoretical uncertainty associated with the 40 PDF sets of the PROSA19 group amounts to ±(20 − 30)% for the (ν τ +ν τ ) number of charged-current (CC) events. Scale uncertainties are much larger, resulting in a range of CC event predictions from ∼ 70% lower to ∼ 90% higher than the central prediction. A comparison of the predictions with those obtained using as input the central PDFs from the 3-flavour NLO PDF sets of the CT14, ABMP16 and NNPDF3.1 collaborations show that far-forward neutrino energy distributions vary by as much as a factor of ∼ 2 − 4 relative to the PROSA19 predictions at TeV neutrino energies. The Forward Physics Facility in the high luminosity LHC era will provide data capable of constraining NLO QCD evaluations with these PDF sets.

Introduction
Forward production of neutrinos at high-energy hadron colliders has been under discussion for several decades [1][2][3][4]. In particular, pp collisions at the Large Hadron Collider (LHC) interaction points can be regarded as sources of beams of neutrinos of different flavours, with an energy spectrum peaked in the hundreds of GeV -TeV energy range, the highest energies ever reached for a neutrino beam in a human-made accelerator. Measurements of the interactions of these neutrinos in detectors placed 100's of meters from the LHC interaction points are on their way to become a reality with recent proposals for LHC experiments, for example, by the FASERν [5,6], XSEN [7,8] and SND@LHC [9] collaborations. The FASERν and SND@LHC experiments were both approved and are expected to take data during Run 3, in the two service tunnels TI12 and TI18 located at a distance of ∼ 480 m from the ATLAS interaction point (IP), on opposite sides. Already a 29 kg pilot detector was placed in TI18 for four weeks in 2018. With 12.2 fb −1 of data from pp collisions at √ s = 13 TeV, neutrino interaction candidates were claimed to be observed [10], with muon neutrino interactions expected to dominate inside the sample. Beyond the pilot and the first phase of experiments, the idea of a forthcoming Forward Physics Facility (FPF), capable of hosting a suite of far-forward experiments active during the HL-LHC phase, has been recently raised, gaining increasing attention [11]. First studies in this direction are ongoing, considering the options of either enhancing the dimensions of one of the two aforementioned tunnels by enlarging it with alcoves, or building a new bigger cavern and access tunnel at approximately ∼ 500-600 m from the ATLAS IP [12]. In that cavern there would be enough space to host various experiments using different technologies at different distances from the beam-axis, significantly larger detectors, and with the completion of the HL-LHC phase, 20 times the integrated luminosity expected during Run 3. All these experiments will be sensitive to neutrinos and antineutrinos arising from the decay of particles and hadrons produced in the ATLAS IP, and propagating along the tangent to the accelerator arc. In this paper, we focus on the tau neutrino plus antineutrino flavour case. Differently from other neutrino flavours, tau neutrinos and antineutrinos are not, or are only very rarely, produced in pion and kaon decays because τ −ν τ final states are kinematically forbidden; K → πν τντ decays require flavour-changing neutral currents, and two body meson decays to ν τντ are suppressed by angular momentum considerations [13,14]. They are instead produced predominantly via D + s (and charge conjugate) production and its leptonic decay D + s → τ + ν τ . The tau neutrino plus antineutrino flux along the beam line will provide the opportunity to make direct tests of lepton flavour universality and to constrain new physics signals in tau neutrino oscillations over the considered baselines of several hundred meters [15]. In what follows, we refer to both neutrinos and antineutrinos generically as "neutrinos." Making a robust prediction of charm hadron production and the corresponding tau neutrino fluxes at very large rapidities can be regarded as a theoretical challenge [15,16]. The initial FASERν experiment is planned for η ν > 8.9 [6]. The XSEN proposal covers the neutrino rapidity ranges 7.4 − 8.1 and 8.0 − 8.6 [17], whereas the SND@LHC experiment explores the rapidity range 7.2−8.6 [9]. These rapidity values follow from space considerations and the limited dimensions of the already available caverns. The FPF could extend the range towards lower rapidities, to an extent that will depend on its dimensions, on the dimension of the experiments included in it and on the distance from the IP. In the aforementioned hypothesis of a new purpose-built cavern, it is foreseen that neutrino rapidities down to at least 6.5 can be covered with sufficient statistics. Lower rapidities down to ∼ 5.6 could be covered as well by some of the experiments, however at the price of reduced statistics, considering an area of ∼ 1 m 2 foreseen for the detectors that will be put off-axis.
On the other hand, charm hadron production has been measured by LHCb in rapidity intervals in the range 2.0 < y < 4.5 for different center-of-mass energies ( √ s = 5, 7 and 13 TeV) [18][19][20]. As we will see in the following, this ensures some overlap with the rapidity range seen by the far-forward experiments, considering that neutrinos at high-rapidity come from the decay of D-mesons from a range of rapidities, including lower rapidities with respect to the neutrino rapidity. Using next-to-leading order (NLO) QCD [21][22][23] and a phenomenologically motivated transverse momentum smearing function, theoretical predictions for open heavy-flavour hadroproduction have been compared with the LHCb data [20], to extract the optimal values of some of the input parameters to be used in the theoretical calculations. Theoretical calculations with input parameters constrained by the LHCb results, have then been used to make predictions in the whole rapidity range explored by far-forward experiments, including higher rapidities [15]. This can be regarded as an approximate extrapolation procedure. Very forward heavy-flavour production evaluated in the QCD-improved parton model probes both small-and large-x parton distribution functions (PDFs), where x is the fraction of longitudinal momentum of the proton that goes to the interacting parton. Theoretical uncertainties in the PDFs are particularly large in both these x regimes, considering that the bulk of HERA ep DIS data that form the backbone of PDF fits are limited to the x range 10 −4 x 10 −1 [24,25]. Some additional data concerning the F L structure functions have also been obtained at HERA. In principle, this would allow to extend the x-coverage down to x ∼ 10 −5 . However, due to the limited statistics of these data, they are either not used in the PDF fits or have a very small constraining power.
In this paper, we evaluate PDF uncertainties in the NLO QCD calculation of pp → D ± s X with decays to ν τ andν τ , using as a basis the 40 PDF sets of the most recent PDF fit delivered by the PROSA collaboration [26], an update including more sets of experimental data and updated theory input with respect to the first PROSA PDF fit [27]. In particular, this fit, besides HERA data, includes LHCb data on open-heavy flavour hadroproduction in pp collisions at center-of-mass energies up to √ s = 13 TeV, which has allowed to extend the x range down to x ∼ 10 −6 . These data have also an impact on the large x range (i.e., 10 −1 < x < 1). Additionally, it includes ALICE D-meson production data, with impact on more central x values, very useful to cross-check the constraints of DIS data + sum rules, playing a role in the same region. The PROSA fit has been first performed in the decoupling factorization and renormalization scheme, considering three active flavours, which is appropriate in the kinematical regime covered by the bulk of the LHCb and HERA charm production crosssections. Besides the fit in the decoupling scheme, a PROSA Variable Flavour Number Scheme PDF set has also been developed and released in Ref. [26]. For the sake of completeness, we observe that a large fraction of the LHCb D-meson production data included in the PROSA PDF fit have also been accounted for in independent PDF sets [28][29][30], built by applying a reweighting procedure on top of the NNPDF3.0 [31] and NNPDF3.1 [32] NLO PDF fits. The PROSA collaboration has explicitly checked the consistency, within uncertainties, of the gluon PDFs in the last version of their fit with the gluon PDFs in the most updated version of this independent PDF set, for factorization scales of the order of Q 2 ∼ 10 GeV 2 and low x values, as relevant for LHCb charm production. The outcome of these checks supports the robustness of these works and of the underlying heterogeneous methodologies. In the following, for the process we are interested in, namely forward ν τ andν τ production through D s production and decay in pp collisions at the LHC, the PDF uncertainties are compared with the uncertainties associated with renormalization and factorization scale variation. The energy distributions of ν τ +ν τ in different rapidity ranges are compared with those obtained using central sets of other NLO PDF fits as well [32][33][34][35], which do not include any constraint from charm hadroproduction data. We also evaluate the associated uncertainties affecting the ν τ N andν τ N charged current DIS cross sections, considering that in the foreseen far-forward experiments ν τ andν τ will be detected thanks to their DIS interactions with appropriate targets. This work aims at the twofold purpose of a) updating the predictions for the number of events in far-forward experiments presented in our previous work [15], by considering more reliable experimental setups, better motivated inputs for the theory calculations, and an improved description of fragmentation; b) providing for the first time an evaluation of PDF uncertainties affecting these predictions, comparing their size to those of other important QCD uncertainties.
The outline of this paper is as follows: in Section 2 the ingredients of the QCD theoretical calculation of charm hadron production are presented, with predictions including NLO radiative corrections compared to the available LHCb data at √ s = 13 TeV. We also make quantitative comparisons of charm production in different rapidity ranges, corresponding to different parton-x ranges. The tau neutrino plus antineutrino rapidity and energy distributions are presented in Section 3, with a focus on PDF uncertainties and the role of scale choice and transverse momentum smearing. Section 4 shows the number of events as a function of energy and in total for different rapidity ranges and detector masses, with an assessment of PDF uncertainties associated with neutrino interactions assuming a tungsten target. Finally, in Section 5 we draw our conclusions and discuss the outlook for future developments. Appendix A describes the procedure for combining PDF uncertainties from various sources. Appendix B reports numerical tables for the tau neutrino plus antineutrino energy distributions for several neutrino rapidity ranges, including QCD scale and PDF uncertainties, and tau neutrino and antineutrino charged-current cross sections per nucleon for a tungsten target. These tables are available in ascii format in the ancillary files of the preprint in the arXiv web archive 1 .
2 Charm hadron production at the LHC

Input factors
A detailed description of the inputs to the NLO QCD evaluation of charm hadron production and decay at the LHC is provided in Ref. [15]. For completeness, we summarize the main features here. Our QCD evaluation of charm production is performed at NLO using the one-particle inclusive results and formulas for the parton-level hard-scattering cross sections first published by Nason, Dawson and Ellis in Ref. [22]. Transverse momentum smearing is implemented phenomenologically by a Gaussian smearing function according to where k 2 T = 4 k T 2 /π is the effective transverse momentum squared that accounts both for intrinsic transverse momentum and as a phenomenological stand-in for higher order QCD effects. Our default is k T = 0.7 GeV [15], a value that makes our predictions for D-meson energy spectra very similar to those of the POWHEG [36]+PYTHIA [37] implementation used in Ref. [15].
We use the Peterson fragmentation functions for the c-quark to charm meson transition [38]. One approach to implement the fragmentation is to perform the reduction of the charm quark 3-momentum p Q (by p H = z p Q ) in the collider center of mass (CM) frame. Such an implementation generates a significant peak at y H = 0 in the hadron rapidity distribution, as can be seen in Fig. 1. This is caused by the larger mass of the heavy hadron compared to that of the heavy quark, m H > m Q , and also by the fragmentation which requires 1/z > 1. However, such a significant peak at y H = 0 is not found in the experimental observations. Alternatively, the fragmentation can be performed with the heavy hadron receiving fraction z of p Q in the CM frame of the initial state partons and then a boost is used to bring the heavy hadron to the collider CM frame. In such an implementation, the unwanted peak at y H = 0 is highly suppressed due to the further integral over y cm , the rapidity of the system of the two initial-state partons. The transverse momentum distributions are identical under these two implementations since the two reference frames are related by a longitudinal boost, and the transverse momentum does not change with such a boost, i.e., p HT = z p QT in both implementations.
A comparison of the energy distributions obtained using the two fragmentation implementations is similar to that of the rapidity distributions, but under the second implementation (in the parton CM frame), the energy of the hadron is a function of not only the energy of the heavy quark and z, but also the heavy quark rapidity and y cm . This makes the energy distribution of the hadron in the parton CM implementation of fragmentation closer to that of the heavy quark energy distribution than in the case of collider CM implementation of fragmentation. Of course, the distributions obtained in the two implementations are all identical in the massless limit. In this work, we use the implementation of fragmentation in the parton CM frame. Implementation in the parton CM frame more closely resembles the frame where fragmentation functions are extracted in e + e − collisions.
The parton level differential cross sections at NLO are convoluted with PDFs for pp collisions. In this work, the default set of NLO PDFs, in a decoupling scheme with three active flavours consistent with the matrix elements we use, is provided by the PROSA collaboration [26] and implemented through the LHAPDF interface [39]. The PROSA19 PDFs (in the following we will rename PROSA19 simply as "PROSA") result from a fit, including, among others, data on open heavy flavour production from HERA, LHCb and ALICE, as already mentioned in the Introduction. Forty PDF eigenvectors in addition to the best fit account for uncertainties associated to the PDFs. They account for i) experimental uncertainties ("fit" uncertainties), ii) model assumptions in performing the fit, including, among others, µ R and µ F variations, as well as variation of the parameters entering fragmentation functions and fractions and of the α s (M Z ) value, iii) parameterization uncertainties related to the functional form of the PDFs at the starting scale for evolution and the parameters appearing there. The prescription for combining the PROSA PDF eigenvectors in order to produce uncertainty bands in kinematic distributions is described in Appendix A. For comparison, we also show the central predictions obtained with other 3-flavour NLO PDFs: the CT14nlo [33], ABMP16 [34] and NNPDF3.1 [32] 3-flavour sets, which have already been used for computing heavy-flavour production in previous works (see e.g. Ref. [40,41]). Differently from the PROSA PDF sets, these fits do not account for any flavour hadroproduction data.
For computations with each PDF set, we use as input the associated α s (M Z ) value, and α s evolution at two-loops. The pole mass m c , also input to the NLO calculation, is related to the MS mass at the mass scale µ = m c (µ) through four loops by the relation [42][43][44] where α s ≡ α s (µ). While converting heavy-quark masses from one mass renormalization scheme to the other, we include the first correction in our NLO evaluation. In particular, for the PROSA PDFs, the MS charm mass is 1.242 GeV, a value that comes from the simultaneous fit of the MS heavy-quark masses and the PDFs performed by the PROSA collaboration in Ref. [26]. This translates into a pole mass of 1.442 GeV. Table 1 shows the charm quark mass values for the various PDFs used in this work. For the renormalization scale µ R and factorization scale µ F , we adopt two different functional forms, the scale used for the PROSA 2019 fit to heavy flavour production and the transverse mass which is more widely used as central µ R and µ F in many computations of heavy-quark pair production. The use of the modified scale m 2 T,2 finds its justification in the fact that it reduces the NNLO/NLO K-factor in the bulk of the phase-space for cc production, with respect to the other option. As we will show in the following Section, we verify that NLO scale uncertainties based on m 2 T,2 are also reduced in the kinematic region relevant for very forward heavy-flavour production compared to those obtained with m 2 T,1 , and that the uncertainty band is more symmetric around the predictions from the central scale. In our characterization of scale uncertainties, we consider variations around the central scale (µ R , µ F ) = (1, 1)m T,2 , related to scale combinations {(0.5, 1), (1, 2), (1, 0.5), (2, 1), (0.5, 0.5) and (2, 2)} m T,2 , by building the envelope of the corresponding predictions [45]. On the other hand, the PROSA PDF uncertainties from the 40 sets are determined by fixing (µ R , µ F ) = (1, 1) m T,2 , assuming the independence of PDF and scale uncertainties, that can then be added in quadrature, as also suggested by the PROSA collaboration in Ref. [26,27]. Details of combining uncertainties associated with the 40 PDF sets are included in Appendix A. These studies are all performed with k T = 0.7 GeV. For comparison, we also show predictions with the scale choice (µ R , µ F ) = (1, 2) m T and k T = 1.2 GeV, parameters that lead to predictions which compare well with the LHCb data at √ s = 13 TeV.
The evaluation of D + s → ν τ decays is straightforward as the relevant process is the twobody decay D + s → τ + ν τ , with branching fraction B = 0.0548 [46]. The neutrino in this two-body decay is denoted as "direct" neutrino. The τ lepton also decays. Its neutrino is called the "chain" neutrino. Details of the implementation of D + s → τ + ν τ and τ decays appear in, for example, Ref. [15].
Detector size and positioning along the beam axis (the z-axis) lead to very small angles θ which translate to neutrino pseudorapidities η ν = y ν = − ln tan(θ/2). Since the neutrino rapidity and pseudorapidity are identical, we use the terms interchangeably for neutrinos. The distinction between rapidity and pseudorapidity is important for charm quark and D ± s distributions, as shown in Fig. 2. The upper panel of Fig. 2, for charm quarks, shows that it is the charm rapidity y c rather than its pseudo-rapidity η c that correlates more closely with the neutrino η ν for the large η ν values most relevant to the forward LHC neutrino experiments. For example, for a charm rapidity window in the range y c = 6.75 − 7.25, the peak of the resulting neutrino distribution (green dotted line) is at η ν = 6.7. By contrast, for a charm pseudorapidity window in the range η c = 6.75−7.25, the corresponding η ν distribution (green solid line) peaks at the lower value η ν =6.3.
The displacements in the location of the η ν peaks with respect to the center of the charm pseudorapidity windows are even larger for angles closer to the beamline. For η c = 8.75−9.25, the peak of the resulting neutrino distribution is η ν = 7.0, while for y c = 8.75 − 9.25, the peak of the neutrino distribution is closer to the charm rapidity window, at η ν = 8.3. The same conclusions apply for the correlations between neutrino rapidities and D s meson rapidity and pseudorapidity windows, as shown in the lower panel of Fig. 2. The neutrino rapidity is better correlated with the meson rapidity than the meson pseudorapidity.
On the other hand, neutrinos with a fixed (pseudo)rapidity receive contributions from decays of charmed mesons with a range of rapidities, including even rapidities y Ds < η ν . For example, in the lower panel of Fig. 2, the red dotted histogram of neutrino η ν distribution extends to larger values than the interval of rapidities y Ds of the parent meson, 8.75-9.25. These kinematics considerations have profound implications for the design of the far-forward neutrino experiments and the possibilities of constraining unknown QCD aspects through their data. In any case, we always use the full kinematics of charm production, fragmentation and decay to direct and chain neutrinos in our evaluation of the neutrino energy distributions for several η ν ranges.

Charm hadron production
We begin with the cc production cross section. Measurements of this cross section have been extracted from LHCb data for D 0 , D + and their charge conjugate mesons produced with 2.0 < y < 4.5 for √ s = 5, 7 and 13 TeV. Accounting for the respective fragmentation fractions for D 0 and D + , the LHCb collaboration extracted the corresponding σ(pp → ccX) in this restricted kinematic regime [18][19][20], shown in Fig. 3 by the black dots with error bars. Our NLO QCD evaluation of the equivalent quantity is performed with charm production and fragmentation to D + , with the D + momenta satisfying the LHC analysis cuts. Normalized by the inverse of the fragmentation fraction to D + , our predictions are shown in Fig. 3 for scales µ R , µ F ∝ m T,2 (orange marker and error bars) and for µ R , µ F ∝ m T,1 = m T (filled green marker and error bar). The figure illustrates the fact that the scale uncertainties are somewhat  [8.75, 9.25]. Lower: the distribution of y ν = η ν for ν τ +ν τ that come from D ± s decays where the D ± s rapidities y D ± s (dotted histograms) and D ± s pseudo-rapidities η Ds (solid histograms) lie in the ranges [2.75, 3.25], [6.75, 7.25] and [8.75, 9.25]. The distributions have unit normalization in each D ± s (pseudo)rapidity interval. They refer to D ± s produced in pp collisions at √ s = 14 TeV, using as input (µ R , µ F ) = (1, 1)m T,2 and k T = 0.7 GeV. smaller for µ R , µ F ∝ m T,2 compared to µ R , µ F ∝ m T,1 = m T . The central predictions for both central scale choices lie below the black data points from the LHCb data from charm mesons in the same transverse momentum and rapidity range. For reference, we also show with the green open marker the corresponding cross sections for (µ R , µ F ) = (1, 2)m T , which lie above the measured cross sections. The feature that central theory predictions with our default central scale choice lie somewhat below the data, although theoretical predictions and experimental data are still compatible within the uncertainties, is also apparent in the LHCb transverse momentum distributions. Figs. 4 and 5 compare our predictions on the double differential cross sections for D ± s and D 0 +D 0 production in p T and y with the corresponding LHCb data [19,20]. The colored bands in the histograms show the uncertainty in the prediction of d 2 σ/dp T dy associated with the 7-point scale variation around (1, 1)m T,2 . The LHCb data generally lie within the scale uncertainty band of theory predictions and have smaller experimental uncertainties than the latter.
For reference in Fig. 6, in addition to predictions with the central scales (µ R , µ F ) = (1, 1)m T,2 and k T = 0.7 GeV (solid histogram), we show d 2 σ/dp T dy with (µ R , µ F ) = (1, 2)m T and k T = 1.2 GeV (dashed histograms). These inputs yield histograms for the D ± s and D 0 +D 0 distributions that lie closer to the LHCb data than our default scale and k T choices. Using (µ R , µ F ) = (1, 2)m T and k T = 1.2 GeV as default input choices, however,  The transverse momentum distribution in the D ± s rapidity ranges indicated in the figure, for production of D + s or its charge-conjugate meson in pp collisions compared with LHCb data for dσ(D + s )/dp T + dσ(D − s )/dp T at √ s = 7 TeV [19] (upper) and for d 2 σ(D + s )/dp T dy + d 2 σ(D − s )/dp T dy at 13 TeV [20] (lower). The ∆y bins are shifted by 10 −m where m = 0, 2, 4, 6 and 8. The central scale is set to (µ R , µ F ) = (1, 1)m T,2 in theory predictions and k T = 0.7 GeV. The colored bands show the 7-point scale variation uncertainty around the solid black histograms for the central scale.   would suffer from some theoretical drawbacks in the context of the current study of scale uncertainties. In particular, the choice of (µ R , µ F ) = (1, 2)m T as central scale leads to the inclusion of scales µ F /µ R = 4. This large ratio of scales leads to negative cross sections at very low p T , indicating a breakdown of perturbative QCD at fixed relatively low (NLO) order, as considered here, in the standard factorization scheme. Furthermore, the larger value of k T = 1.2 GeV is more difficult to interpret as an intrinsic k T , considering that we expect the latter to be less than the scale of the proton mass. However, we note that large values of k T are suggested in other processes [47][48][49] and may reflect the importance of higher-order QCD corrections, somehow mimicked by the use of high k T values. With our default scale and k T values, Fig. 7 shows the uncertainty band associated with the variation of the PROSA PDFs, on p T distributions for D ± s mesons for the two LHCb more extreme rapidity ranges. 2.0 < y < 2.5 and 4.0 < y < 4.5. The PROSA PDFs include 40 different sets representative of fit, model, and parameterization uncertainties, in addition to the central set, corresponding to the best-fit. The 40 variations included in the PROSA PDF fit bring in uncertainties in the differential cross section of D ± s production within about ± 15%, with the largest deviation in the low p T region, especially in the forward direction. The dominant contribution to the uncertainties come from the so-called PROSA "model" uncertainties, which involve, among others, the uncertainty on various theory inputs for heavy quark production used in the fit (see Ref. [26] for more detail).
Other PDF fits in the 3-flavour scheme are available. We compare the PROSA PDF sets with the NLO sets of the CT14 [33], ABMP16 [34] and NNPDF3.1 [32] collaborations. Particularly important for this work and heavy-flavour production in general, are gluon PDFs. The small-x and large-x behavior of the gluon PDFs is shown in the upper and lower panels of Fig. 8, respectively, where the PROSA PDF best-fit results (black curves) are plotted together with the related uncertainty bands (orange). Additionally the best-fit gluon distributions from the other aforementioned PDF fits are shown. Beneath each panel ratios with respect to the central PROSA PDF are also shown. The upper panel is limited to the x range 10 −8 < x < 0.3, whereas the lower panel focuses on x > 0.3. Both ranges are interesting because far-forward production of charm quarks at large center-of-mass energies involves the product of small-x and large-x PDFs. The lower panel shows large deviations of the various PDF best-fits among each other for the largest x values (note the different scale in the y axis, when comparing the ratios of the upper and the lower panels).
To illustrate the regions of the partonic longitudinal momentum fractions (x 1 , x 2 ) that contribute to charm production in pp collisions, the upper panel of Fig. 9 shows the range of involved (x 1 , x 2 ) values, depending on the rapidity of the produced charm quark in pp → cX. For y c > 0, the charm quark has a momentum component in the +z-axis direction, so x 1 > x 2 where x 1 is the parton momentum fraction of the involved parton in the proton in the beam traveling in the +z-axis direction. The blue region in the upper panel and the colored regions to its right, show that combinations of partons with x 1 ∼ > 0.04 and x 2 ∼ > 4 × 10 −8 contribute to the production of charm with y c > 6. The (x 1 , x 2 ) ranges are peaked on more extreme values when the minimum charm rapidity increases. For example, for y c > 9, x 1 ∼ > 0.8 and x 2 ∼ > 10 −7 . On the other hand, as we saw in Fig. 8, the PDF vary widely for large values of x.
The lower panel of Fig. 9 shows the region of (x 1 , x 2 ) that contribute to energy distributions above a given minimum charm energy. Since there is not a rapidity restriction in the lower panel, there is not a condition on which parton momentum fraction is greater.   Figure 6: The double-differential distribution d 2 σ/dp T dy for D + s + D − s (upper) and D 0 +D 0 (lower) production at √ s = 13 TeV [20], with input (µ R , µ F ) = (1, 1)m T,2 and k T = 0.7 GeV (solid) and (µ R , µ F ) = (1, 2)m T and k T = 1.2 GeV (dashed), compared with LHCb data. The distributions are scaled as in Fig. 4.
The dotted boxes in both panels of Fig. 9 show the range 10 −6 < x < 0.3 in which the considered PDFs either are best constrained or better agree among each other. At large x, the gluon PDF is constrained by inclusive jet, dijet and tt distributions [32,[50][51][52][53][54][55]. However, gluon PDF uncertainties remain large, in part due to tensions between different data sets. For example, for the tt case, this issue is discussed in, e.g., refs. [53,56]. In any case, using the PROSA19 PDFs, y(tt) data from CMS at √ s = 8 TeV [57] are reasonably well reproduced for a range of tt invariant masses (χ 2 = 23/15), slightly better than with HERAPDF2.0 (χ 2 = 24/15) [58]. Different PDF sets will most diverge from each other for charm produced at large rapidities and/or at high energies.
To uncover the effect of large x PDFs, we set k T = 0 and integrate over all rapidities. The impacts of large-x differences in the PROSA, ABMP16, CT14 and NNPDF3.1 PDFs on the evaluations of the charm quark energy distributions with these fits are shown in Fig. 10. We note that predictions for the various PDF sets are evaluated with the corresponding charm quark mass values from Table 1. When x 1 , x 2 < 0.3, the kinematic limit to the charm quark energy is ∼ 2.1 TeV. As shown in the upper panel of Fig. 10, the charm energy distributions from the four central PDF sets differ by at most 50% at large E c , and less at lower energies. The dashed purple histograms show the charm energy distribution using PROSA PDFs and m c = 1.3 GeV rather than 1.442 GeV, showing an increase up to ∼ 25% at low energy, and much lower at the highest energies, in the energy distribution. The distribution based on the NNPDF3.1 PDFs, computed using the charm mass value accompanying this PDF fit, i.e., m c = 1.51 GeV, turns out to be the lowest.
The lower panel of Fig. 10 shows the charm energy distribution with integration over the full ranges of x 1 and x 2 . As anticipated, the energy distributions differ significantly for E c greater than a few TeV. The ratios of the deviations show a qualitatively similar behavior to the ratios of the large-x gluon PDFs shown in Fig. 8. Figs. 9 and 10 show that predictions for very high charm quark rapidities (y c ∼ > 8) have significantly larger deviations derived from different PDF sets than predictions for lower charm rapidities. This is reflected in the wider range of predictions for FASERν than for SND@LHC, for example, as discussed in Section 4.2. We note that in Fig. 10 the energy distribution of the charm quark is considered. When the charm quark is fragmented and the meson decays to tau neutrinos, the predictions from different PDFs start to differ already for lower neutrino energy values, below 2 TeV, as we show in the next section.

Tau neutrino and antineutrino production from charm
We now turn to tau neutrino and antineutrino production from D ± s production and decay, including both the direct neutrinos D s → ν τ and chain decay neutrinos D s → τ → ν τ . We observe that data on ν τ andν τ are not available in standard LHC experiments, for which all neutrinos produced in an event just contribute to the event as missing energy. We assume a branching fraction B(D s → τ ν τ ) = 0.0548 [46]. Our perturbative evaluation produces equal numbers of D + s and D − s . The energy distributions of left-handed ν τ 's are equal to the energy     distributions of right-handedν τ 's. All of the predictions shown in this section are for the sum ν τ +ν τ , and are generically referred to as "neutrinos." We begin with the neutrino rapidity distributions. Fig. 11 shows the PROSA PDF uncertainties for the η ν distribution of the tau neutrinos from direct (D s → ν τ , upper left panel) and chain decay (D s → τ → ν τ , upper right panel). The two contributions to the neutrino rapidity distributions show similar PDF uncertainties, in the range of ±20 − 30% in the ratio to the central PDF. The PROSA model uncertainty dominates. The lower left plot of Fig. 11 shows the ratio of the η ν distributions of the direct and chain neutrinos. Their contributions are nearly equal across the considered rapidity range η ν > 6.9, with a maximum difference of very few percent at the largest η ν .
The lower right plot of Fig. 11 shows the η ν distributions from the sum of direct plus chain contributions and the scale uncertainty. For reference, the dashed line is just the direct contribution. The scale uncertainties, again, are significantly larger than the total PROSA PDF uncertainties. The scale uncertainty envelope for the rapidity distributions range between the lower edge of the envelope that is ∼ 25% of the central value (a correction of ∼ −75% to the central value) to close to a factor of 2 times the central scale result for the upper edge of the scale uncertainty envelope, nearly independent of energy.
At high rapidity, the lower right plot shows that the rapidity distributions for different scale combinations have a common behavior (i.e., the same shape) as a function of rapidity, and they only differ among each other for the normalization. The dotted line corresponds to the analytical formula dσ This represents the central scale histogram values to within ±5% for η ν > 8.3 in Fig. 11, but eq. (3.1) lies above the histograms in the lower right panel for η ν < 8.3. The scaling behaviour at large η ν appears to be a universal feature, independent of N R and N F for (N R , N F )m T,2 scales when η ν > 8.3. We probe the functional behavior of dσ/dη ν ∼ e −2ην for η ν > 8.3 with the following considerations. We consider a cylindrical detector at a distance D d from the interaction point, aligned and centered along the beam axis (z-axis), "on-axis," for which the minimum detectable η ν is labelled as η 1 . Given the relation between the pseudorapidity and angle θ relative to the z-axis, the cross sectional area of the detector is where Ω ν (η 1 ) ≡ 4πe −2η 1 is the approximate solid angle enclosed by a circle around the z-axis for which the angle θ relative to the z-axis corresponds to η 1 and is a distance D d from the interaction point. The functional form of dσ/dη ν in eq. (3.1) is therefore related to Ω ν (η ν ). We will come back to this point in our evaluation of the number of events per ton of target (see Section 4.2), which depends on the transverse area and the depth of the target. The scaling of dσ/dη ν with A d (η ν ) for large η ν has already been noted in Ref. [59]. For a detector on-axis of radius 1 m at a distance of 480 m from the interaction point, η ν > 6.87. We show results for η ν > 6.87 and other higher neutrino rapidity ranges corresponding to SND@LHC (7.2 < η ν < 8.6), FASERν (η ν > 8.9) and the range of 8.0 -9.2. For the energy distributions of tau neutrinos whose rapidity is restricted to be higher than η min = 6.87 6.9, the PROSA PDF uncertainties gradually increase from low tau neutrino energies up to the 2 TeV energy range, as illustrated in Fig. 12. The black histogram and yellow band reflect the central PROSA PDF predictions and PDF uncertainty bands, respectively. Relative to the central result, the PDF uncertainty band has approximately the same shape as a function of      energy for all the η ν ranges shown. As also in case of the rapidity distributions, the ABMP16 and NNPDF3.1 predictions lie within the PROSA PDF uncertainty band. However, in the E ν range of 1-2 TeV, they tend towards the upper edge of the PDF uncertainty band. Again, the CT14 predictions are higher. The energy distributions obtained with the CT14 PDFs   are larger by almost a factor of ∼ 4 for η ν > 8.9 at high energy. Only the PROSA PDF uncertainty bands are shown in Fig. 12. One expects that the uncertainty bands for the other PDFs will partly overlap with the PROSA bands. While the considered PDFs are compatible for a range of partonic x values, Fig. 12 illustrates that for tau neutrino energy distributions at high rapidities, predictions from different PDFs may vary widely. Fig. 13 shows that the scale uncertainties for the same rapidity ranges as in Fig. 12 are largely independent of energy and much larger than the PROSA PDF uncertainties. Consistent with the neutrino rapidity distribution scale uncertainties, neutrino energy distributions corresponding to different scale combinations as obtained in the scale variation procedure range between ∼ 25% to a factor of ∼ 2 compared to predictions with the central scale.
In addition to our default scale dependence on m T,2 , Fig. 13 shows the neutrino energy predictions using as input (1.0, 2.0)m T,1 with k T = 1.2 GeV, a choice that better matches the LHCb experimental charm meson transverse momentum distributions. With these inputs, the upper two panels show that the corresponding histogram lies at the upper edge of the scale uncertainty band below ∼ 1 TeV, and, for rapidities below ∼ 8, is above the scale uncertainty band for tau neutrino energies above 1 TeV, as we discuss below.
Tables 5-7 in Appendix B list, for various rapidity ranges, the central values and uncertainties of the predictions for energy distribution obtained with input (1, 1)m T,2 , k T =0.7 GeV and the PROSA PDFs, corresponding to the plots of Fig. 12 and 13. Tables 8-10 list the predictions for scales (1, 2)m T with k T = 1.2 GeV and the PROSA PDFs, and the NNPDF3.1, CT14 and ABMP16 predictions with (1, 1)m T,2 and k T = 0.7 GeV. The Tables are available as ancillary files on the arXiv link to this paper. The upper panel of Fig. 14 shows the ratio of energy distributions for k T = 1.2 GeV to k T = 0.7 GeV, for the scales fixed at (µ F , µ R ) = (1, 1)m T,2 . The k T smearing moves very low p T mesons to higher p T , thereby increasing the distribution in the high-p T tails up to 10 GeV and larger at LHCb, as shown in Ref. [15]. At high energies, At very forward rapidity, the larger of the two k T values considered here, namely k T = 1.2 GeV, pushes a fraction of the differential cross section outside of the allowed rapidity region, thereby decreasing the neutrino energy distribution relative to the distribution evaluated with a smaller k T . For η ν > 8.9, the ratio of the energy distributions dσ( kT = 1.2 GeV)/dE ν to dσ( kT = 0.7 GeV)/dE ν shown in the upper panel of Fig. 14 is about 0.75 for E ν = 2 TeV. For η ν > 6.9 and 7.2 < η ν < 8.6, neutrinos with higher p T are detectable. The ratio   Figure 14: Ratio of the tau neutrino plus antineutrino distribution dσ/dE ν for η ν > 6.9, 7.2 < η ν < 8.6 and η ν > 8.9. The upper plot has the ratio with m T,2 scale dependence and k T = 1.2 GeV (numerator) and 0.7 GeV (denominator). The middle plot has k T = 1.2 GeV fixed with scale dependence (1, 2)m T (numerator) and (1, 1)m T,2 (denominator). The lower plot has both k T and scale differences in the ratio. of energy distributions increases by ∼ 25% for E ν = 2 TeV for k T = 1.2 GeV compared to k T = 0.7 GeV.
A larger effect on the neutrino energy distribution is generated by the difference in scale choice, for a fixed k T . The middle panel of Fig. 14 shows the ratio of the predictions for the neutrino energy distribution obtained with the (1, 2)m T,1 and (1, 1)m T,2 scale choices, for a fixed k T = 1.2 GeV. The ratio of dσ(µ R = m T , µ F = 2m T )/dE ν to dσ(µ R = µ F = m T,2 )/dE ν in the three different rapidity ranges largely overlap. This means that the effect of changing the renormalization scales is largely rapidity independent. At low energy, dσ/dE ν evaluated with (1, 2)m T,1 is ∼ 1.6 time larger than with (1, 1)m T,2 . The ratio between the energy distributions in the middle panel of Fig. 14 increases to ∼ 1.9 as E ν → 2 TeV. Specifically for the two scale choices: (1, 2)m T,1 and (1, 1)m T,2 , for the charm transverse momentum near p T → 0, m T,2 2m T,1 , so the factorization scales are approximately equal and the renormalization scales differ by a factor of 2. For higher p T the renormalization scales are nearly the same, whereas the PDFs are evaluated at factorization scales that differ by a factor of ∼ 2.
The lower panel of Fig. 14 shows the ratio of the neutrino energy distributions for different rapidity ranges obtained with (1, 2)m T,1 and (1, 1)m T,2 , in association with k T = 1.2 GeV and 0.7 GeV, respectively. The trends of the ratios follow the dependence on the ratio at fixed scale for different k T , scaled by the roughly energy independent ratio of differential cross sections evaluated at two scales and fixed k T value. For η ν > 6.9 and 7.2 < η ν < 8.6, the ratio increases with neutrino energy, whereas for η ν > 8.9, the ratio is nearly energy independent.
4 Tau neutrino and antineutrino charged current events

Interaction cross sections
PDF uncertainties also apply to the tau neutrino and antineutrino charged current (CC) cross sections. However, their size is much smaller than in case of charm production due to the fact that in the characteristic ranges of x and Q 2 m 2 c values relevant for the calculation of the CC cross sections, the PDF fits are very well constrained by the already available experimental data. Fig. 15 shows the tau neutrino and antineutrino cross sections per nucleon for charged current interaction in a tungsten target, scaled by incident (anti-)neutrino energy. For the interaction cross sections, we use the nCTEQ15 [60] as a default PDF set as in Ref. [15], here for tungsten, and evaluate the deep inelastic scattering (DIS) cross sections at NLO in QCD including target mass and tau lepton mass corrections [61][62][63]. For Q 2 < 2 GeV 2 , we extrapolate the neutrino and antineutrino structure functions with a prescription outlined in Ref. [64] that is based on the Capella et al. parameterization [65]. This prescription has low-Q 2 behavior that is similar to the Bodek-Yang PDF dependent prescription [66][67][68][69][70].
The target mass corrections are small, even for tau antineutrino scattering, for energies above 25 GeV [62]. The tau mass correction accounts for a ∼ 25% suppression of the tau neutrino cross section relative to the muon neutrino cross section for 100 GeV incident energies, reducing to ∼ 5% tau mass effect at 1 TeV [63]. The tau neutrino and antineutrino CC cross section for incident neutrino energies below 10 GeV do not make significant contributions to the total number of events over all energy spectrum. We therefore neglect quasi-elastic and resonant scattering contributions, which are smaller than the DIS contribution for E ν ∼ > 10 GeV. For neutrino energies above 100 GeV, the resonance region accounts for only a few percent of the total cross section [12,71]. We find that the contribution from Q 2 < 1 GeV 2 to the cross section at 100 GeV is less than 3% (5%) for ν τ (ν τ ) scattering with our low-Q 2 extrapolation of the structure functions [12,71]. The fact that NLO QCD corrections to the neutrino interaction cross section for most of the energy range of interest are small [63] means that our neglect of factorization and renormalization scale dependence of the neutrino cross section is reasonable. The bands in the figure display the uncertainty due to the 32 different PDF sets of nCTEQ15 for tungsten, which are obtained according to eq. (A.6). Deviation of the charged current cross sections due to the nCTEQ15 PDF uncertainty from the one with the central PDF set is about 2-3 % for E ν ∼ > 100 GeV for both tau neutrinos and antineutrinos. Table 11 in Appendix B lists the ν τ andν τ CC cross sections per nucleon, including PDF uncertainties from the nCTEQ15 PDFs, for scattering with a tungsten target.
For comparison, we also present the cross sections evaluated with the PROSA PDF set in the variable flavour number scheme (VFNS) [26]. We use the VFNS version to be consistent with the nCTEQ15 PDFs, extracted in the same scheme, and considering the kinematical range of Q 2 values we are mostly interested in. For example, for E ν = 100 GeV, the ν τ charged current DIS cross section on nucleons has Q 2 = 23 GeV 2 , while forν τ at the same incident energy, Q 2 = 13 GeV 2 [72]. Most of the CC events come from higher energies where Q 2 is even higher. The PROSA collaboration so far only fitted proton PDFs, and not yet nuclear PDFs. In that case, we obtained the neutron PDFs according to isospin symmetry. Only shown, with the dashed curves, are the predictions for the tau neutrino and antineutrino CC cross sections per nucleon for interactions with tungsten with the best fit of the PDF since other PROSA PDF sets necessary for evaluating uncertainty are not provided in the VFNS. Table 12 in Appendix B lists the ν τ andν τ CC cross sections per nucleon for the PROSA VFNS PDFs. The resulting cross sections with the PROSA VFNS PDFs are smaller than the central results with nCTEQ15 by ∼ 3% for tau neutrinos and ∼ 7-11% for antineutrinos in the energy range of interest, 100 ∼ < E ντ (ντ ) /GeV ∼ < 2000, hence not in the uncertainty band of the nCTEQ15, although one could expect that, in case an uncertainty band would accompany the PROSA VFNS PDF fit as well, this would overlap with the nCTEQ15 uncertainty band. The discrepancy between the central predictions with the two PDF best-fits is due to the combined effect of the differences between the PDFs and the way nuclear effects are incorporated. The largest impact comes from the difference between the nCTEQ15 and PROSA PDFs in the large x region, where the PDFs are not well constrained by the data. In Fig. 16, we present the fractions of the CC cross section according to the minimum values of parton-x in the evaluation. Fig. 16 illustrates that the x range extends to lower x values as the neutrino energy increases. We find that for E ν = 2 TeV, 68 (58)% of the contribution to the CC tau neutrino (antineutrino) cross sections comes from phase-space configurations characterized by x > 0.1, and 98 (96)% is from x > 0.01.

Charged-current event distributions in the forward detectors
In this section, we investigate the expected number of the tau neutrinos and antineutrinos events for the next stages of the LHC, Run 3 and High Luminosity LHC (HL-LHC). In our evaluation here, we consider contributions only from D ± s meson production and decays, which yield most of the tau neutrino and antineutrino events. The contribution from B mesons was investigated in Ref. [15]; these amount to ∼ 4 % of total events.
As in Ref. [15], we calculate the number of neutrinos that can be detected at the LHC using with the integrated luminosity L = 150 fb −1 and 3000 fb −1 for Run 3 and HL-LHC, respectively. All of our results presented above are integrated over 2π azimuthal angle. Here, we use the interaction probability in the detector given by assuming that detectors are made of tungsten, thus the mass number is A W = 184, and the density ρ W = 19.3 g/cm 3 . We show results for a cylindrical detector aligned with the beam axis and capable of detecting neutrinos with η ν > η 1 , constituted by a fixed mass of tungsten equal to 1 ton (M 0 = 10 6 g). The mass can be written in terms of the solid angle Ω ν (η) defined in Eq. (3.2), the distance D d from the interaction point and the length L d of the cylinder, according to the following relation, L d depends on η 1 , since η 1 is related to the radius of the detector placed at a distance D d from the IP. When looking in the rapidity range η 1 < η ν < η 2 , the relation is: Here, the length L d (η 1 , η 2 ) is related to the inner and outer radii of the cylindrical detector through the corresponding detectable rapidity range. For D d = 480 m and detection of η ν > η 1 = 8.9 in the full 2π azimuthal angle around the z-axis, L d = 1.06 m for one ton of tungsten. For η ν > 6.9, the corresponding detector depth is L d = 0.016 m, whereas for 7.2 < η ν < 8.6, L d = 0.034 m. For a fractional azimuthal coverage ∆φ, the detector depth for 1 ton increases, L d → L d · 2π/∆φ. By including a factor 1/M 0 , the number of events per unit energy per detector mass at a distance D d from the interaction point, under the neutrino pseudorapidity cuts incorporated in dσ/dE ν to select forward neutrinos, can be written as where we have suppressed the dependence of dσ/dE ν and Ω ν on η 1 (or (η 1 , η 2 ) where relevant). For large η 1 (η 1 8.3), the approximate solid angle scaling of the rapidity dependent differential cross section shown in eq. (3.1) approximately cancels the rapidity dependence of the solid angle to yield similar numbers of events per unit detector mass for each rapidity range at large enough rapidity [59]. The number of events per unit mass for η ν > 8.9 is larger than the number of events per unit mass for η ν > 6.9, since for η ν 8.3, the rapidity dependent differential cross section lies below the rapidity scaling curve (the dotted line in the right panel of Fig. 11).
On the other hand, in the case of two detectors of equal thickness L d of the same material, looking respectively to pseudorapidities η > η 1 and η > η 2 , with η 1 η 2 , the total number of events will be much smaller in detector 1 (η > η 1 ) with respect to detector 2 (η > η 2 )(see Table 2) due to both the much smaller dσ/dE ν , e.g., as shown in Fig. 12, and the smaller mass of detector 1. In designing a far-forward experiment, both the transverse and the longitudinal size of the spaces to host them need to be considered.

During the Run 3 stage
The LHC is scheduled to operate Run 3 from 2022 to 2024. During this stage, two new experiments, FASERν and SND@LHC will be conducted to detect neutrinos in the large rapidity region. The FASERν detector is made of tungsten and emulsion with a size of 25 cm × 25 cm × 1.3 m, and the target mass of 1.2 ton [73]. Considering that the baseline is 480 m, the FASERν experiment can probe the pseudorapidity range of η 8.9. SND@LHC is designed to cover the pseudorapidity range 7.2 η 8.6 and uses a tungsten, emulsion, and scintillating fiber detector with target mass of 830 kg [9]. While the rapidity ranges are not precise due to the fact that the detector shapes might differ from simple cylinders, we use these representative ranges for all the results presented here. Fig. 17 shows the number of tau neutrino and antineutrino events per energy per ton of tungsten detector for the two aforementioned rapidity cuts, considered for the plots in the upper and lower panel, respectively. The presented predictions were obtained with our default scale choice, (µ R , µ F ) = (1, 1) m T,2 and k T = 0.7 GeV. The shape of the energy distribution of events reflects those of contributions from direct neutrinos at low energies and from chain decay neutrinos at high energies. In our perturbative QCD evaluation of charm pair production, using a fragmentation function c → D + s andc → D − s , the ν τ andν τ energy distributions are identical. However, the resulting number of events is larger for tau neutrinos than for tau antineutrinos, since the tau neutrino charged-current deep-inelastic-scattering (DIS) cross sections are about two times larger than the antineutrino cross sections [63]. The fact that the number of events per unit energy per ton is larger for η > 8.9 than for 7.2 < η ν < 8.6, as discussed above, follows from the deviation from solid angle scaling of the neutrino rapidity distribution for η ν < 8.3. In the lower right panel of Fig. 11, for the bin centered at η ν = 7.3, dσ/dη ν is ∼ 2/3 of the solid angle scaling result that applies to larger η ν , the dotted line in the lower right panel of Fig. 11. This factor of 2/3 is the approximate ratio of the M −1 0 dN/dE ν peaks in the upper and lower panels in Fig. 17.  In Fig. 18, we show the uncertainties in the energy distribution of ν τ +ν τ events due to the conventional QCD seven-point scale variation procedure (green) and the 40 different PROSA PDF sets (yellow). These uncertainties affect the charm production cross sections, which, in turn, yield the neutrino flux. We show as well the uncertainty band related to the 32 nCTEQ15 PDF sets, affecting the cross section for neutrino CC DIS inside the detector. The total uncertainties from these three factors combined in quadrature are shown with red. The upper boundary of the total uncertainty band is about 2 times larger than the    Figure 18: Energy distribution of tau neutrino and antineutrino CC interaction events per ton of detector at a distance D d = 480 m from the collider interaction point, in the pseudorapidity ranges of 7.2 < η ν < 8.6 (upper) and η ν > 8.9 (lower). Central predictions refer to (µ R , µ F ) = (1, 1) m T,2 and k T = 0.7 GeV. The presented uncertainties are due to seven variations of the QCD scales (green) around the central value, the PDF uncertainty sets of the PROSA FFNS fit [26] (yellow) and of the nCTEQ15 fit for tungsten [60] (orange) used for production and interaction, respectively, as well as the combination of the QCD and PDF uncertainties. Central predictions with (µ R , µ F ) = (1, 2) m T and k T = 1.2 GeV are also shown (magenta triangles) for comparison. The corresponding ratios to the central predictions of energy distributions of charged-current event numbers with error bands and with (µ R , µ F ) = (1, 2) m T and k T = 1.2 GeV are shown in the ratio plots.
central predictions, whereas the lower boundary is about 80% lower, as one can infer from the lower insets of the panels where ratios of the uncertainty bands to the central predictions for the energy distributions of the CC event numbers are shown. As shown in the figure, the largest contribution to the total uncertainty comes from the QCD scale variation, and the smallest contribution comes from the uncertainty in the interaction cross section arising from the variants of nCTEQ15 PDFs for tungsten, which amounts to about ±2% for E ν ∼ > 100 GeV. The different sets of the PROSA PDF fit, used in charm production, lead to a further uncertainty on the energy distribution of the events amounting to about +(20 − 30)% and −(30 − 40)% of the central prediction, similar in the two rapidity ranges shown in Fig. 18. We also present the predictions evaluated with our parameter set, (µ R , µ F ) = (1, 2) m T and k T = 1.2 GeV, (magenta triangles) for comparison. As discussed above, the predictions obtained with this set are larger than those with the default parameter set. As can be expected from Fig. 14, the alternative parameter set leads to about a 60 -70% increase in the event rate with respect to our default parameter set for E ν ∼ < 500 GeV and 7.2 < η ν < 8.6. The discrepancy increases with energy, so that the alternative parameter set prediction is a factor of ∼ 2 (3) relative to the default set at E ν ∼ 1000 (2000) GeV. On the other hand, for η > 8.9, the discrepancy between the predictions evaluated with the two parameter sets is about 60% for all the presented energies.
For the two neutrino rapidity ranges also considered in previous plots, in Fig. 19 we compare the energy distributions of the ν τ +ν τ CC events obtained using as input for charm production the PROSA PDFs with those computed by using as input the central NLO PDF sets provided by other groups, i.e. ABMP16 [34], CT14 [33] and NNPDF3.1 [32]. The NNPDF3.1 PDFs yield results consistent with PROSA , with differences between central predictions within 10% for E ν ∼ < 1000 GeV. At higher energies, the discrepancy increases with energy, so that the NNPDF predictions reach the upper edge of the PROSA PDF uncertainty band at E ν 1500 GeV. The ABMP16 PDF predictions are almost on the edge of upper boundary of the PROSA PDF uncertainty band, while those with CT14 PDFs are out of the uncertainty ranges. They are larger than the central prediction with the PROSA PDFs by a factor of about 1.5 -2 for E ν 200 -1000 GeV, and the discrepancies become even larger at higher energies, reaching a factor ∼ 3 -3.5 at E ν = 2000 GeV. The difference in the charm mass value used in the simulations with PROSA and CT14 NLO PDFs (see Table 1) plays a decreasing role with energy and is definitely not enough to explain such a large discrepancy. Fig. 8 shows that the gluon PDF for CT14 is more than a factor of ∼ 2 larger than the gluon PDF for PROSA for x = 0.45, with an increasing factor as x gets larger, while at small x, the ratio of the CT14 gluon PDF to PROSA gluon PDF for Q 2 = 10 GeV 2 is closer to ∼ 1.1. Considering Fig. 8 together with the lower panel of Fig. 9, we can infer that it is the large x behavior of the PDFs that is responsible for the discrepancy between the CT14 and PROSA results for neutrino energies above ∼ 1 TeV.
In Table 2, the expected numbers of CC events induced by tau neutrinos, tau antineutrinos, as well as the sum of these two components, are shown, respectively for a 1 m length of tungsten detector. For the given rapidity ranges, the target mass is 29.56 ton for 7.2 < η < 8.6  than those with the first set (default) for 7.2 < η ν < 8.6 (η ν > 8.9). For the results with the default set, we also present the total uncertainty in the number of (ν τ +ν τ ) induced CC events, as well as the uncertainty component due to QCD scale variation, the PDF uncertainties related to neutrino production, using the PROSA PDFs and those related to neutrino interactions, using the nCTEQ15 PDFs. The total CC event numbers range between values that are larger by a factor of ∼ 2 and about 80% smaller compared to the central prediction.
As seen in Fig. 13, the total uncertainty is dominated by the QCD scale variations, which yield 90% higher and 70% lower values than the central CC event numbers, for the upper and lower predictions. The 40 different sets of tbe PROSA PDF fit in neutrino production impact the total events number by +20% and −30% while the 32 variants of nCTEQ15 in neutrino-nucleus interaction make a difference of only ±2%.
In Table 3, we estimate the numbers of ν τ ,ν τ and ν τ +ν τ induced CC interaction events for the two experiments, FASERν (η ν > 8.9) and SND@LHC (7.2 < η ν < 8.6) which will be carried out during the Run 3 stage of the LHC by taking into account the different target masses of the two detectors. The expected total event numbers are in the range of 0.9 − 8.0 for SND@LHC and in the range of 2.3 − 23.7 for FASERν, according to the estimate with our default parameter set. The number of events is dominated by neutrinos with hundreds of GeV rather than TeV energies. The results from different central PDF sets, although not always lying within the PROSA PDF uncertainty band, turn out to always lie within the scale uncertainty range in the evaluation. With the m T dependent parameter set, the central prediction is 7.5 and 19.9 for SND@LHC and FASERν, respectively.

During the High Luminosity-LHC stage
The High-Luminosity (HL)-LHC stage is planned to start from 2027, and to lead to an overall integrated luminosity L = 3000 fb −1 . Possible forward experiments for the HL-LHC era are under study with upgrades of FASERν and SND@LHC, and with more additional experiments. The detector specifics are not yet determined for the experiments at this stage,  Table 3: The numbers of events induced in the FASERν (1.2 tons of tungsten, η ν > 8.9) and SND@LHC (830 kg of tungsten, 7.2 < η ν < 8.6) detectors, by the CC interactions of ν τ andν τ arising from the decay of D ± s mesons produced in pp collisions at √ s = 14 TeV for an integrated luminosity L = 150 fb −1 . The actual shape of the detectors is taken into account in an approximate way, by assuming that they are cylinders or a portion of a cylindrical shell as described in Sec. 4.2.
and it is also possible that the baseline will be changed. Thus, in order to estimate the event numbers, we use a hypothetical detector by assuming that both the radius and the length of detector is 1 m as in Ref. [15]. The pseudorapidity range covered by this setup corresponds to η ∼ > 6.9 for 480 m of the baseline. As in Fig. 17 related to different rapidity ranges, the energy distribution of the number of ν τ ,ν τ and ν τ +ν τ induced CC events per energy per ton of tungsten detector for the pseudorapidity range η ν > 6.9 is shown in Fig. 20 for the default parameter set, (µ R , µ F ) = (1, 1) m T,2 with k T = 0.7 GeV. Even in this case, as visible in the figure, the transition between direct and chain neutrino contributions to the number of events, here occurring at E ν ∼ 300 GeV, is visible in the shape of the distribution. Fig. 21 shows the energy spectrum of the total number of ν τ +ν τ events for the different rapidity range probed in this work for the default input parameter set, (µ R , µ F ) = (1, 1) m T,2 with k T = 0.7 GeV (upper) and the alternative set (µ R , µ F ) = (1, 2) m T and k T = 1.2 GeV (lower). In order to compare the shape of the spectra, we present the results per unit luminosity. As can be inferred from Fig. 9, the neutrino events at high energies are mostly from charm production at large rapidity. This feature combined with the scaling behavior of the neutrino rapidity distributions makes the energy spectrum of the event number for η > 8.9 harder at high energies compared to those s = 14 TeV η ≳ 6.9, L tot = 3000 fb -1  Figure 20: The same as Fig. 17, except for the pseudorapidity range and the integrated luminosity, which in this plot are set to η > 6.9 and L = 3000 fb −1 , respectively.
for the lower rapidity ranges. As discussed above, one can also see that the predictions for the neutrino energy distributions for η ν > 6.9 and 7.2 < η ν < 8.6 are lower than for η ν > 8.9 due to deviation from the scaling behaviour with Ω ν (η ν ) for the neutrino rapidity distribution when η ν < 8.3 (see Fig. 11). The figure also shows that the transition from the dominance of events from direct neutrinos to chain neutrinos somewhat increases with energy as the minimum rapidity increases. Fig. 22 shows the uncertainties in the energy distribution of ν τ +ν τ events due to the QCD scale variations (green) and PROSA PDF eigenvectors (yellow) involved in charm production, the uncertainty of the interaction cross section from the nCTEQ15 PDF sets (orange), as well as the total combined uncertainties (red) using the quadrature formula. The percentage size of uncertainties with respect to the central predictions are similar to Fig. 18. For the difference in the central predictions with the alternative parameter set (µ R , µ F ) = (1, 2) m T with k T = 1.2 GeV and the default parameter set (µ R , µ F ) = (1, 1) m T,2 with k T = 0.7 GeV, the trend is similar to the case of 7.2 < η ν < 8.6. The result evaluated with the first one is 60 -70% larger, for E ν ∼ < 500 GeV, and the difference further increases at higher energies, approaching a factor of ∼ 2 (2.5) at E ν ∼ 1000 (2000) GeV.
As in Fig. 19, Fig. 23 presents the comparison of the energy distributions of the events for ν τ +ν τ for the PDFs by different groups. While the results with the NNPDF3.1 and ABMP16 PDFs can be considered as being consistent with the predictions with the PROSA PDFs, being within the uncertainties of the latter, the central CT14 prediction is obviously out of range of the PROSA PDFs. The size of the difference is similar to the case of the 7.2 < η ν < 8.6 range. Unfortunately, the CT14 PDFs with three active flavours at all scales s = 14 TeV,  do not come with a group of variations out of which an uncertainty band can be computed. Table 4 presents the prediction for the number of CC interaction events due to tau neutrinos, antineutrinos and their sum, measured in a 1 m long tungsten detector, in the pseudorapidity range η ∼ > 6.9, which corresponds to a target mass of 60.63 ton. We present the results with our default set of input parameters, (µ R , µ F ) = (1, 1)m T,2 with k T = 0.7 GeV, together with their scale and PDF uncertainties affecting neutrino production and  Ratio Figure 22: The same as Fig. 18, except for the pseudorapidity range and the integrated luminosity, which in this plot are set to η ν > 6.9 and L = 3000 fb −1 . Ratio Figure 23: The same as Fig. 19 except for the pseudorapidity range and the integrated luminosity, which are set to η > 6.9 and L = 3000 fb −1 , respectively.
the PDF uncertainties affecting neutrino interactions, as well as central results with different PDFs and those from the set of parameters (µ R , µ F ) = (1, 2)m T with k T = 1.2 GeV.  TeV for an integrated luminosity L = 3000 fb −1 and η ν 6.9.
The total (ν τ +ν τ )-induced CC interaction event number with the default parameter set is predicted to be about 4800 with a variation in the range of 1000 -9000 due to the considered uncertainties. The central prediction evaluated with the m T -related parameter set leads to about 8600 events, that is 80% larger than the results from the default set, but still in its uncertainty bands.

Comparisons with previous computations
In Ref. [15], we evaluated the number of events for η ν > 6.87 in a lead detector of length 1 m. A volume of tungsten has ∼ 1.9× more nucleons than the same volume of lead. Depending on the renormalization and factorization scales and k T , the number of ν τ +ν τ CC events for L = 3000 fb −1 , approximately converted to tungsten, range between 3600−7300, with a wider uncertainty band associated with the scale uncertainties. Differences between our evaluations in the work presented here and in our prior work [15] are three-fold: we used different scale, PDF and k T inputs, charm quark fragmentation was implemented in the hadron center-ofmass frame in our earlier work, as compared to implementation in the parton center-of-mass frame in the present work, and the neutrino interaction cross section was evaluated in Ref. [15] using the nCTEQ PDFs for lead. In our comparisons below, we convert our results from lead to tungsten by multiplying by the ratio of the number of targets in tungsten relative to lead and neglect nuclear effects that may arise in the different PDFs for tungsten and lead. The scale dependence considered in Ref. [15] was based on m T rather than m T,2 used here. For the central scales (1, 1)m T and k T , Table 2 of Ref. [15], converted to a tungsten detector with the same volume, yields a central prediction of 4541 CC events from ν τ +ν τ from D ± s decay. This agrees well with the central value shown here in Table 4, 4775 for (1, 1)m T,2 and k T = 0.7 GeV. However, as noted in Sec. 2.2, the error bands from the scale dependence are more asymmetric for scales that depend on m T rather than m T,2 . Ref. [15] shows that the upper edge of the envelope for the scale dependence around the number of events as a function of energy for (1, 1)m T and k T = 0.7 GeV is of order a factor of 3 larger than the central value, whereas here with m T,2 scale dependence, the upper edge of the uncertainty envelope for the scale dependence is close to a factor of 2 larger than the central value. The lower edge of the scale uncertainty envelope in our previous work deviates less from the central value than what we find here for m T,w dependent scales. In Ref. [15], we found that the event number on the lower edge of the scale uncertainty envelope is ∼ 75% of the central value for the event number.
We also considered as a central scale (1, 1.5)m T with k T = 0.7 GeV, and, for better agreement with the LHCb data [20], k T = 2.2 GeV. A seven-point scale variation around (1, 1.5)m T with k T = 0.7 GeV yielded an event number [15] which, when converted to tungsten, amount to 7284 +14143 −3698 accounting for the scale variation. Again, this shows an asymmetric scale uncertainty band. We found that for (1, 1.5)m T with k T = 2.2 GeV, the predicted number of events, converted to tungsten, is 5251 CC events, which can be compared with our predictions shown here in Table 4, 8616 CC events for (1, 2)m T with k T = 1.2 GeV. Overall, however, the scale uncertainty band in Ref. [15] leads to a higher minimum number of CC events than what we have presented here in Table 4.
An estimate of the number of tau neutrino and antineutrino induced interaction events computed with various standalone event generators initially developed for simulations of cosmic ray-induced extended air showers appears in Ref. [59]. Event generators, used in Ref. [59] with their default options, include SIBYLL2.3c [74][75][76] and DPMJET (version from 2017) [77,78]. As an alternative to SIBYLL and DPMJET, the PYTHIA8 [79] code, very popular for LHC phenomenology, was also adopted in Ref. [59], using as the seed a tree-level description of cc quark pair production, on top of which parton shower, hadronization and other soft physics effects plus hadron decay are included.
The theory frameworks implemented in the SIBYLL and DPMJET generators are oriented especially to the description of soft physics effects. Their use in the description of charm production, a process characterized by a hard scale even at small p T , presents a number of challenges, the first one being related to the lack of radiative corrections in the description of the hard-scattering events. Radiative corrections affect not only the normalization but also the shape of differential distributions. It is worth mentioning that, notwithstanding this lack, extensive and systematic comparisons of D-meson energy distributions obtained using a QCD based approach including NLO QCD radiative corrections, parton shower and nonperturbative QCD effects, with those from SIBYLL, were carried out by one of us in close collaboration with the authors of this code, finding reasonable agreement on a very wide set of center-of-mass energies, at least as far as the p−Air collisions that matter most for cosmic ray interactions in the atmosphere are concerned. Unfortunately these comparisons, although of practical interest for those experiments still relying on these Monte Carlo generators, less CPU-intensive than more complex calculations, do not allow for deep QCD insights, considering that the uncertainties associated with leading-order evaluations of charm production are so large that leading-order calculations can provide, at most, an order of magnitude estimate for this process. By no means can they lead to an accurate evaluation and, therefore, should not be misunderstood as such.
The Monte Carlo evaluations in Ref. [59] that account for specific detector geometries (in contrast to our use of approximate rapidity ranges) yield estimates of the number of CC events in the range 10.1-22.4 for SND@LHC and 21.2-131 for FASERν, where the range of values arises from the use of different generators, and is far from being an estimate of the uncertainties within each event generator. The lower numbers are compatible with our central predictions with the µ 0 = (1, 2)m T and k T = 1.2 GeV parameter set, which yield a good match to cross sections for D ± s production at LHCb. On the other hand, our central predictions with our default parameter set yield a smaller number of tau neutrinos, which enables a very cautious perspective concerning the possibility to detect a significant number of tau neutrinos in the Run-3 far-forward experiments. In any case, considering the large uncertainty due to scale variations affecting our predictions, and the possible enhancing effect of additional radiative corrections (in particular the NNLO ones are positive for inclusive cross sections), one can still hope for a more optimistic scenario. Furthermore, considering the even larger uncertainty affecting calculations as approximate as those in Ref. [59], which, as explained above, can provide at most just an order-of-magnitude estimate, one can still claim for consistency between their predictions, the predictions of other authors using tools with similar accuracy (see, e.g., the evaluations conducted in Ref. [17] and in Ref. [9]) and ours.

Conclusions
We have performed a QCD evaluation of the (ν τ +ν τ ) rapidity and energy distributions, as well as distributions of CC events induced by them, of interest to experiments in the forward region of the LHC. NLO QCD radiative corrections are included in our estimates of charm production and neutrino-induced DIS. Our focus is on PDF uncertainties. As for charm production, we use as a basis the 40 sets of the most recent PROSA NLO PDF fit, and, as an alternative, the central sets of other 3-flavour NLO PDFs. In particular, we consider those provided by the ABMP16, CT14, and NNPDF3.1 collaborations, together with their associated charm mass values, converted to the on-shell renormalization scheme for those PDFs accompanied by charm mass values in the MS mass renormalization scheme, and their associated α s (M Z ) values. All these sets have been widely used for LHC predictions at small and mid-rapidities and their effects have also already been studied in a more inclusive rapidity range in the context of prompt neutrino production in the atmosphere (see e.g. [26,28,40,80]). We investigate their effects in the computation of forward ν τ andν τ spectra at the LHC for the first time in this work.
Within the PROSA PDFs, the largest uncertainties turned out to come from the scale dependence in the theoretical predictions for charm production used for the fit, rather than from the fit procedure itself, the experimental data, and the parameterization ansatz specifying the form of the PDFs at the initial scale or the variations of this scale. These PDF features are reflected in our predictions, where the largest component of the PDF uncertainty turned out to be the so-called "model" uncertainty from the PROSA fit. Our computation of charm hadroproduction has the same accuracy as in the PROSA fit, and, as default, we adopt the same renormalization and factorization scales (m T,2 ). While the scale uncertainty is reduced when setting the scale equal to m T,2 , rather than to the m T functional form adopted more often in the literature, the total number of events ranges between ∼ 0.2 − 2 times the central prediction, depending on the scale. The percentage size of scale uncertainty with respect to central predictions is approximately uniform with rapidity. This large scale dependence is a signal that higher-order corrections to charm pair production are needed to improve our current estimates of D-meson production at the LHC. This applies to all rapidity ranges discussed in this paper.
For fixed input renormalization and factorization scales and k T , the PROSA PDF uncertainty band does not include all the other PDF central predictions for the tau neutrino energy distribution and the corresponding number of events, although it is expected that the uncertainty bands for the other PDFs overlap those of the 40 PROSA sets. The central CT14 NLO prediction is the largest. It yields a number of events that is a factor of ∼ 1.6 (∼ 1.9) times higher than the central PROSA prediction for 7.2 < η ν < 8.6 (η ν > 8.9) as compared to the PROSA PDF uncertainty of ∼ ±20 − 30%. However, the scale uncertainty for the PROSA evaluation still covers the range of event number predictions from the other central PDF sets.
If we restrict our attention to E ν ∼ > 1 TeV, the CT14 PDF predicted number of events are more than a factor of 2 larger than the PROSA ones. Predictions corresponding to high neutrino energies are sensitive to gluon PDFs in the x ∼ > 0.3 region, where the constraints from data are not strong. Besides the invaluable role that will be played by future collider experiments, like the foreseen electron-ion collider (EIC) [81], even future measurements of tau neutrino and antineutrino interactions in this energy regime may help pin down the behavior of large-x PDFs. At such high energies at forward rapidities, ν e +ν e fluxes from charm meson decays dominate with respect to the corresponding contributions from kaon decays [59], thereby enhancing the statistics for CC events from high energy neutrinos from charm decays. Kaon decay contributions to ν µ +ν µ would have to be well modeled in order to exploit high energy fluxes of all three flavours of neutrinos to better understand large-x PDFs at a Forward Physics Facility.
We have considered different renormalization and factorization scale dependencies and different k T values. The choice of the functional form of the scale dependence (m T or m T,2 ) causes an offset to the energy distribution of ν τ +ν τ that changes slowly with energy above a few hundred GeV. The offset is largely independent of rapidity range. On the other hand, the value of k T impacts different rapidity ranges in a different way for neutrino energies above ∼ 500 GeV. In particular, the ratio of the energy distribution with k T = 1.2 GeV to the energy distribution with k T = 0.7 GeV is larger, and grows with energy, for η ν > 6.9 as compared to η ν > 8.9.
The energy distributions of ν τ +ν τ , shown in Figs. 12 and 13, are reported in Tables  5-10 in Appendix B and in ancillary files on the arXiv, distributed together with the preprint version of this manuscript.
Our ν τ +ν τ energy distributions differ somewhat from those estimated with Monte Carlo generators, as in, e.g., Ref. [59]. Detailed comparisons are difficult as the detector geometries are handled in different ways. SND@LHC reports that they expect 25 ν τ +ν τ CC events for their 850 kg tungsten detector for 150 fb −1 [9], much more than our prediction of 1 − 8 events for 830 kg of tungsten for 7.2 < η ν < 8.6. The estimate of 1-8 events accounts for our uncertainty band in Table 3, whereas the value of 25 events is not accompanied by an uncertainty. However, if it would, the corresponding band would be bigger than our one. Compared to the SIBYLL 2.3c result for FASERν in Ref. [59], our NLO neutrino energy distribution for η ν > 8.9 is harder. For (µ R , µ F ) = (1, 2)m T and k T = 1.2 GeV, our results for the neutrino energy distribution are larger by a factor as large as ∼ 2 − 6 for E ν = 10 − 300 GeV, increasing to a factor of ∼ 10 at 2 TeV neutrino energy, nearly matching the DPMJET3.2017 histogram in Ref. [59] at TeV neutrino energies. Our neutrino energy distribution evaluated using default parameter set with PROSA PDFs most resembles the PYTHIA8.2 predictions for E ν ∼ > a few hundred GeV. However, our distributions at low energies are higher than the PYTHIA8.2 ones. That the central predictions using perturbative QCD at NLO are larger than Monte Carlo evaluations is not surprising, since the Monte Carlos simulations use as a seed leading order evaluations of charm quark pair production. While performing all these comparisons, we should recall that leading order predictions for heavy-quark hadroproduction can, at most, provide an order-of-magnitude estimate, being accompanied by a huge scale uncertainty band, much larger than the corresponding NLO scale uncertainty band. A NNLO/NLO comparison will be interesting, and will become possible as soon as NNLO predictions for differential distributions for charm pair and Dmeson production become available. For the total number of events, the uncertainties due to the scale dependence and the PROSA PDFs in the evaluation of charm production are similar in all the rapidity ranges considered here, namely rapidities larger than 6.9. With our defaults scales proportional to m T,2 and k T value of 0.7 GeV, the scale variations lead to approximately +90% and -70% uncertainty around the central value. The PROSA PDF uncertainties are between +20% and -30%. The NNPDF3.1, CT14 and ABMP16 predictions for the number of events with the default scale and k T value lie within these large error bands. The CT14 set yields the largest prediction for the number of events given the default inputs, depending on the rapidity interval, which is a factor of 1.5-1.8 larger than the prediction using the PROSA PDF central set with (µ R , µ F ) = (1, 1)m T and k T = 0.7 GeV, as shown in Tables 3 and 4. For an integrated luminosity of 150 fb −1 and nominal rapidity and mass values for SND@LHC (7.2 < η ν < 8.6, 830 kg) and FASERν (η ν > 8.9, 1.2 ton), the predicted total number of ν τ +ν τ -induced charged-current events ranges in the intervals of 0.9-8.0 and 2.3-23.7 events, respectively.
The high luminosity era and larger neutrino detectors will bring a significant increase in the number of ν τ +ν τ charged current events. Our predictions with central PDF set range in the interval ∼ 4, 500−8, 600 events for a tungsten detector with a 1 m radius and a 1 m length. The high statistics of neutrino events from charm in experiments in a Forward Physics Facility will enable better constraints on the small-x and large-x PDFs, as LHCb measurements of charm hadrons for y = 2 − 4.5 has already done for gluon distributions with x as small as 10 −6 . [26,27,29]. Meanwhile, it will be important to develop differential calculations including NNLO and higher-order corrections (see e.g. the recent work of Ref. [82]), to refine the theoretical predictions by reducing the large scale uncertainties that are inherent to NLO, regardless of the PDF set.
Note added: Additional data tables with the double-differential neutrino energy and pseudorapidity distributions of ν τ +ν τ from D ± s production and decay, and of ν e +ν e from D ± production and decay, appear on the arXiv as ancillary files for Ref. [83].
The total PROSA PDF uncertainty comes from For the nCTEQ15 PDFs, the prescription for the PDF uncertainty comes from n = 32 sets. The PDF uncertainty bars, symmetric around the central PDF set, come from evaluating the differential cross sections [60] Σ ± = ± 1 2 k=1,16 B Tables for the energy distributions and charged-current cross sections of tau neutrinos plus antineutrinos Tables 5-10 show the energy distributions and errors that correspond to Fig. 12 and 13. Table  11 shows the tau neutrino and antineutrino charged-current cross section per nucleon with PDF uncertainties for the nCTEQ15 PDF sets for a tungsten target and Table 12 shows the tau neutrino charged-current cross section per nucleon for the VFNS PROSA PDF set, as in Fig. 15. The data tables are also available at https://arxiv.org/src/2112.11605v1/anc.   Table 5: Sum of ν τ andν τ energy distributions dσ/dE ν for √ s = 14 TeV and η ν > 6.9. The predictions are shown for the PROSA NLO PDF set with scales (µ R , µ F ) = (1, 1)m T,2 with k T = 0.7 GeV. This table with 20 GeV energy bins for the whole energy range is available as an ancillary file on the arXiv link to this paper.   Table 6: Sum of ν τ andν τ energy distributions dσ/dE ν for √ s = 14 TeV, 7.2 < η ν < 8.6.
The predictions are shown for the PROSA PDF set with scales (µ R , µ F ) = (1, 1)m T,2 with k T = 0.7 GeV. This table with 20 GeV energy bins for the whole energy range is available as an ancillary file on the arXiv link to this paper.  The predictions are shown for the PROSA PDF set with scales (µ R , µ F ) = (1, 1)m T,2 with k T = 0.7 GeV. This table with 20 GeV energy bins for the whole energy range is available as an ancillary file on the arXiv link to this paper. Table 8: Sum of ν τ andν τ energy distributions dσ/dE ν from D ± s for √ s = 14 TeV, η ν > 6.9.
The predictions are shown for the PROSA PDF set with scales (µ R , µ F ) = (1, 2)m T with k T = 1.2 GeV, and the CT14, ABPM16 and NNPDF3.1 with (µ R , µ F ) = (1, 1)m T,2 and k T = 0.7 GeV. This table with 20 GeV energy bins for the whole energy range is available as an ancillary file on the arXiv link to this paper. Table 9: Sum of ν τ andν τ energy distributions dσ/dE ν from D ± s for √ s = 14 TeV, 7.2 < η ν < 8.6. The predictions are shown for the PROSA NLO PDF set with scales (µ R , µ F ) = (1, 2)m T with k T = 1.2 GeV, and the CT14, ABPM16 and NNPDF3.1 NLO PDFs with (µ R , µ F ) = (1, 1)m T,2 and k T = 0.7 GeV. This table with 20 GeV energy bins for the whole energy range is available as an ancillary file on the arXiv link to this paper.  Table 11: The ν τ andν τ charged current cross section per nucleon for interactions with tungsten, evaluated using the nCTEQ PDFs. The PDF error is also shown, where (σ CC ± ∆)/A represents the PDF error band for a tungsten target. Tables for incident tau neutrino and antineutrino energies between 5 GeV and 5 TeV are available as ancillary files on the arXiv link to this paper.  Table 12: The ν τ andν τ charged current cross section per nucleon for interactions with tungsten, evaluated using the PROSA VFNS PDFs. Isospin symmetry is used for the neutron PDFs. Tables for incident tau neutrino and antineutrino energies between 5 GeV and 5 TeV are available as ancillary files on the arXiv link to this paper.