Far-forward neutrinos at the Large Hadron Collider

We present a new calculation of the energy distribution of high-energy neutrinos from the decay of charm and bottom hadrons produced at the Large Hadron Collider (LHC). In the kinematical region of very forward rapidities, heavy-flavor production and decay is a source of tau neutrinos that leads to thousands of charged-current tau neutrino events in a 1 m long, 1 m radius lead neutrino detector at a distance of 480 m from the interaction region. In our computation, next-to-leading order QCD radiative corrections are accounted for in the production cross-sections. Non-perturbative intrinsic-kT effects are approximated by a simple phenomenological model introducing a Gaussian kT -smearing of the parton distribution functions, which might also mimic perturbative effects due to multiple initial-state soft-gluon emissions. The transition from partonic to hadronic states is described by phenomenological fragmentation functions. To study the effect of various input parameters, theoretical predictions for Ds±\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {D}_s^{\pm } $$\end{document} production are compared with LHCb data on double-differential cross-sections in transverse momentum and rapidity. The uncertain- ties related to the choice of the input parameter values, ultimately affecting the predictions of the tau neutrino event distributions, are discussed. We consider a 3+1 neutrino mixing scenario to illustrate the potential for a neutrino experiment to constrain the 3+1 parameter space using tau neutrinos and antineutrinos. We find large theoretical uncertainties in the predictions of the neutrino fluxes in the far-forward region. Untangling the effects of tau neutrino oscillations into sterile neutrinos and distinguishing a 3+1 scenario from the standard scenario with three active neutrino flavours, will be challenging due to the large theoretical uncertainties from QCD.

JHEP06(2020)032 produce large fluxes of tau neutrinos in the forward direction [12][13][14][15]. Electron and muon neutrinos and antineutrinos from D and B meson decays will also be produced. Charm and bottom quark production in the Standard Model mostly occurs in quark-antiquark pairs, leading to an approximately equal number of hadrons and antihadrons, 1 so there will be approximately equal fluxes of neutrinos and antineutrinos from heavy flavor. In our discussion below, we will refer to both particles and antiparticles as neutrinos. As described in 1984 by De Rujula and Ruckl [12], using Feynman scaling arguments, a quark-gluon string model and empirical extrapolations based on collider data available at the time, a few thousand tau neutrino events in future pp and pp colliders could be detected with a 2.4 ton detector placed 100 m distant from the interaction point along the tangent to the accelerator arc. An estimation of the neutrino flux using Pythia [17] tuned to Tevatron data and rescaled to a √ s = 14 TeV center-of-mass energy gives qualitatively consistent event rates [5].
In the last few years, there is renewed interest in far-forward neutrino production and detection [6-11, 18, 19]. The ForwArd Search ExpeRiment at the LHC (Faser), primarily dedicated to searches for light, extremely weakly interacting particles [6][7][8], with phase 1 approved and under construction, will be sensitive to tau neutrinos if the detector mass is sufficiently large. The location of the detector is projected to be 480 m from the ATLAS interaction point along the colliding beam axis. The Faser collaboration uses a pseudorapidity cut η > 6.87 on particle momenta (here, neutrinos) in its second phase to determine if they enter a detector of radius 1.0 m. Other evaluations for this baseline use η > 6.7 [9]. For a half-cylinder, 2-meter long lead detector covering pseudorapidities η > 6.7, Beni et al. [10] used Pythia 8 [17,20] to find ∼ 8, 700 tau neutrino events for a 3, 000 fb −1 integrated luminosity at the LHC. Other configurations for detectors, for example, in η ranges of 8-9.5 and 7.4-8.2 for XSEN [11] and in the range 7.2-8.7 for SND@LHC [21] are also under consideration. The prototype Faser-ν has η > 8.6 [18,19]. We adopt a minimum pseudorapidity η > 6.87 in this work for definiteness. Our results for η > 6.7 are qualitatively similar.
A source of high-energy tau neutrinos opens the possibility of new tests of the Standard Model that can not be done at DUNE. Measurements of LHC forward tau neutrino interactions can be used for direct tests of lepton flavor universality in charged current interactions with much higher statistics than achieved by Donut [22,23] and Opera [24,25]. Both Super-Kamiokande and IceCube have reported signs of tau neutrino appearance in their datasets, but the statistics is still very limited [26,27]. Measurements of the interaction cross-sections of muon neutrinos from heavy-flavor decays at the LHC with the nucleons/nuclei of the target will help to close the gap between direct neutrino cross-section 1 Differences in the forward region between the total number of hadrons and antihadrons from quarkantiquark pair production arise from the recombination of forward final-state heavy quarks with partons from the initial-state protons non participating to the hard scattering, considering the fact that the parton distribution functions of the up and down quarks in the nucleon differ from those of the antiquarks. Baryon-antibaryon asymmetries, in principle, would show a larger effect than meson-antimeson asymmetries. We show below that the Λc gives a small contribution to neutrino production. The D ± s has no valence component, and experimental studies by LHCb show a D + s − D − s asymmetry below 1% [16].

JHEP06(2020)032
measurements for E ν < 370 GeV [1] and the IceCube Collaboration's determination of the averaged cross-section for neutrino plus antineutrino deep-inelastic scattering with nucleons for E ν = 6.3-980 TeV [28] (see also, e.g., ref. [29]). LHC forward neutrinos will provide the first opportunity for direct neutrino and antineutrino cross-section measurements for neutrino energies up to E ν 2 TeV. In this paper, we perform a new evaluation of the D and B meson contributions to the pp differential cross section as a function of tau neutrino and muon neutrino energy for η > 6.87. Differently from previous evaluations which are limited to leading order/leading logarithmic accuracy, our evaluation accounts for the effects of next-to-leading order (NLO) QCD radiative corrections to the heavy-quark hadroproduction [30][31][32] and neutrino deepinelastic-scattering cross-sections. The effects due to the intrinsic transverse momentum of initial state partons confined in the protons are accounted for by a simple model adding k T -smearing effects to the standard collinear parton distribution functions (PDF) as an input for the calculation. The same model might also mimic the effects of the resummation of logarithms related to initial-state soft gluon emissions in an approximate and purely phenomenological way. Furthermore, we include a description of the fragmentation of partons into heavy mesons, relying on widely used phenomenological fragmentation functions.
Our main focus is on tau neutrinos which come predominantly from D ± s decays. To study the effects of different choices of the values of various parameters entering our computation, we compare our theoretical predictions with the LHCb data on D ± s production [33] in the rapidity range of 2 < y < 4.5. These same parameters give predictions in similar agreement with experimental data for other charm hadron distributions at LHCb.
As discussed below, the effect of the transverse momenta of initial state partons can significantly impact the predictions for forward tau neutrino event rates. In general, the LHCb and other charm data give hints of the need for higher-order effects in the description of single-inclusive open D-hadron production, beyond NLO and the limited logarithmic accuracy of the parton shower implementations and of the analytical resummations of various kinds of logarithms presently available. Power suppressed non-perturbative terms might also play a relevant role, considering the smallness of the charm quark mass, still larger but not too large with respect to Λ QCD . At present, it is not clear if the discrepancies between theoretical predictions and experimental data at low transverse momenta can be completely cured within the collinear factorization framework, or if it is necessary to go beyond this scenario. Considering the unavailability of rigorous perturbative and non-perturbative QCD theoretical calculations accurate enough to reproduce the shape of the experimental transverse-momentum distributions, we incorporate initial-state partonic k T -smearing effects in a purely phenomenological way in the QCD calculation of D ± s production for the phase space covered by LHCb and assess their impact on the predicted number of neutrino charged-current interaction events.
The elements of the PMNS mixing matrix are least constrained in the tau sector, compared to the other flavor sectors. We demonstrate how sterile neutrino mass and mixing parameters can begin to be constrained by measuring oscillations of LHC tau neutrinos, and we discuss the challenges to pushing these constraints to sterile neutrino masses of order ∼ 20 eV.

