Lepton fluxes from atmospheric charm revisited

We update predictions for lepton fluxes from the hadroproduction of charm quarks in the scattering of primary cosmic rays with the Earth's atmosphere. The calculation of charm-pair hadroproduction applies the latest results from perturbative QCD through next-to-next-to-leading order and modern parton distributions, together with estimates on various sources of uncertainties. Our predictions for the lepton fluxes turn out to be compatible, within the uncertainty band, with recent results in the literature. However, by taking into account contributions neglected in previous works, our total uncertainties are much larger. The predictions are crucial for the interpretation of results from neutrino experiments like IceCube, when disentangling signals of neutrinos of astrophysical origin from the atmospheric background.


Introduction
Atmospheric lepton fluxes are important backgrounds in the search of neutrinos of astrophysical origin [1]. In particular recent claims from the IceCube experiment, which has detected a statistically significant sample of leptonic events at very high energies [2,3], whose interpretation is still under debate [4,5], require an estimate of the background as accurately as possible. One of the most uncertain components of this background is the prompt contribution due to the hadroproduction of charm quarks in the hard scattering of primary cosmic rays with the Earth's atmosphere, the so-called atmospheric charm. In this paper we will concentrate on the contribution to the lepton fluxes that can be ascribed to atmospheric charm.
After initial studies on the basis of phenomenological models (see e.g. Ref. [6,7] and references therein), previous predictions for lepton fluxes from atmospheric charm have been obtained within the framework of perturbative Quantum Chromodynamics (QCD) for proton-proton collisions according to the standard QCD collinear factorization formalism, with hard-scattering evaluated at leading order (LO) in Ref. [8] and including radiative corrections at next-to-leading order (NLO) in Ref. [9], respectively. As an alternative description motivated by the high collision energies of the underlying hard scattering, Ref. [10] has used the so-called color dipole approach as an effective model for the production of colored particles at high energies, in order to compute the production rates for atmospheric charm.
All these predictions, however, are subject to very large theoretical uncertainties. While the results of Ref. [10] are very sensitive to the parameters of the phenomenological model for the color dipole, which are poorly constrained by experimental data, also the standard perturbative QCD predictions for the hadroproduction of charm quarks acquire big uncertainties, of the order of several ten percents, in the kinematical regions of interest for astrophysical applications. In the latter case, these uncertainties are due to estimates of the missing radiative corrections at higher orders, the knowledge of parton distribution functions (PDFs), especially the gluon PDF at small fractions x of the momenta of the colliding protons, as well as the precision on the charm quark mass.
Since the start of the Large Hadron Collider (LHC) and thanks to both theoretical and experimental progress, our understanding of charm-pair hadroproduction at high energies has significantly improved. It is, therefore, the aim of the present paper to provide new predictions for atmospheric charm and its contribution to lepton fluxes, on the basis of standard perturbative QCD, taking into account the most recent developments in this field.
For inclusive charm-pair hadroproduction we use QCD predictions up to nextto-next-to-leading order (NNLO) in order to establish the apparent convergence of the perturbative expansion and the stability under variation of the renormalization and factorization scales, together with recent determinations of PDFs compatible with constraints from LHC measurements. We also investigate the dependence on the renormalization scheme and the value for the charm quark mass and discuss differences between the running mass and the pole mass schemes. A consistency check is performed by a comparison of the theory predictions to available LHC data from ALICE [11], ATLAS [12] and the LHCb [13] experiments obtained in the runs at √ S = 7 and 8 TeV center-of-mass energy, because the data are within the kinematic region of interest for atmospheric charm. The differential distributions for charmed hadron production which are necessary in order to compute the lepton fluxes are obtained in perturbative QCD with a consistent matching between NLO QCD corrections and parton showers, as the respective predictions at NNLO are currently not available. The interface to the PYTHIA event generator [14] accounts for the full effect of parton showers and the hadronization. Our study features a detailed discussion of the different sources of theoretical uncertainties affecting predictions, which are propagated to the computation of the lepton fluxes. In this way, a total uncertainty band for the final predictions of the prompt lepton fluxes is established and compared to previous results in the literature. The effects on the fluxes due to modifications in primary cosmic ray spectra are also shown by making use of the latest spectra available in an analytic form.
Recently, the authors of Ref. [9] have proposed an update of that work in Ref. [15].
Our computation is independent and differs from Ref. [15] because of the up-to-date perturbative QCD results and methods used in the computation of charm and Dhadron production cross-sections and because of the choice and the variation of the input parameters. In summary, this leads to a more comprehensive estimate of the related uncertainties for the prompt lepton fluxes. The paper is organized as follows: In Sec. 2 we present the method, the input, and the tools we have used for the calculation, together with examples of results from its intermediate steps. Sec. 3 contains our predictions for the lepton fluxes along with a discussion of the related uncertainties. In Sec. 4 we sketch the astrophysical implications of these predictions, after comparing them to those used so far by the astrophysical community, and discuss their potential implications for the IceCube experiment. Finally, in Sec. 5 we draw our conclusions, with reference to future theory progress and measurements which could help to decrease the uncertainties on the predictions presented.

Method: cascade equations and their solution
The particle evolution through an air column of depth X in the Earth's atmosphere can be obtained by solving a set of coupled differential equations, so-called cascade equations. Following Ref. [9,16] one has S decay (k → j) + S reg (j → j) .

