Prompt neutrino fluxes in the atmosphere with PROSA parton distribution functions

Effects on atmospheric prompt neutrino fluxes of present uncertainties affecting the nucleon composition are studied by using the PROSA fit to parton distribution functions (PDFs). The PROSA fit extends the precision of the PDFs to low x, which is the kinematic region of relevance for high-energy neutrino production, by taking into account LHCb data on charm and bottom hadroproduction. In the range of neutrino energies explored by present Very Large Volume Neutrino Telescopes, it is found that PDF uncertainties are far smaller with respect to those due to renormalization and factorization scale variation and to assumptions on the cosmic ray composition, which at present dominate and limit our knowledge of prompt neutrino fluxes. A discussion is presented on how these uncertainties affect the expected number of atmospheric prompt neutrino events in the analysis of high-energy events characterized by interaction vertices fully contained within the instrumented volume of the detector, performed by the IceCube collaboration.


Introduction
The production of leptons due to the interaction of cosmic rays with the Earth's atmosphere introduces a background in the searches for astrophysical neutrinos performed by means of Very Large Volume Neutrino Telescopes (VLVνTs), like IceCube and KM3NeT [1,2]. Although at least part of this background can be eliminated by using veto experimental techniques, accurate predictions for atmospheric lepton fluxes are of crucial importance both to refine or further develop these techniques and to get a precise estimate of the actual spectrum of astrophysical neutrinos (the signal), in terms of its slope and normalization, after subtraction of the background. This may allow to understand the origin of astrophysical neutrinos, distinguishing between different potential galactic and extragalactic sources. The importance of this issue relies on the fact that neutrinos may travel to the Earth from very distant sources, subject only to weak interactions, undeflected by cosmic magnetic fields, thus providing direct information on the sources and regions of the sky where they are produced and where cosmic rays are probably accelerated. In this respect neutrinos are unique carriers of information, complementing those obtained by the detection of other messengers (charged cosmic rays, gamma rays, gravitational waves. . . ) in the multimessenger approach to astroparticle physics, allowing to reach a deeper understanding of our Universe. Furthermore, an accurate knowledge of the atmospheric neutrino spectrum, at least for energies up to a few TeV, is fundamental for the precise determination of the parameters governing the oscillations of atmospheric neutrinos and for understanding the neutrino mass hierarchy [3]. Nowadays this research can be performed by making use of extensions of VLVνTs, like the DeepCore instrumentation and its extension PINGU in JHEP05(2017)004 the inner core of IceCube [4], or the ORCA project in KM3NeT [5], complementing more traditional searches for atmospheric neutrino oscillations.
Atmospheric leptons are typically separated into two components, according to the decaying hadrons from which they originate: the conventional component is created by the decay of light mesons, i.e. charged pions and kaons, whereas the prompt component is generated by the decay of heavier mesons and baryons, including charm or bottom valence quarks. In particular, the prompt component is expected to become larger than the conventional component at energies larger than ∼ 6 +12 −3 × 10 5 GeV [6]. At present, uncertainties affect the shape and normalization of both the conventional and prompt contributions. As for the conventional component, uncertainties of about 10%-20% are quoted in the literature [7]. Uncertainties on the prompt component are larger, due to the poor knowledge of the charm hadroproduction process [6,8].
In this paper, we investigate how the uncertainties on nucleon composition in terms of parton distribution functions (PDFs) affect prompt neutrino fluxes, by making use of the PROSA PDF fit [9], and we compare this uncertainty to that from QCD and astrophysical origin. Current PDF fits are based on data collected at the HERA ep collider, in fixed-target Deep Inelastic Scattering (DIS) and Drell-Yan experiments and in some cases also make use of pp data from Tevatron and pp data from the Large Hadron Collider (LHC) at central rapidities. Due to the limited phase-space coverage in Bjorken-x by these experiments all the data allow to constrain PDFs only in the range x 10 −4 . The PROSA fit was the first exploiting the potential of the LHCb data in constraining gluon PDFs at low Bjorken-x values, encouraging global fits to include these data as well. A first study in this direction was subsequently done in ref. [10], including, a-posteriori, information from the LHCb charm data at √ s = 7 TeV, on top of the NNPDF3.0 NLO fit, by using the Bayesian reweighting method [11,12], leading to the so-called NNPDF3.0 + LHCb set. 1 The method applied in the extraction of these PDFs is detailed in ref. [10], together with the approximations inherent to the heavy flavour scheme adopted. The PROSA analysis [9] contains a full PDF fit in the self-consistent fixed flavour number scheme employed, and also includes LHCb data on D ± s , D * ± and on beauty hadroproduction [16] at √ s = 7 TeV. These developments have allowed to constrain gluon PDFs for x 10 −4 , a regime that is important to study in view of high-energy applications such as discussed in this work. In fact, VLVνTs have so far detected neutrinos with deposited energies up to a few PeV in the laboratory frame, i.e. E ν, lab ∼ (few) × 10 6 GeV. In this context it is important to note that the production of neutrinos at a given laboratory energy E ν, lab is affected by collisions of cosmic rays with nuclei in the Earth's atmosphere at energies E p, lab even larger by a factor of order O(10-1000). Further, E p, lab ∼ 10 8 GeV corresponds to E pp, cm ∼ 13.7 TeV, close to the present LHC center-of-mass energy. Plans exist to extend the accessible neutrino energy range even further, by increasing the active volume of ice/water probed by 1 Very recently, during the completion of our paper, the study of ref. [10] was extended further in ref. [13] by using additional selected LHCb charm data at √ s = 13 and 5 TeV [14,15]. In particular, charm data in the transverse momentum range limited to 0 < pT < 8 GeV were considered. The data on D ± s at 7 TeV were not considered. In the case of (D 0 +D 0 ), data in the rapidity bins far from the rapidity reference bin 3 < y < 3.5 were discarded.