JHEP06(2020)032
Heavy-flavor decays are not the only sources of forward neutrinos, but they dominate the forward tau-neutrino flux. The contributions from pp → W, Z production followed by W, Z leptonic decays are negligible for η 6.5 [10]. On the other hand, for muon and electron neutrinos, charged pion (π ± ) and kaon (K ± ) decays to ν µ and K L → ν e decays (K 0 e3 ) are most important. For our evaluation of muon and electron neutrino oscillations, we use parametrizations of light meson distributions based on Pythia distributions [34].
We begin in section 2 with an overview of the far-forward geometry used for our discussion here. In section 3, we present our D and B hadron production results. Energy distributions of charged-current interaction events generated in a forward detector by neutrinos from heavy-flavor production and decay are shown in section 4 for both tau neutrinos and muon neutrinos. In our estimation, we also account both for muon neutrinos from the decays of charged pions and kaons and for electron neutrinos from kaon decays. Section 5 shows an application to the study of tau neutrino oscillations, considering a 3+1 oscillation framework with three active neutrinos and one sterile neutrino. We conclude in section 6. Appendix A collects formulas for the decay distributions to neutrinos.

Overview of forward neutrino detection geometry
A forward detector along a line tangent to the LHC beam line necessitates calculations in high-pseudorapidity regimes. A detector with radius of 1.2 m placed at 480 m from the LHC interaction point is used for the evaluations in ref. [10]. This corresponds to a neutrino pseudorapidity of η > 6.7 for detection. The Faser2 proposal [8] has R = 1.0 m, corresponding to η > 6.87, which we use for the results shown below. Other smaller detectors like XSEN [11] and SND@LHC [21] are discussed in the recent literature. A prototype Faser-ν at even higher η, 2 with a 25 cm × 25 cm cross sectional area and length of 1.35 m of tungsten interleaved with emulsion detectors, will have a few tens of tau neutrino events with an integrated luminosity of 150 fb −1 delivered during Run 3 of the LHC [18]. A detector radius of 1.0 meter, for a comparable target mass, increases the number of events by a factor of ∼ 50 when scaling by cross-sectional area alone.
The LHC interaction region is a very compact source of tau neutrinos. Most of the tau neutrinos come from D s → ν τ τ decay. The τ → ν τ X decay is also a prompt process. The characteristic size of the region where tau neutrinos are produced by D s decays is of order γcτ Ds E Ds /m Ds · 150 µm, so of order of 1.5-15 cm for E Ds = 200 GeV-2 TeV. The tau decay length cτ = 87.11 µm, multiplied by the γ-factor for the same energy range, gives a size of 0.98-9.8 cm. Similarly, γcτ B + = E B + /m B + · 496 µm produces a size of 1.9-19 cm for the same energy range. Thus, for tau neutrinos produced along the beam pipe, the longitudinal distance is a few to 20 cm. The transverse size is 17 µm [35] from the proton bunch size. The compact source means the assumed detector radius of 1 m and distance of 480 m from the interaction point translates to a maximum angle relative to the beam axis for the tau neutrino three-momentum of θ max = 2.1 mrad (η min = 6.87). This same constraint applies to the momenta of muon and electron neutrinos from heavy-flavor decays.

JHEP06(2020)032
While the focus of this paper is on heavy-flavor production of neutrinos in the forward region, consideration of oscillations in a 3+1 mixing framework also requires an estimate of the number of electron and muon neutrinos from both heavy-flavor decays and light-meson decays. In section 4.2, we make an estimation of the production of electron neutrinos and muon neutrinos from light-meson decays. The light-meson decay lengths are long compared to heavy-meson decay lengths, so the detector and magnets near the interaction point play a role. In ref. [18], an evaluation of the number of ν µ +ν µ events in a detector of 25 cm × 25 cm cross sectional area finds that most of the events below 1 TeV come from charged pion and kaon decays that occur within 55 m of the interaction point and stay within the opening of the front quadrupole absorber with inner radius of 17 mm. This corresponds to light-meson momenta lying within 1 mrad from the beam axis. A more detailed discussion of this point appears in section 4.2.
Heavy-meson production at small angles with respect to the beam axis receives dominant contributions from the transverse momentum region below few GeV. For example, for E Ds = 1 TeV, approximating the neutrino direction by the D s direction means that the p T of the D s meson must be smaller than 2.1 GeV, and even smaller for lower energies. Non-perturbative effects related to the intrinsic k T of the partons confined in the initial state nucleons are important at such low transverse momenta. Additionally, perturbative effects related to the appearance of large high-energy logarithms together with those due to initital-state multiple soft-gluon emissions, are potentially relevant in the p T range [0, 15] GeV of the LHCb data considered in this paper, covering low to intermediate p T values.
In our evaluation of differential cross sections for open charm and bottom production, we include Gaussian k T -smearing to better match LHCb data, using a purely phenomenological approach. In this approach the non-perturbative effects are reabsorbed in the k T smearing model, which may also mimic part of the all-order perturbative effects in a rough way. We found that a better description of the LHCb data, given the renormalization and factorization scale choices discussed below, is provided in our approach by a k T value higher than the naive estimate due to Fermi motion and even larger than the upper end of the range of typical non-perturbative k T values of ∼ 1-2 GeV reported in the literature [36][37][38][39][40]. This gives hints, on the one hand, that non-perturbative physics aspects not yet well understood might be particularly relevant for the process we are studying. On the other hand, it hints that the contribution from the Sudakov resummation of double logarithmic perturbative terms related to the emission of an arbitrarily large number of soft gluons, missing in fixed-order calculations, may be large. Since our focus is on ν τ +ν τ production, we use the LHCb data for D s production at √ s = 13 TeV for rapidities in the range of 2.0-4.5 [33].
3 Forward heavy-flavor production and decay at the LHC We evaluate the single-particle inclusive heavy-flavor energy and angular distributions. Our approach relies on perturbation theory (pQCD) in the collinear factorization framework. In particular, we include NLO QCD corrections to the heavy-quark production cross sections [30][31][32].
-5 -JHEP06(2020)032 As noted above, at very high rapidities, the effect of relatively small transverse momenta of the initial state partons can affect the acceptance of neutrino events. In shower Monte Carlo event generators like Pythia, the transverse momentum distribution of the produced heavy quarks is affected, on the one hand, by the effects of multiple soft and collinear gluon emissions, accounted for with a limited logarithmic accuracy by the introduction of Sudakov form factors resumming the relative main logarithmic contributions, and, on the other hand, by the inclusion of a small intrinsic transverse momentum, related to the confinement of partons in finite-size nucleons and the uncertainty principle. An alternative to the collinear factorization approach is to use unintegrated parton distribution functions that have a transverse momentum k T dependence in addition to the usual longitudinal momentum fraction x dependence [41,42]. In both collinear and k T factorization, further higher-order effects can also modify the transverse momentum distribution of heavy-flavor hadroproduction in the low p T region. In particular, the effect of the resummation of high-energy logarithms, not yet implemented into a publicly available code, and the joint resummation of these and other logarithms could also play a relevant role, which deserves future dedicated investigations, but is beyond the scope of this paper.
In the present absence of calculations capable of fully reproducing the experimental shape of the transverse momentum distributions of charm mesons at small transverse momenta at the LHC, our approach here is phenomenological and uses Gaussian smearing of the outgoing charm quark. With the NLO QCD calculation of g NLO (q T , y) = d 2 σ(NLO)/dq T dy for the quark, we evaluate the Gaussian smeared one-particle inclusive charm-quark distribution g(p T , y) = d 2 σ/dp T dy according to where f ( k T , k 2 T ) is a two-dimensional Gaussian function normalized to unity, The smearing generates a shift of the outgoing heavy-quark momentum vector q T by k T after the hard scattering, before fragmentation. The single Gaussian f ( k T , k 2 T ) is equivalent to starting with two Gaussian functions, one for each incoming parton (k iT ), making a change of variables to k T = ( k 1T + k 2T )/2 for one of the integrals, and integrating over the other parton's k T , This definition of the effective k T variable distributes the transverse momentum smearing equally to the one-particle transverse momentum and to the recoil hadrons [36]. The quantity k 2 T is related to the average magnitude of k T , k T , by As we discuss below, we set the k 2 T value based on comparisons with double-differential distributions in rapidity and transverse momentum for forward D s production measured by LHCb at √ s = 13 TeV [33].