(2.1)
A dependence on the energy E j is understood in all terms of eq. (2.1), j labels a particle species, λ j,int and λ j,dec its interaction and decay lengths, respectively, while S prod and S decay denote the generation functions for production and decay: Here, φ k (E k , X) is the flux of particle k, σ k is the total inelastic cross-section for the interaction of particle k in the atmosphere, dσ k→j /dE j is the energy distribution of particle j produced by k, Γ j is the total decay width of particle j and dΓ j→l /dE l is the energy distribution of particle l produced by the decay of j. Regeneration functions in eq. (2.1), i.e., S reg (j → j), can be viewed as a particular case of S prod (k → j) when k = j. According to the nature of particle j (nucleon, heavy-hadron, neutrino), some of the terms in eq. (2.1) may be absent.
The right hand sides of eqs. (2.2) and (2.3), defining the so-called Z-moments for the production and decay of particle j, respectively, are obtained after noticing that the X dependence of fluxes approximately factorizes from their E dependence. In this approximation, analytic solutions exist for eq. (2.1) in the limit where the energy of intermediate particles leading to final leptons is either very small or very large with respect to their critical energy, the latter being proportional to the particle mass m and to the inverse of its proper life-time τ 0 . In the vertical direction, E crit = m c 2 h 0 /(cτ 0 ), where h 0 is the vertical depth of the atmosphere, for which an isothermal model is assumed with the density of the atmosphere evolving with depth as ρ(h) = ρ 0 exp(−h/h 0 ).
In fact, the competition between hadron interaction and hadron decay is crucial in determining the final lepton fluxes and E crit represents an approximate energy above which the hadron decay probabilities are suppressed with respect to their interaction probabilities. In particular, one can distinguish between the conventional neutrino flux and the prompt neutrino flux, according to the nature of the intermediate hadrons. The conventional flux is produced by the decays of charged kaons and pions which dominate over their interaction rates at relatively low energies, as the critical energies for these particles are smaller than 1 TeV. On the other hand, for larger energies, the probability of secondary interactions overcomes the probability that these mesons decay, thereby progressively suppressing the flux of neutrinos from them. At energies above 10 5 − 10 6 GeV, neutrinos are thus mainly produced by the decay of other particles. In the framework of the Standard Model, these are, in particular, charmed and bottomed heavy-hadrons, which are characterized by a larger critical energy (E crit > 10 7 GeV) than pions and kaons 1 . These immediately decaying particles (τ ∼ 10 −12 s) give rise to the so called prompt flux, that is the object of study of this paper.
In case of hadrons decaying into leptons, the flux of leptons coming from low energy hadrons, i.e., from hadrons with E E crit , can be approximated by whereas the flux of leptons from hadrons with E E crit is approximated by fluxes φ, all Z-moments and all attenuation lengths Λ. Note that the low energy lepton flux is isotropic, whereas the high energy lepton flux is characterized by an angular dependence f (θ) ∼ 1/ cos(θ) for θ < 60 o , where θ is the angle with respect to the zenith, and by a more complex angular dependence close to the horizon. The solution in the intermediate energy range E ∼ E crit is obtained by the geometric approximation , (2.6) whose quality and validity depend on the particular shapes of φ low h→l and φ high h→l as a function of the lepton energy, see, e.g., Fig. 10 in Sec. 3 below, for an example of this interpolation. In this way one gets the contribution φ h→l to the lepton flux from each hadron species h. Summation over all hadron species finally provides the total lepton flux φ l for each lepton species l, that is φ l = h φ h→l . Alternatively, the system of differential equations in eq. (2.1) can also be solved numerically.

Input: cosmic ray primary spectrum
The primary cosmic ray (CR) spectrum is an important input of our calculation as it enters the solution of the cascade equations (2.4) and (2.5) both implicitly through Z-moments and explicitly.
CR spectra in the upper layer of our atmosphere are very well constrained at low energies by many direct measurements. Balloon-borne and space experiments, like AMS and CREAM, are able to discriminate with high precision the individual elements included in the cosmic ray composition up to energies around E lab = 10 5 GeV [17]. On the other hand, the high energy tail of cosmic rays is subject to significant uncertainties, in particular related to the different possible CR compositions (protons or heavier ions, up to iron) and the CR origin (galactic or extra-galactic). The high energy region is investigated by ground-based experiments, like KASCADE, KASCADE-Grande and the Pierre Auger Observatory, which globally cover the energy range between E lab = 10 6 GeV and up to several 10 11 GeV. In this context it is important to note that the lepton fluxes at a given energy are affected by CR spectra at energies even larger by a factor of order O(100 − 1000), due to the integration over primary energies in the expressions of the generation functions, eqs. (2.2) and (2.3). Therefore, in order to parametrize the uncertainty on our knowledge of cosmic ray spectra at high energies we consider the following possibilities, i.e., we evaluate lepton fluxes separately for each of the following primary cosmic ray spectra, which are available in literature 2 : 1) Power-law spectrum, composed by two parts: This is a reference spectrum used in earlier works on prompt lepton fluxes, and, although recent measurements have shown that it basically overestimates nucleon fluxes at the highest energies, we consider it for reference and comparison with older works [9,10,15].
2) Gaisser 2012 (variant 1 and 2) [18]: The first variant of the Gaisser spectrum, fitting available experimental data of different origin to an analytic expression with a number of parameters, is based on the hypothesis that three populations, one including CR particles accelerated by SuperNova remnants in our galaxy, a second one still of galactic origin but with an higher energy, and a third one of particles accelerated at extra-galactic sources, contribute to the measured CR spectrum. The three populations are characterized by different rigidities 3 , they all include protons and nuclear groups (He, CNO, Mg-Si, Fe) with different spectral indices. The second variant of the Gaisser spectrum provides a special treatment of the third population, which is assumed to be composed of protons only, with large rigidity.
3) Gaisser 2014 (variant 1 and 2) [19,20]: This uses the same functional form as in Gaisser 2012, but with updated parameters for an alternative fit of experimental data. In particular, the first variant of the spectrum involves three populations, two of galactic and one of extra-galactic origin, involving the p, He, C, O, Fe nuclear groups, with different rigidities with respect to the Gaisser 2012 case. The second variant differs from the first one because it includes an additional component from heavier nuclei, plus a fourth population, characterized by extra-galactic protons only, with large rigidity. This affects the ultra-high-energy part of the spectrum and improves the agreement with Auger data on cosmic ray composition at high-energy [21].
The all nucleon spectra corresponding to these cases are shown in Fig. 1. The effect of the different options on the shape of lepton fluxes is extensively discussed in Sec. 3.1.  Figure 1. The all-nucleon primary cosmic ray spectra used as input in this work. See text for more detail.