JHEP05(2017)004
the optical module instrumentation. With the IceCube-Gen2 program [17], an extension from an instrumented fiducial volume of ∼ 1 km 3 to a larger one, of ∼ 10 km 3 , has been proposed, which claims the potential to deliver significant samples of neutrinos with laboratory energies in the EeV range, i.e. increased by a factor of ∼ 1000. In nucleon-nucleon collisions such high-energy neutrinos may arise from collision energies E pp, cm up to a few hundred TeV.
So far IceCube has detected several tens of leptonic events with a characteristic topology of tracks or showers, with interaction vertices fully contained within the instrumented volume of the detector. This feature, in addition to a minimum charge/energy requirement, define a category of so-called high-energy-starting events (HESE), in relation to their specific signatures. 2 Three of those HESE events populate the high-energy tail in the O(PeV) region of neutrino deposited energies in the laboratory frame. As a practical application of our study we consider the case of the IceCube HESE analyses [1,18,19], and show how the theoretical uncertainties affect the number of expected events from prompt neutrinos as a function of the energy deposited in the detector. We also compare our predictions to the upper limit on the prompt (ν µ +ν µ ) flux recently obtained by the IceCube collaboration in a complementary analysis, restricted to ∼ 350000 up-going events from the Northern hemisphere, with neutrino interaction vertices either outside or inside the instrumented volume [20].
The paper is organized as follows: a short review of the PROSA PDF fit and of the details of the QCD framework we adopt is presented in section 2. Predictions for prompt neutrino fluxes and their uncertainties are discussed in section 3, whereas comparisons with the data collected in the aforementioned IceCube analyses are provided in section 4, followed by our conclusions in section 5. Appendix A contains details of the QCD predictions for charm hadroproduction compared with the LHCb experimental data taken in collisions at center-of-mass energy √ s = 7 and 13 TeV.