JHEP06(2020)032
Charm quark fragmentation to mesons is accounted for by using fragmentation functions of the Peterson form [43]. For c → D 0 , D + and D + s , we take fragmentation fractions 0.6086, 0.2404 and 0.0802, respectively [44]. The parameter in the Peterson fragmentation function that describes the hardness of the meson spectrum relative to the heavy quark is taken to be = 0.028, 0.039 and 0.008 for D 0 , D + and D + s , respectively, and = 0.003 for B production [45]. The fragmentation fraction for B's is b → B + = b → B 0 = 0.362 from ref. [46].
We use the NLO nCTEQ15 parton distribution function (PDF) grids [47], available for free proton and nuclear targets, for our evaluations here. We use their best fit set for free proton targets in our evaluation of charm production in pp collisions. As noted below, we use the NLO nCTEQ15 nuclear PDF set for lead targets in our neutrino cross section calculation. We take a charm pole mass value m c = 1.3 GeV consistent with the choice of PDFs. A b-quark pole mass of 4.5 GeV is used here, also consistent with the PDFs. We use renormalization and factorization scales (µ R , µ F ) which are factors (N R , N F ) of the transverse mass m T = m 2 Q + p 2 T , where p T is the magnitude of the transverse momentum of the heavy quark Q = c, b. Conventionally, renormalization and factorization scales are chosen in a range with factors N R , N F = 0.5-2 in (µ R , µ F ) = (N R , N F )m T , with N R = N F = 1 set as conventional scale factors [48]. The values N R = 1.0 and N F = 1.5 lie in the standard scale factor uncertainty range and are used as our default parameters as discussed below, but we also show predictions obtained with the standard conventional choice N R = N F = 1.0 multiplying m T for (µ R , µ F ) used in most, (see, e.g., ref. [48]), though not all (see, e.g., refs. [49,50]), of the literature on heavy-flavor production. In particular, in the following we will discuss predictions obtained with the following configurations: • N R = 1, N F = 1 (conventional central scale choice), with k T = 0.7 GeV, • N R = 1, N F = 1.5, (alternative central scale choice, used as default in this paper, as better motivated in the following), with k T = 0.7 GeV, k T = 0 GeV and k T = 2.2 GeV. The scale input, the PDFs, the fragmentation functions and the non-perturbative transverse momenta all influence the predicted heavy-flavor energy and rapidity distributions. Since our focus is on tau neutrino production, LHCb data on forward D s production [33] are used to set k T for selected (µ R , µ F ) combinations, after fixing the PDF and fragmentation function details. The D s data cover the range 0 < p T < 14 GeV 3 and 2.0 < y < 4.5. While varying the input parameters, a χ 2 is computed, which combines the double-differential predictions, the binned data for 71 data points and the experimental uncertainties. As additional constraint, the integrated cross section evaluated with our Monte Carlo integration program is also required to be within the experimental error bars of the measured cross section for the same kinematic region. Fitting the LHCb data for D s production with p T < 14 GeV, by minimizing the χ 2 with µ R = m T , µ F = 1.5 m T and varying k T , leads to k T = 2.2 ± 0.7 GeV ( k 2 T = 6.2 GeV 2 ). This represents a JHEP06(2020)032 reasonably good fit, with χ 2 /DOF= 2.8 and the cross section within 10% of the experimentally measured cross section in this kinematic region. With this parameter choice, the theoretical total cross section for σ cc is within 2% of the central value of the LHCb estimate reported in ref. [33].
The predictions with µ R = m T , µ F = 1.5 m T and k T = 2.2 GeV, together with those for a k T = 0 evaluation using the sames scales and the LHCb data [33], are shown in the upper left panel of figure 1 with the solid and dashed histograms, respectively. From top to bottom, the panel shows data and theoretical evaluations in five different rapidity bins of ∆y = 0.5 width for 2.0 < y < 4.