Input: p-Air total inelastic cross-section
The total inelastic proton-Air cross-section as a function of the laboratory energy is an input in the denominator of the integrand of the generation function S prod in eq. (2.2), and as a consequence, for the Z-moments for heavy-hadron hadroproduction. Several measurements exist for this quantity, performed by different experiments (for a collection of results see Ref. [22] and references therein), together with theoretical predictions, on the basis of phenomenological models. In this paper, for compatibility with previous works, we consider both the analytical formula [23], already used in old estimates of prompt neutrino fluxes (see, e.g., [9]), σ inel p−Air (E) = 290 − 8.7 ln(E/GeV) + 1.14 ln 2 (E/GeV) mb , (2.8) and predictions from the SIBYLL2.1 [24] and the QGSJet0.1c [25] models for hadronic interactions included in the CORSIKA package [26]. Cross-sections corresponding to  The p-Air total inelastic cross-section according to different models (QGSJet0.1c, SIBYLL2.1, analytical) as compared to measurements from Auger and older experiments. The experimental data are taken from Ref. [22] and references therein.
those possible options are shown in Fig. 2 together with presently available experimental data. We also point out that predictions from other CORSIKA models, like EPOS 1.99 [27] lie within the band that one can draw from these two choices, as shown in Fig. 2 of Ref. [22], so that we consider them as upper and lower limits. The recent measurement from the Auger collaboration at √ s = 57 TeV, reported in Ref. [22] turns out to be in agreement, within the error bands, with the predictions from QGSJet0.1c.
Discussions on the effects of the different assumptions for the p-Air cross-section on our final results of lepton fluxes are reported in Sec. 3.2.