Charm hadroproduction in QCD and parton distribution functions
The hadroproduction of heavy-quark pairs has received a lot of attention since the time when early colliders were built. The attention focused on the production of charmanticharm and bottom-antibottom pairs [21,22] already well before the advent of the high-energy colliders Tevatron and LHC. The charm and bottom masses m c and m b are large enough with respect to Λ QCD , so that the use of perturbative QCD is justified to describe the hard partonic scatterings giving rise to these quarks. However, these masses are far smaller than present day collider energies. Thus, processes involving charm and bottom hadroproduction at high energies are typically multiscale processes for which it is worth to investigate the role of logarithms of the ratios of the different scales involved. The effect of resummation of different kinds of logarithms has been studied in some specific cases (e.g. resummation up to a certain accuracy of threshold logarithms and of transverse JHEP05(2017)004 momentum dependent logarithms (p T /m q ) at high p T ), while there are other resummations which have not yet been performed. NLO QCD corrections to charm hadroproduction are available since the eighties [23,24], and have been further implemented in Monte Carlo event generators matching NLO QCD corrections to parton shower (starting from pioneering work of [25] and [26]), which are routinely used in the analysis of LHC data. NLO electroweak corrections are available in the case of top-antitop quark hadroproduction and, consistently with the naive expectations, they have been found to be far smaller than the QCD ones, see e.g., [27,28]. In case of cc and bb hadroproduction the electroweak corrections are expected to lie well within the NLO QCD scale uncertainty band. Thus we neglect electroweak corrections in the present work. On the other hand, NNLO QCD predictions are already available at least for the total cross-section [6,8,29]. Extensions of the first differential studies on tt pairs at NNLO [30,31] to the case of cc pairs still require further effort, particularly to take into account the small value of the charm mass, m c < 2 GeV, which may lead to stability issues in the fully numerical methods used so far for investigating heavy-quark hadroproduction at NNLO. The calculation of charm-pair hadroproduction in this paper follows ref. [6] and uses recent results from perturbative QCD together with estimates on various sources of uncertainties. In particular, we consider NLO QCD corrections matched to parton showers (with their standard logarithmic accuracy) as implemented in the POWHEGBOX [32] + PYTHIA6 [33]/PYTHIA8 [34] frameworks. We used the last available PYTHIA6 version, PYTHIA 6.4.28, and the most recent PYTHIA8 version (PYTHIA 8.223, published in January 2017). Hadronization is also taken into account by means of PYTHIA6/PYTHIA8, with non-perturbative parameters regulated according to one of the most recent Perugia tunes [35] in case of PYTHIA6 and the A14 [36] or the Monash 2013 [37] tune in case of PYTHIA8. For all predictions shown in the plots of this paper, the PYTHIA8 generator and A14 tune are used by default, unless stated otherwise. The POWHEGBOX framework includes hard-scattering matrix-elements for charm hadroproduction according to refs. [24,26]. The charm is considered massive throughout the calculation of the scattering amplitudes performed in the fixed flavour number scheme (n f = 3). Consistently, we adopt 3-flavour PDFs. POWHEGBOX is capable of producing events in the Les Houches format, including up to a first resolved radiation emission. These events are subsequently further showered by the p T -ordered parton shower algorithms included in PYTHIA6/PYTHIA8. The POWHEGBOX framework was also used in ref. [38] in association with PYTHIA8 and the Monash 2013 tune.
Alternatively, it is possible to use approaches where predictions at the parton level are directly convoluted with phenomenological non-perturbative fragmentation functions (FFs) for the transformation of partons into charmed hadrons. The latter choice was adopted in the computation of prompt neutrino fluxes by refs. [39,40]. It was also adopted when performing the original PROSA fit for PDFs, in the computation of theoretical predictions for charmed-and bottomed-meson hadroproduction which were compared with LHCb experimental data. In particular, data for D ± , D ± s , Λ ± c , D * ± , D 0 andD 0 hadroproduction [41] in pp collisions at √ s = 7 TeV in different transverse momentum and rapidity bins in the interval 0 < p T < 8 GeV and 2 < y < 4.5 and data for B + , B 0 and B 0 s meson production [16] with 0 < p T < 40 GeV in the same rapidity range and at the JHEP05(2017)004 same center-of-mass energy were used in the fit, as detailed in ref. [9]. The performances of the PROSA PDF fit in association with our POWHEGBOX + PYTHIA8/PYTHIA6 setup in comparison with these data at both √ s = 7 and 13 TeV are shown in appendix A. We emphasize here that the PROSA PDFs used in this paper were derived by fitting the ratios of LHCb rapidity distributions in each fixed p T bin. It turns out that we obtain absolute p T distributions reasonably consistent with LHCb experimental data within QCD uncertainties. We consider this fact as a confirmation of the robustness of the PROSA fit. This consideration remains valid both if we produce theory predictions for differential cross-sections by the NLO QCD computation interfaced with FFs adopted in the PROSA fit, with fragmentation fractions in agreement with the most recent measurements [42], and if we use for the same purpose the approach of this paper, matching NLO QCD matrix elements with parton shower and hadronization, pointing to the consistency of the two theoretical approaches.