JHEP06(2020)032
The lower left panel of figure 1 shows the central results and uncertainty band for the same scale choices, now with k T = 0.7 GeV, the transverse momentum smearing that approximates the Powheg + Pythia results (see figure 2). Finally, the lower right panel shows results with k T = 0.7 GeV and the conventional central choice of QCD scales, namely, N R,F = 1.0.
The uncertainty bands in figure 1 are consistent across the rapidity bins. We expect that the large scale uncertainties here dominate the uncertainties at even higher rapidity. PDF uncertainties are typically smaller than the scale uncertainties for forward charm production at these energies [50]. Results for the p T and y distributions, in similar agreement with the LHCb data, are obtained for D ± and D 0 ,D 0 production, used below as sources of ν µ +ν µ and ν e +ν e .
In general varying (N F , N R ) changes both the rate and the shape of the distributions. The shape is particularly sensitive to N F , whereas the normalization is particularly sensitive to N R . At fixed N F , the larger N R gives a smaller cross section. The variation of the shape with N F depends on the √ s. At LHC energies, the p T distribution is steeper for larger N F than for smaller N F .
Here, B-meson production follows from the same heavy-quark dynamics, with cc replaced by bb, so it is generally expected that the renormalization and factorization scaling factors are the same for both heavy-quark flavors. To the extent that k T describes the intrinsic parton transverse momentum in the initial state protons, this parameter also is process independent for pp → QQX for Q = c, b. For the same central scale and k T choices, the p T distribution of B ± mesons at √ s = 13 TeV for y = 2.0-4.5 lies ∼ 10% below the LHCb data [51]. The B-decay contribution to the total number of tau neutrino events amounts to less than 10%, as shown in the following, so to approximate bottom production, we use the same scale and k T choices as for charm production, with the replacement m c → m b . Refinements in the parameters used for the description of the B-mesons will not change the conclusions of this paper.
Alternatively, one could also consider varying N F and k T , with N R fixed, to find the two-parameter combination which provides the best fit to the LHCb 13 TeV D s data on double-differential cross sections in p T and y at √ s = 13 TeV. We find that N F = 1.44 and k T = 2.23 GeV is the best fit in that case when N R = 1, with χ 2 /DOF=2.68 and with corresponding predicted σ cc for 1 GeV< p T < 8 GeV and 2.0 < y < 4.5 amounting to 87% of the central value of the experimental result by the LHCb collaboration, which they extrapolate from D s data. The 3σ allowed interval of k T ( k 2 T ) values turns out to be an ellipse that spans the range 2.02-2.44 GeV (5.20-7.58 GeV 2 ) for N F = 1.26-1.62, favoring a lower N F and higher k T (and vice versa) for excursions from the best fit when N R = 1.0. An examination of the differences between the full charm meson production data from LHCb, their total charm-anticharm pair production cross sections, and theoretical predictions varying simultaneously the three parameters N F , N R and k T , would be interesting to better understand the roles of µ R , µ F and k T in theoretical predictions of charm production [52]. 4 Since our focus here is on ν τ +ν τ production, in the following JHEP06(2020)032 sections we use N F = 1.5, N R = 1.0 and k T = 2.2 GeV as representatives values of the range of parameter choices that lead to a reasonable description of the experimental LHCb data, together with more widely adopted choices of the same input parameters. As we will see, k T = 2.2 GeV leads to a suppression, relative to predictions with smaller k T , of neutrino production in the forward region.
The large value of k T = 2.2 GeV is difficult to reconcile with theoretical expectations for the strong interaction effects that we approximately model with the Gaussian factor. Large k T values have been used in some analyses [36,38,39]. For example, the NLO evaluation of direct photon production in pp collisions at the Tevatron, without resummation effects but including k T -smearing with a Gaussian function, shows good agreement with CDF and D0 data when k T = 2.5 GeV ( k 2 T = 8.0 GeV 2 ) for the photon [36]. A smaller value of k 2 T = 1 GeV 2 for the gluon PDF in the unpolarized proton is used to describe polarized proton-unpolarized proton scattering to J/ψX and DX in ref. [53]. On the other hand, in an analysis of di-J/ψ production at LHCb [54], the unpolarized transverse momentum gluon distribution, factorized in terms of the usual collinear PDF and a Gaussian as in eq. (3.2), yields k 2 T = 3.3 ± 0.8 GeV 2 [55]. Non-perturbative Sudakov factors are introduced on top of resummation procedures to describe the transverse momentum distributions of Drell Yan and W, Z production. For example, non-perturbative parameters in impact parameter space that translate to values of k 2 T in the range of ∼ 1-2 GeV 2 are obtained in ref. [56] to augment the corrections from resummations to better fit the experimental data. One interpretation of the need for a large value of k T is that it compensates for missing higher-order perturbative QCD corrections.
The effect of a more conservative value of k T = 0.7 GeV ( k 2 T = 0.6 GeV 2 ) is shown in the lower left panel of figure 1 along with the scale uncertainty band, again with our default central scale choices (µ R , µ F ) = (1.0, 1.5)m T . Figure 2 shows predictions for energy distributions for D s production in the very forward region, with η > 6.87, using different ap--10 -JHEP06(2020)032 proaches, at two different center-of-mass energies. Predictions of a computation relying on NLO QCD matrix-elements matched to the parton shower and hadronization algorithms implemented in Pythia, with matching performed according to the Powheg [57,58] method, are compared to those obtained by combining NLO QCD results with the Gaussian transverse momentum smearing model using k T = 0, 0.7 and 1.4 GeV, followed by fragmentation. All the distributions shown use the same default central scales and charm quark mass. In the NLO + shower Monte Carlo computation, the parameters related to fragmentation and intrinsic transverse momenta are tuned to pre-existing experimental data. The parton shower algorithm accounts for the effect of multiple soft and collinear partonic emissions on top of the hard-scattering process, affecting the kinematics and the dynamics of the event, inducing modifications of the distributions at the fixed-order level. The energy distributions of the Powheg + Pythia computation for D s production in the forward region agree well with our Gaussian smeared NLO predictions when k T = 0.7 GeV, as shown in both panels of figure 2. The lower left panel of figure 1 shows that the lower value of k T agrees less well with the data at higher transverse momentum but agrees well in case of lower p T of the meson, the kinematic region that dominates in the total cross section and in the production of a forward neutrino beam.
As already mentioned, it is more conventional to use scales equal to (µ R , µ F ) = (1.0, 1.0)m T . We show in the lower right panel of figure 1 the double differential distributions d 2 σ/dp T dy for D s production in different rapidity bins, rescaled as above, obtained using as input of our computation these scales and the same value of k T = 0.7 GeV as in the lower left panel. The conventional scale combination (µ R , µ F ) = (1.0, 1.0)m T typically gives rise to predicted cross sections lower than the LHCb data, at least considering the present limited accuracy of the theoretical calculations, but the data are still included within the large theoretical scale uncertainty bands.
To allow for comparisons with what could be considered more conventional k T , we show results for both k T = 0.7 GeV, a value consistent with a fully non-perturbative interpretation of k T effects that most closely follows the Powheg + Pythia result, and for k T = 2.2 GeV, the value corresponding to our best fit of LHCb experimental data. We also show some results with k T = 0-1.4 GeV. Also for comparisons with other QCD evaluations of heavy-flavor production, we show results for both central scales (µ R , µ F ) = (1.0, 1.5)m T and for conventional central scales (µ R , µ F ) = (1.0, 1.0)m T , with their respective uncertainty bands that cover the seven point scale variation of 0.5-2 times the central scale. With an interpretation of k T = 2.2 GeV, rather than with a value ∼ m p , as compensating for missing higher order QCD corrections, the seven point scale variations encompass k T uncertainties, at least to a degree.
Finally, we illustrate the effect of fragmentation in figure 3. In the left panel, the charm and D s energy distributions are shown for η > 4.5 at √ s = 14 TeV, while on the right, the energy distribution is shown when η > 6.87. For both panels, the fragmentation fraction for c → D s is set to unity to allow a direct comparison of the charm and meson distributions. The right panel of figure 3 shows that the impact of the fragmentation function extends to very high energy in the very forward region.  [45]. The decay of the τ itself is also prompt and produces a tau neutrino, which is the dominant source of high-energy tau neutrinos in the forward region since the tau carries most of the energy of the D s , given m Ds = 1.97 GeV and m τ = 1.78 GeV. The energy distribution of the ν τ that comes directly from the D s leptonic decay (called the "direct" neutrino) is straightforward to calculate from the isotropic decay of the D s in its rest frame, followed by a boost to the collider frame where the D s has a four-momentum p D . Keeping the polarization of the tau [59], the energy distribution of the ν τ from the tau decay (here called the "chain" decay neutrino) can also be obtained. The latter distribution, after integrating over angles relative to the tau momentum direction, is discussed in refs. [60,61] in the context of the atmospheric tau neutrino flux. Here, we use the full energy and angular distribution in the collider frame to apply the requirement that the neutrino pseudorapidity fulfills the η > 6.87 constraint. Details of the tau neutrino energy and angular distribution from D s → τ → ν τ appear in ref. [62]. The distributions of tau neutrinos and antineutrinos from the D s decays are identical because of the zero-spin of the meson (direct neutrinos) and the compensation of particle/antiparticle and left-handed/right-handed effects in the chain decays of the tau [59].
B-meson production and decay also contributes to the number of tau neutrinos, but at a level of less than 10% of the contribution from D s production and decay. The branching fractions to taus from charged and neutral B's are B(B + →D 0 τ + ν τ ) = (7.7 ± 2.5) × 10  all the branching fractions. The energy and angular distributions of tau neutrinos from B meson decays are discussed in appendix A, as are the energy and angular distributions of muon and electron neutrinos from heavy-meson decays. Figure 4 shows the energy distributions for tau neutrinos and antineutrinos from both the D s and B meson decays with neutrino pseudorapidity η > 6.87. The contribution of the direct and the chain decays are shown separately, as well as their sum. The left panel shows the distributions with our default scale combination (µ R , µ F ) = (1.0, 1.5) m T , while the right panel shows the same distributions, but using as input (µ R , µ F ) = (1.0, 1.0) m T . Qualitatively, the distributions are similar, although, as expected from the discussion in section 3, the differential distributions using the conventional scale are lower than with our default scale choice.
The number of neutrino events per unit energy can be written as where the interaction probability in the detector is Here, we use an integrated luminosity L = 3000 fb −1 and a lead density ρ P b = 11.34 g/cm 3 . The nucleon number of lead is A P b = 208 and l det is the length of the lead neutrino target in the detector which is also characterized by a cross sectional area of radius 1 m (thus η > 6.87) for our discussion here. For reference, a detector of lead with radius of 1 m and length l det = 1 m has a mass of ∼ 35.6 ton. For the number of events, we quote the number of events per ton of lead (per 2.8 cm depth of a lead disk with radius 1 m). As discussed in more detail in ref. [62], we evaluate the neutrino and antineutrino charged-current cross sections at NLO, including mass effects [63][64][65][66][67], using the nCTEQ15 PDFs for lead [47]. For -13 -JHEP06(2020)032 neutrino energies above ∼ 10 GeV, deep inelastic scattering (DIS) dominates quasi-elastic scattering and few-pion production [68,69]. At the energies of interest, the neutrino DIS cross section is roughly a factor of ∼ 2 larger than the antineutrino cross section. Kinematic corrections due to the tau mass reduce the charged-current cross section by ∼ 25% for a neutrino energy E ν = 100 GeV, and by ∼ 5% for E ν = 1000 GeV [66]. NLO corrections are small, less than ∼ 5% for neutrino energies of 100 GeV, as shown e.g. in ref. [65].
In figure 5, we show the number of events per unit neutrino energy per ton of detector lead for tau neutrinos and antineutrinos from the D ± s (left panels) and B meson (right panels) decays, evaluated for a total integrated luminosity L = 3000 fb −1 . The upper panels include results with our default scales (µ R , µ F ) = (1.0, 1.5) m T , while the lower panels show the same quantities using as input the conventional scale combination (µ R , µ F ) = (1.0, 1.0) m T . The central solid histograms are obtained for k T = 0.7 GeV, and the bands reflect the uncertainty range due to the variation of k T in the range 0-1.4 GeV. We also present predictions for k T = 2.2 GeV, that are shown with the dashed histograms.