Charm hadroproduction cross-section
Heavy-quark hadroproduction has been extensively studied in perturbative QCD. The QCD corrections at NLO have first been obtained in Refs. [28][29][30] and are available in public tools, like hvqmnr [31], MCFM [32] or HELAC-NLO [33] for the automatic computation of fully differential observables. For the inclusive cross-section, the QCD corrections are complete to NNLO [34][35][36][37] and, thus far, have been applied to top-quark pair production. All these theory predictions have adopted the on-shell renormalization scheme for the heavy-quark mass. The conversion to the MS scheme for the heavy-quark mass has been discussed in Refs. [38][39][40].
Beyond the perturbative expansion at fixed order, the resummation of large logarithms features important improvements, cf., e.g., the review in Ref. [41] for charm and bottom production at the LHC. In dynamical regimes where the transverse momentum p T of the heavy quark is much larger than its mass m the semi-analytical resummmation of logarithms in p T /m has been performed in Ref. [42] in the socalled FONLL approach. On the other hand, when the NLO corrections are consistently matched with parton showers (PS) as in the POWHEG [43,44] or MC@NLO [45] approaches using the frameworks of POWHEG-BOX [46] or (a)MC@NLO [47], respectively, the leading logarithms are effectively resummed through the Monte Carlo (NLO + PS) event generators.
In summary, there exists a robust theoretical framework with a set of well tested tools for the computation of top, bottom and charm-pair hadroproduction at high energies, which has been developed for and used extensively in the LHC environment. For estimating lepton fluxes from atmospheric charm, the core of our calculation is an updated estimate of the N N → cc production cross-section in perturbative QCD, where N indicates a nucleon 4 . We use the QCD predictions for the inclusive crosssection pp → cc at NNLO, and their comparison to those at NLO, to study the stability of the perturbative expansion as a guide for fixing parameters and inputs to be adopted in the application to atmospheric charm. In detail, this includes the PDF dependence, the choice of the central values for renormalization and factorization scales µ R and µ F and the charm mass, as well as plausible intervals for their variations, given the fact that the global uncertainty bands at NLO due to scale, PDF and charm mass variation are large. Our findings are summarized in the sequel. Fig. 3 displays the dependence of the total cross-section σ pp→cc on the laboratory energy E lab . The computation is performed in the fixed flavor number scheme (FFNS) with the number of flavors n f = 3, implying that charm is considered as a heavy state consistently included in the matrix elements with its mass different from zero and its presence excluded from the PDFs. The computation is performed in the theoretical framework as implemented in the HATHOR code [39]. Fig. 3 applies two different schemes for the heavy quark mass renormalization, the commonly chosen on-shell scheme with the pole mass m pole c and the MS scheme with the running mass m c (µ R ), where the renormalization scale for the evaluation of the mass has been fixed at µ R = m c .
The experimental data in Fig. 3 for the fixed target experiments with energies up to E lab = 10 3 GeV are taken from Ref. [48] and for HERA-B from Ref. [49] (purple points in Fig. 3). RHIC data from PHENIX and STAR have been published in Refs. [50,51] (black points in Fig. 3) and LHC data are available from ALICE [11], ATLAS [12] and LHCb [13] (blue points in Fig. 3), see also Ref. [15].  strates the stability of the perturbative expansion of the σ pp→cc cross-section through NNLO up to very high energies and good consistency of the predictions with the experimental data.    Figure 6. Dependence of the total cross-section for pp → cc on the PDF choice at LO (dotted), NLO (dashed), NNLO (solid) QCD accuracy in the MS mass scheme. The charm mass and the scales were fixed as in Fig. 3. The upper and the lower lines at NLO and NNLO indicate the total 1σ PDF uncertainty band for ABM11 (left) and NNPDF3.0 PDF set (right) with n f = 3. Experimental data are the same as in Fig. 3.
Related to the heavy quark mass renormalization is the choice of the numerical value for the charm quark mass. The Particle Data Group (PDG) [52] reports a very precise value of m c (m c ) = 1.275 ± 0.025 GeV in the MS scheme. In case of charm, the conversion of the MS to the pole mass suffers from well-known convergence problems, see, e.g., Ref. [53]. In addition, the definition of the pole mass is based on the unphysical idea of quarks as asymptotic states of the S-matrix, so the accuracy in the pole mass is limited to be of the order O(Λ QCD ) by the renormalon ambiguity. The comparison of the σ pp→cc cross-sections in the two mass renormalization schemes at the nominal scales µ R = µ F equal to 2m pole c and 2m c (m c ) and taking the result with the precise PDG value as a reference, motivates our choice for the charm pole mass m pole c = 1.40 ± 0.15 GeV, as illustrated in Fig. 4. The behavior of the total cross-section as a function of the factorization scale µ F is further explored in Fig. 5 by considering three different renormalization scales, µ R = µ F /2, µ F and 2µ F , at LO, NLO and NNLO QCD. Fig. 5 demonstrates, that the choice of the mass renormalization scheme is important, because the MS scheme leads to predictions with slightly improved convergence. Scale stability of the perturbative expansion at NNLO is reached in both schemes for scales µ R ∼ µ F > ∼ 2 GeV. As shown in Fig. 5 the use of the running mass scheme leads to a somewhat reduced scale uncertainty band at NNLO for the independent variation of µ R and µ F in the standard range µ R /m c (m c ) and µ F /m c (m c ) ∈ [1/2, 2] and restricting the ratio 1/2 ≤ µ R /µ F ≤ 2, as compared to the pole mass scheme. Similar features, although much more pronounced, have been found already for the tt hadroproduction crosssections and differential distributions in Ref. [40]. Both in the pole and in the running mass scheme the point of minimal sensitivity, i.e., the region where the cross-sections predictions at NLO and NNLO approximately coincide, turned out to be around scales µ R = µ F ∼ 2m c and larger. This justifies scale choice adopted in Fig. 3. Translating this value into a dynamical scale, more suitable to describe the dynamics of heavy-quarks in differential distributions, we will use in the following a dynamical central scale for our calculation fixed to µ F = µ R = p 2 T + 4m 2 c , where p T is the transverse momentum of the emitted charm quark.
Another important factor contributing to the theoretical uncertainties of the cc hadroproduction cross-section originates from the choice of the PDF set. We have taken predictions at NNLO accuracy as the basis of our central PDF choice. Among the different possibilities, currently available in the LHAPDF interface [54], we have chosen the ABM11 one [55], together with the respective value for the strong coupling constant α s (M Z ), as a default. In the FFNS with n f = 3, this PDF features a central set complemented by 28 variations, allowing to estimate a PDF uncertainty band at the 1σ level. The ABM11 PDFs are compatible with ABM12 [56], where the latter set has been tuned to LHC data. Moreover, the predictions of the ABM11 and ABM12 family for gluon PDF at low Bjorken-x values are in complete compatibility with the only PDF fit available in literature so far including LHCb data on cc and bb hadroproduction, that has recently been performed by the PROSA collaboration [57,58] As an alternative, we have used the 3-flavor central PDF set at NLO available from CT10 [59], which also provides results characterized by partial compatibility with the PROSA fit (differences lie within about 2σ) and with ABM11. At NNLO, both the CT10 and the ABM11 PDF sets give positive results for the total cross-section for cc hadroproduction in the highest energy range ranging up to E lab ∼ 10 10 GeV, together with an uncertainty band for the PDFs which always stays positive as well, see Fig. 6 (left) for predictions from ABM11. On the other hand, PDF sets with unconstrained gluons at small x lead to very different results. In particular, the NNPDF3.0 set [60], characterized by a different parameterization, leads in case of the highest energies, to a huge uncertainty band, even covering a range with negative cross-section values. As shown in Fig. 6 (right), the cross-sections obtained with NNPDF3.0 at NNLO do not remain positive anymore already for E lab > ∼ 5 · 10 7 GeV, an energy well below the one so far covered in run 1 by the LHC 6 . Thus, in the remainder of this paper we only consider the central value of the NNPDF3.0 set, with the purpose of quantifying differences with respect to the central values of the other families.
The PDFs from MMHT [62], and in particular their central fit, lead to negative cc hadroproduction cross-sections at NNLO for energies above E lab > ∼ 5 × 10 8 GeV.
Such an unphysical feature also affects predictions for the longitudinal structure function F L in deep inelastic scattering as noted, e.g., in Ref. [63]. While the MMHT set seems to be valid for the description of the production of heavier particles at LHC energies, its extrapolation to higher energies characterizing several astroparticle physics problems is quite questionable. Thus, we neglect this PDF family in the present study.
Finally, the small x region is not only of importance for the gluon PDF, but also when considering the behavior of the perturbative hard parton scattering crosssection. The high-energy factorization of the cross-section [64,65] in the limit when the center-of-mass energy S is much larger than the heavy-quark mass provides an effective theory for the description of the high-energy logarithms in S/m 2 . These behave as ln 0 (S/m 2 ) const. at NLO, as ln 1 (S/m 2 ) at NNLO, and so on, see, e.g., Ref. [66] and studies of operator matrix elements in deep-inelastic scattering at three loop order in the small-x limit [67,68]. At the energies currently considered, even up to E lab ∼ 10 10 GeV, their numerical importance is, however, strongly suppressed in the convolution integral of the hard partonic cross section at small x with the large x part of the gluon PDF (and vice versa). The apparent convergence of perturbative expansion for the σ pp→cc cross-section through NNLO observed in Fig. 3 underpins this fact.
The cc differential distributions which we use in this paper are at parton level exact to NLO in QCD, because differential predictions for cc hadroproduction are not yet available at NNLO. For generating these distributions we use the POWHEG-BOX [46], complemented by the event generator PYTHIA-6.4.28 [14], in a p T -ordered tune belonging to the family of Perugia tunes [69], for describing parton shower and hadronization. This provides us with differential distributions of D-hadrons at NLO accuracy in QCD with NLO matching to parton showers according to the POWHEG method. Beyond the resummation to leading logarithmic accuracy provided by the parton showers approaches, next-to-leading logarithmic (NLL) contributions of p T /m, as obtainable by an approach like FONLL, are not considered here. This is justified because the computation of Z-moments requires an integration over the whole kinematically accessible range in p T and thus exhibits less sensitivity to the shape of the p T distribution at large p T , which is mainly influenced by the NLL corrections provided by the FONLL approach.
Moreover, following Ref. [9], we derive the total charm cross-section in the inelastic proton-Air (pA) collisions, from the pN cross-section, by using the formula where A = 14.5 for a nucleus of air. Here, we take γ = 1 assuming a linear superposition since, for light nuclei, the effect of nuclear shadowing is expected to be small [15,70].

