CERN-TH-2020-007, BNL-213623-2020-FORE

aainstitutetext: Department of Physics and Astronomy, University of Iowa, Iowa
City, IA 52242, USAbbinstitutetext: Brookhaven National Laboratory,
USAccinstitutetext: Università degli Studi di Firenze, Dipartimento di Fisica
e Astronomia & INFN, Firenze, Italyddinstitutetext: Theoretical Physics
Department, CERN, 1211 Geneva 23, Switzerland

# Far-forward neutrinos at the Large Hadron Collider

Weidong Bai b Milind Diwan c Maria Vittoria Garzelli d Yu Seon Jeong a Mary
Hall Reno<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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-$k_{T}$
effects are approximated by a simple phenomenological model introducing a
Gaussian $k_{T}$-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 $D_{s}^{\pm}$ production are compared
with LHCb data on double-differential cross-sections in transverse momentum
and rapidity. The uncertainties 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.

##  1 Introduction

Since the discovery of oscillation properties of neutrinos, their fundamental
roles in shaping the universe have become an important line of inquiry
pdg:2016 . A high-luminosity neutrino program with a flux of neutrinos
produced at an energy scale of few GeV is the focus of the large experimental
neutrino physics community working on the Deep Underground Neutrino Experiment
(DUNE) Acciarri:2015uup ; Acciarri:2016crz ; Abi:2018dnh . With a proton beam
energy of 80-120 GeV, pions are produced with the highest multiplicity in
proton-nucleon interactions. Thus the beams at DUNE are predominantly muon
neutrino (and muon antineutrino) beams, generated by charged pion decays.
Future measurements of muon neutrino disappearance and electron neutrino
appearance, and of the differences between neutrino and antineutrino rates,
will allow the extraction of elements of the Pontecorvo-Maki-Nakagawa-Sakata
(PMNS) mixing matrix with better precision than they are currently known. The
mixing of tau neutrinos with muon and electron neutrinos will be determined
indirectly through disappearance and appearance probabilities since at these
energies, tau neutrinos are not directly produced in the beam and the tau mass
threshold severely suppresses the number of tau neutrino charged-current
events at the far detector.

Recently, attention has turned to opportunities of measuring the interactions
of highly energetic tau neutrinos produced by $pp$ collisions at the LHC
Park:2011gh ; Feng:2017uoz ; Ariga:2018pin ; Ariga:2019ufm ; Buontempo:2018gta
; Beni:2019gxv ; Beni:2019pyp . Tau neutrino beams would allow for direct
tests of lepton universality and to explore tau neutrino PMNS mixing in the
traditional three-neutrino mixing paradigm and in scenarios including
additional exotic neutrinos. It was already recognized several decades ago
that through the production and prompt decays of the $D_{s}^{\pm}$ and $B$
mesons to taus and tau neutrinos, hadron colliders produce large fluxes of tau
neutrinos in the forward direction DeRujula:1984pg ; Winter:1990ry ;
DeRujula:1992sn ; Vannucci:1993ud . 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
antihadrons111Differences in the forward region between the total number of
hadrons and antihadrons from quark-antiquark 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 for meson-antimeson
asymmetries. We show below that the $\Lambda_{c}$ gives a small contribution
to neutrino production. The $D_{s}^{\pm}$ has no valence component, and
experimental studies by LHCb show a $D_{s}^{+}-D_{s}^{-}$ asymmetry below 1%
Aaij:2018afd , 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 DeRujula:1984pg , 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
$p\bar{p}$ 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 Sjostrand:2019zhc  tuned to
Tevatron data and rescaled to a $\sqrt{s}=14$ TeV center-of-mass energy gives
qualitatively consistent event rates Park:2011gh .

In the last few years, there is renewed interest in far-forward neutrino
production and detection Feng:2017uoz ; Ariga:2018pin ; Ariga:2019ufm ;
Beni:2019gxv ; Beni:2019pyp ; Buontempo:2018gta ; Abreu:2019yak ;
Abreu:2020ddv . The ForwArd Search ExpeRiment at the LHC (Faser), primarily
dedicated to searches for light, extremely weakly interacting particles
Ariga:2019ufm ; Ariga:2018pin ; Feng:2017uoz , 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 $\eta>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 $\eta>6.7$ Buontempo:2018gta .
For a half-cylinder, 2-meter long lead detector covering pseudorapidities
$\eta>6.7$, Beni et al. Beni:2019gxv  used Pythia 8 Sjostrand:2019zhc ;
Sjostrand:2014zea  to find $\sim 8,700$ tau neutrino events for a $3,000$ fb-1
integrated luminosity at the LHC. Other configurations for detectors, for
example, in $\eta$ ranges of $8-9.5$ and $7.4-8.2$ for XSEN Beni:2019pyp  and
in the range $7.2-8.7$ for SND$@$LHC Ahdida:2020evc  are also under
consideration. The prototype Faser-$\nu$ has $\eta>8.6$ Abreu:2019yak ;
Abreu:2020ddv . We adopt a minimum pseudorapidity $\eta>6.87$ in this work for
definiteness. Our results for $\eta>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 Kodama:2000mp ; Kodama:2007aa  and Opera Pupilli:2016zrb ;
Agafonova:2018auq . Both Super-Kamiokande and IceCube have reported signs of
tau neutrino appearance in their datasets, but the statistics is still very
limited Li:2017dbe ; Aartsen:2019tjl . 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 measurements for $E_{\nu}<370$ GeV pdg:2016  and the
IceCube Collaboration’s determination of the averaged cross-section for
neutrino plus antineutrino deep-inelastic scattering with nucleons for
$E_{\nu}=6.3-980$ TeV Aartsen:2017kpd  (see also, e.g., ref.
Bustamante:2017xuy ). LHC forward neutrinos will provide the first opportunity
for direct neutrino and antineutrino cross-section measurements for neutrino
energies up to $E_{\nu}\lesssim 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 $\eta>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 Nason:1987xz ;
Nason:1989zy ; Mangano:1991jk  and neutrino deep-inelastic-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}^{\pm}$
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}^{\pm}$ production Aaij:2015bpa  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 $\Lambda_{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}^{\pm}$ 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 $\sim 20$ eV.

Heavy-flavor decays are not the only sources of forward neutrinos, but they
dominate the forward tau-neutrino flux. The contributions from $pp\to W,Z$
production followed by $W,Z$ leptonic decays are negligible for
$\eta\mathrel{\hbox{\hbox to0.0pt{\lower 3.22916pt\hbox{$\sim$}\hss}
\kern-3.00003pt\raise 1.72218pt\hbox{$>$}}}6.5$ Beni:2019gxv . On the other
hand, for muon and electron neutrinos, charged pion ($\pi^{\pm}$) and kaon
($K^{\pm}$) decays to $\nu_{\mu}$ and $K_{L}\to\nu_{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 Koers:2006dd .

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.

##  2 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. Beni:2019gxv . This corresponds to a neutrino pseudorapidity of
$\eta>6.7$ for detection. The Faser2 proposal Ariga:2019ufm  has $R=1.0$ m,
corresponding to $\eta>6.87$, which we use for the results shown below. Other
smaller detectors like XSEN Beni:2019pyp  and SND$@$LHC Ahdida:2020evc  are
discussed in the recent literature. A prototype Faser-$\nu$ at even higher
$\eta$,222The Faser-$\nu$ experiment has recently been approved. The
Faser-$\nu$ detector will be installed in front of Faser. It is expected to be
fully operational and collect data during the LHC Run-III. with a 25 cm
$\times$ 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 Abreu:2019yak . A detector radius of $1.0$ meters, for a comparable target
mass, increases the number of events by a factor of $\sim 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}\to\nu_{\tau}\tau$ decay. The
$\tau\to\nu_{\tau}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
$\gamma c\tau_{D_{s}}\simeq E_{D_{s}}/m_{D_{s}}\cdot 150\ \mu$m, so of order
of $1.5-15$ cm for $E_{D_{s}}=200$ GeV$-2$ TeV. The tau decay length
$c\tau=87.11\ \mu$m, multiplied by the $\gamma$-factor for the same energy
range, gives a size of $0.98-9.8$ cm. Similarly, $\gamma
c\tau_{B^{+}}=E_{B^{+}}/m_{B^{+}}\cdot 496\ \mu$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\
\mu$m Bruning:2004ej  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 $\theta_{\rm max}=2.1$ mrad ($\eta_{\rm
min}=6.87$). This same constraint applies to the momenta of muon and electron
neutrinos from heavy-flavor decays.

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 sec. 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. Abreu:2019yak , an evaluation of the number of
$\nu_{\mu}+\bar{\nu}_{\mu}$ events in a detector of 25 cm $\times$ 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_{D_{s}}=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 $\langle k_{T}\rangle$ 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 $\langle k_{T}\rangle$ value higher than the naive
estimate due to Fermi motion and even larger than the upper end of the range
of typical non-perturbative $\langle k_{T}\rangle$ values of $\sim 1-2$ GeV
reported in the literature Apanasevich:1998ki ; Mangano:1998oia ; Miu:1998ju ;
Balazs:2000sz ; Skands:2010ak . 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
$\nu_{\tau}+\bar{\nu}_{\tau}$ production, we use the LHCb data for $D_{s}$
production at $\sqrt{s}=13$ TeV for rapidities in the range of $2.0-4.5$
Aaij:2015bpa .

##  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 Nason:1987xz ;
Nason:1989zy ; Mangano:1991jk .

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 Catani:1990eg ; Collins:1991ty . 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_{\rm NLO}(q_{T},y)=d^{2}\sigma({\rm NLO})/dq_{T}\,dy$ for
the quark, we evaluate the Gaussian smeared one-particle inclusive charm-quark
distribution $g(p_{T},y)=d^{2}\sigma/dp_{T}\,dy$ according to

$g(p_{T},y)=\int d^{2}\vec{k}_{T}\,f(\vec{k}_{T},\langle
k_{T}^{2}\rangle)g_{\rm NLO}(q_{T}=|\vec{p}_{T}-\vec{k}_{T}|,y)\ ,$ (1)

where $f(\vec{k}_{T},\langle k_{T}^{2}\rangle)$ is a two-dimensional Gaussian
function normalized to unity,

$\int\,d^{2}\vec{k}_{T}\,f(\vec{k}_{T},\langle
k_{T}^{2}\rangle)=\int\,d^{2}\vec{k}_{T}\,\frac{1}{\pi\langle
k_{T}^{2}\rangle}\exp[{-k_{T}^{2}/\langle k_{T}^{2}\rangle}]=1\ .$ (2)

The smearing generates a shift of the outgoing heavy-quark momentum vector
$\vec{q}_{T}$ by $\vec{k}_{T}$ after the hard scattering, before
fragmentation. The single Gaussian $f(\vec{k}_{T},\langle k_{T}^{2}\rangle)$
is equivalent to starting with two Gaussian functions, one for each incoming
parton ($k_{iT}$), making a change of variables to
$\vec{k}_{T}=(\vec{k}_{1T}+\vec{k}_{2T})/2$ for one of the integrals, and
integrating over the other parton’s $k_{T}$,

$d^{2}\vec{k}_{1T}\,f(\vec{k}_{1T},2\langle
k_{T}^{2}\rangle)\,d^{2}\vec{k}_{2T}\,f(\vec{k}_{2T},2\langle
k_{T}^{2}\rangle)\to d^{2}\vec{k}_{T}\,f(\vec{k}_{T},\langle
k_{T}^{2}\rangle)\ .$ (3)

This definition of the effective $\vec{k}_{T}$ variable distributes the
transverse momentum smearing equally to the one-particle transverse momentum
and to the recoil hadrons Apanasevich:1998ki . The quantity $\langle
k_{T}^{2}\rangle$ is related to the average magnitude of $\vec{k}_{T}$,
$\langle k_{T}\rangle$, by

$\langle k_{T}\rangle^{2}=\langle k_{T}^{2}\rangle\pi/4\ .$ (4)

As we discuss below, we set the $\langle k_{T}^{2}\rangle$ value based on
comparisons with double-differential distributions in rapidity and transverse
momentum for forward $D_{s}$ production measured by LHCb at $\sqrt{s}=13$ TeV
Aaij:2015bpa .