JHEP06(2020)032
Incorporating k T effects has a large impact on the predictions at low energies as expected from figure 2. In particular, computing charm production with k T = 0 GeV, enhances the number of events per unit energy by up to ∼ 8% at E ν ∼ 100 GeV with respect to the case with the default k T = 0.7 GeV, whereas the differences are less than 1% for E ν 1000 GeV. The lower limit on the number of events per unit energy, corresponding to the case k T = 1.4 GeV, is lower by ∼ 16 (5)% at E ν ∼ 100 (1000) GeV with respect to the case with k T = 0.7 GeV, whereas the difference reduces to 1% or even less for E ν 1500 GeV. On the other hand, if k T = 2.2 GeV, the number of events from D ± s in the peak of the distribution is a factor of ∼ 2/3 lower than in the peak for k T = 0.7 GeV. The k T sensitivity of the predictions of the number of (ν τ +ν τ ) neutrinos from B meson decays is much smaller than for D s meson decays. A comparison of upper and lower panels shows the stronger impact of the factorization scale dependence on the predictions for ν τ +ν τ from charm mesons than from B mesons. Figure 6 shows the uncertainty bands associated with the QCD scale variation in the range considered in figure 1 For neutrinos from the D s decay, the upper boundary of the band is larger than the central prediction by a factor of ∼ 3-4, while the lower edge of the band is 40 − 60% smaller than the central prediction for E ν 1500 GeV. Thus, the QCD scale uncertainty band in our evaluation with k T = 0.7 GeV has overlap with the prediction using k T = 2.2 GeV, as follows from comparing figures 5 and 6. For neutrinos from B meson decays, the scale uncertainty bands are smaller than for neutrinos from D mesons. In particular, the edge of the upper uncertainty band is a factor 1.5-2 larger than the central prediction, whereas the lower uncertainty band extends to about 20% below the central prediction for E ν 1000 GeV. The scale uncertainty bands around central predictions with the conventional scale combination (µ R , µ F ) = (1.0, 1.0)m T , as shown in figure 7, have similar sizes to those in figure 6. However, the conventional scale choice leads to an overall lower central prediction for the number of events from both D s and B → ν τ andν τ than the case (µ R , µ F ) = (1.0, 1.5)m T , as follows by comparing figures 6 and 7.
In tables 1 and 2, we present the total event numbers for our preferred scale choice (µ R , µ F ) = (1.0, 1.5) m T and the conventional QCD scale choice (µ R , µ F ) = (1.0, 1.0) m T , respectively, for 1 m of lead (35.6 ton). Assuming an integrated luminosity of 3,000 fb −1 , each table shows separately the number of ν τ andν τ events from D ± s and B meson decays, respectively, and the total number, each for selected (µ R , µ F ) scale and k T choices. Table 1 shows that when using (µ R , µ F ) = (1.0, 1.5)m T and the k T value that yields the best match to the Powheg + PYTHIA D ± s energy distribution in the far forward η > 6.87 region ( k T = 0.7 GeV), the total number of tau neutrino plus antineutrino events is predicted to be ∼ 3, 800. k T = 0-2.2 GeV, produces a smaller uncertainty in the number of events, which span the range ∼ 3, 900-2, 900, with the lower end of the range of number of events corresponding to the (µ R , µ F and k T )-combination that is favored by LHCb data. Given our expectation that, for k T = 2.2 GeV, some effects of missing higher-order perturbative corrections are compensated by the large value of this parameter, an uncertainty in k T cannot be added -16 -JHEP06(2020)032 to the scale uncertainty to get a total uncertainty for the number of events. In any case, the scale uncertainty yields a range of number of events that is significantly larger than the range found by varying k T from 0 to 2.2 GeV.
We set (µ R , µ F ) = (1.0, 1.0) m T , the conventional scales used in most QCD theory evaluations of heavy-flavor production, for the predictions in table 2. This table shows as well a broad range of predicted number of events, with a central value of ∼ 2, 500 events for k T = 0.7 GeV. This number amounts to a factor of ∼ 2/3 of the one obtained in the evaluation with (µ R , µ F ) = (1, 1.5) m T for k T = 0.7 GeV, however, it is comparable to the 2, 900 events for k T = 2.2 GeV with (µ R , µ F ) = (1, 1.5) m T .