Z-moments: Z
The inclusive cc cross-section is an essential ingredient for the estimate of the charm contribution to lepton fluxes. The latter ones, in fact, depend on the Z production moments which can be expressed as integrals over the differential distribution dσ p→charm (E/x E )/dx E , through the formula, with the ratio x E = E/E k . Here, E k is the nucleon energy in the laboratory frame and E the energy of the produced particle (charm). The primary CR nucleon flux is φ p and we have assumed that charm is all produced in cc pairs from primary CR nucleons (denoted by p) interacting with the Earth atmosphere, i.e., and that σ p (E) in eq. (2.10) coincides with the total inelastic proton-Air crosssection σ inel p−Air (E). The lower integration limit would ideally correspond to the case of E k → + ∞. We thus replace it with , as we compute the cross-sections for E k limited to a finite value, with small enough that the results for Z-moments are almost independent of its variations var < . As discussed in Sec. 2.3, for the differential cross-sections dσ pA→cc /dx E and dσ pp→cc /dx E no predictions at NNLO are available at present. Thus we compute it at NLO through POWHEG-BOX, using PDFs and parameters such as m c and scale choices as described above, so that we are close to the point of minimal sensitivity for the total cross-sections, where the differences between NLO and NNLO QCD predictions are small.
The hadronic moments Z p h with h = D 0 ,D 0 , D ± , D ± s and Λ ± c were calculated in Ref. [9] from the partonic moment Z p,charm , through the relation Z p h = f charm,h Z p,charm , including conversion factors f charm,h which represent the fraction of charm quarks emerging as specific hadrons after fragmentation. In Ref. [15] a more refined approach was used: differential distributions for hadrons were obtained from those for quarks after convoluting the latter ones with fragmentation functions that were assumed to depend on the energies through the ratio E h /E charm and to be independent of the beam energy. On the other hand, the use of POWHEG-BOX allows us to follow a different path, i.e., we take into account parton shower and fragmentation effects by means of the Monte Carlo PYTHIA event generator applied to the Les Houches events at first radiation emission level obtained by running POWHEG-BOX. This allows us to directly extract differential distributions dσ pp→h+X /dx E for D-hadrons (x E = E h /E k ), whose shape may, in general, differ from those of the charm quarks dσ pp→c+X /dx E (x E = E charm /E k ), implying that a global rescaling factor is too naive an approximation for the translation of quark distributions at parton level into the corresponding ones at hadron level.
As an example, the differential distribution dσ pp→h+X /dx E is shown for the case of D 0 hadrons in Fig. 7 for a pp → cc collision characterized by E lab = 10 7 GeV. The uncertainties due to scale and mass variation around the central predictions are shown in the left and right panels of the figure, respectively. They are, in general, large. The scale uncertainties turn out to be almost constant with x E , whereas the uncertainties due to the variation of the mass increase with increasing x E on kinematical grounds. Differential distributions for different E lab show a qualitatively similar shape, a scaling property already pointed out in the literature.
Uncertainties due to the PDF variation are shown in the left panel of Fig. 8. For illustration, the predictions for the central fit and the 28 additional variations which parametrize the PDF and α s uncertainties of the ABM11 PDFs with n f = 3 in the FFNS, are displayed individually 7 . and PDFs) and mass variation (at fixed µ R , µ F and PDFs) are shown in the left and right panels, respectively. The lower subpanels display the band for the relative uncertainties when normalized with respect to the central prediction.
As is evident from Fig. 8, the differences with the central values of other PDF sets (CT10 and NNPDF3.0 at NLO in the n f = 3 FFNS with their respective default value for α s (M Z )) turn out to be larger than the combined PDF and α s uncertainty coming from the 28+1 ABM individual sets. This reflects a feature already observed at NLO for several examples of LHC cross-sections, where differences arising from the use of different PDF families turned out to be larger than those from the variation of α S and the PDFs through the sets belonging to a same family, see, e.g., Refs. [55,71].
These differences propagate to the computation of the Z-moments, although the latter quantities are integrals over all possible x E values. As an example, in the right panel of Fig. 8, the Z-moments for D 0 hadroproduction are shown as a function of the D 0 energy in the laboratory frame E lab, D 0 for the different PDF choices just discussed above (left panel, Fig. 8). Here, the power-law spectrum as been chosen as input for the CR flux. Thus, the change of shape visible around 5 ·10 6 GeV is due to the change of the spectral index in the power-law spectrum around the knee. The largest differences between different PDF sets appear at the lowest and at the highest D 0 energies.
It is worth noting that pp collisions with E lab > E lab,D 0 contribute to the Zand scale in the whole integrand, in the separate generation of each set of events. moment at any given fixed energy E lab,D 0 . Although, in line of principle, E lab can be very large, in practice it turns out that the largest contribution comes from values of E lab within the range E lab,D 0 < E lab < (100 − 1000) × E lab,D 0 , due to the fact that the distribution in x E = E D 0 /E lab is rapidly suppressed for large x E . As a consequence, for energies as those probed by IceCube, the contributions to the Z-moments come mainly from regions with a center-of-mass energy √ S non too high with respect to the energy range reached and probed so far at the LHC, where perturbative QCD in the standard formalism of collinear factorization has been tested to work. Any deviations from this formalism which may exist at the highest energies, e.g., in the form of non-linear effects (like gluon recombination as opposed to gluon splittings) or due to the dominance of large logarithms ln(S/m 2 ) subject to resummation on the basis of a different factorization formalism (k T factorization) [65], are, thus, expected to have only a small impact on the Z-moments we are interested in for the aim of understanding the IceCube results.
The relative importance of Z-moments of different D-hadron species is shown in Fig. 9, where the Z-moments of positively charged D-hadrons and D 0 are shown. The D 0 contribution is the dominant one at all hadron energies. The different shape of the Λ + c contribution with respect to those of other hadrons is partly related to the shape of the differential distribution dσ/dx E of this hadron, which turns out to be less steep at large x E than those of the D-meson distributions.
The expressions for Z p h enter directly into those for the fluxes eqs. (2.4) and (2.5), obtained after solving the system of coupled differential equations describing the linear development of the hadronic cascade in the atmosphere, under the approximations outlined in Sec. 2.