Predictions for prompt neutrino fluxes and their uncertainties
The collision of primary cosmic rays with the Earth's atmosphere lead to secondary particles, which during their propagation towards the Earth surface may in turn reinteract or decay. This mechanism, characterized by an interplay between interactions and decays, leads to the production of fluxes of different kinds of particles, including neutrinos, and is described by cascade equations, which define a set of coupled differential equations regulating particle flux evolution through an air column of slant depth X, representing the amount of matter traversed downward along the direction of the particle that initiated the cascade. These equations admit approximate solutions with the help of the so-called Z-moment approach initially proposed in refs. [43,44] (see also refs. [6,39]).
The differential distributions for the hadroproduction of cc quark pairs are one of the essential ingredients to the solution of those equations, yielding the prompt neutrino fluxes. The necessary Z-moments, e.g., Z p hc for charm hadroproduction in the atmosphere, are then obtained from the differential distributions dσ/dx E for the process N N → cc → h c + X, for all lowest-lying charmed hadrons h c = (D ± , D 0 ,D 0 , D ± s and Λ ± c ), with x E = E h, lab /E N, lab , by an integration over the entire kinematically accessible range. To that end we assume that the interaction of cosmic rays with nuclei in the atmosphere can be approximated by the superposition of nucleon-nucleon (N N ) interactions.
Those moments Z p hc are combined together with the Z-moments for h c decays into neutrinos Z hc ν , those for h c regeneration, Z hc hc , and those for nucleon regeneration, Z N N , see for instance refs. [6,44,45]. In this work, all moments Z p hc were computed using as input the normalized variant of the PROSA PDF fit [9], with central values for the factorization and renormalization scales as µ R = µ F = µ 0 = p 2 T,c + m 2 c . The scales µ R and µ F were allowed to vary independently in the range µ 0 /2 < µ R , µ F < 2µ 0 , with exclusion of the combinations (1/2, 2)µ 0 and (2, 1/2)µ 0 , see, e.g. ref. [46], and the charm pole mass was fixed to m c = 1.40 ± 0.15 GeV, as in ref. [6]. All other moments Z hc ν , Z hc hc and Z N N were computed as in ref. [6] (see also refs. [39,40]).
The resulting predictions for the (ν µ +ν µ ) fluxes are presented in figures 1, 2 and 3. In particular, the central predictions were computed using as input the central set of the PROSA PDFs, whereas PDF uncertainties were computed using all different PROSA PDF variations. The latter consist of three components, as described in detail in ref. [9]: fit uncertainties originating from experimental uncertainties of the measurements, model uncertainties and parametrization uncertainties, cf. table 1. Model uncertainties arise from the variations of the model assumptions, such as the value of the strong coupling constant α s (M Z ), the strangeness fraction f s in the PDF fit, the minimum virtuality cut Q 2 min on the ep DIS data entering the fit, the choices of the renormalization and factorization scales, µ R and µ F , and the parameters α K of the charm and beauty fragmentation functions [50]. Parametrization uncertainties are assessed by varying the functional form of the PDFs using additional parameters D uv , DŪ , DD as described in refs. [9,51] at the starting scale Q 2 0 of the QCD evolution, as well as the value of the starting scale.
The fit uncertainties (referred to as experimental uncertainties in ref. [9]) were provided as 30 eigenvectors arising from 13 fitted PDF parameters and charm and beauty masses left free in the fit. The model uncertainties were determined as positive and negative variations of each model parameter and were technically provided as 16 eigenvectors. The parametrization uncertainties were assessed with 4 individual variations, and the total parametrization uncertainty was built as an envelope of the maximal differences between these variations and the central value. The total PDF uncertainties were obtained by adding fit, model and parametrization uncertainties in quadrature. In figure 4 the PROSA gluon distribution is compared to the results of other PDF groups [52][53][54][55][56][57][58] in the relevant kinematic region of low Bjorken-x, at a scale of 10 GeV 2 . The use of the LHAPDF 6.1.6 framework [59] for figure 4 ensures the extrapolation of the PDFs to the low-x region beyond the kinematic range provided and fitted by the individual groups [52][53][54][55][56][57][58]. The differences in the behaviour of different PDFs can be traced to a variety of potential origins [8]: among the most important ones we wish to mention the use of different sets of experimental data, different theory assumptions in fitting the data, and the different PDF parameterizations adopted.
The separate contribution of each of the uncertainty sources in the PROSA PDFs listed in table 1 is shown in figure 1 together with the total PROSA PDF uncertainty and using a broken power-law all-nucleon spectrum as input for the cosmic ray flux, cf. ref. [6]. The total PDF uncertainty turns out to be dominated by the model uncertainties, which in turn for the gluon distribution in the region x < 10 −4 are influenced by theoretical uncertainties on heavy-flavor hadroproduction, mostly arising from renormalization and factorization scale variations. However, due to the use of the normalized heavy-flavor LHCb cross sections in the PROSA fit, these theoretical uncertainties are strongly reduced, since variations of the renormalization and factorization scales as well as of the fragmentation parameters do (ν µ + anti-ν µ ) flux total PDF uncertainty parameterization PDF uncertainty Figure 1. PROSA PDF uncertainties on the prompt (ν µ +ν µ ) atmospheric flux as a function of the neutrino energy E ν,lab : the contribution due to (fit + experimental), model and parameterization PDF uncertainties are shown in separate panels, respectively, and compared to the total PDF uncertainty (blue band). A broken power-law all-nucleon spectrum for the cosmic ray flux impinging on the Earth atmosphere is used as input, cf. ref. [6].  Left: central prediction for the prompt (ν µ +ν µ ) flux together with its QCD uncertainties, computed by means of POWHEGBOX + PYTHIA8, as a function of the neutrino energy E ν,lab . The uncertainty contributions due to µ R and µ F scale variation around µ 0 , m c and the PDF eigenvalues within the PROSA fit, are shown separately by bands of different styles and colors, together with their combination in quadrature. The same broken power-law all-nucleon spectrum for the cosmic ray flux as in figure 1 is used as input. Right: comparison between central predictions and (scale + mass + PDF) uncertainty bands obtained by using our setups for POWHEGBOX + PYTHIA8 and POWHEGBOX + PYTHIA6. See text for more detail. The predictions using PYTHIA8 have been obtained by using the A14 tune. In case of PYTHIA8 in association with the Monash 2013 tune, predictions (not plotted here) turn out to be even closer to those with PYTHIA6.
not significantly affect the shape of the rapidity distributions for heavy-flavor production, while this shape remains sensitive to PDFs. For this reason, in the calculation of prompt neutrino flux the scale variations in the matrix elements were performed independently of the scale variations in the PDF fit. The same applies for the value of the charm mass, which is well constrained in the PROSA fit by the data on DIS charm production. The choice and uncertainity for the charm pole mass m c = 1.4±0.15 GeV adopted here are compatible and have been motivated by detailed studies of the inclusive cross sections in ref. [6] computed with the very precise value for the charm mass reported by the Particle Data Group [62] in the MS scheme. In addition, since the scale variations have a much larger effect on the predictions than the m c value, possible correlations between the two can be neglected. The fit and parameterization PDF uncertainties turn out to have a minor role in the kinematic region of interest. The role of PDF uncertainties with respect to other sources of QCD uncertainties affecting (ν µ +ν µ ) fluxes is shown in figure 2 again for the broken power-law primary cosmic-ray input spectrum as in figure   renormalization scale variations. On the other hand, the role of charm mass and PDF uncertainties is complementary: the former dominate over the latter at lower energies, whereas at higher energies the behavior is the opposite, with PDF uncertainties increasing due to the increasing number of collisions occurring in an asymmetric situation when one of the partons inside the nucleon is characterized by low x and the other one by high x. In this respect we would like to note that the present PROSA fit extends down to x ∼ 10 −6 , while below this value the gluon distribution is not directly constrained by any data and should be considered as an extrapolation which relies on assumptions for the parametrization of the PDFs. 4 In particular, for neutrino energies E lab, ν around 1 and 2 PeV, corresponding to the JHEP05(2017)004 leptonic events with highest energy observed so far by IceCube, the PDF uncertainties look to be already far better constrained than the scale ones, and amount to (+42%, −13.5%) and (+52%, −13.5%), respectively, for the power-law cosmic ray spectrum. These values decrease by few percent when considering more realistic input cosmic ray spectra. They can be compared to the charm mass uncertainties amounting to about (+26%, −22%) and (+25.8%, −20%) at 1 and 2 PeV, respectively.
It is plausible that, with the increasing availability of LHC experimental data even at higher energies, PDF uncertainties will decrease. Therefore, we can already conclude that, as for present day investigations of prompt neutrinos at neutrino telescopes, the uncertainty related to PDFs does not form a bottleneck. Of course, if neutrino telescopes will be able to discover events at higher energies, either thanks to the extension of their fiducial volume, or simply by accumulating much more statistics over the years, the uncertainty due to the PDF variation may become a more critical issue.
Another important input to cascade equations is represented by primary cosmic ray fluxes, i.e. the energy spectra of cosmic rays on top of the Earth atmosphere, as a function of their mass number. At energies above a PeV, an important aspect, at least presently, is related to our superficial knowledge of the composition of cosmic ray fluxes [63]. Here the problem is that at the highest energies, cosmic ray spectra cannot be measured directly by satellites or balloon-born experiments, because the flux of cosmic rays decreases too rapidly. It is presently impossible to build an instrument capable of measuring a small, but JHEP05(2017)004 non negligible, flux because the detector surface and exposure time necessary are too big for current capabilities. Therefore balloon-born instruments or those in satellites are used at present days only for measuring cosmic ray spectra at lower energies, i.e. below the knee [64].
On the other hand, at higher energies, cosmic ray spectra are investigated indirectly, through extended air shower (EAS) experiments [65,66]. EAS experiments count the rate of leptons of different origin (e and µ) reaching the array of detectors on the Earth's surface and compare it with the data on the maximum development of the EAS electromagnetic component seen by fluorescence telescopes pointing to the upper layer of the Earth's atmosphere. Last studies in this direction, exploiting correlations between different EAS observables, seem to point to a CR spectrum characterized by a mixed composition in the energy region around the dip/ankle [67], which tends to become heavier at E lab ∼ a few 10 19 eV. However, not all questions are solved, in particular in the comparison of the experimental data with the expectations from the Monte Carlo generators used for simulating the formation and development of EAS, characterized by large uncertainties in their hadronic interaction models [68]. As a consequence, uncertainties on the composition of cosmic ray spectra above E lab = 10 16 -10 17 eV are still large.
In order to have an idea of the effect of these uncertainties on prompt (ν µ +ν µ ) fluxes, we show in figure 3 prompt neutrino spectra with their QCD uncertainties, for four different recent all-nucleon cosmic ray spectra [47][48][49], provided by the group of T. Gaisser, which were recently created in order to fit the measured cosmic ray all-particle spectrum data. The transformation of an all-particle spectrum into an all-nucleon spectrum requires additional knowledge (or assumption) on the composition. We have used these same spectra in ref. [6], where we have provided a more extensive discussion of their content. The resulting neutrino spectra at energies 1 PeV are larger if a proton component or, more generally, a light component in the cosmic rays is assumed and dominates over the nuclear ones at the highest energies. The QCD uncertainties behave in a similar way in all cases, comparable to the case of figure 2 already discussed, but at the highest E lab,ν energies the uncertainties due to our poor knowledge of cosmic ray composition are becoming large, as follows from comparing one with each other the various panels of figure 3.
It is worth to compare our results on cosmic ray fluxes using the PROSA PDFs to those obtained in 2015 with the ABM PDFs in ref.  choice adopted in this paper, µ 0 = p 2 T,c + m 2 c , leads to uncertainty bands larger than for GMS 2015. Note that the scale in the GMS 2015 study was explicitly chosen to fulfill the principle of fastest convergence, i.e. in order to reduce the difference between NLO and NNLO predictions for the total cross section, cf. ref. [6]. In the present study we have used the scale which retains full consistency with the PROSA PDF fit.
Interestingly earlier predictions, evaluated at the central scale, even the leading-order ones presented in ref. [69] several years ago, turn out to lie within the uncertainty band of the PROSA predictions, as shown in figure 6. In this plot, besides these old predictions, we also show more recent ones obtained by different groups in the last few years. Even those predictions which are obtained in approaches quite different from ours, e.g., on the basis of the so-called dipole model, or those obtained by a recent update of the SYBILL Monte-Carlo [72], turn out to be compatible with our prediction, cf. figure 6. This is an important result because dipole model predictions have been extensively used by IceCube as a theory reference for prompt neutrino fluxes, while the SYBILL Monte-Carlo is one of the event generators extensively used in EAS simulations. It includes a perturbative QCD leading-order core and a soft phenomenological component for the description of hadronic interactions.
Finally, we briefly address the superposition model assumption used so far, in which the collision of cosmic rays with air nuclei, mostly nitrogen, are approximated by the mere superposition of nucleon interactions, i.e. the N N hard scattering described within perturbative QCD. Initial work to test the superposition approximation has recently been presented in ref. [40], where also two nuclear PDFs were used in the cross section computation of cosmic ray interactions with atmospheric nuclei. In figure 7 we compare the predictions of ref. [40] obtained with nuclear PDFs EPS09 [73] and nCTEQ15 [74] with  . Comparison of our predictions for (ν µ +ν µ ) fluxes with the PROSA proton PDFs and the superposition approximation, to those of ref. [40] using nuclear PDFs, with their respective uncertainty bands. Predictions of ref. [40] using the EPS09 nuclear PDFs are shown on the left, those based on the nCTEQ15 nuclear PDFs on the right. Uncertainties affecting nuclear PDFs are not accounted for in these plots. For consistency with ref. [40], and differently from figure 5 and 6, the cosmic ray primary flux H3p is used in all these predictions, cf. figure 3. our own predictions. Figure 7 shows that the predictions of ref. [40] using nuclear PDFs, although being systematically lower than ours, are still close to the lower limit of, or stay within, our total QCD uncertainty band. At present, any uncertainties on nuclear PDFs are not included in these plots. However, these uncertainties are actually quite large in the full x range (i.e. also for x > 10 −4 ), due to the limited amount of nuclear data that can be used in the nuclear PDF fits and the limited kinematic coverage of those data, cf. [73,74]. In addition, model assumptions for the description of hadronic interactions involving nuclei JHEP05(2017)004 are ingredients entering nuclear PDF fits. The uncertainty in this modellization introduces further uncertainties on the resulting PDFs. Thus, in order to fully address the superposition approximation, a systematic experimental and theoretical study of collisions of nuclei at high energies, including nuclear and possibly quark-gluon plasma effects, is required. This is clearly beyond the scope of the present study and we leave these improvements and new developments on this issue for the future.
4 Uncertainties on prompt neutrino expected events in the IceCube HESE analysis and comparison of the prompt (ν µ +ν µ ) flux with the present IceCube upper limit One of the most intriguing results of IceCube has been reported in the HESE analysis. Over the years this analysis has collected leptonic events also known as contained-vertex events, where the incoming lepton has its first interaction inside the active volume of the detector, with a deposited energy in the range from tens of TeV to a few PeV. In this analysis, an excess with respect to the expected atmospheric background considered so far by the IceCube collaboration has been registered with increasing statistical evidence along the years. The experimental data include events with shower and muon track topologies, coming from both the Northern and Southern hemisphere. In particular, results were reported for a 662-day analysis with 28 candidates in the energy range [50 TeV-2 PeV], corresponding to a 4.1 σ excess in the year 2013 [1]. These results were updated in the year 2014, thanks to the 988-day analysis [18], with 37 events in the energy range [30 TeV-2 PeV] (5.7 σ excess over the background), and featuring also an "empty" window corresponding to the [400 TeV-1 PeV] interval. A further update was presented in 2015, in the 1347day analysis [19], which has collected 54 events, corresponding to a ∼ 7 σ rejection of the atmospheric-only hypothesis, and with the empty window partially filled and thus reduced to the [∼600 PeV-1 PeV] bin. In these analyses, the experimental data were fitted considering two possible kinds of sources of neutrinos: the atmosphere and astrophysical ones. For the shape of the astrophysical signal a power-law was assumed. The best-fit parameter values turned out to change slightly from one analysis to the other, although they always remained compatible with each other within the quoted uncertainties. For the atmospheric component, the possibilities of both a conventional and a prompt contribution were considered. In particular, the IceCube collaboration used the Honda predictions [7], extended to higher energies and modified for taking into account more recent cosmic ray primary spectra with a knee component, as theoretical input for the modellization of conventional neutrino fluxes in these analyses. A prior was used for the normalization of this contribution in the fit to IceCube experimental data, performed under the assumption that leptons of both astrophysical and atmospheric origin contribute to the total signal seen by IceCube.
In the most recent versions of the analysis, part of the atmospheric background was vetoed and subtracted from the signal, by using information on the atmospheric muons detected in coincidence with neutrino events and the techniques of ref. [75] in the modellization of this case. As for prompt neutrinos, IceCube has used the ERS 2008 predictions [71], reweighted to a cosmic ray spectrum with a knee component, as a basis for modelling prompt neutrino fluxes. This component was included in the fit as well, and it happened its normalization fits to zero, although with a large uncertainty. As a consequence of this big uncertainty, an upper limit on the total atmospheric neutrino flux, as derived from other IceCube analyses (discussed in the second part of this section), was adapted as well to the HESE analyses, in order to show how big a potential (prompt + conventional) component could look like.
Given the model dependence understood in the description of atmospheric neutrinos by IceCube, it is instructive to investigate the effect of predictions for prompt neutrino fluxes considered in this study, on the expected number of events in the IceCube analyses. In particular, we consider the case of the 988-day analysis, discussed previously also in ref. [39]. We note that even considering the most recent IceCube 1347-day analysis gives rise to qualitatively similar results and does not alter our discussion and conclusions. The number of prompt neutrino events in the different energy bins reported by IceCube, computed on the basis of PROSA predictions for prompt neutrino fluxes, is shown in figure 8 for the various cosmic ray spectra already presented in section 3. It is evident that modern CR spectra give a suppressed number of events with respect to the broken power-law spectrum (red band). Furthermore, for those events seen by IceCube so far with the maximum energies deposited, limited to ∼ 2 PeV, the difference between the variants of Gaisser CR spectra considered in this work do not introduce any dramatic changes in the number of observed events. The bands reported in the plot refer to QCD theoretical uncertainties stemming from the combination of scale, charm mass and PDF uncertainties, computed as described previously in section 3. Even considering these uncertainty bands, our theoretical predictions for prompt fluxes turn out to lie below the data, thus confirming a different origin for those most energetic IceCube events, also shown in the plot.  Figure 9. Predictions for the number of prompt neutrino events as a function of the deposited energy in the detector for the IceCube 988-day HESE analysis. The H3a cosmic ray primary flux is used as input. Results from the PROSA calculation of this paper, shown in blue, with their QCD uncertainty band, are compared with the central predictions of GMS 2015 [6] and with the BERSS one [39], including the uncertainty band of the latter. IceCube experimental data on the total number of leptons detected are also shown. These data include leptons of both atmospheric and astrophysical origin.
A comparison between our predictions for prompt neutrino events, and predictions computed by using different models is shown in figure 9, considering as a basis for the cosmic ray primary flux the H3a model. In particular, predictions computed by using the GMS 2015 central flux [6] and those reported in ref. [39] are shown. These earlier predictions are compatible with the PROSA predictions. Furthermore, it turns out that the uncertainty band of ref. [39] is completely contained, in all the event bins, inside the PROSA uncertainty band. The latter is much larger because it includes a broader range of scale variations (with µ R = µ F ) and also the PDF uncertainties, not considered in [39].
The total atmospheric flux, which can be computed by summing the conventional and the prompt flux, is plotted in figure 10. As for the conventional flux, the Honda spectrum [7] was adopted, extended to the highest energies and reweighted to the H3a cosmic ray flux as described in ref. [39]. Uncertainties on the total flux, due to uncertainties on the prompt component, are shown, under the assumption that the conventional flux does not contribute any additional uncertainty. 5 The IceCube upper limit on the total neutrino flux at 90% confidence level [76] is also shown. At the highest energies, the IceCube upper limit lies well inside our uncertainty band. The latter limit corresponds to the result of a separate analysis of muonic events from the Northern hemisphere, characterized by neutrino interaction vertices which can lie both inside and outside the instrumented volume, with the Earth acting as an efficient shielding for atmospheric muons. IceCube has progressively updated IceCube data Figure 10. Predictions for the number of prompt, conventional and total expected atmospheric neutrino events for the IceCube 988-day HESE analysis, as compared to the IceCube lepton data. The H3a cosmic ray primary flux is used as input. On the top the central predictions are shown and also compared to the total central predictions by BERSS [39]. On the bottom the effect of uncertainties on prompt neutrino contribution is propagated to an uncertainty on the total atmospheric neutrino contribution under the assumption that the conventional component does not contribute any additional uncertainty. IceCube experimental data are shown in black, and the IceCube 90% confidence level upper limit on the total atmospheric neutrino flux obtained in ref. [76] and reproduced in ref. [18] is shown in magenta.
this analysis in refs. [20,76,77]. As for the prompt neutrino component, the experimental results of this kind of analysis are presented in the form of upper limits on the prompt (ν µ +ν µ ) flux. 6 The comparison of the most recent IceCube estimate of this limit [20] JHEP05(2017)004  Figure 11. Comparison of the prompt (ν µ +ν µ ) flux using the PROSA PDFs with the present upper limit on prompt neutrino flux at 90% confidence level obtained by the IceCube experiment [20] (solid red line) and its extrapolation (dotted red line), which adopted the ERS model [71] as a basis for modelling prompt neutrinos. Central predictions using the scale µ R = µ F = p 2 T + 4m 2 c and PROSA PDFs and ABM PDFs (GMS 2015) are also shown. The limit and all predictions refer to the H3p CR flux.
with the fluxes computed in the present study and the central prediction of GMS 2015 is shown in figure 11. Again, the published IceCube upper limit, although being larger than our central prompt predictions, is well inside the uncertainty band of the prompt (ν µ +ν µ ) flux, both over the limited range of neutrino energies E lab probed by present IceCube data [8 · 10 3 -8 · 10 4 ] GeV and even when considering the extrapolation to a larger energy range [10 3 -10 7 ] GeV, also presented by the IceCube collaboration. This might point to a need of revising the model assumptions in the aforementioned IceCube analyses.
In summary, figures 10 and 11 show the future potential of astrophysical measurements at VLVνTs, especially when a higher statistical accuracy will be reached, in complementing accelerator-based measurements by putting constraints and providing complementary information on the physics related to charm hadroproduction.