Muon neutrinos
At the interaction region, the muon neutrino and electron neutrino fluxes from heavyflavor production and decay will be nearly the same, coming primarily from the neutral and charged D semileptonic decays. Figure   factor of ∼ 1/20 at high energy. As in figure 4 for ν τ +ν τ , figure 8 shows that the predicted energy distribution of ν µ +ν µ from charm is ∼ 1.5 larger for (µ R , µ F ) = (1.0, 1.5) m T than for (µ R , µ F ) = (1.0, 1.0) m T , while the smaller contributions from B hadrons are much less sensitive to the scale choice.
The corresponding predictions for the number of muon neutrino and antineutrino charged-current events per unit energy for a 1 ton lead target for η > 6.87 are shown in figures 9 and 10. Figure 9 shows the number of muon neutrino and muon antineutrino events per unit energy from heavy flavor, including uncertainty bands from k T variation in the range [0, 1.4] GeV for (µ R , µ F ) = (1.0, 1.5) m T . Again, the dashed histograms correspond to results with k T = 2.2 GeV. The left panels in figure 10 show, from top to bottom, the sum of ν µ +ν µ charged current events per ton of lead, for ν µ and forν µ , all for (µ R , µ F ) = (1.0, 1.5) m T . The right panels show the same, but with (µ R , µ F ) = (1.0, 1.0) m T . Each panel includes a wide uncertainty band reflecting the theoretical uncertainties associated with scale variation. Our evaluation using NLO perturba--19 -JHEP06(2020)032  tive QCD gives a factor of ∼ 13 in the ratio of ν µ +ν µ to ν τ +ν τ events based on heavy flavor alone, with our default input parameters.
Heavy-flavor hadron decays, however, are not the only sources of ν µ +ν µ . In principle, pion, kaon and weak gauge boson decays can also contribute. In ref. [10], the contributions from W and Z production and decay to neutrinos are studied. They show that neutrinos from weak gauge boson decays populate mostly the pseudorapidity range of |η| < 4.5, but are negligible in the region η > 6.5 [10].
At first glance, one might expect that pion and kaon decays will not be important sources of neutrinos. Pions with E π > 9 GeV have γcτ > 500 m. Charged kaons have γcτ > 500 m for E K > 67 GeV. Pion and kaon contributions to the number of ν µ +ν µ at high energy were neglected in the earlier work of ref. [14]. Park, in ref. [5], used Pythia to account for pions and kaons, requiring the decay to occur within 50 m of the interaction point to guarantee that particles decay inside the beam pipe, leading to a factor of ∼ 100 times more ν µ +ν µ events than ν τ +ν τ events.
In ref. [18], a more detailed evaluation of ν µ +ν µ events is performed for a detector with 25 × 25 cm 2 cross-sectional area at 480 meters from the Atlas interaction point. They find a factor of ∼ 1000 more interactions by ν µ +ν µ than by ν τ +ν τ when pions and kaons are included, although with a different energy distribution. Their evaluation of heavy-flavor contributions is done using Pythia [17], while their light hadron production is estimated using the Crmc [70] simulation package with Epos-LHC [71], Qgsjet-II-04 [72] and Sibyll 2.3c [73][74][75]. They find that most of the neutrinos with energies above 1 TeV from charged light hadron decays come from pion and kaon decays in a region within ∼ 55 m from the Atlas interaction point. They find that magnetic fields sweep lower-energy (below ∼ 100 GeV) charged particles away if the particles have travelled 20 m downstream from the interaction point and have passed through the front quadrupole absorber of inner -20 -JHEP06(2020)032 -21 -JHEP06(2020)032 radius 17 mm [18]. They also find that two-body decays of charged pions and charged kaons are the dominant sources of ν µ +ν µ production.
Our primary focus here is on the heavy-flavor contributions to the number of events from ν µ +ν µ and ν τ +ν τ . Instead of simulating light charged hadron trajectories in magnetic fields, guided by the results of ref. [18], we approximately evaluate the number of ν µ +ν µ that come from pions and kaons as follows. We evaluate the π ± and K ± two-body decay contributions to the flux of ν µ +ν µ with the requirement that the mesons decay within 55 m of the interaction point and the decaying hadron's momentum lies within an angle of θ < 1 mrad relative to the beam axis to stay within the opening of the quadrupole absorber. While the 2-body light meson decays dominate, a more complete calculation of the light meson contributions to ν µ +ν µ would include K L semileptonic decays.
We use the parametrization of Koers et al. [34], based on fits to Pythia distributions for charged pions and kaons as a function of energy and rapidity, to generate the pion and kaon distributions. For reference, pion and kaon charged-particle multiplicities per interaction in pp collisions at √ s = 14 TeV are ∼ 50 and ∼ 6, respectively [34]. Two-body decays are implemented, as for the D ± s → ν τ τ case, with the requisite changes to the initial and final particle masses. Compared to figure 8, charged pions and kaons contribute a factor of ∼ 100 more ν µ +ν µ than heavy flavors do, as shown by the blue and red curves in figure 11. We will turn to the issue of the oscillations of neutrinos of different flavors to tau neutrinos. To estimate the number of ν e +ν e , we consider the dominant contribution which is from K L → πeν e . With the same geometry requirements, we show with the green curve in figure 11 the contribution from K L production and decay into ν e +ν e . We use the Koers distributions for K + + K − , then divide by two for K L . The three-body semileptonic kaon decay is evaluated following ref. [76]. The peak of the electron neutrino distribution from K L decays is about a factor of ∼ 2 larger than the peak of the electron neutrino distribution from heavy-flavor decays. While K L semileptonic decays to ν e are dominant, a more complete calculation would include the K + → ν e semileptonic decays as well.
Three-flavor neutrino oscillations of the much larger number of ν µ +ν µ from charged pions and kaons to ν τ +ν τ could, in principle, overwhelm the number of ν τ +ν τ from heavyflavor decays. However, the baseline of 480 m is not at all optimal for ν µ → ν τ oscillations in the standard scenario with 3 active flavors. In two flavor approximation, for a mass squared difference of ∆m 2 32 2.5 × 10 −3 eV 2 and a mixing angle θ 23 π/4, the oscillation probability is P (ν µ → ν τ ) 2.3 × 10 −6 /(E 2 /GeV 2 ), so for energies above E νµ = 10 GeV, ν µ → ν τ oscillations give a negligible contribution to the number of ν τ +ν τ , even given the large theoretical uncertainties in the number of ν τ +ν τ events per unit energy. The ν e → ν τ oscillation probability is even smaller.
In the next section, we use the approximate numbers of electron and muon neutrinos from light mesons, together with our heavy flavor results for all three neutrino flavors and turn to tau neutrino oscillations in a 3+1 neutrino mixing framework. This illustrates an example of a signal of new physics that could be probed by a forward neutrino detector at the LHC when heavy-flavor uncertainties are under better theoretical control. While not necessary for an analysis of ν τ +ν τ because tau neutrinos come from D s and B decays, for an analysis of all three neutrino flavors, a full accounting of light meson production and π ± →ν µ /ν µ K ± →ν µ /ν µ K L →ν e /ν e Figure 11. The ν µ +ν µ energy distributions from charged pions (blue) and kaons (red) and the ν e +ν e energy distributions from K L semileptonic decays (green) for forward neutrinos (η > 6.87) from the decay of mesons produced in pp interactions at √ s = 14 TeV, considering those mesons whose decay occurs within 55 m of the interaction point and whose momentum lies within an angle of θ < 1 mrad from the beam axis. The meson energy and rapidity distributions are evaluated using the parametrization of Koers et al. [34]. decay to neutrinos, using all available data on forward production (e.g., from LHCf [77], TOTEM [78] and CMS [79]) along with magnet and detector configurations in the interaction region, will be necessary. A forward tune of Pythia that is underway [19] would also guide future work.