Z h l , Z p p and Z h h
In the following we briefly summarize our treatment of the other Z-moments Z h l , Z p p and Z h h entering eqs. (2.4) and/or (2.5).
For Z h l , our treatment of the semileptonic decay of D-hadrons follows closely Ref. [15]. Form factors for analytical decay distributions h → µν µ X were extracted from Ref. [7] and for the decay branching ratios the most recent values reported by the PDG [52] were taken.
In order to evaluate the proton regeneration Z-moment, Z p p , we have approximated the inelastic x E distribution for the process pA → pX by a scaling form dσ/dx E ∼ σ inel pA (E lab )(1 − x E ) n (1 + n) with n = 0.51, as already done in Ref. [15]. Here, for σ inel pA we have considered the three different models already described in Sec. 2.2, which also enter the generation of the production moments Z p h .
Due to the obvious difficulties in measuring hA cross-sections with h being a D or B hadron, caused by the short lifetime of these particles, the moments Z h h are approximated by considering the available estimates for KA cross-sections, on the basis of analogies between K and D mesons. Both include quarks belonging to the same flavor family (charm quarks in case of D's are replaced by strange quarks in case of K's). In particular, following Ref. [15], the inelastic x E distribution for the reaction hA → hX is estimated by the ansatz dσ/dx E ∼ A 0.75 σ inel KN (E lab ) (1 − x E ) n (1 + n), where n = 1 and σ inel K ± N is the total inelastic cross-section for K ± -nucleon interactions. To estimate the latter one, total and elastic cross-sections were extracted from the latest version of the PDG. However, the behavior of the K ± p elastic cross-section at high energies is uncertain because no data above E lab ∼ 300 GeV exist. We have thus assumed that the slope of the K ± p elastic cross-section at high energies is similar to the one of the pp elastic cross-section, which was recently constrained at LHC energies by TOTEM data [72].
The regeneration Z-moments Z p p and Z h h enter the expression for the attenuation lengths Λ p and Λ h , respectively, as defined below eq. (2.5). It is worth noting that estimates of the Z hh -moments and of the related uncertainties only affect the high-energy approximated solution of the cascade equations, eq. (2.5).

Neutrino fluxes and their uncertainties
The IceCube experiment is looking for a diffuse flux of neutrinos at high-energies, including both downward and upward going neutrinos, by trying to establish the nature of the observed events as due to astrophysical signals or to the atmospheric (conventional + prompt) background. In this Section we focus on (ν µ +ν µ )-fluxes, by taking into account that the largest contribution to the atmospheric conventional neutrino flux (to which we will compare our prompt flux) comes from this flavor. Predictions for other leptons can be obtained with the same method as well. The qualitative/quantitative difference between the results for (ν µ +ν µ ) fluxes and those for other leptons depends on the specific decay modes and branching fractions of D hadrons in each species.
The lepton fluxes are derived after evaluating all quantities entering eqs. (2.4) and (2.5), already described in previous Sections, and by interpolating between the high energy and the low energy solutions according to eq. (2.6). An example of the typical behavior of the two solutions and of their interpolation is shown in Fig. 10 for the case of a power-law primary CR spectrum as input of the whole calculation.  In the following we will present the central values of our fluxes, together with the uncertainty bands arising from the different source of uncertainties, both of QCD and of astrophysical origin.

Main uncertainties from QCD and astrophysics
Uncertainties on the fluxes whose origin can be ascribed to perturbative QCD mainly reflect those uncertainties already found in the differential distributions dσ/dx E and in the Z-moments. In particular, we discuss in the following the scale, charm mass and PDF variation, as well as matching uncertainties, related to the NLO matching to the parton shower.
Uncertainties in the (ν µ +ν µ )-fluxes due to µ R and µ F scale variation for a power-law CR spectrum as input, are reported in Fig. 11, while the corresponding uncertainties for other CR spectra are shown in Fig. 12. The scale variation turns out to be the largest source of uncertainties. Including cases with µ R = µ F , i.e., the independent variation of µ R and µ F , leads to an uncertainty band which is almost uniform on a wide interval of energies E lab . In this respect our findings in Figs. 11  Figure 11. The (ν µ +ν µ )-fluxes as a function of the neutrino energy E lab, ν with uncertainties due to renormalization and factorization scale variation. Charm mass, PDF and (µ R ,µ F ) scales were fixed as in the left panel of Fig. 7. and 12 are different from the result of Ref. [15], where the non-diagonal choices with µ R = µ F were neglected and the scale uncertainty is underestimated, especially at low energies. Uncertainties arising from the variation of the charm mass m c within the range motivated in Sec. 2.3 are illustrated in Fig. 13 for a power-law CR spectrum and in Fig. 14 for other CR spectra. The mass variation turns out to be the second largest source of QCD uncertainties, with an uncertainty band slightly decreasing for increasing laboratory energies E lab .
Uncertainties in the neutrino fluxes related to the PDF variation are displayed in Figs. 15 and 16 in case of a power-law primary CR spectrum and for the different variants of Gaisser spectra, respectively. As already discussed for the case of the Z p h -moments, the difference of the predictions with the ABM11 set (3-flavor FFNS) and the central value of other PDF families (CT10 and NNPDF3.0 at NLO with n f = 3) turn out to be larger than those coming from the 28 sets in the ABM11 fit for the combined PDF and α s uncertainty. While the neutrino fluxes from the different PDF fits look quite consistent among each other for energies in the interval 10 2 < E lab, ν < 4 · 10 4 GeV, visible differences between the fluxes from different PDF families start to appear at higher energies. In this region, the predictions based on the ABM11 PDFs are the smallest ones, which is related to differences in the shape of the gluon PDF, the nominal values of the strong coupling α s (M Z ) at NLO being largely the same among the ABM11, CT10 and NNPDF3.0 sets used in this work.
Finally, we provide a first estimate of the uncertainties in the NLO matching to the parton shower (NLO + PS) by varying the h damp value in the POWHEG-BOX, which parameterizes the freedom in choosing the form of the separation of the NLO real contribution R into a singular piece plus a piece damped in the singular region and thus treatable as a finite remainder [73], [74]. Only R s enters the exponent of the Sudakov form factor and h damp = +∞ corresponds to the default choice in POWHEG-BOX so that R = R s , whereas the limit h damp → 0 allows to decrease the amount of radiation that is exponentiated and to recover the α 3 s dependence (pure NLO) in the high-p T limit. In this work we use variations of h damp in the interval {m c , 2m c , 4m c , +∞}. This is inspired by similar choices performed in experimental studies of tt hadroproduction,  Figure 13. The (ν µ +ν µ )-fluxes as a function of the neutrino energy E lab, ν with uncertainties due to the variation of the pole mass m pole c = 1.40 ± 0.15. PDF, (µ R ,µ F ) scales, charm mass were fixed as in the right panel of Fig. 7. see, e.g., the ATLAS note [75]. The uncertainty is estimated as the envelope of the predictions corresponding to the different choices above, as shown in Fig. 17 in case of a power-law primary CR spectrum and in Fig. 18 in case of the Gaisser spectra. Although several discussions are still on-going about the most meaningful way of providing NLO + PS matching uncertainties which is why we consider our result as a first rough estimate, we would like to point out that the uncertainty we got is quite small (less than 10%) with respect to other uncertainties of QCD origin. This is related to the fact that the key quantities in perturbative QCD to compute Zmoments, the differential cross sections dσ/dx E , are integrated over the entire range of transverse momenta p T . We thus believe that this conclusion is quite robust, i.e., it does not depend on very specific details of the way the matching uncertainty is estimated.
A summary of the main QCD uncertainties (mass, scale, PDF variation) in relation to uncertainties of astrophysical origin, in particular those arising from the variations of the primary CR flux used as input in the (ν µ +ν µ )-flux calculation, is provided in Fig. 19. In the first three panels of this figure we show separately the uncertainties due to scale, mass and PDF variation (by restricting ourselves to the ABM11 PDF set) for all five primary CR fluxes considered as input in this paper.  The uncertainties due to scale variation are the dominant component. Apart from that, it is evident that at energies > ∼ 10 6 GeV uncertainties due to variations in the CR fluxes dominate over those from mass and PDFs. On the other hand, uncertainties related to QCD effects always dominate at energies < ∼ 10 5 GeV, where the primary CR fluxes are well constrained by several measurements (see Sec. 2.1). Finally, in the last panel of Fig. 19 we show the quadratic combination of the uncertainties above, assumed as independent, i.e., ∆ QCD = ∆ 2 mc + ∆ 2 (µ R ,µ F ) + ∆ 2 P DF . For E lab, ν = 10 6 GeV, −72% ≤ ∆ QCD ≤ +84%, i.e., the uncertainty is slightly asymmetric, and it slightly changes (a few percent) at higher energies.
We also observe that with a restricted scale variation interval, neglecting the combinations (µ R , µ F ) = (0.5, 1) and (1, 0.5) µ 0 for µ 0 = p 2 T,c + 4m 2 c , scale uncertainties and, as a consequence, the total ones, are reduced, as shown in Fig. 20. In particular, for E lab, ν = 10 6 GeV, the combined uncertainty amounts to −48% ≤ ∆ QCD ≤ +63%, and it changes by a few percent at higher energies.  Figure 15. The (ν µ +ν µ )-fluxes as a function of the neutrino energy E lab, ν with uncertainties due to PDF variation in the 3-flavor ABM11 PDF set at NLO (red band) and predictions for the central set of the 3-flavor NNPDF3.0 (light-blue line) and CT10 (solid blue line) PDFs at NLO, respectively. Charm mass and (µ R ,µ F ) scales were fixed as in Fig. 8. The power-law cosmic ray flux has been used as input in the calculation of Z-moments.