Charm quark fragmentation to mesons is accounted for by using fragmentation
functions of the Peterson form Peterson:1982ak . For $c\to D^{0}$, $D^{+}$ and
$D_{s}^{+}$, we take fragmentation fractions 0.6086, 0.2404 and 0.0802,
respectively Lisovyi:2015uqa . The parameter $\epsilon$ in the Peterson
fragmentation function that describes the hardness of the meson spectrum
relative to the heavy quark is taken to be $\epsilon$ = 0.028, 0.039 and 0.008
for $D^{0}$, $D^{+}$ and $D_{s}^{+}$, respectively, and $\epsilon=0.003$ for
$B$ production pdg:2019 . The fragmentation fraction for $B$’s is $b\to
B^{+}=b\to B^{0}=0.362$ from ref. Aaij:2019pqz .

We use the NLO nCTEQ15 parton distribution function (PDF) grids
Kovarik:2015cma , 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
($\mu_{R},\mu_{F}$) which are factors $(N_{R},N_{F})$ of the transverse mass
$m_{T}=\sqrt{m_{Q}^{2}+p_{T}^{2}}$, 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 $(\mu_{R},\mu_{F})=(N_{R},N_{F})m_{T}$, with
$N_{R}=N_{F}=1$ set as conventional scale factors Cacciari:2012ny . 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 $(\mu_{R},\mu_{F})$ used in most, (see, e.g., ref.
Cacciari:2012ny ), though not all (see, e.g., ref. Benzke:2017yjn ;
Zenaiev:2019ktw ), 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 $\langle
k_{T}\rangle$ = 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 $\langle
k_{T}\rangle$ = 0.7 GeV, $\langle k_{T}\rangle$ = 0 GeV and $\langle
k_{T}\rangle$ = 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 Aaij:2015bpa  are used to
set $\langle k_{T}\rangle$ for selected $(\mu_{R},\mu_{F})$ combinations,
after fixing the PDF and fragmentation function details. The $D_{s}$ data
cover the range $0<p_{T}<14$ GeV 333The LHCb measurements reported in Ref.
Aaij:2015bpa  extend to a $p_{T}$ value of 15 GeV for $D^{\pm}$, $D^{0}$, and
$\bar{D}^{0}$, while, for the less abundant $D_{s}^{\pm}$, the experimental
results in the largest $p_{T}$ bin are not reported. and $2.0<y<4.5$. While
varying the input parameters, a $\chi^{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 $\chi^{2}$ with $\mu_{R}=m_{T}$,
$\mu_{F}=1.5\,m_{T}$ and varying $\langle k_{T}\rangle$, leads to $\langle
k_{T}\rangle=2.2\pm 0.7$ GeV ($\langle k_{T}^{2}\rangle=6.2$ GeV2). This
represents a reasonably good fit, with $\chi^{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 $\sigma_{c\bar{c}}$ is within 2% of the central value of the LHCb
estimate reported in ref. Aaij:2015bpa .

Figure 1: Comparison between LHCb experimental data on double-differential
distributions in the meson $p_{T}$ and $y$ for $D_{s}^{\pm}$ production in
$pp$ collisions Aaij:2015bpa  and our QCD predictions. Data (and predictions)
for different $\Delta y$ bins are shifted by $10^{-m}$ where values of $m=0$,
2, 4, 6 and 8. The upper left panel refers to the case $N_{R}=1.0$ and
$N_{F}=1.5$, where $(\mu_{R},\mu_{F})=(N_{R},N_{F})\,m_{T}$, both fixing
$\langle k_{T}\rangle=2.2$ GeV and in the collinear approximation ($\langle
k_{T}\rangle=$ 0 GeV). The upper right and lower left panels refer to the same
scale configuration with the shaded band showing the uncertainty built by
seven-point scale variation in the range $N_{R,F}=0.5-2$ for $\langle
k_{T}\rangle=2.2$ GeV (best fit to data) and $\langle k_{T}\rangle=0.7$ GeV
(our default value), respectively. The lower right panel shows the same as
lower left, but with $(N_{R},N_{F})$ = (1,1) (conventional scale choice).

The predictions with $\mu_{R}=m_{T},\ \mu_{F}=1.5\,m_{T}$ and $\langle
k_{T}\rangle=2.2$ GeV, together with those for a $\langle k_{T}\rangle=0$
evaluation using the sames scales and the LHCb data Aaij:2015bpa , 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 $\Delta y=0.5$ width for
$2.0<y<4.5$ and normalization shifted by $10^{-m}$ where $m=0$, 2, 4, 6, 8.
The upper right panel of figure 1 shows, for $\langle k_{T}\rangle=2.2$ GeV,
the scale uncertainty band obtained as an envelope of seven combinations of
scales between factors of 0.5 and 2.0 of the central scales
$(\mu_{R},\mu_{F})=(1.0,1.5)\,m_{T}$, namely with the combinations of
$(N_{R},N_{F})$ equal to factors of $m_{T}$ of $(0.5,0.75)$, $(2,3)$,
$(1.0,0.75)$, $(0.5,1.5)$, $(1,3)$, $(2.0,1.5)$ and $(1.0,1.5)$
Cacciari:2012ny . The lower left panel of figure 1 shows the central results
and uncertainty band for the same scale choices, now with $\langle
k_{T}\rangle=0.7$ GeV, the transverse momentum smearing that approximates the
Powheg + Pythia results (see figure 2). Finally, the lower right panel shows
results with $\langle k_{T}\rangle=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 Zenaiev:2019ktw .
Results for the $p_{T}$ and $y$ distributions, in similar agreement with the
LHCb data, are obtained for $D^{\pm}$ and $D^{0},\ \bar{D}^{0}$ production,
used below as sources of $\nu_{\mu}+\bar{\nu}_{\mu}$ and
$\nu_{e}+\bar{\nu}_{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 $\sqrt{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
$c\bar{c}$ replaced by $b\bar{b}$, so it is generally expected that the
renormalization and factorization scaling factors are the same for both heavy-
quark flavors. To the extent that $\langle k_{T}\rangle$ describes the
intrinsic parton transverse momentum in the initial state protons, this
parameter also is process independent for $pp\to Q\bar{Q}X$ for $Q=c,\ b$. For
the same central scale and $\langle k_{T}\rangle$ choices, the $p_{T}$
distribution of $B^{\pm}$ mesons at $\sqrt{s}=13$ TeV for $y=2.0-4.5$ lies
$\sim 10\%$ below the LHCb data Aaij:2017qml . 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
$\langle k_{T}\rangle$ choices as for charm production, with the replacement
$m_{c}\to 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 $\langle
k_{T}\rangle$, 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 $\sqrt{s}$ = 13 TeV. We find that
$N_{F}=1.44$ and $\langle k_{T}\rangle=2.23$ GeV is the best fit in that case
when $N_{R}=1$, with $\chi^{2}$/DOF=2.68 and with corresponding predicted
$\sigma_{c\bar{c}}$ 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$\sigma$ allowed interval of
$\langle k_{T}\rangle\ (\langle k_{T}^{2}\rangle)$ values turns out to be an
ellipse that spans the range $2.02-2.44$ GeV ($5.20-7.58$ GeV2) for
$N_{F}=1.26-1.62$, favoring a lower $N_{F}$ and higher $\langle k_{T}\rangle$
(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 cross sections, and theoretical
predictions varying simultaneously the three parameters $N_{F}$, $N_{R}$ and
$\langle k_{T}\rangle$, would be interesting to better understand the roles of
$\mu_{R}$, $\mu_{F}$ and $\langle k_{T}\rangle$ in theoretical predictions of
charm production underway .444The parameter $\langle k_{T}\rangle$ is also
tied to the choice of fragmentation function since both influence the heavy-
meson rapidity and transverse momentum in the collider frame underway . Since
our focus here is on $\nu_{\tau}+\bar{\nu}_{\tau}$ production, in the
following sections we use $N_{F}=1.5$, $N_{R}=1.0$ and $\langle
k_{T}\rangle=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, $\langle k_{T}\rangle=2.2$ GeV leads to a suppression, relative to
predictions with smaller $\langle k_{T}\rangle$, of neutrino production in the
forward region.

The large value of $\langle k_{T}\rangle=2.2$ GeV is difficult to reconcile
with theoretical expectations for the strong interaction effects that we
approximately model with the Gaussian factor. Large $\langle k_{T}\rangle$
values have been used in some analyses Apanasevich:1998ki ; Miu:1998ju ;
Balazs:2000sz . For example, the NLO evaluation of direct photon production in
$p\bar{p}$ 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 $\langle k_{T}\rangle=2.5$ GeV ($\langle
k_{T}^{2}\rangle=8.0$ GeV2) for the photon Apanasevich:1998ki . A smaller
value of $\langle k_{T}^{2}\rangle=1$ GeV2 for the gluon PDF in the
unpolarized proton is used to describe polarized proton-unpolarized proton
scattering to $J/\psi X$ and $DX$ in ref. DAlesio:2017rzj . On the other hand,
in an analysis of di-$J/\psi$ production at LHCb Aaij:2016bqq , the
unpolarized transverse momentum gluon distribution, factorized in terms of the
usual collinear PDF and a Gaussian as in eq. (2), yields $\langle
k_{T}^{2}\rangle=3.3\pm 0.8$ GeV2 Lansberg:2017dzg . 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 $\langle k_{T}^{2}\rangle$ in the range of $\sim 1-2$ GeV2 are
obtained in ref. Bacchetta:2018lna  to augment the corrections from
resummations to better fit the experimental data. One interpretation of the
need for a large value of $\langle k_{T}\rangle$ is that it compensates for
missing higher-order perturbative QCD corrections.

Figure 2: Differential distributions of the energy $E$ of the $D_{s}$ produced
in $pp$ collision at $\sqrt{s}=$ 13 (left) and 14 (right) TeV, for
$D_{s}^{\pm}$ pseudorapidity $\eta>6.87$. Predictions obtained within our
framework using $(\mu_{R},\mu_{F})=(1.0,1.5)m_{T}$ and $\left\langle
k_{T}\right\rangle=0,\ 0.7,\ 1.4$ GeV, are compared to those of a computation
based on NLO hard-scattering matrix-element matched to parton shower followed
by hadronization, according to the Powheg + Pythia framework. See text for
more detail.

The effect of a more conservative value of $\langle k_{T}\rangle=0.7$ GeV
($\langle k_{T}^{2}\rangle=0.6$ GeV2) is shown in the lower left panel of
figure 1 along with the scale uncertainty band, again with our default central
scale choices $(\mu_{R},\mu_{F})=(1.0,1.5)m_{T}$. Figure 2 shows predictions
for energy distributions for $D_{s}$ production in the very forward region,
with $\eta>6.87$, using different approaches, 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 Frixione:2007vw ;
Frixione:2007nw  method, are compared to those obtained by combining NLO QCD
results with the Gaussian transverse momentum smearing model using $\langle
k_{T}\rangle=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
$\langle k_{T}\rangle=0.7$ GeV, as shown in both panels of figure 2. The lower
left panel of figure 1 shows that the lower value of $\langle k_{T}\rangle$
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.

Figure 3: The effect of fragmentation of charm quark to $D_{s}$ which shifts
the energy distribution to lower energy for both $\eta>4.5$ (left) and
$\eta>6.87$ (right). Shown are the differential cross sections of charm quarks
and $D_{s}$ mesons produced in $pp$ collision at $\sqrt{s}=14$ TeV for
$(\mu_{R},\mu_{F})=(1.0,1.5)m_{T}$ and $\left\langle k_{T}\right\rangle=0.7$
GeV, as a function of the energy $E$. Note that the $D_{s}$ fragmentation
fraction has been set to unity to demonstrate this effect clearly.

As already mentioned, it is more conventional to use scales equal to
$(\mu_{R},\mu_{F})=(1.0,1.0)m_{T}$. We show in the lower right panel of figure
1 the double differential distributions $d^{2}\sigma/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 $\langle
k_{T}\rangle=0.7$ GeV as in the lower left panel. The conventional scale
combination $(\mu_{R},\mu_{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
$\langle k_{T}\rangle$, we show results for both the $\langle
k_{T}\rangle=0.7$ GeV, a value consistent with a fully non-perturbative
interpretation of $\langle k_{T}\rangle$ effects that most closely follows the
Powheg + Pythia result, and for $\langle k_{T}\rangle=2.2$ GeV, the value
corresponding to our best fit of LHCb experimental data. We also show some
results with $\langle k_{T}\rangle=0-1.4$ GeV. Also for comparisons with other
QCD evaluations of heavy-flavor production, we show results for both central
scales $(\mu_{R},\mu_{F})=(1.0,1.5)m_{T}$ and for conventional central scales
$(\mu_{R},\mu_{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 $\langle k_{T}\rangle=2.2$ GeV, rather than with a
value $\sim m_{p}$, as compensating for missing higher order QCD corrections,
the seven point scale variations encompass $\langle k_{T}\rangle$
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 $\eta>4.5$ at
$\sqrt{s}=14$ TeV, while on the right, the energy distribution is shown when
$\eta>6.87$. For both panels, the fragmentation fraction for $c\to 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.

##  4 Neutrinos from the heavy-flavor hadrons

###  4.1 Tau neutrinos

Tau neutrinos and antineutrinos arise primarily from the prompt decays of
$D_{s}$ mesons into $\tau+\nu_{\tau}$, where
$B(D_{s}\to\tau\nu_{\tau})=(5.48\pm 0.23)\times 10^{-2}$ pdg:2019 . The decay
of the $\tau$ 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_{D_{s}}=1.97\ {\rm
GeV}$ and $m_{\tau}=1.78\ {\rm GeV}$. The energy distribution of the
$\nu_{\tau}$ 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
Barr:1988rb , the energy distribution of the $\nu_{\tau}$ 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. Pasquali:1998xf ; Bhattacharya:2016jce  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 $\eta>6.87$ constraint. Details of
the tau neutrino energy and angular distribution from
$D_{s}\to\tau\to\nu_{\tau}$ appear in ref. Bai:2018xum . 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 Barr:1988rb .

$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^{+}\to\bar{D}^{0}\tau^{+}\nu_{\tau})=(7.7\pm 2.5)\times
10^{-3}$, $B(B^{+}\to\bar{D}^{*}(2007)^{0}\tau^{+}\nu_{\tau})=(1.88\pm
0.20)\times 10^{-2}$, $B(B^{0}\to{D}^{-}\tau^{+}\nu_{\tau}=(1.08\pm
0.23)\times 10^{-2}$, and
$B(B^{0}\to{D}^{*}(2010)^{-}\tau^{+}\nu_{\tau})=(1.57\pm 0.10)\times 10^{-2}$
pdg:2019 . We use the central values of 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:  The neutrino energy distributions for tau neutrinos and
antineutrinos from the direct decay
$D_{s}^{\pm}/B^{0,\pm}\rightarrow\overset{{\scriptscriptstyle(-)}}{\nu}_{\tau}$
(green) and the chain decay
$D_{s}^{\pm}/B^{0,\pm}\rightarrow\tau^{\pm}\rightarrow\overset{{\scriptscriptstyle(-)}}{\nu}_{\tau}$
(orange) and their sum (blue) for neutrinos with pseudorapidity $\eta>6.87$,
produced in $pp$ collisions at $\sqrt{s}=14$ TeV. Predictions are obtained
using as input $\langle k_{T}\rangle=0.7$ GeV with our default scale
combination $(\mu_{R},\mu_{F})=(1.0,1.5)\,m_{T}$ (left) and with the
conventional scale combination $(\mu_{R},\mu_{F})=(1.0,1.0)\,m_{T}$ (right).
The contributions from both $D_{s}$ and $B$ mesons are shown separately.

Figure 4 shows the energy distributions for tau neutrinos and antineutrinos
from both the $D_{s}$ and $B$ meson decays with neutrino pseudorapidity
$\eta>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 $(\mu_{R},\mu_{F})=(1.0,1.5)\,m_{T}$, while the
right panel shows the same distributions, but using as input
$(\mu_{R},\mu_{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

$\frac{dN}{dE}=\frac{d\sigma(pp\rightarrow\nu X)}{dE}\times{\cal P}_{\rm
int}\times{\mathcal{L}}\ ,$ (5)

where the interaction probability in the detector is

${\cal P}_{\rm int}=\bigl{(}\rho_{\rm Pb}\times l_{\rm det}\times N_{\rm
avo}\bigr{)}\,\frac{\sigma_{\nu_{\tau}Pb}}{A_{\rm Pb}}\ .$ (6)

Here, we use an integrated luminosity ${\mathcal{L}}=3000{\ \rm fb}^{-1}$ and
a lead density $\rho_{Pb}=11.34{\ \rm g/cm^{3}}$. The nucleon number of lead
is $A_{Pb}=208$ and $l_{\rm 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 $\eta>6.87$) for our discussion here. For reference, a detector of
lead with radius of 1 m and length $l_{\rm det}$ = 1 m has a mass of $\sim$
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. Bai:2018xum , we evaluate the neutrino and antineutrino
charged-current cross sections at NLO, including mass effects Kretzer:2002fr ;
Kretzer:2003iu ; Jeong:2010za ; Jeong:2010nt ; Reno:2006hj , using the nCTEQ15
PDFs for lead Kovarik:2015cma . For neutrino energies above $\sim 10$ GeV,
deep inelastic scattering (DIS) dominates quasi-elastic scattering and few-
pion production Lipari:1994pz ; Kretzer:2004wk . At the energies of interest,
the neutrino DIS cross section is roughly a factor of $\sim 2$ larger than the
antineutrino cross section. Kinematic corrections due to the tau mass reduce
the charged-current cross section by $\sim 25\%$ for a neutrino energy
$E_{\nu}$ = 100 GeV, and by $\sim 5\%$ for $E_{\nu}$ = 1000 GeV Jeong:2010nt .
NLO corrections are small, less than $\sim 5\%$ for neutrino energies of 100
GeV, as shown e.g. in ref. Jeong:2010za .

Figure 5:  Our predictions for the tau neutrino and/or antineutrino number of
charged-current events per GeV per ton as a function of the incident neutrino
energy for neutrino pseudorapidity $\eta>6.87$ for $pp$ collisions with
$\sqrt{s}=14$ TeV and integrated luminosity ${\cal L}=3000$ fb-1. The central
predictions (solid histogram) refer to the $\left\langle k_{T}\right\rangle$ =
0.7 GeV value and the bands arise from the $\left\langle k_{T}\right\rangle$
variation in the range $0<\left\langle k_{T}\right\rangle<1.4$ GeV. The
results obtained with $\left\langle k_{T}\right\rangle$ = 2.2 GeV are also
shown (dashed histogram). The upper plots are obtained by setting the QCD
scales to $(\mu_{R},\mu_{F})=(1.0,1.5)m_{T}$, while the lower plots refer to
the conventional scale combination $(\mu_{R},\mu_{F})=(1.0,1.0)m_{T}$. NLO
corrections are included in the DIS cross sections.

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}^{\pm}$ (left
panels) and $B$ meson (right panels) decays, evaluated for a total integrated
luminosity ${\cal L}=3000\,{\rm fb^{-1}}$. The upper panels include results
with our default scales $(\mu_{R},\mu_{F})=(1.0,1.5)\,m_{T}$, while the lower
panels show the same quantities using as input the conventional scale
combination $(\mu_{R},\mu_{F})=(1.0,1.0)\,m_{T}$. The central solid histograms
are obtained for $\left\langle k_{T}\right\rangle$ = 0.7 GeV, and the bands
reflect the uncertainty range due to the variation of $\left\langle
k_{T}\right\rangle$ in the range $0-1.4$ GeV. We also present predictions for
$\left\langle k_{T}\right\rangle$ = 2.2 GeV, that are shown with the dashed
histograms.

Figure 6: Uncertainty range due to the QCD scale variation in the tau neutrino
and antineutrino number of charged-current events per GeV per ton as a
function of the incident neutrino energy for neutrino pseudorapidity
$\eta>6.87$ and for the $pp$ collision with $\sqrt{s}=14$ TeV and integrated
luminosity ${\cal L}=3000$ fb-1. The central predictions are obtained using as
input $(\mu_{R},\mu_{F})=(1,1.5)m_{T}$ and $\left\langle k_{T}\right\rangle$ =
0.7 GeV. The upper and lower limits arise from the QCD 7-point scale variation
range in figure 1. NLO QCD corrections are accounted for in the DIS pQCD cross
section.

Incorporating $\langle k_{T}\rangle$ effects has a large impact on the
predictions at low energies as expected from figure 2. In particular,
computing charm production with $\langle k_{T}\rangle=0$ GeV, enhances the
number of events per unit energy by up to $\sim 8\%$ at $E_{\nu}\sim$ 100 GeV
with respect to the case with the default $\langle k_{T}\rangle=0.7$ GeV,
whereas the differences are less than 1% for $E_{\nu}\mathrel{\hbox{\hbox
to0.0pt{\lower 3.22916pt\hbox{$\sim$}\hss} \kern-3.00003pt\raise
1.72218pt\hbox{$>$}}}$ 1000 GeV. The lower limit on the number of events per
unit energy, corresponding to the case $\langle k_{T}\rangle=$ 1.4 GeV, is
lower by $\sim$ 16 (5)% at $E_{\nu}\sim$ 100 (1000) GeV with respect to the
case with $\langle k_{T}\rangle=0.7$ GeV, whereas the difference reduces to 1%
or even less for $E_{\nu}\mathrel{\hbox{\hbox to0.0pt{\lower
3.22916pt\hbox{$\sim$}\hss} \kern-3.00003pt\raise 1.72218pt\hbox{$>$}}}$ 1500
GeV. On the other hand, if $\left\langle k_{T}\right\rangle$ = 2.2 GeV, the
number of events from $D_{s}^{\pm}$ in the peak of the distribution is a
factor of $\sim 2/3$ lower than in the peak for $\left\langle
k_{T}\right\rangle$ = 0.7 GeV. The $\left\langle k_{T}\right\rangle$
sensitivity of the predictions of the number of ($\nu_{\tau}$ \+
$\bar{\nu}_{\tau}$) 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
$\nu_{\tau}+\bar{\nu}_{\tau}$ 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 a central scale choice
$(\mu_{R},\mu_{F})=(1,1.5)\,m_{T}$. These bands are computed as envelopes of
seven combinations of $(N_{R},N_{F})$, equal to factors of $m_{T}$ of
$(0.5,0.75)$, $(2.0,3.0)$, $(1.0,0.75)$, $(0.5,1.5)$, $(1.0,3.0)$, $(2.0,1.5)$
and $(1.0,1.5)$. For neutrinos from the $D_{s}$ decay, the upper boundary of
the band is larger than the central prediction by a factor of $\sim 3-4$,
while the lower edge of the band is 40 $-$ 60% smaller than the central
prediction for $E_{\nu}\lesssim 1500$ GeV. Thus, the QCD scale uncertainty
band in our evaluation with $\left\langle k_{T}\right\rangle$ = 0.7 GeV has
overlap with the prediction using $\left\langle k_{T}\right\rangle$ = 2.2 GeV,
as follows from comparing figure 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_{\nu}\lesssim 1000$ GeV. The
scale uncertainty bands around central predictions with the conventional scale
combination $(\mu_{R},\mu_{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\to\nu_{\tau}$ and $\bar{\nu}_{\tau}$ than the case
$(\mu_{R},\mu_{F})=(1.0,1.5)m_{T}$, as follows by comparing figures 6 and 7.

Figure 7: The same as figure 6, but adopting the central scale conventional
choice $(\mu_{R},\mu_{F})=(1.0,1.0)\,m_{T}$ and seven-point scale variation
around it, again with $\left\langle k_{T}\right\rangle$ = 0.7 GeV.
| $\nu_{\tau}$ | $\bar{\nu}_{\tau}$ | $\nu_{\tau}+\bar{\nu}_{\tau}$ | $\nu_{\tau}+\bar{\nu}_{\tau}$
---|---|---|---|---
$(\mu_{R},\ \mu_{F})$ | (1, 1.5) $m_{T}$ | (1, 1.5) $m_{T}$ | (0.5, 1.5) $m_{T}$ | (1, 0.75) $m_{T}$
$\langle k_{T}\rangle$ | 0.7 GeV | 0 GeV | 1.4 GeV | 2.2 GeV | 0.7 GeV
$D_{s}$ | 2451 | 1191 | 3642 | 3799 | 3261 | 2735 | 11008 | 1716
$B^{\pm,0}$ | 96 | 46 | 142 | 144 | 137 | 127 | 214 | 115
Total | 2547 | 1237 | 3784 | 3943 | 3398 | 2862 | 11222 | 1831
Table 1: The charged-current event numbers for tau neutrinos and antineutrinos
in 1 m length of the lead detector (equivalent to $M_{\rm pb}\simeq$ 35.6 ton)
assuming central scales ($\mu_{R},\mu_{F})=(1.0,1.5)\,m_{T}$ in the
computation of heavy-meson production in $pp$ collisions at $\sqrt{s}$ = 14
TeV and an integrated luminosity ${\cal L}=3000$ fb-1. The scale combinations
$(0.5,1.5)m_{T}$ and $(1,0.75)m_{T}$ give the upper and lower limits of the
scale variation envelope obtained from the seven combinations of
$(N_{R},N_{F})$ shown in figure 6.
| $\nu_{\tau}$ | $\bar{\nu}_{\tau}$ | $\nu_{\tau}+\bar{\nu}_{\tau}$ | $\nu_{\tau}+\bar{\nu}_{\tau}$
---|---|---|---|---
$(\mu_{R},\ \mu_{F})$ | (1, 1) $m_{T}$ | (1, 1) $m_{T}$ | (0.5, 1) $m_{T}$ | (1, 0.5) $m_{T}$
$\langle k_{T}\rangle$ | 0.7 GeV | 0 GeV | 1.4 GeV | 2.2 GeV | 0.7 GeV
$D_{s}$ | 1591 | 774 | 2365 | 2455 | 2143 | 1822 | 7834 | 1779
$B^{\pm,0}$ | 87 | 42 | 129 | 131 | 124 | 115 | 202 | 91
Total | 1678 | 816 | 2494 | 2586 | 2267 | 1937 | 8036 | 1870
Table 2: Same as Table 1, but adopting ($\mu_{R},\mu_{F})=(1.0,1.0)m_{T}$ as
central scales in the computation of heavy-meson hadroproduction (conventional
scale choice). The scale combinations $(0.5,1)m_{T}$ and $(1,0.5)m_{T}$ give
the upper and lower limits of the scale variation envelope from the seven
combinations of $(N_{R},N_{F})$ shown in figure 7.

In tables 1 and 2, we present the total event numbers for our preferred scale
choice $(\mu_{R},\mu_{F})=(1.0,1.5)\,m_{T}$ and the conventional QCD scale
choice $(\mu_{R},\mu_{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 $\nu_{\tau}$ and $\bar{\nu}_{\tau}$ events from
$D_{s}^{\pm}$ and $B$ meson decays, respectively, and the total number, each
for selected ($\mu_{R}$, $\mu_{F}$) scale and $\langle k_{T}\rangle$ choices.
Table 1 shows that when using $(\mu_{R},\mu_{F})=(1.0,1.5)m_{T}$ and the
$\langle k_{T}\rangle$ value that yields the best match to the Powheg + PYTHIA
$D_{s}^{\pm}$ energy distribution in the far forward $\eta>6.87$ region
($\langle k_{T}\rangle=0.7$ GeV), the total number of tau neutrino plus
antineutrino events is predicted to be $\sim 3,800$. Factorization and
renormalization scale variations around the QCD central scales yield a very
broad uncertainty band on the number of events, varying in the interval $\sim
1,800-11,200$. The large scale uncertainties indicate that contributions from
missing higher orders in the perturbative calculation in collinear
factorization are relatively large. On the other hand, the central scale
choice, with $\langle k_{T}\rangle$ variation in the range $\langle
k_{T}\rangle=0-2.2$ GeV, produces a smaller uncertainty in the number of
events, which span the range $\sim 3,900-2,900$, with the lower end of the
range of number of events corresponding to the ($\mu_{R}$, $\mu_{F}$ and
$\langle k_{T}\rangle$)-combination that is favored by LHCb data. Given our
expectation that, for $\langle k_{T}\rangle=2.2$ GeV, some effects of missing
higher-order perturbative corrections are compensated by the large value of
this parameter, an uncertainty in $\langle k_{T}\rangle$ cannot be added 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 $\langle k_{T}\rangle$
from 0 to 2.2 GeV.

We set $(\mu_{R},\mu_{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 $\sim 2,500$ events for $\langle k_{T}\rangle=0.7$
GeV. This number amounts to a factor of $\sim 2/3$ of the one obtained in the
evaluation with $(\mu_{R},\mu_{F})=(1,1.5)\,m_{T}$ for $\langle
k_{T}\rangle=0.7$ GeV, however, it is comparable to the $2,900$ events for
$\langle k_{T}\rangle=2.2$ GeV with $(\mu_{R},\mu_{F})=(1,1.5)\,m_{T}$.

###  4.2 Muon neutrinos

Figure 8: Our predictions for the muon neutrinos plus antineutrinos energy
distribution for $pp$ collisions at $\sqrt{s}=14$ TeV, with neutrino
pseudorapidity $\eta>6.87$. Left: Shown are the contributions from each charm
hadron $D^{+}$, $D^{0}$, $D_{s}^{+}$ and $\Lambda_{c}$ and their
antiparticles, together with their sum. Right: Total contributions from the
charm hadrons and bottom hadrons are presented, respectively. The total for
$B$-hadrons accounts for the contributions from $B^{+}$, $B^{0}$, $B_{s}^{+}$
and $\Lambda_{b}$ and their antiparticles. The upper plots refer to the scale
choice $(\mu_{R},\mu_{F})=(1.0,1.5)\,m_{T}$, while the lower plots correspond
to the conventional choice $(\mu_{R},\mu_{F})=(1.0,1.0)\,m_{T}$. The value
$\left\langle k_{T}\right\rangle=0.7$ GeV is used as input in all cases.

At the interaction region, the muon neutrino and electron neutrino fluxes from
heavy-flavor production and decay will be nearly the same, coming primarily
from the neutral and charged $D$ semileptonic decays. Figure 8 shows the
energy distributions of the sum of muon neutrinos and antineutrinos. The upper
two panels are for $(\mu_{R},\mu_{F})=(1.0,1.5)\,m_{T}$ and the lower two
panels are for $(\mu_{R},\mu_{F})=(1.0,1.0)\,m_{T}$. The left plots show the
contributions from the charm hadrons, $D^{+}$, $D^{0}$, $D_{s}^{+}$ and
$\Lambda_{c}$ and their sum. In the right plots, the total contributions from
all charm hadrons and the bottom hadrons are shown. As can be seen in the left
plots, the decays of $D^{\pm}$, $D^{0}$ and $\bar{D}^{0}$ dominate the muon
neutrino and antineutrino distributions. This is due to larger fragmentation
fractions and decay branching fractions to muon neutrinos compared to those of
$D_{s}^{\pm}$ and $\Lambda_{c}$. Similarly, the muon neutrinos and
antineutrinos from $B$ meson decays are mainly from $B^{\pm}$, $B^{0}$ and
$\bar{B}^{0}$. The bottom hadron contributions, compared to the charm hadron
contributions, to the inclusive muon neutrino plus antineutrino energy
distribution at $\sqrt{s}=14$ TeV with $\eta>6.87$ is about a factor $\sim
1/60$ smaller than the distributions from charm at low energy, and a factor of
$\sim 1/20$ at high energy. As in figure 4 for $\nu_{\tau}+\bar{\nu}_{\tau}$,
figure 8 shows that the predicted energy distribution of
$\nu_{\mu}+\bar{\nu}_{\mu}$ from charm is $\sim 1.5$ larger for
$(\mu_{R},\mu_{F})=(1.0,1.5)\,m_{T}$ than for
$(\mu_{R},\mu_{F})=(1.0,1.0)\,m_{T}$, while the smaller contributions from $B$
hadrons are much less sensitive to the scale choice.

Figure 9: Our predictions for the muon neutrino and antineutrino number of
charged current events per GeV for 1 ton of lead target as a function of the
incident neutrino energy for neutrinos with a pseudorapidity $\eta>6.87$
generated by heavy-flavor decays in $pp$ collisions with $\sqrt{s}=14$ TeV. A
value of $\left\langle k_{T}\right\rangle$ = 0.7 GeV is adopted for the
central predictions, while the uncertainty bands accounts for the effect of
$\left\langle k_{T}\right\rangle$ variation in the interval $0-1.4$ GeV. NLO
QCD corrections are accounted for in the DIS cross section. The integrated
luminosity amounts to ${\cal L}=3000$ fb-1.

The corresponding predictions for the number of muon neutrino and antineutrino
charged-current events per unit energy for a 1 ton lead target for $\eta>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 $\langle k_{T}\rangle$ variation in the range [0, 1.4]
GeV for $(\mu_{R},\mu_{F})=(1.0,1.5)\,m_{T}$. Again, the dashed histograms
correspond to results with $\langle k_{T}\rangle=2.2$ GeV. The left panels in
figure 10 show, from top to bottom, the sum of $\nu_{\mu}+\bar{\nu}_{\mu}$
charged current events per ton of lead, for $\nu_{\mu}$ and for
$\bar{\nu}_{\mu}$, all for $(\mu_{R},\mu_{F})=(1.0,1.5)\,m_{T}$. The right
panels show the same, but with $(\mu_{R},\mu_{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
perturbative QCD gives a factor of $\sim 13$ in the ratio of
$\nu_{\mu}+\bar{\nu}_{\mu}$ to $\nu_{\tau}+\bar{\nu}_{\tau}$ events based on
heavy flavor alone, with our default input parameters.

Figure 10: Uncertainty range due to the QCD scale variation in the muon
neutrino and antineutrino number of charged current events per GeV per ton as
a function of the incident neutrino energy for neutrinos with pseudorapidity
$\eta>6.87$ generated by heavy-flavor decays in $pp$ collisions with
$\sqrt{s}=14$ TeV and integrated luminosity ${\cal L}=3000$ fb-1. The central
predictions are obtained using as input $(\mu_{R},\mu_{F})=(1,1.5)\,m_{T}$
(left) and $(\mu_{R},\mu_{F})=(1.0,1.0)\,m_{T}$ (right), respectively. The
upper and lower limits in each panel arise from 7-point scale variation in the
same range as in figure 1.

Heavy-flavor hadron decays, however, are not the only sources of
$\nu_{\mu}+\bar{\nu}_{\mu}$. In principle, pion, kaon and weak gauge boson
decays can also contribute. In ref. Beni:2019gxv , the contributions from $W$
and $Z$ production and decay to neutrinos are studied. They show that
neutrinos from weak gauge boson decays populate mostly in the pseudorapidity
range of $|\eta|<4.5$, but are negligible in the region $\eta>6.5$
Beni:2019gxv .

At first glance, one might expect that pion and kaon decays will not be
important sources of neutrinos. Pions with $E_{\pi}>9$ GeV have $\gamma
c\tau>500$ m. Charged kaons have $\gamma c\tau>500$ m for $E_{K}>67$ GeV. Pion
and kaon contributions to the number of $\nu_{\mu}+\bar{\nu}_{\mu}$ at high
energy were neglected in the earlier work of ref. DeRujula:1992sn . Park, in
ref. Park:2011gh , 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 $\sim 100$ times
more $\nu_{\mu}+\bar{\nu}_{\mu}$ events than $\nu_{\tau}+\bar{\nu}_{\tau}$
events.

In ref. Abreu:2019yak , a more detailed evaluation of
$\nu_{\mu}+\bar{\nu}_{\mu}$ events is performed for a detector with $25\times
25$ cm2 cross-sectional area at 480 meters from the Atlas interaction point.
They find a factor of $\sim 1000$ more interactions by
$\nu_{\mu}+\bar{\nu}_{\mu}$ than by $\nu_{\tau}+\bar{\nu}_{\tau}$ when pions
and kaons are included, although with a different energy distribution. Their
evaluation of heavy-flavor contributions is done using Pythia
Sjostrand:2019zhc , while their light hadron production is estimated using the
Crmc crmc  simulation package with Epos-LHC Pierog:2013ria , Qgsjet-II-04
Ostapchenko:2010vb  and Sibyll 2.3c Engel:2019dsg ; Fedynitch:2018cbl ;
Riehn:2017mfm . 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 $\sim 55$ m from the Atlas interaction point. They find that magnetic
fields sweep lower-energy (below $\sim 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 radius 17 mm
Abreu:2019yak . They also find that two-body decays of charged pions and
charged kaons are the dominant sources of $\nu_{\mu}+\bar{\nu}_{\mu}$
production.

Our primary focus here is on the heavy-flavor contributions to the number of
events from $\nu_{\mu}+\bar{\nu}_{\mu}$ and $\nu_{\tau}+\bar{\nu}_{\tau}$.
Instead of simulating light charged hadron trajectories in magnetic fields,
guided by the results of ref. Abreu:2019yak , we approximately evaluate the
number of $\nu_{\mu}+\bar{\nu}_{\mu}$ that come from pions and kaons as
follows. We evaluate the $\pi^{\pm}$ and $K^{\pm}$ two-body decay
contributions to the flux of $\nu_{\mu}+\bar{\nu}_{\mu}$ 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 $\theta<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 $\nu_{\mu}+\bar{\nu}_{\mu}$ would include $K_{L}$
semileptonic decays.

We use the parametrization of Koers et al. Koers:2006dd , 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
$\sqrt{s}=14$ TeV are $\sim 50$ and $\sim 6$, respectively Koers:2006dd . Two-
body decays are implemented, as for the $D_{s}^{\pm}\to\nu_{\tau}\tau$ case,
with the requisite changes to the initial and final particle masses. Compared
to figure 8, charged pions and kaons contribute a factor of $\sim 100$ more
$\nu_{\mu}+\bar{\nu}_{\mu}$ 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
$\nu_{e}+\bar{\nu}_{e}$, we consider the dominant contribution which is from
$K_{L}\to\pi e\nu_{e}$. With the same geometry requirements, we show with the
green curve in figure 11 the contribution from $K_{L}$ production and decay
into $\nu_{e}+\bar{\nu}_{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. Cirigliano:2001mk . The peak of the
electron neutrino distribution from $K_{L}$ decays is about a factor of $\sim
2$ larger than the peak of the electron neutrino distribution from heavy-
flavor decays. While $K_{L}$ semileptonic decays to $\nu_{e}$ are dominant, a
more complete calculation would include the $K^{+}\to\nu_{e}$ semileptonic
decays as well.

Figure 11: The $\nu_{\mu}+\bar{\nu}_{\mu}$ energy distributions from charged
pions (blue) and kaons (red) and the $\nu_{e}+\bar{\nu}_{e}$ energy
distributions from $K_{L}$ semileptonic decays (green) for forward neutrinos
($\eta>6.87$) from the decay of mesons produced in $pp$ interactions at
$\sqrt{s}=14$ TeV, considering those mesons whose decay occurs within 55 m of
the interaction point and whose momentum lies within an angle of $\theta<1$
mrad from the beam axis. The meson energy and rapidity distributions are
evaluated using the parametrization of Koers et al. Koers:2006dd .

Three-flavor neutrino oscillations of the much larger number of
$\nu_{\mu}+\bar{\nu}_{\mu}$ from charged pions and kaons to
$\nu_{\tau}+\bar{\nu}_{\tau}$ could, in principle, overwhelm the number of
$\nu_{\tau}+\bar{\nu}_{\tau}$ from heavy-flavor decays. However, the baseline
of 480 m is not at all optimal for $\nu_{\mu}\to\nu_{\tau}$ oscillations in
the standard scenario with 3 active flavors. In two flavor approximation, for
a mass squared difference of $\Delta m_{32}^{2}\simeq 2.5\times 10^{-3}$ eV2
and a mixing angle $\theta_{23}\simeq\pi/4$, the oscillation probability is
$P(\nu_{\mu}\to\nu_{\tau})\simeq 2.3\times 10^{-6}/(E^{2}/{\rm GeV^{2}}$), so
for energies above $E_{\nu_{\mu}}=10$ GeV, $\nu_{\mu}\to\nu_{\tau}$
oscillations give a negligible contribution to the number of
$\nu_{\tau}+\bar{\nu}_{\tau}$, even given the large theoretical uncertainties
in the number of $\nu_{\tau}+\bar{\nu}_{\tau}$ events per unit energy. The
$\nu_{e}\to\nu_{\tau}$ 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 $\nu_{\tau}+\bar{\nu}_{\tau}$ 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 decay to neutrinos, using all
available data on forward production (e.g., from LHCf Adriani:2015iwv , TOTEM
Antchev:2014lez  and CMS Sirunyan:2017nsj ) along with magnet and detector
configurations in the interaction region, will be necessary. A forward tune of
Pythia that is underway Abreu:2020ddv  would also guide future work.

##  5 New Physics

The detection of a large number of identifiable $\nu_{\tau}+\bar{\nu}_{\tau}$
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 $\nu_{l}$ can be expressed as a superposition of the
mass eigenstates $\nu_{j}$ according to the formula

$\nu_{l}=\sum_{j=1}^{n_{\nu}}U_{lj}\nu_{j}\ ,$ (7)

where $n_{\nu}$ is the total number of neutrinos subject to oscillations: for
the standard oscillation scenario with three active neutrinos, $n_{\nu}=3$,
while for an oscillation scenario with three active neutrinos plus one sterile
neutrino, $n_{\nu}=4$. The full transition probability appears in, for
example, ref. pdg:2016 . The blue lines in the two panels of figure 12 show
that tau neutrino survival probability $P(\nu_{\tau}\to\nu_{\tau})$ in the
three-flavor scenario approaches unity in the energy range
$E_{\nu}\in[1,1000]$ GeV. The values of the parameters used here for three-
flavor mixing are $\sin^{2}(\theta_{12})=0.310$,
$\sin^{2}(\theta_{13})=0.02241$, $\sin^{2}(\theta_{23})=0.580$, $\Delta
m_{21}^{2}=7.5\times 10^{-5}$ eV2 and $\Delta m_{31}^{2}=2.457\times 10^{-3}$
eV2 Esteban:2018azc . The values of all CP phases are taken to be zero. For a
baseline of 480 m, the tau neutrino survival probability has significant
features in the few MeV energy range and below, but not in the energy range of
interest (see the position of the peaks of the energy distributions in section
4).

Figure 12: The survival probability $P(\nu_{\tau}\to\nu_{\tau})$ as a function
of the tau neutrino energy $E_{\nu}$, for the standard 3-active-flavor
oscillation framework pdg:2016  and in a $3+1$ oscillation scenario,
considering a baseline of $L=480$ m. Besides results in the standard
framework, the upper panel shows results for the 3+1 scenario with $m_{4}=20$
eV, for $|U_{\tau 4}|^{2}=0.08$ and 0.15, whereas the lower panel shows
results for the 3+1 scenario with $|U_{\tau 4}|^{2}=0.15$ and $m_{4}=1$ eV, 20
eV and 1 keV.

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}\gg m_{1,2,3}$, can be simplified to Giunti:2019aiy

$P(\nu_{\alpha}\rightarrow\nu_{\beta})\simeq\delta_{\alpha\beta}-4(\delta_{\alpha\beta}-\mid
U_{\beta n_{\nu}}\mid^{2})\mid U_{\alpha
n_{\nu}}\mid^{2}\sin^{2}\Biggl{(}\frac{\Delta m^{2}L}{4E_{\nu}}\Biggr{)}\ ,$
(8)

where $\Delta m^{2}=m_{4}^{2}-(m_{1}^{2}+m_{2}^{2}+m_{3}^{2})/3\simeq
m_{4}^{2}$. For $m_{4}\sim 20$ eV, oscillation effects in the tau neutrino
survival probability can be pronounced. The highest energy oscillation node
with this hierarchy occurs at

$E_{\nu\textrm{-max [GeV]}}=\frac{\Delta m^{2}L}{2\pi}=0.807\,\Delta
m_{\textrm{[eV]}}^{2}L_{\textrm{[km]}}\ .$ (9)

Introducing a heavy sterile neutrino extends the region of pronounced
oscillation dips in the tau neutrino survival probability to a higher energy
with examples shown in figure 12. When $m_{4}=20$ eV, a survival probability
dip occurs at $E_{\nu_{\tau}}=155$ GeV, an energy near that of the peak of the
unoscillated tau neutrino number of events per unit energy. Using eq. (8) and
$|U_{\tau 4}|^{2}=0.08$ and 0.15, the $\nu_{\tau}\to\nu_{\tau}$ survival
probability is shown for $m_{4}=20$ eV (upper panel) and $m_{4}=1$ eV, 20 eV
and 1 keV (lower panel), all for a baseline of $L=480$ m. The values of
$|U_{\tau 4}|^{2}$ are acceptable according to current IceCube constraints
Jones:2019nix ; Blennow:2018hto . Matter effects can be neglected with the
density of the Earth’s crust $\rho\approx 2.6$ g/cm3 and the electron fraction
$Y_{e}\approx 0.5$ Ohlsson:1999xb ; Li:2018ezt .

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 $\nu_{\tau}$ survival probability.
Given the uncertainties in the absolute scale of the
$\nu_{\tau}+\bar{\nu}_{\tau}$ 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
three-active-flavor oscillation framework (black histogram) and in the 3+1
oscillation framework with $m_{4}=20$ eV and $|U_{\tau 4}|^{2}=0.15$,
considering our default heavy-flavor QCD input parameter set and $\langle
k_{T}\rangle=0.7$ GeV (left) and 2.2 GeV (right). The orange-dashed histogram
shows the effect of $\nu_{\tau}$ disappearance due to oscillations in a 3+1
scenario with $|U_{e4}|^{2}=|U_{\mu 4}|^{2}=0$, where a dip is visible at
neutrino energies $\sim$ 150 GeV.

Figure 13: Predictions of the number of $\nu_{\tau}+\overline{\nu}_{\tau}$
charged-current events as a function of neutrino energy in absence of
oscillations and in presence of oscillations in a 3+1 mixing framework, for
various choices of the oscillation parameters in the same experimental setup
already used in section 4. Numbers of events are reported for a ton of lead
detector, for $\langle k_{T}\rangle$ = 0.7 GeV (left) and 2.2 GeV (right).

Figure 14: The ratio of the number of $\nu_{\tau}+\overline{\nu}_{\tau}$
charged-current interaction events and the number of events at 1 TeV in each
case, as a function of neutrino energy, without oscillations (upper left plot)
and in a 3+1 oscillation scenarios with different mixing parameters. The
numbers of events are evaluated for $\langle k_{T}\rangle$ = 0.7 GeV for three
different ($\mu_{R}$, $\mu_{F}$) scale choices that bracket the scale
uncertainties.

Figure 15: The same as figure 14, but considering different values of $\langle
k_{T}\rangle$, for our default ($\mu_{R}$, $\mu_{F}$) scale choice.

The oscillation dip may be partially filled in by $\nu_{\mu}\to\nu_{\tau}$ and
$\nu_{e}\to\nu_{\tau}$ oscillations, where the $\nu_{\mu}$ and $\nu_{e}$ come
from heavy-flavor decays and from $\pi^{\pm},\ K^{\pm}$ and $K_{L}$ decays.
The NOMAD experiment set the most stringent limits on effective mixing angles
with $\nu_{\tau}$ for 3+1 scenarios with $\Delta m^{2}\simeq
m_{4}^{2}\mathrel{\hbox{\hbox to0.0pt{\lower 3.22916pt\hbox{$\sim$}\hss}
\kern-3.00003pt\raise 1.72218pt\hbox{$>$}}}30$ eV2 Astier:2001yj ;
Astier:2003gs ,

$\displaystyle\sin^{2}2\theta_{e\tau}$ $\displaystyle\simeq$ $\displaystyle
4|U_{e4}|^{2}|U_{\tau 4}|^{2}<1.5\times 10^{-2}$ (10)
$\displaystyle\sin^{2}2\theta_{\mu\tau}$ $\displaystyle\simeq$ $\displaystyle
4|U_{\mu 4}|^{2}|U_{\tau 4}|^{2}<3.3\times 10^{-4}$ (11)
$\displaystyle\sin^{2}2\theta_{e\mu}$ $\displaystyle\simeq$ $\displaystyle
4|U_{e4}|^{2}|U_{\mu 4}|^{2}<1.4\times 10^{-3}\ .$ (12)

We use $|U_{\tau 4}|^{2}=0.15$ to illustrate sterile neutrino mixing
effects.555 This is consistent with the 99%CL limit from ref. Dentler:2018sju
, where results are shown only for $\Delta m_{41}^{2}\mathrel{\hbox{\hbox
to0.0pt{\lower 3.22916pt\hbox{$\sim$}\hss} \kern-3.00003pt\raise
1.72218pt\hbox{$<$}}}10$ eV2. For $|U_{\mu 4}|^{2}$ and $|U_{e4}|^{2}$, eqs.
(10-12) must be satisfied. For illustration purposes, we take $|U_{\mu
4}|^{2}=5\times 10^{-4}$, a value close to the maximum mixing consistent with
eq. (11). The Troisk tritium beta decay data can be used to set limits on
$|U_{e4}|^{2}$ for sterile neutrino masses in the range of $1-100$ eV
Belesev:2012hx ; Belesev:2013cba . For $m_{4}=20$ eV,
$|U_{e4}|^{2}\mathrel{\hbox{\hbox to0.0pt{\lower 3.22916pt\hbox{$\sim$}\hss}
\kern-3.00003pt\raise 1.72218pt\hbox{$<$}}}6.7\times 10^{-3}$ comes from the
Troisk upper bound. The other constraints, eqs. (10) and (12), are satisfied
by this value of $|U_{e4}|^{2}$, considering our aforementioned choice for
$|U_{\tau 4}|^{2}$. We use parameters $|U_{\mu 4}|^{2}=5\times 10^{-4}$ and a
slightly lower value of $|U_{e4}|^{2}=6\times 10^{-3}$, still satisfying the
constraints discussed above, in the blue histogram in figure 13. In the
figure, $\nu_{e}$ and $\nu_{\mu}$ (and anti-neutrinos) from heavy-flavor
decays and from $\pi^{\pm},\ K^{\pm}$ and $K_{L}$ decays are included.

The spectral shape of the number of $\nu_{\tau}+\bar{\nu}_{\tau}$ events
changes with different choices of mixing parameters. If $|U_{e4}|^{2}=|U_{\mu
4}|^{2}=0$, spectral distortions will be conspicuous, even if $\langle
k_{T}\rangle=2.2$ GeV (corresponding to a total smaller number of events than
our default case $\langle k_{T}\rangle=0.7$ GeV), as can be seen in the right
panel of figure 13. Keeping the same value of $|U_{e4}^{2}|$ and increasing
the $|U_{\mu 4}|^{2}$ value to $|U_{\mu 4}|^{2}<10^{-3}$, the
$\nu_{\tau}\to\nu_{\tau}$ 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_{\mu 4}|^{2}=5\times
10^{-4}$. When $|U_{e4}|^{2}\neq 0$, electron neutrinos coming primarily from
$K^{0}_{e3}$ that oscillate to $\nu_{\tau}$ because of sterile neutrino mixing
additionally fill in the spectral distortion of the
$\nu_{\tau}+\bar{\nu}_{\tau}$ event number as shown in figure 13 with the
solid blue histogram for $|U_{e4}|^{2}=6\times 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
$\langle k_{T}\rangle$ smearing parameter, will be difficult if the mixing
parameters $|U_{e4}|^{2}$ and $|U_{\mu 4}|^{2}$ are close to $6\times 10^{-3}$
and $5\times 10^{-4}$, respectively. Although the the values of $|U_{e4}|^{2}$
and $|U_{\mu 4}|^{2}$ are small, the large number of $\nu_{e}$ and $\nu_{\mu}$
from light-meson production and decay means that $\nu_{e}$ and $\nu_{\mu}$
mixing with sterile neutrinos can have an non-negligible impact. The challenge
is illustrated in figs. 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 $\langle k_{T}\rangle$ on the shape of the distribution. In figure
14, $\langle k_{T}\rangle=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_{\ell 4}|^{2}\neq 0$ in the upper right panel presents a slightly
shallower dip in the oscillated $3+1$ spectrum as compared to the case with
$|U_{e4}|^{2}=0$ and $|U_{\mu 4}|^{2}=0$ shown in the lower right panel. The
normalization of each distribution in figure 14 to its own distribution at
$E=1$ TeV collapses the scale uncertainty band of figure 6. Figure 15 shows
that a large transverse momentum smearing can somewhat obscure the spectral
features of $3+1$ oscillations as well. The effects of $\nu_{e}$ and
$\nu_{\mu}$ mixing with sterile neutrinos are more important in each panel of
figures 14 and 15, namely, for $(\mu_{R},\,\mu_{F})=(1.0,\,0.75)m_{T}$ in
figure 14 and $\langle k_{T}\rangle$ = 2.2 GeV in figure 15. The dip in the
spectrum is least pronounced in the case of
$(\mu_{R},\,\mu_{F})=(1.0,\,0.75)m_{T}$ and $\langle k_{T}\rangle$ = 2.2 GeV.
A dedicated study of $3+1$ oscillation scenarios that includes detector
resolution effects would be important to understand in detail the reach of a
forward tau neutrino experiment along the LHC beamline.

##  6 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 $\nu_{\tau}+\bar{\nu}_{\tau}$ 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 $\eta>6.87$,
consistently with the geometry of a 1 m radius cylindrical detector at a
distance of 480 m from the interaction point Feng:2017uoz ; Ariga:2018pin ;
Ariga:2019ufm  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. We show results for both
$(\mu_{R},\mu_{F})=(1.0,1.5)~{}m_{T}$ and
$(\mu_{R},\mu_{F})=(1.0,1.0)~{}m_{T}$, to better match the LHCb data and for
comparisons with the usual scale conventions adopted in QCD studies of heavy-
flavor production, respectively. Seven-point scale variation around the
central scale combinations $(\mu_{R},\mu_{F})=(1.0,1.0)~{}m_{T}$ and
$(1.0,1.5)~{}m_{T}$ causes differences amounting to a factor of $\sim 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 $\langle k_{T}\rangle$ in the Gaussian
smearing factor $f(\vec{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 Bai:2018xum . 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(\vec{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 $\sqrt{s}=13$ TeV, for $p_{T}$ $\in$
[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 heavy-flavor physics in the very forward region (e.g.,
ref. Aoki:2019jry ) 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. Bacchetta:2018ivt ; DAlesio:2019qpk ). In another arena, predictions of
the prompt atmospheric neutrino flux Bhattacharya:2015jpa ;
Bhattacharya:2016jce ; Gauld:2015kvh ; Garzelli:2016xmx ; Benzke:2017yjn ;
Zenaiev:2019ktw  will also be constrained by measurements of the tau neutrino
energy spectrum in the far-forward region at the LHC.

For a baseline of $\sim$ 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 ${\cal 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 $\nu_{\tau}$ events due to
$\nu_{\mu}\to\nu_{\tau}$ oscillations is not important due to the small value
of $|U_{\mu 4}|^{2}$. This is the case even though there are $\sim 100$ times
more $\nu_{\mu}+\bar{\nu}_{\mu}$ produced from charged pion and kaon decays
than from heavy-flavor decays. The number of $\nu_{\mu}+\bar{\nu}_{\mu}$ from
heavy flavor decays is larger than the number of $\nu_{\tau}+\bar{\nu}_{\tau}$
by a factor of $\sim$ 10\. The $\nu_{e}\to\nu_{\tau}$ oscillation effects are
also small.

Therefore a large number of $\nu_{\tau}+\bar{\nu}_{\tau}$ 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
$\nu_{\tau}$ 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_{\tau 4}|^{2}$ is currently poorly constrained. We showed that, as long as
$\nu_{e}\to\nu_{\tau}$ appearance is suppressed, for $|U_{\tau 4}|^{2}=0.15$
in principle the oscillation effect would be unambiguous. However,
uncertainties in heavy-flavor 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.

##  Appendix A Decay Distributions

The two-body decays of the $D_{s}$ in the $D_{s}$ rest frame come from energy
and momentum conservation. Keeping the polarization of the $\tau$, the energy
and angular distribution of the tau neutrino can also be obtained Barr:1988rb
; Barr:1989 ; Gaisser:1990vg . Detailed formulas for the direct $\nu_{\tau}$
and chain decay $\nu_{\tau}$ distributions appear in Appendix B of ref.
Bai:2018xum . Decays in the $D_{s}$ and $\tau$ 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)\to K(D)l\nu_{l}\,(l=e,\mu)$. The distribution of
neutrinos from the decay in the rest frame is respectively given by

$\displaystyle\frac{{\rm d}\Gamma(h_{i}\to\nu_{l})}{{\rm d}x_{\nu}}$
$\displaystyle\sim$
$\displaystyle\begin{dcases}\frac{m_{D}^{5}\,x_{\nu}^{2}(1-r_{h}-x_{\nu})^{2}[3+r_{h}(3-x_{\nu})-5x_{\nu}+2x_{\nu}^{2}]}{(1-x_{\nu})^{3}}\quad{\rm
for}\quad h_{i}=D\\\
\frac{m_{B}^{5}\,x_{\nu}^{2}(1-r_{h}-x_{\nu})^{2}}{(1-x_{\nu})}\hskip
135.0002pt{\rm for}\quad h_{i}=B\end{dcases}$ (13)

where the fraction of energy transferred to the neutrino is
$x_{\nu}=2E_{\nu}/m_{i}$ and the hadron mass fraction is
$r_{h}=m_{f}^{2}/m_{i}^{2}$, with $m_{i}$ and $m_{f}$ being the hadron masses
in the initial and final states, respectively. The corresponding cumulative
distribution functions (CDF) are

$\displaystyle F(x_{\nu})$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma}\int_{0}^{x_{\nu}}dx^{\prime}\frac{\textrm{d}\Gamma}{\textrm{d}x^{\prime}}$
$\displaystyle=$
$\displaystyle\begin{dcases}\frac{1}{D(r_{h})(1-x_{\nu})^{2}}\biggl{[}x_{\nu}\Bigl{(}2r_{h}^{3}x_{\nu}^{2}-6\,r_{h}^{2}(2-x_{\nu})(1-x_{\nu})+(2-x_{\nu})(1-x_{\nu})^{2}x_{\nu}^{2}\\\
\hskip
87.00009pt-2r_{h}(1-x_{\nu})^{2}x_{\nu}^{2}\Bigr{)}-12r_{h}^{2}(1-x_{\nu})^{2}\ln[1-x_{\nu}]\biggr{]}\quad{\rm
for}\quad h_{i}=D\\\
\frac{1}{D(r_{h})}\biggl{[}x_{\nu}\Bigl{(}x_{\nu}^{2}(4-3x_{\nu})-8r_{h}x_{\nu}^{2}-6r_{h}^{2}(2+x_{\nu})\Bigr{)}-12r_{h}^{2}\ln[1-x_{\nu}]\biggr{]}\hskip
11.99998pt{\rm for}\quad h_{i}=B\end{dcases}$

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\to D\tau\nu_{\tau}$, and the
$\tau$ lepton decays subsequently to $\nu_{\tau}$. In this case, there is an
additional effect due to the non-negligible mass of $\tau$. With the
additional mass term $r_{\tau}=m_{\tau}^{2}/m_{B}^{2}$, one obtains the
distributions of $\tau$ and $\nu_{\tau}$ from $B$ decays:

$\displaystyle\frac{{\rm d}\Gamma(B\to\tau)}{{\rm d}x_{\tau}}$
$\displaystyle\sim$
$\displaystyle\frac{m_{B}^{5}}{(1+r_{\tau}-x_{\tau})^{3}}(1-x_{\tau}+r_{\tau}-r_{h})^{2}(x_{\tau}^{2}-4r_{\tau})^{1/2}\biggl{[}r_{\tau}^{2}(3x_{\tau}-4)$
$\displaystyle+~{}x\Bigl{(}3+2x_{\tau}^{2}-5x_{\tau}+r_{h}(3-x_{\tau})\Bigr{)}-r_{\tau}\Bigl{(}4+5x_{\tau}^{2}-10x_{\tau}+r_{h}(8-3x_{\tau})\Bigr{)}\biggr{]},$
$\displaystyle\frac{{\rm d}\Gamma(B\to\nu_{\tau})}{{\rm d}x_{\nu}}$
$\displaystyle\sim$
$\displaystyle\frac{m_{B}^{5}}{1-x_{\nu}}\,x_{\nu}^{2}(1-r_{h}-r_{\tau}-x_{\nu})\Bigl{(}r_{h}^{2}-2r_{h}(1+r_{\tau}-x_{\nu})+(1-r_{\tau}-x_{\nu})^{2}\Bigr{)}^{1/2}.$

Here, $x_{\tau}=2E_{\tau}/m_{i}$ is the fraction of energy transferred to tau.
The expression of the CDF for $\tau$ and $\nu_{\tau}$ 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 $\nu_{\mu}$, $\nu_{\tau}$ and $\tau$ 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:

$\displaystyle 0<$ $\displaystyle x_{\nu_{l}}$
$\displaystyle<1-(\sqrt{r_{l}}+\sqrt{r_{h}})^{2}$ $\displaystyle
2\sqrt{r_{\tau}}<$ $\displaystyle x_{\tau}$ $\displaystyle<1+r_{\tau}-r_{h}\
.$ (16)

For electron neutrino ($\nu_{e}$) and muon neutrino ($\nu_{\mu}$), the maximum
value of $x_{\nu}$ is simply $x_{\nu}^{\rm max}=1-r_{h}$.

Figure 16: Left panel: The fractional energy distributions $x_{i}=E_{i}/(m/2)$
in three-body semi-leptonic decays in the rest frame of the decaying heavy
particle with $m=m_{D_{s}}$ and $m=m_{B}$ and $i=\nu_{\mu},\ \nu_{\tau}$ and
$\tau$. Right panel: The cumulative distribution function for each of the
three-body decays shown in the left panel.

###### Acknowledgements.

This work is supported in part by Department of Energy grants DE-SC-0010113,
DE-SC-0012704, and the Korean Research Foundation (KRF) through the CERN-Korea
Fellowship program. The authors would like to express a special thanks to the
Mainz Institute for Theoretical Physics (MITP) of the Cluster of Excellence
PRISMA+ (Project ID 39083149) for its hospitality and support. We thank the
Rencontres du Vietnam/VietNus 2017 workshop at the International Center for
Interdisciplinary Science and Education (ICISE), QuiNhon, Vietnam for support
and discussions in the early stage of this work, and we thank C. Giunti and
M.L. Mangano for comments and discussions. M.V.G. is grateful to the II
Institute for Theoretical Physics of the University of Hamburg for hospitality
during the completion of this work.

## References

  * (1) Particle Data Group collaboration, C. Patrignani et al., _Review of Particle Physics_ , _Chin. Phys._ C40 (2016) 100001.
  * (2) DUNE collaboration, R. Acciarri et al., _Long-Baseline Neutrino Facility (LBNF) and Deep Underground Neutrino Experiment (DUNE)_ , 1512.06148.
  * (3) DUNE collaboration, R. Acciarri et al., _Long-Baseline Neutrino Facility (LBNF) and Deep Underground Neutrino Experiment (DUNE)_ , 1601.05471.
  * (4) DUNE collaboration, B. Abi et al., _The DUNE Far Detector Interim Design Report Volume 1: Physics, Technology and Strategies_ , 1807.10334.
  * (5) H. Park, _The estimation of neutrino fluxes produced by proton-proton collisions at $\sqrt{s}=14$ TeV of the LHC_, _JHEP_ 10 (2011) 092, [1110.1971].
  * (6) J. L. Feng, I. Galon, F. Kling and S. Trojanowski, _ForwArd Search ExpeRiment at the LHC_ , _Phys. Rev._ D97 (2018) 035001, [1708.09389].
  * (7) FASER collaboration, A. Ariga et al., _Technical Proposal for FASER: ForwArd Search ExpeRiment at the LHC_ , 1812.09139.
  * (8) FASER collaboration, A. Ariga et al., _FASER: ForwArd Search ExpeRiment at the LHC_ , 1901.04468.
  * (9) S. Buontempo, G. M. Dallavalle, G. De Lellis, D. Lazic and F. L. Navarria, _CMS-XSEN: LHC Neutrinos at CMS. Experiment Feasibility Study_ , 1804.04413.
  * (10) N. Beni et al., _Physics Potential of an Experiment using LHC Neutrinos_ , _J. Phys._ G46 (2019) 115008, [1903.06564].
  * (11) XSEN collaboration, N. Beni et al., _XSEN: a $\nu$N Cross Section Measurement using High Energy Neutrinos from pp collisions at the LHC_, 1910.11340.
  * (12) A. De Rujula and R. Ruckl, _Neutrino and muon physics in the collider mode of future accelerators_ , in _ECFA-CERN Workshop on large hadron collider in the LEP tunnel, Lausanne and CERN, Geneva, Switzerland, 21-27 Mar 1984: Proceedings. 1._ , pp. 571–596, 1984.  DOI.
  * (13) K. Winter, _Detection of the tau-neutrino at the LHC_ , in _ECFA Large Hadron Collider Workshop, Aachen, Germany, 4-9 Oct 1990: Proceedings.2._ , pp. 37–49, 1990.
  * (14) A. De Rujula, E. Fernandez and J. J. Gomez-Cadenas, _Neutrino fluxes at future hadron colliders_ , _Nucl. Phys._ B405 (1993) 80–108.
  * (15) F. Vannucci, _Neutrino physics at LHC / SSC_ , in _5th International Symposium on Neutrino Telescopes Venice, Italy, March 2-4, 1993_ , pp. 57–68, 1993.
  * (16) LHCb collaboration, R. Aaij et al., _Measurement of $D_{s}^{\pm}$ production asymmetry in $pp$ collisions at $\sqrt{s}=7$ and 8 TeV_, _JHEP_ 08 (2018) 008, [1805.09869].
  * (17) T. Sjöstrand, _The PYTHIA Event Generator: Past, Present and Future_ , _Comput. Phys. Commun._ 246 (2020) 106910, [1907.09874].
  * (18) FASER collaboration, H. Abreu et al., _Detecting and Studying High-Energy Collider Neutrinos with FASER at the LHC_ , 1908.02310.
  * (19) FASER collaboration, H. Abreu et al., _Technical Proposal: FASERnu_ , 2001.03073.
  * (20) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An Introduction to PYTHIA 8.2_ , _Comput. Phys. Commun._ 191 (2015) 159–177, [1410.3012].
  * (21) SHiP collaboration, C. Ahdida et al., _SND@LHC_ , 2002.08722.
  * (22) DONUT collaboration, K. Kodama et al., _Observation of tau neutrino interactions_ , _Phys. Lett._ B504 (2001) 218–224, [hep-ex/0012035].
  * (23) DONuT collaboration, K. Kodama et al., _Final tau-neutrino results from the DONuT experiment_ , _Phys. Rev._ D78 (2008) 052002, [0711.0728].
  * (24) OPERA collaboration, F. Pupilli, _Recent results of the OPERA experiment_ , _AIP Conf. Proc._ 1743 (2016) 060004.
  * (25) OPERA collaboration, N. Agafonova et al., _Final Results of the OPERA Experiment on $\nu_{\tau}$ Appearance in the CNGS Neutrino Beam_, _Phys. Rev. Lett._ 120 (2018) 211801, [1804.04912].
  * (26) Super-Kamiokande collaboration, Z. Li et al., _Measurement of the tau neutrino cross section in atmospheric neutrino oscillations with Super-Kamiokande_ , _Phys. Rev._ D98 (2018) 052006, [1711.09436].
  * (27) IceCube collaboration, M. G. Aartsen et al., _Measurement of Atmospheric Tau Neutrino Appearance with IceCube DeepCore_ , _Phys. Rev._ D99 (2019) 032007, [1901.05366].
  * (28) IceCube collaboration, M. G. Aartsen et al., _Measurement of the multi-TeV neutrino cross section with IceCube using Earth absorption_ , _Nature_ 551 (2017) 596–600, [1711.08119].
  * (29) M. Bustamante and A. Connolly, _Extracting the Energy-Dependent Neutrino-Nucleon Cross Section above 10 TeV Using IceCube Showers_ , _Phys. Rev. Lett._ 122 (2019) 041101, [1711.11043].
  * (30) P. Nason, S. Dawson and R. K. Ellis, _The Total Cross-Section for the Production of Heavy Quarks in Hadronic Collisions_ , _Nucl. Phys._ B303 (1988) 607–633.
  * (31) P. Nason, S. Dawson and R. K. Ellis, _The One Particle Inclusive Differential Cross-Section for Heavy Quark Production in Hadronic Collisions_ , _Nucl. Phys._ B327 (1989) 49–92.
  * (32) M. L. Mangano, P. Nason and G. Ridolfi, _Heavy quark correlations in hadron collisions at next-to-leading order_ , _Nucl. Phys._ B373 (1992) 295–345.
  * (33) LHCb collaboration, R. Aaij et al., _Measurements of prompt charm production cross-sections in $pp$ collisions at $\sqrt{s}=13$ TeV_, _JHEP_ 03 (2016) 159, [1510.01707].
  * (34) H. B. J. Koers, A. Pe’er and R. A. M. J. Wijers, _Parameterization of the energy and rapidity distributions of secondary pions and kaons produced in energetic proton-proton collisions_ , _Submitted to: Phys. Rev. D_ (2006) , [hep-ph/0611219].
  * (35) O. S. Bruning, P. Collier, P. Lebrun, S. Myers, R. Ostojic, J. Poole et al., _LHC Design Report Vol.1: The LHC Main Ring_ , _CERN-2004-003-V1_ (2004) .
  * (36) L. Apanasevich et al., _$k_{T}$ effects in direct photon production_, _Phys. Rev._ D59 (1999) 074007, [hep-ph/9808467].
  * (37) M. L. Mangano, _Two lectures on heavy quark production in hadronic collisions_ , _Proc. Int. Sch. Phys. Fermi_ 137 (1998) 95–137, [hep-ph/9711337].
  * (38) G. Miu and T. Sjostrand, _$W$ production in an improved parton shower approach_, _Phys. Lett._ B449 (1999) 313–320, [hep-ph/9812455].
  * (39) C. Balazs, J. Huston and I. Puljak, _Higgs production: A Comparison of parton showers and resummation_ , _Phys. Rev. D_ 63 (2001) 014021, [hep-ph/0002032].
  * (40) P. Z. Skands, _Tuning Monte Carlo Generators: The Perugia Tunes_ , _Phys. Rev._ D82 (2010) 074018, [1005.3457].
  * (41) S. Catani, M. Ciafaloni and F. Hautmann, _High-energy factorization and small x heavy flavor production_ , _Nucl. Phys._ B366 (1991) 135–188.
  * (42) J. C. Collins and R. K. Ellis, _Heavy quark production in very high-energy hadron collisions_ , _Nucl. Phys._ B360 (1991) 3–30.
  * (43) C. Peterson, D. Schlatter, I. Schmitt and P. M. Zerwas, _Scaling Violations in Inclusive $e^{+}e^{-}$ Annihilation Spectra_, _Phys. Rev._ D27 (1983) 105.
  * (44) M. Lisovyi, A. Verbytskyi and O. Zenaiev, _Combined analysis of charm-quark fragmentation-fraction measurements_ , _Eur. Phys. J._ C76 (2016) 397, [1509.01061].
  * (45) Particle Data Group collaboration, M. Tanabashi et al., _Review of Particle Physics_ , _Phys. Rev._ D98 (2018) 030001.
  * (46) LHCb collaboration, R. Aaij et al., _Measurement of $b$-hadron fractions in 13 TeV $pp$ collisions_, 1902.06794.
  * (47) K. Kovarik et al., _nCTEQ15 - Global analysis of nuclear parton distributions with uncertainties in the CTEQ framework_ , _Phys. Rev._ D93 (2016) 085037, [1509.00792].
  * (48) M. Cacciari, S. Frixione, N. Houdeau, M. L. Mangano, P. Nason and G. Ridolfi, _Theoretical predictions for charm and bottom production at the LHC_ , _JHEP_ 10 (2012) 137, [1205.6344].
  * (49) M. Benzke, M. V. Garzelli, B. Kniehl, G. Kramer, S. Moch and G. Sigl, _Prompt neutrinos from atmospheric charm in the general-mass variable-flavor-number scheme_ , _JHEP_ 12 (2017) 021, [1705.10386].
  * (50) PROSA collaboration, O. Zenaiev, M. V. Garzelli, K. Lipka, S. O. Moch, A. Cooper-Sarkar, F. Olness et al., _Improved constraints on parton distributions using LHCb, ALICE and HERA heavy-flavour measurements and implications for the predictions for prompt atmospheric-neutrino fluxes_ , 1911.13164.
  * (51) LHCb collaboration, R. Aaij et al., _Measurement of the $B^{\pm}$ production cross-section in pp collisions at $\sqrt{s}=$ 7 and 13 TeV_, _JHEP_ 12 (2017) 026, [1710.04921].
  * (52) W. Bai, M. Diwan, M. V. Garzelli, Y. S. Jeong and M. H. Reno _in preparation_ .
  * (53) U. D’Alesio, F. Murgia, C. Pisano and P. Taels, _Probing the gluon Sivers function in $p^{\uparrow}p\to J/\psi\,X$ and $p^{\uparrow}p\to D\,X$_, _Phys. Rev._ D96 (2017) 036011, [1705.04169].
  * (54) LHCb collaboration, R. Aaij et al., _Measurement of the J/ $\psi$ pair production cross-section in pp collisions at $\sqrt{s}=13$ TeV_, _JHEP_ 06 (2017) 047, [1612.07451].
  * (55) J.-P. Lansberg, C. Pisano, F. Scarpa and M. Schlegel, _Pinning down the linearly-polarised gluons inside unpolarised protons using quarkonium-pair production at the LHC_ , _Phys. Lett._ B784 (2018) 217–222, [1710.01684].
  * (56) A. Bacchetta, G. Bozzi, M. Radici, M. Ritzmann and A. Signori, _Effect of Flavor-Dependent Partonic Transverse Momentum on the Determination of the $W$ Boson Mass in Hadronic Collisions_, _Phys. Lett._ B788 (2019) 542–545, [1807.02101].
  * (57) S. Frixione, P. Nason and C. Oleari, _Matching NLO QCD computations with Parton Shower simulations: the POWHEG method_ , _JHEP_ 11 (2007) 070, [0709.2092].
  * (58) S. Frixione, P. Nason and G. Ridolfi, _A Positive-weight next-to-leading-order Monte Carlo for heavy flavour hadroproduction_ , _JHEP_ 09 (2007) 126, [0707.3088].
  * (59) S. M. Barr, T. K. Gaisser, P. Lipari and S. Tilav, _Ratio of $\nu_{e}/\nu_{\mu}$ in Atmospheric Neutrinos_, _Phys. Lett._ B214 (1988) 147–150.
  * (60) L. Pasquali and M. H. Reno, _Tau-neutrino fluxes from atmospheric charm_ , _Phys. Rev._ D59 (1999) 093003, [hep-ph/9811268].
  * (61) A. Bhattacharya, R. Enberg, Y. S. Jeong, C. S. Kim, M. H. Reno, I. Sarcevic et al., _Prompt atmospheric neutrino fluxes: perturbative QCD models and nuclear effects_ , _JHEP_ 11 (2016) 167, [1607.00193].
  * (62) W. Bai and M. H. Reno, _Prompt neutrinos and intrinsic charm at SHiP_ , _JHEP_ 02 (2019) 077, [1807.02746].
  * (63) S. Kretzer and M. H. Reno, _Tau neutrino deep inelastic charged current interactions_ , _Phys. Rev._ D66 (2002) 113007, [hep-ph/0208187].
  * (64) S. Kretzer and M. H. Reno, _Target mass corrections to electroweak structure functions and perturbative neutrino cross-sections_ , _Phys. Rev._ D69 (2004) 034002, [hep-ph/0307023].
  * (65) Y. S. Jeong and M. H. Reno, _Quark mass effects in high energy neutrino nucleon scattering_ , _Phys. Rev._ D81 (2010) 114012, [1001.4175].
  * (66) Y. S. Jeong and M. H. Reno, _Tau neutrino and antineutrino cross sections_ , _Phys. Rev._ D82 (2010) 033010, [1007.1966].
  * (67) M. H. Reno, _Electromagnetic structure functions and neutrino nucleon scattering_ , _Phys. Rev._ D74 (2006) 033001, [hep-ph/0605295].
  * (68) P. Lipari, M. Lusignoli and F. Sartogo, _The Neutrino cross-section and upward going muons_ , _Phys. Rev. Lett._ 74 (1995) 4384–4387, [hep-ph/9411341].
  * (69) S. Kretzer and M. H. Reno, _sigma DIS (nu N), NLO perturbative QCD and O(1 GeV) mass corrections_ , _Nucl. Phys. Proc. Suppl._ 139 (2005) 134–139, [hep-ph/0410184].
  * (70) “Cosmic ray monte carlo package.” https://web.ikp.kit.edu/rulrich/crmc.html.
  * (71) T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko and K. Werner, _EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider_ , _Phys. Rev._ C92 (2015) 034906, [1306.0121].
  * (72) S. Ostapchenko, _Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: I. QGSJET-II model_ , _Phys. Rev._ D83 (2011) 014018, [1010.1869].
  * (73) R. Engel, A. Fedynitch, T. K. Gaisser, F. Riehn and T. Stanev, _The hadronic interaction model Sibyll 2.3c and extensive air showers_ , 1912.03300.
  * (74) A. Fedynitch, F. Riehn, R. Engel, T. K. Gaisser and T. Stanev, _Hadronic interaction model sibyll 2.3c and inclusive lepton fluxes_ , _Phys. Rev._ D100 (2019) 103018, [1806.04140].
  * (75) F. Riehn, H. P. Dembinski, R. Engel, A. Fedynitch, T. K. Gaisser and T. Stanev, _The hadronic interaction model SIBYLL 2.3c and Feynman scaling_ , _PoS_ ICRC2017 (2018) 301, [1709.07227].
  * (76) V. Cirigliano, M. Knecht, H. Neufeld, H. Rupertsberger and P. Talavera, _Radiative corrections to K(l3) decays_ , _Eur. Phys. J._ C23 (2002) 121–133, [hep-ph/0110153].
  * (77) LHCf collaboration, O. Adriani et al., _Measurements of longitudinal and transverse momentum distributions for neutral pions in the forward-rapidity region with the LHCf detector_ , _Phys. Rev._ D94 (2016) 032007, [1507.08764].
  * (78) TOTEM collaboration, G. Antchev et al., _Measurement of the forward charged particle pseudorapidity density in pp collisions at $\sqrt{s}=8$ TeV using a displaced interaction point_, _Eur. Phys. J._ C75 (2015) 126, [1411.4963].
  * (79) CMS collaboration, A. M. Sirunyan et al., _Measurement of the inclusive energy spectrum in the very forward direction in proton-proton collisions at $\sqrt{s}=13$ TeV_, _JHEP_ 08 (2017) 046, [1701.08695].
  * (80) I. Esteban, M. C. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni and T. Schwetz, _Global analysis of three-flavour neutrino oscillations: synergies and tensions in the determination of $\theta_{23}$, $\delta_{CP}$, and the mass ordering_, _JHEP_ 01 (2019) 106, [1811.05487].
  * (81) C. Giunti and T. Lasserre, _eV-scale Sterile Neutrinos_ , _Ann. Rev. Nucl. Part. Sci._ 69 (2019) 163–190, [1901.08330].
  * (82) IceCube collaboration, B. J. P. Jones, _IceCube Sterile Neutrino Searches_ , _EPJ Web Conf._ 207 (2019) 04005, [1902.06185].
  * (83) M. Blennow, E. Fernandez-Martinez, J. Gehrlein, J. Hernandez-Garcia and J. Salvado, _IceCube bounds on sterile neutrinos above 10 eV_ , _Eur. Phys. J._ C78 (2018) 807, [1803.02362].
  * (84) T. Ohlsson and H. Snellman, _Three flavor neutrino oscillations in matter_ , _J. Math. Phys._ 41 (2000) 2768–2788, [hep-ph/9910546].
  * (85) W. Li, J. Ling, F. Xu and B. Yue, _Matter Effect of Light Sterile Neutrino: An Exact Analytical Approach_ , _JHEP_ 10 (2018) 021, [1808.03985].
  * (86) NOMAD collaboration, P. Astier et al., _Final NOMAD results on muon-neutrino to tau-neutrino and electron-neutrino to tau-neutrino oscillations including a new search for tau-neutrino appearance using hadronic tau decays_ , _Nucl. Phys._ B611 (2001) 3–39, [hep-ex/0106102].
  * (87) NOMAD collaboration, P. Astier et al., _Search for $\nu_{\mu}\to\nu_{e}$ oscillations in the NOMAD experiment_, _Phys. Lett._ B570 (2003) 19–31, [hep-ex/0306037].
  * (88) M. Dentler, A. Hernández-Cabezudo, J. Kopp, P. A. N. Machado, M. Maltoni, I. Martinez-Soler et al., _Updated Global Analysis of Neutrino Oscillations in the Presence of eV-Scale Sterile Neutrinos_ , _JHEP_ 08 (2018) 010, [1803.10661].
  * (89) A. I. Belesev, A. I. Berlev, E. V. Geraskin, A. A. Golubev, N. A. Likhovid, A. A. Nozik et al., _An upper limit on additional neutrino mass eigenstate in 2 to 100 eV region from ’Troitsk nu-mass’ data_ , _JETP Lett._ 97 (2013) 67–69, [1211.7193].
  * (90) A. I. Belesev, A. I. Berlev, E. V. Geraskin, A. A. Golubev, N. A. Likhovid, A. A. Nozik et al., _The search for an additional neutrino mass eigenstate in the 2–100 eV region from ‘Troitsk nu-mass’ data: a detailed analysis_ , _J. Phys._ G41 (2014) 015001, [1307.5687].
  * (91) DsTau collaboration, S. Aoki et al., _DsTau: Study of tau neutrino production with 400 GeV protons from the CERN-SPS_ , 1906.03487.
  * (92) A. Bacchetta, D. Boer, C. Pisano and P. Taels, _Gluon TMDs and NRQCD matrix elements in $J/\psi$ production at an EIC_, _Eur. Phys. J._ C80 (2020) 72, [1809.02056].
  * (93) U. D’Alesio, F. Murgia, C. Pisano and P. Taels, _Azimuthal asymmetries in semi-inclusive $J/\psi\,+\,\mathrm{jet}$ production at an EIC_, _Phys. Rev._ D100 (2019) 094016, [1908.00446].
  * (94) A. Bhattacharya, R. Enberg, M. H. Reno, I. Sarcevic and A. Stasto, _Perturbative charm production and the prompt atmospheric neutrino flux in light of RHIC and LHC_ , _JHEP_ 06 (2015) 110, [1502.01076].
  * (95) R. Gauld, J. Rojo, L. Rottoli, S. Sarkar and J. Talbert, _The prompt atmospheric neutrino flux in the light of LHCb_ , _JHEP_ 02 (2016) 130, [1511.06346].
  * (96) PROSA collaboration, M. V. Garzelli, S. Moch, O. Zenaiev, A. Cooper-Sarkar, A. Geiser, K. Lipka et al., _Prompt neutrino fluxes in the atmosphere with PROSA parton distribution functions_ , _JHEP_ 05 (2017) 004, [1611.03815].
  * (97) G. Barr, T. K. Gaisser and T. Stanev, _Flux of atmospheric neutrinos_ , _Phys. Rev. D_ 39 (Jun, 1989) 3532–3534.
  * (98) T. Gaisser, _Cosmic rays and particle physics_.  Cambridge University Press, 1990.