New physics
The detection of a large number of identifiable ν τ +ν τ events would offer opportunities to explore a new corner of parameter space for neutrino oscillations into a fourth "sterile" neutrino. The baseline of 480 m is most sensitive to a fourth neutrino mass m 4 of the order of tens of eV.
The flavor eigenstates ν l can be expressed as a superposition of the mass eigenstates ν j according to the formula where n ν is the total number of neutrinos subject to oscillations: for the standard oscillation scenario with three active neutrinos, n ν = 3, while for an oscillation scenario with three active neutrinos plus one sterile neutrino, n ν = 4. The full transition probability appears in, for example, ref. [1]. The blue lines in the two panels of figure 12 show that tau neutrino survival probability P (ν τ → ν τ ) in the three-flavor scenario approaches unity in the energy range The transition probability in the scenario with three active and one sterile neutrino flavor, in the case when the mass eigenstates fulfill the hierarchy m 4 m 1,2,3 , can be simplified to [81]  are acceptable according to current IceCube constraints [82,83]. Matter effects can be neglected with the density of the Earth's crust ρ ≈ 2.6 g/cm 3 and the electron fraction Y e ≈ 0.5 [84,85].
As the lower panel of figure 12 shows, for increasing m 4 masses, oscillations become rapid over the full neutrino energy range, even at low energies, so their average determines the ν τ survival probability. Given the uncertainties in the absolute scale of the ν τ +ν τ flux in the very forward direction, an average decrease in the number of events due to oscillations into sterile neutrinos would be difficult to extract with measurements of forward LHC neutrinos. For this reason, we focus on our example of the case m 4 = 20 eV. Figure 13 shows the number of events as a function of energy, in the standard threeactive-flavor oscillation framework (black histogram) and in the 3+1 oscillation framework with m 4 = 20 eV and |U τ 4 | 2 = 0.15, considering our default heavy-flavor QCD input parameter set and k T = 0.7 GeV (left) and 2.2 GeV (right). The orange-dashed histogram shows the effect of ν τ disappearance due to oscillations in a 3+1 scenario with |U e4 | 2 = |U µ4 | 2 = 0, where a dip is visible at neutrino energies ∼ 150 GeV.
The spectral shape of the number of ν τ +ν τ events changes with different choices of mixing parameters. If |U e4 | 2 = |U µ4 | 2 = 0, spectral distortions will be conspicuous, even if k T = 2.2 GeV (corresponding to a total smaller number of events than our default case k T = 0.7 GeV), as can be seen in the right panel of figure 13. Keeping the same value of |U 2 e4 | and increasing the |U µ4 | 2 value to |U µ4 | 2 < 10 −3 , the ν τ → ν τ oscillation dip in the energy distribution of the events is not significantly modified because muon mixing parameters are quite constrained as discussed above. This effect is shown in the green histogram of figure 13, obtained by setting |U e4 | 2 = 0 and |U µ4 | 2 = 5 × 10 −4 . When |U e4 | 2 = 0, electron neutrinos coming primarily from K 0 e3 that oscillate to ν τ because of sterile neutrino mixing additionally fill in the spectral distortion of the ν τ +ν τ event number as shown in figure 13 with the solid blue histogram for |U e4 | 2 = 6 × 10 −3 .
Untangling the physics of a 3 + 1 flavor oscillation scenario, from the standard 3-flavor oscillation scenario, considering the QCD theoretical uncertainties related to the choice of the scales and of the phenomenological k T smearing parameter, will be difficult if the mixing parameters |U e4 | 2 and |U µ4 | 2 are close to 6 × 10 −3 and 5 × 10 −4 , respectively.
Although the values of |U e4 | 2 and |U µ4 | 2 are small, the large number of ν e and ν µ from light-meson production and decay means that ν e and ν µ mixing with sterile neutrinos can have a non-negligible impact. The challenge is illustrated in figures 14 and 15. These figures show the numbers of tau neutrino plus antineutrino charged-current events, with each distribution normalized to the corresponding number of events at 1 TeV for that distribution, to highlight the effect of the scale dependence and of various values of k T on the shape of the distribution. In figure 14, k T = 0.7 GeV is fixed and the scales are varied. In particular, although a peak is still visible, the case with all three elements |U 4 | 2 = 0 in