Other uncertainties
In the previous Sec. 3.1 we have provided a minimal estimate of the combined QCD and astrophysical uncertainties which affect our predictions for (ν µ +ν µ )-flux. In the following we shortly describe other sources of uncertainties which could be added to the previous ones. A further QCD contribution arises from heavier hadrons, in particular B-hadrons, which are also a source of neutrinos, but whose effect has been neglected here. Given, that the cross-section for bb hadroproduction with respect to the one for cc hadroproduction is smaller by a factor of order 20 at LHC energies and still suppressed by a factor of order 10 at E lab = 100 TeV, we expect that the bottom-quark contribution can be neglected with respect to the charm one at the energies of interest for IceCube. However, bb hadroproduction may play a larger role at ultra-high-energies.
Other uncertainties can be attributed to the approximate description of the decay of heavy hadrons. In particular, a component of secondary neutrinos coming from the decay of the lighter mesons (baryons) produced as decay products of Dmesons (baryons), is missing in our computation as well as in many previous ones (see e.g. Ref. [15]). Furthermore, from the QCD point of view, non-perturbative effects, suppressed by powers of Λ QCD /m, are increasingly important for smaller quark masses m. In this respect, one should account for a contribution to the uncertainties due to fragmentation, arising from both fragmentation fractions and fragmentation functions [76]. The latter can be estimated, for instance, by varying the choice of the functional form of fragmentation functions for heavy flavors in PYTHIA together with the parameters involved. The corresponding uncertainty estimates have been discussed in the literature, see, e.g., Ref. [41]. Potential uncertainties related to the variation of the partonic intrinsic transverse momentum k T ∼ Λ QCD are indeed smaller since they mostly affect the p T distributions to which our work is not particularly sensitive.
From an astrophysical point-of-view, further uncertainties to be included encompass those related to a change in the p-Air cross-section, which affect both the Z p h production moments and the Z p p regeneration moments, and, as a consequence, the lepton fluxes. To that end, we consider the theoretical predictions coming from the three different models described in Sec. 2.2 (QGSJet0.1c, SYBILL2.1 and the analytical model eq. (2.8)). The larger the inelastic cross section σ inel (p-Air) is, the smaller are the predictions for our fluxes, as is evident when comparing Fig. 21 with Fig. 2. However, these global effects on neutrino fluxes turn out to be not too relevant, i.e., the uncertainties coming from the use of different models, amount to less than 10% over the whole E lab, ν energy range considered. They are therefore much smaller than those from QCD effects discussed previously and those from the choice of different primary CR spectra, as shown in Fig. 22. energies of interest for the IceCube experiment) than the central values in Ref. [15], that lie in any case within our quoted uncertainty band, for the various primary CR spectra already considered in that work, as shown in Figs. 23 and 24. Note, that Ref. [15] is based on a completely independent QCD calculation and on different inputs and methods.
On the other hand, differences with older calculations, like the one in Ref. [8] on the basis of PYTHIA and including QCD hard-scattering effects at leading order only, are obvious, especially regarding the shape of the distributions. This is the case not only for the lepton fluxes, but already for the Z-moments for D-hadron production.
Interestingly, our results are very well compatible with those from the dipole model of Ref. [10], altough the latter were computed with older sets of PDFs: as shown in Fig. 23, for E lab,ν > 10 3 GeV, our central predictions are included in the uncertainty band of the latter, whereas the central predictions from the dipole model are included in our uncertainty band.  Figure 19. Summary of the main QCD and astrophysical uncertainties affecting our central predictions for (ν µ +ν µ )-fluxes. Uncertainties due to scale, mass and PDF variation (considering the ABM11 PDF and α s uncertainty band), are shown separately and combined for each of the five primary CR spectra, cf. Sec. 2.1.
In order to infer a value for the transition energy E trans where the prompt neutrino flux overcomes the conventional one, in Fig. 24 we compare our prompt lepton flux with the conventional neutrino flux originally computed in Ref. [77] for a powerlaw CR primary spectrum and rescaled to one variant of the Gaisser spectra in Ref. [15]. We obtain E trans = 6.0 +12 −3 ·10 5 GeV. Interestingly, the central value lies well within the interval (4 ·10 5 − 10 6 ) GeV where the IceCube experiment did not observe any event after the full 988-day analysis [3]. In fact, the IceCube collaboration has reported an excess of neutrinos in the diffuse flux, all lying in the neutrino energy regions [0.3 − 4] 10 5 GeV and [1 − 2] 10 6 GeV. According to our predictions, the "empty" region of IceCube corresponds to the "conventional-prompt" transition region, i.e., the region where the contributions of conventional neutrinos and prompt neutrinos to the total neutrino flux become of the same order of magnitude. We thus believe that the "empty" region seen by IceCube so far, should not be empty, but actually dominated by prompt neutrinos. However, the IceCube error bars in the "empty" region are still quite large, and we stress that the accumulation of more statistics is necessary before judgment can be made, whether this lack of signal is just an artifact due to poor statistics or due to some other technical issue, or instead has a real physical interpretation.
At higher energies, on the other hand, the total observed neutrino flux E 2 ν φ(E ν ) for E ν in the [1 -2] 10 6 GeV energy interval looks to be slightly suppressed with respect to that in the [2 -3] 10 5 GeV bins. However, looking at our central prompt flux distributions and summing them with the distributions for the conventional flux, as a first rough estimate it turns out that we would expect a much larger suppression in the [1 -2] 10 6 GeV region with respect to the [2 -3] 10 5 GeV one, disfavoring the interpretation that the events seen by IceCube in the [1 -2] 10 6 GeV window are just due to a prompt neutrino component 8 . The difference between the IceCube (signal + background) observed total yield at high-energy and the yield for prompt neutrinos as predicted by our calculation is slightly reduced if we observe that our predictions have a sizable uncertainty band, meaning that even the shapes of the distributions can change in a non-negligible way when a higher-order calculation in QCD will be available, if we consider neutrino flux values corresponding to the upper value of our uncertainty band and if we use as input primary CR fluxes including a population of extra-galactic protons with very-high rigidity (i.e. variants 2 of Gaisser spectra, instead of variants 1 which have a mixed extra-galactic component with a lower global rigidity for the extra-galactic population). In fact, the latter give rise to neutrino spectra which are less severely suppressed at the highest energies than those from models with extragalactic mixed components, as is evident when comparing, e.g., the left and right panels of Fig. 24, obtained with the variants 1 and 2 of the Gaisser 2014 spectrum, respectively. In order to go beyond these purely qualitative considerations and to draw more definite quantitative conclusions, one should definitely wait for more experimental statistics and, also insert our fluxes into the specific experimental analysis software.
In any case, we would like to emphasize that the transition region for the prompt (ν µ +ν µ )-flux in our calculation turns out to be also a transition region for uncertainties, i.e., the QCD uncertainties dominate the total uncertainties at energies below the transition region whereas the astrophysical ones start to give a progressively sizable contribution above it, pointing out the importance and necessity of pursuing further studies of cosmic ray composition at the highest energies [78], and, possibly, future measurements independent of Monte Carlo simulations of hadronic-interactions at the highest energies.