Conclusions
We have used the PROSA PDF fit to provide predictions for the flux of prompt neutrinos in the atmosphere. The PROSA PDFs are the first fit to include LHCb open charm and beauty data in order to constrain the gluon distribution inside the proton in regions of Bjorken-x not previously covered by any other experiment. We have shown that present PDF uncertainties on the prompt neutrino flux increase with increasing neutrino energies. Moreover, in the region of interest for present day neutrino telescopes, which have so far JHEP05(2017)004 detected neutrinos up to a few PeV, PDF uncertainties are already quite well constrained and are subdominant with respect to the dominating QCD uncertainties related to the renormalization and factorization scale variation. Our flux turns out to be compatible with several previous computations of prompt fluxes, obtained in a variety of approaches, which may or may not involve a perturbative QCD description of the hard scattering.
As a practical application, we have studied the uncertainties on the number of expected prompt neutrino events in the IceCube HESE analysis and we have compared the theoretical predictions on (ν µ +ν µ ) fluxes with the IceCube upper limit published in complementary analyses, including also non-contained events. We have found that the adoption of different assumptions for the composition of the cosmic ray primary flux has a small effect on the shape of the distribution of prompt neutrino events in the HESE analysis, at least when considering the energy range tested so far by IceCube, and that the high energy tail of the atmospheric neutrino flux has a steeper slope than the slope of IceCube events, even when including in the analysis prompt neutrino uncertainties of QCD origin. This confirms that, to explain IceCube HESE events, is necessary to add the existence of at least one additional neutrino component of non-atmospheric origin. Furthermore, we have found that the upper limit on prompt neutrino fluxes at 90% confidence level, published by IceCube, although being model dependent, is just slightly above the central predictions for the prompt (ν µ +ν µ ) flux obtained in this study, but well inside our global QCD uncertainty band. This holds over the entire range of relevant neutrino energies and challenges the model assumptions on atmospheric fluxes at high-energies adopted in the IceCube analyses.
In summary, this paper has presented the first application of the PROSA PDFs to an astrophysical problem. This opens up the possible use of these PDFs in many other problems arising in the description of microscopic interactions at low Bjorken-x at present and future high-energy pp and ep colliders, and in cosmic ray hadronic interactions occurring in their astrophysical sources, during propagation and in the atmosphere.
In future, the LHCb measurements of charm hadroproduction at √ s = 13 TeV [14] could be used for checking the extrapolation of the present PROSA fit and further reduce the PDF uncertainties at low Bjorken-x. On the other hand, incorporating into the fit recently appeared LHCb charm data at √ s = 5 TeV [15], could provide insights on the self-consistency of the data and of their theoretical description. Likewise, the inclusion of nuclear effects in the theory description of the hadroproduction of D-mesons in collisions of ultra-high-energy cosmic ray nuclei with the nuclei of our atmosphere will be left for future work.
Our predictions for prompt neutrino fluxes with PROSA PDFs are publicly available at https://prosa.desy.de. The contribution arising from the feeddown from D * 's and from other excited charmed states is accounted for in the theory predictions for the lowest-lying charmed mesons shown in the plots.
While the data at 7 TeV are succesfully described in all rapidity and p T bins by the QCD computation adopted in this paper, the data at 13 TeV in the low p T regime turn out to lie above the central theoretical predictions. 7 This tendency remains essentially unchanged when using POWHEGBOX + PYTHIA6 instead of POWHEGBOX + PYTHIA8, with the former giving slightly harder distributions. In any case, for (D 0 +D 0 ) and D ± , whose contribution actually dominate the total prompt neutrino flux at the neutrino energies explored by IceCube, giving rise to more than 80% of it, the LHCb experimental data turn out to lie within our scale uncertainty limits in all rapidity bins, with the exception of (D 0 +D 0 ) data in the 1.5 < p T < 7 GeV interval, in case of the most central bin 2 < y < 2.5. We expect that even in case the theoretical description of these data will be improved in the future, the resulting central predictions for prompt neutrino fluxes will likely lie within our present uncertainty bands.  [41] in different rapidity bins. Theoretical predictions are accompanied by their uncertainty bands, due to µ R and µ F scale variation (green), to m c (magenta) and to PROSA PDF (light-blue hatched) variation, as described in the text. The LHCb experimental data [41] are shown together with their statistical and systematic uncertainties, added in quadrature. In the lower panel, ratios of the uncertainties with respect to the theoretical central predictions are shown, together with the ratio of the central theoretical predictions obtained with POWHEGBOX + PYTHIA6 with respect to those with POWHEGBOX + PYTHIA8, and the ratio of experimental data with respect to theoretical predictions by POWHEGBOX + PYTHIA8. dσ / dp T ( pb / GeV )

JHEP05(2017)004
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.