Conclusions
Theoretical proposals to exploit collider production of tau neutrinos and antineutrinos via heavy-flavor decays in the far-forward region have a long history. Recent proposals of experiments to detect BSM particles with feeble interactions have made our evaluation timely. This work provides a first evaluation of the number of ν τ +ν τ charged-current interaction events in the very forward region, including QCD effects beyond the leading order/leading logarithmic accuracy considered in previous estimates and studying the effects of different sources of QCD uncertainties previously neglected. We focus on neutrinos with pseudorapidities η > 6.87, consistently with the geometry of a 1 m radius cylindrical detector at a distance of 480 m from the interaction point [6][7][8] to illustrate a number of effects.
Thousands of tau neutrino plus antineutrino events are predicted for 1 meter of lead (35.6 ton) target. Theoretical uncertainties are quite large due to, in particular, the renormalization and factorization scale variation in the heavy-flavour production cross-sections. amounting to a factor of ∼ 4-6 between the upper and lower limit of the predicted event numbers, for the intrinsic momentum parameter with which our NLO predictions convoluted with fragmentation functions approximately reproduce Powheg + Pythia results.
The introduction of the parameter k T in the Gaussian smearing factor f ( k T ) in eq. (3.1) distorts the far-forward tau neutrino and antineutrino spectra. This is also the case in fixed-target experiments like SHiP [62]. The collinear parton model is sufficient for central collisions, but for forward production, as we have shown, non-collinear effects have a significant impact on the number of events for forward neutrino production at the LHC. Although the Gaussian form of f ( k T ) is imperfect, it can approximate non-perturbative QCD effects and mimic part of the perturbative effects beyond fixed-order. A comparison with LHCb double-differential cross sections in p T and rapidity for D s production in pp collisions at √ s = 13 TeV, for p T ∈ [0, 14] GeV in five rapidity bins in the range y = 2.0-4.5, shows that the experimental data are reasonably reproduced by theoretical predictions in a framework combining NLO pQCD corrections to the hard-scattering, k T -smearing effects and phenomenological fragmentation functions. Experiments that probe heavyflavor physics in the very forward region (e.g., ref. [91]) will help to better quantify and disentangle the need of accurate procedures jointly resumming different kinds of logarithms and that of more rigorous descriptions of non-perturbative QCD effects, as well as the limits of the collinear approximation. The inclusion of both small x and k T effects in theoretical evaluations of heavy-flavor production at the Electron-Ion Collider (EIC) have implications for measurements (see, e.g., refs. [92,93]). In another arena, predictions of the prompt atmospheric neutrino flux [49,50,61,[94][95][96] will also be constrained by measurements of the tau neutrino energy spectrum in the far-forward region at the LHC.
For a baseline of ∼ 500 m, oscillations in the standard 3-flavor scenario are negligible at the energy scales of the tau neutrino beams in our setup. On the other hand, in a 3+1 oscillation scenario, including a sterile neutrino in the O(20 eV) mass range with three active flavor eigenstates, oscillations could give rise to visible signals in the energy dependence of the tau neutrino events. The contribution to ν τ events due to ν µ → ν τ oscillations is not important due to the small value of |U µ4 | 2 . This is the case even though there are ∼ 100 times more ν µ +ν µ produced from charged pion and kaon decays than from heavy-flavor decays. The number of ν µ +ν µ from heavy flavor decays is larger than the number of ν τ +ν τ by a factor of ∼ 10. The ν e → ν τ oscillation effects are also small. Therefore a large number of ν τ +ν τ events in this energy range would provide an opportunity to constrain 3 + 1 oscillation models with a sterile neutrino in the 10's of eV mass range. In the case of ν τ disappearance, the location of a dip in the charged-current event distribution as a function of tau neutrino energy will constrain the mass of the fourth mass eigenstate (mostly sterile neutrino) m 4 . The quantity |U τ 4 | 2 is currently poorly constrained. We showed that, as long as ν e → ν τ appearance is suppressed, for |U τ 4 | 2 = 0.15 in principle the oscillation effect would be unambiguous. However, uncertainties in heavyflavor production present challenges to precision constraints on the 3+1 sterile neutrino parameter space accessible to a far-forward neutrino experiment at the LHC. Practical aspects of tau neutrino detection in the high-luminosity environment will also be a challenge because of the muon neutrino background.

A Decay distributions
The two-body decays of the D s meson in its rest frame come from energy and momentum conservation. Keeping the polarization of the τ , the energy and angular distribution of the tau neutrino can also be obtained [59,97,98]. Detailed formulas for the direct ν τ and chain decay ν τ distributions appear in appendix B of ref. [62]. Decays in the D s and τ rest frames are appropriatedly boosted to the collider frame, including full angular dependence.
Neutrinos can be produced through the three-body semileptonic decay process of D and B mesons, D(B) → K(D)lν l (l = e, µ). The distribution of neutrinos from the decay in the rest frame is respectively given by where the fraction of energy transferred to the neutrino is x ν = 2E ν /m i and the hadron mass fraction is r h = m 2 f /m 2 i , with m i and m f being the hadron masses in the initial and final states, respectively. The corresponding cumulative distribution functions (CDF) are where D(r) = 1 − 8r − 12r 2 ln(r) + 8r 3 − r 4 . The CDF is used to determine the neutrino energy in the heavy meson rest frame, and along with an isotropic decay distribution, the neutrino four-momentum in the rest frame is boosted to the collider frame, where the heavy meson D(B) has four-momentum p D(B) . The B meson also decays to tau neutrinos via B → Dτ ν τ , and the τ lepton decays subsequently to ν τ . In this case, there is an additional effect due to the non-negligible mass of τ . With the additional mass term r τ = m 2 τ /m 2 B , one obtains the distributions of τ and ν τ from B decays: Here, x τ = 2E τ /m i is the fraction of energy transferred to tau. The expression of the CDF for τ and ν τ from B decays is too complicated to be presented here. Again, both energy and angular distributions are determined in the rest frame, then boosted to the collider frame. For reference, we show the fractional energy distribution (left) and CDFs (right) for ν µ , ν τ and τ in the rest frame of the decaying D and B in figure 16 for the three-body decays.
The energy fractions x span the following range: For electron neutrino (ν e ) and muon neutrino (ν µ ), the maximum value of x ν is simply