Conclusions
We have computed the prompt neutrino fluxes from atmospheric charm using up-todate theoretical results and tools for charm hadroproduction in perturbative QCD. Our results for the neutrino fluxes are several tens percent larger over a wide range of neutrino energies than predictions in the recent literature making use of Z-moments computed with the standard QCD hard-scattering formalism, that we also adopt. At the energies of interest for the IceCube experiment the increase of our prompt (ν µ +ν µ )-flux amounts to some 40%. However, our uncertainties on the fluxes both of QCD and astrophysical origin are dramatically larger. Partly as an effect of this fact, even predictions obtained by making use of the dipole picture, representing an alternative description to the undelying hard-scattering, lie within our uncertainty band over a wide energy range. is shown by open magenta squares, the ERS central flux (Ref. [10]) and its uncertainty is shown in yellow, whereas the more recent BERSS flux (Ref. [15]) is shown by filled light-blue squares.
We have discussed extensively the different sources of uncertainties which affect the fluxes. The main sources come from (i) the renormalization and factorization scale variation allowing for independent variations of µ R = µ F , (ii) the charm mass uncertainties for the pole mass choice, and (iii) PDF uncertainties evaluated for the ABM11 set and studied by comparing its predictions to the central predictions of different PDF sets (CT10, ABM11, NNPDF3.0) at NLO. Further uncertainties due to hadronization and hadron decay have been discussed as well. In particular (i) and (ii) had not been included in a systematic way in studies in literature before, so we conclude that previous uncertainties on prompt neutrino fluxes are underestimated.
The uncertainties of QCD origin dominate at low neutrino energies, whereas for increasing energies E lab, ν > ∼ 10 5 − 10 6 GeV the uncertainties in the astrophysical input, in particular the primary CR flux and its composition in terms of different populations, turn out to add a progressively important contribution to those from QCD.
The results presented may benefit from a number of future developments. On the QCD side, a fully differential NNLO computation of charm hadroproduction, when available, will be of great help in reducing the theoretical uncertainties from scale, mass and PDF variation. In this respect, the role of resummation of different kinds of logarithms deserves further exploration as well. Furthermore, a dedicated systematic survey of the uncertainties related to both the fragmentation functions in the Monte Carlo parton shower matched to NLO predictions and the fragmentation fractions, would allow to quantify those effects. This could be a step towards the optimization of Monte Carlo tunes, to make them especially tailored to studies like those performed in this work. This optimization also concerns the search for the best parameter values for the description of dual and multiple particle interactions. On the experimental side, measurements of the cc and bb production cross-section at the LHC, looking not only at central rapidities but also in the forward rapidity regions, can be of importance especially at the highest energies, where the contribution of low x events becomes increasingly important.
Finally, from the astrophysical point of view, one could obtain a substantial reduction of the uncertainties on prompt neutrino fluxes at the highest energies once issues related to the transition between a galactic and an extragalactic component in the CR primary spectrum and to the composition of the latter will be understood better.
Our lepton fluxes will be made available as numerical tables for download at http://www.desy.de/∼promptfluxes. Further predictions can be requested to the authors of this paper by e-mail.