Prompt atmospheric neutrino fluxes: perturbative QCD models and nuclear effects

We evaluate the prompt atmospheric neutrino flux at high energies using three different frameworks for calculating the heavy quark production cross section in QCD: NLO perturbative QCD, $k_T$ factorization including low-$x$ resummation, and the dipole model including parton saturation. We use QCD parameters, the value for the charm quark mass and the range for the factorization and renormalization scales that provide the best description of the total charm cross section measured at fixed target experiments, at RHIC and at LHC. Using these parameters we calculate differential cross sections for charm and bottom production and compare with the latest data on forward charm meson production from LHCb at $7$ TeV and at $13$ TeV, finding good agreement with the data. In addition, we investigate the role of nuclear shadowing by including nuclear parton distribution functions (PDF) for the target air nucleus using two different nuclear PDF schemes. Depending on the scheme used, we find the reduction of the flux due to nuclear effects varies from $10\%$ to $50 \%$ at the highest energies. Finally, we compare our results with the IceCube limit on the prompt neutrino flux, which is already providing valuable information about some of the QCD models.


Introduction
Measurements of high-energy extraterrestrial neutrinos by the IceCube Collaboration [1,2] have heightened interest in other sources of high-energy neutrinos. A background to neutrinos from astrophysical sources are neutrinos produced in high energy cosmic ray interactions with nuclei in the Earth's atmosphere. While pion and kaon production and decay dominate the low energy "conventional" neutrino flux [3][4][5], short-lived charmed hadron decays to neutrinos dominate the "prompt" neutrino flux [6][7][8][9][10][11][12][13][14] at high energies. The precise cross-over energy where the prompt flux dominates the conventional flux depends on the zenith angle and is somewhat obscured by the large uncertainties in the prompt flux. The astrophysical flux appears to dominate the atmospheric flux at an energy of E ν ∼ 1 PeV. Atmospheric neutrinos come from hadronic interactions which occur at much higher energy. With the prompt neutrino carrying about a third of the parent charm energy E c , which in turn carries about 10% of the incident cosmic ray nucleon energy E CR , the relevant center of mass energy for the pN collision that produces E ν ∼ 1 PeV is √ s ∼ 7.5 TeV, making a connection to LHC experiments, e.g., [15,16].
There are multiple approaches to evaluating the prompt neutrino flux. The standard approach is to use NLO perturbative QCD (pQCD) in the collinear approximation with the integrated parton distribution functions (PDFs) and to evaluate the heavy quark pair production which is dominated by the elementary gluon fusion process [17][18][19]. Such calculations were performed in [7][8][9] (see also [6]). Recent work to update these predictions using modern PDFs and models of the incident cosmic ray energy spectrum and composition appears in [11], and including accelerator physics Monte Carlo interfaces, in [12][13][14]. Using x c = E c /E CR ∼ 0.1 for charm production, one can show that high energies require gluon PDF with longitudinal momentum fractions x 1 ∼ x c and x 2 ∼ 4m 2 c /(x c s) x 1 . For a factorization scale M F ∼ 0.5 − 4m c , this leads to large uncertainties. In addition, due to the small x of the gluon PDFs in the target one may need to address the resummation of large logarithms at low x.
In particular, comparisons with LHCb data at 7 TeV [15] were used in ref. [14] to reduce uncertainties in pQCD calculation (see also ref. [20]). Using FONLL [21][22][23][24] predictions for the p T distribution of charm mesons obtained with different PDFs, they have shown that LHCb data for D mesons and B mesons can reduce the theoretical uncertainty due to the choice of scales to 10% and the uncertainty due to the PDF by as much as a factor of 2 at high energies in the region of large rapidity and small p T . Still, the uncertainty due to the low x gluon PDF remains relatively large.
Given the fact that the gluon PDF is probed at very small values of x, it is important to investigate approaches that resum large logarithms ln(1/x) and that can incorporate other novel effects in this regime, such as parton saturation. Such effects are naturally incorporated in the so-called dipole model approaches  and within the k T (or high energy) factorization framework [46][47][48][49].
There is another major source of uncertainty in the low x region. The target air nuclei have an average nucleon number of A = 14.5. Traditionally in the perturbative approach, the nuclear effects are entirely neglected and a linear scaling with A is used for the cross section. Nuclear shadowing effects, however, may be not negligible at very low x and low factorization scale.
In the present paper, we expand our previous work (BERSS) [11] to include nuclear effects in the target and analyze the impact of the low x resummation and saturation effects on the prompt neutrino flux.
We incorporate nuclear effects in the target PDFs by using in our perturbative calculation two different sets of nuclear parton distribution functions: nCTEQ15 [50] and EPS09 [51]. As there is no nuclear data in the relevant energy regime, these nuclear PDFs are largely unconstrained in the low x region (x < 0.01) and there is a substantial uncertainty associated with nuclear effects. Nevertheless, for charm production, the net effect is a suppression of the cross section and the corresponding neutrino flux. At E ν = 10 6 GeV, the central values of the nCTEQ PDF yields a flux as low as ∼ 73% of the flux evaluated with free nucleons in the target, while the corresponding reduction from the EPS09 PDF is at the level of 10%.
We also show our results using the dipole approach, with significant theoretical improvements with respect to our previous work (ERS) [10]. These include models of the dipole cross sections that are updated to include more precise experimental data. Furthermore, we calculate the prompt neutrino flux in the k T factorization approach, using unintegrated gluon distribution functions with low x resummation and also with saturation effects. We compare these calculations to the dipole and NLO pQCD results.
Overall we find that for all calculations, there is a consistent description of the total charm cross section at high energies, for pp and pN production of cc. We also evaluate the bb cross section and the contribution of beauty hadrons to the atmospheric lepton flux. For each approach we find that our choice for theoretical parameters is in agreement with the latest LHCb data [15,16] on charm transverse momentum and rapidity distributions in the forward region, and the total cross sections at 7 TeV and at 13 TeV.
In addition to including nuclear and low x effects, we also consider four different cosmic ray fluxes [52][53][54] and show how the prompt neutrino flux strongly depends on the choice of the primary cosmic ray flux.
The present paper is organized as follows. In the next section we present calculations of the total and differential charm cross section. We present comparisons of all three approaches, pQCD, dipole model and k T factorization, and we show the impact of nuclear effects on the total charm cross sections. We show comparisons of our theoretical results with the rapidity distributions measured at LHCb energies. In sec. 3 we compute neutrino fluxes for muon and tau neutrinos and compare them with the IceCube limit. Finally, in sec. 4 we state our conclusions. Detailed formulas concerning the fragmentation functions and meson decays are collected in the Appendix.

Heavy quark cross sections 2.1 NLO perturbation theory
We start by expanding our recent work on heavy quark cross section with NLO perturbation theory [11] to constrain QCD parameters by comparison with RHIC and LHC data, and by including nuclear effects. In particular, we shall compare the results of the calculation with the latest LHCb data on forward production of charm mesons [15,16]. Gauld et al. [13,14] have evaluated charm forward production to constrain gluon PDFs and to compute the prompt atmospheric lepton flux. Garzelli et al. [12] have recently evaluated the total charm production cross section at NNLO and used the NLO differential charm cross section to evaluate the prompt atmospheric lepton flux. Below, we discuss differences between our approaches to evaluating first charm production, then the prompt fluxes.
We use the HVQ computer program to evaluate the energy distribution of the charm quark at NLO in pQCD [17][18][19]. The resummation of logarithms associated with large transverse momentum p T as incorporated by the FONLL calculation [21][22][23] is not necessary for this application since the low p T kinematic region dominates the cross section.
For heavy quark production, one important parameter is the charm quark mass. In ref. [12], neutrino fluxes were evaluated using NLO QCD on free nucleon targets with m c = 1.40 GeV taken as the central choice of charm quark mass, based on the pole mass value of m c = 1.40 ± 0.15 GeV. Values of m c = 1.5 ± 0.2 GeV are used in refs. [13,14]. In our work, we use the running charm quark mass of m c = 1.27 GeV, which is consistent with the average value quoted in [55], m c (m c ) = 1.275 ± 0.025 GeV. A direct translation between the pole mass and running mass is not possible because of poor convergence of the perturbative series, as discussed in, e.g., ref. [56]. By using m c = 1.27 GeV, we can make use of the data-constrained analysis of the factorization and renormalization scale dependence discussed in ref. [57].
The mass dependence enters through the renormalization and factorization scale dependence as well as through the kinematic threshold. By keeping the values of the factorization and renormalization scales fixed and only varying the charm mass dependence in the matrix element and phase space integration, one can show that there is a strong dependence on mass at low incident beam energies, but at higher energies, the mass dependence is much weaker. For example, keeping the renormalization and factorization scales fixed at M R = M F = 2.8 GeV, the cross section σ(pp → ccX) with m c = 1.27 GeV is a factor of only 1.26-1.16 larger than the cross section with m c = 1.4 GeV for incident proton beam energies of 10 6 − 10 10 GeV. The uncertainties due to the choice of scales are larger than those due to the mass variation. We discuss below the impact of the scale variations on both the cross section and prompt fluxes.
For the NLO pQCD bb contribution to the prompt flux, we use a fixed value of m b = 4.5 GeV and consider the same range of scale factors as for cc production.
In the perturbative calculation of the heavy quark pair production cross section in cosmic ray interactions with air nuclei with A = 14.5, one has to take into account the fact that the nucleons are bound in nuclei, as opposed to free nucleons. Nuclear effects can result in both suppression and enhancement of the nuclear parton distribution functions (nPDF) relative to the free nucleon PDF, depending on the kinematic variables (x, Q). The extraction of nPDF at NLO has been done by several groups, among them Eskola, Paukkunen and Salgado (EPS09) [51] and Kovarik et al. (nCTEQ15) [50]. The nuclear PDFs in the EPS framework [51] are defined by a nuclear modification factor multiplying the free proton PDFs. For example, the up quark PDF for the quark in a proton bound in nucleus A is is the nuclear modification factor to the free proton PDF. For our calculations, we use the central CT14 NLO [58] PDF for free protons and to approximate nitrogen targets, the EPS09 NLO results for oxygen. The recent nCTEQ15 PDF sets [50] instead provide directly the parton distribution functions for partons in protons bound in the nucleus, e.g., u A (x, Q). As usual, one uses isospin symmetry to account for neutrons bound in nuclei. For the calculation in this work we take as our standard free proton PDFs those of [50], labeled here as nCTEQ15-01, and PDFs for nucleons in nitrogen, labeled here as nCTEQ15-14.
In figs. 1 and 2 we show the impact of nuclear modification on the gluon distribution in the small x region using nCTEQ15 and EPS PDFs respectively. In the standard distribution of the nCTEQ15 grids, low x extrapolations must be used to avoid the unphysical behavior shown by the dashed lines in fig. 1. The dotted lines show a power law extrapolation xg(x, Q) ∼ x −λ(Q) below x min = 10 −6.5 . The solid lines in fig. 1 show the nCTEQ15 results with grids extended to low x [59]. The shaded band shows the range of nuclear PDF uncertainties in the 32 sets provided. We use the corresponding lower and upper curves (sets 27 and 28) to quantify the nuclear PDF uncertainty, which is likely underestimating the uncertainty given the lack of data in this kinematic regime. Similarly, for the EPS09 gluon distribution in fig. 2, we also show the uncertainty band, which is now computed as the maximal deviation from the central band due to a combination of uncertainties from the 57 different members of the base proton CT14NLO PDF set and those from the different members of EPS09 modification factors themselves. As a result of incorporating PDF uncertainties from both the proton PDF and nuclear modification factors, the net uncertainty bands at low x in results obtained using the EPS09 scheme are generally larger than those from the nCTEQ15 PDF's. Overall, we find that within the nuclear PDF sets used here, the uncertainty is rather modest, which is due to the constraints stemming from the parametrization. We note that the real uncertainty for the nuclear PDFs can be much larger in the low x region. 1 Depending on the observable, the nuclear effects in the nCTEQ15 and EPS09 frameworks can be sizeable. For the total cross section, the dominant contribution comes from the symmetric configuration of partons' longitudinal momenta, i.e., x 1,2 ∼ 2m c / √ s. On the other hand, the differential distribution in outgoing charm energy fraction x c = E c /E CR , for the forward production is dominated by asymmetric configurations x c ∼ x 1 x 2 , and thus probes deeper into the shadowing region of the target nucleus. We will show below that the impact of shadowing on the total charm cross section is less significant than it is on the neutrino flux, which is dominated by the forward charm production.
In table 1, we show the cross sections σ(pp → ccX) and σ(pA → ccX)/A using the nCTEQ15 PDFs for our central set of factorization and renormalization scale factors For an incident beam energy E p = 10 6 GeV, σ(M F,R ∝ m c ) is larger than σ(M F,R ∝ m T ) by a factor of 1.16 − 1.17, while at 10 8 GeV, the cross sections are nearly equal. These choices of factorization and renormalization scales proportional to m c are the central values constrained by the data in an analysis using NLO pQCD charm cross section calculation in ref. [57] and used in ref. [11]. As noted, we find similar results for the cross sections for scales proportional to m T . Scale variations of (M F , M R ) = (1.25, 1.48)m T and (M F , M R ) = (4.65, 1.71)m T bracket the results of ref. [57], and we use this range here as well. While the total cross section requires extrapolations of the fiducial to inclusive phase space for data comparisons with theory, we show below that our choices of scales are consistent with forward charm measurements at LHCb.
The total charm and bottom cross sections per nucleon in pp and pA collisions as functions of incident proton energy are shown in the left panel of fig. 3 for nCTEQ15 PDFs for free nucleons (dashed-magenta curves) and for the case when nucleons are bound in nitrogen (solid blue curves). The range of curves reflects the uncertainty in the cross section due to the scale dependence. The dependence of the cross section on 1 Set number 55 from the CT14NLO PDF leads to total cross-sections that significantly exceed experimental upper limits from ALICE and LHCb results at √ s = 7 TeV, even when using our central values of the factorization and renormalization scales. Consequently, it has been excluded when computing the uncertainty bands throughout this work.
The gluon distribution functions for free protons (upper, magenta) and isoscalar nucleons bound in nitrogen (lower, blue) in the nCTEQ15 PDF sets [50] with Q = 2m c . The standard distribution of the PDF sets are shown with dashed lines. Small-x extrapolations with xg(x, Q) ∼ x −λ(Q) for x < 10 −6.5 are shown with dotted lines. The solid lines show PDFs with grids extended to treat the small-x regime [59], with a shaded band to show the range of predictions for the 32 sets for nitrogen, likely an underestimate of the uncertainty since the fits were made for x > 0.01. the nuclear PDFs is on the order of a few percent at the highest energies when one uses the 32 sets of nCTEQ15-14 PDFs. The right panel of fig. 3 shows with the solid blue curve the total charm cross section per nucleon, σ(pA → ccX)/A, for nitrogen with the EPS09 nuclear modification factor. For each fixed set of scales, the maximal deviation from the central cross-section due to uncertainties from the different members of EPS09 and CT14NLO PDFSets is at the level of 30% at energies of 10 10 GeV. The cross section with nitrogen (per nucleon) falls within the data constrained QCD scale uncertainties (shaded blue area) evaluated for the isoscalar nucleon cross sections in ref. [11]. In fig. 3, we vary the factorization scale from M F = 1.25m c to 4.65m c and the renormalization scale from M R = 1.48m c to 1.71m c . The data points for the total charm cross section in proton-proton collisions at RHIC and LHC energies in the figures are from [15,[60][61][62][63][64][65][66][67][68][69], while the lower energy data are from a compilation of fixed target data in [70].
The nCTEQ15-01 free nucleon sets yield slightly larger isoscalar nucleon cross section for charm production than the CT10 evaluation of BERSS [11] which are shown  [51] with Q = 2m c with CT14 PDFs [58]. The uncertainty band (blue shaded) around the central nuclear gluon distribution is obtained by combining the maximal uncertainties from the proton CT14NLO PDFs sets and those from the different EPS09 nuclear modification factors. Set 55 of CT14NLO PDFs is not included here.
by the black dotted lines in fig. 3. The nuclear corrections to the CTEQ15-01 set decrease the cross section relative to the BERSS evaluation using CT10, with a net decrease relative to CT10 of 10% at the highest energies, where the differences in the small x distribution of the PDFs are most important. The EPS09 parametrizations incorporate less nuclear shadowing at small x than the nCTEQ15 nuclear corrected PDFs. Fig. 4 shows the cross section ratio for (σ(pA → QQX)/A)/σ(pN → QQX)) for Q = c (solid lines) and Q = b (dashed lines) for isoscalar target N and A = 14. The ratio of the cross section per nucleon for partons in nitrogen and free nucleons for (M F , M R ) = (2.1, 1.6)m c using nCTEQ15 PDFs are shown in blue curves in fig. 4, and for EPS09 with CT14 free proton PDFs using the magenta curves. At low energies, where the cross section is quite small due to threshold effects, the anti-shadowing dominates, however for the energy range of interest, shadowing is more important, resulting in a 20% (10%) decrease in the cross section at high energies for cc production with the nCTEQ15 (EPS09) PDF. For bb production, the cross section is decreased by ∼ 6%-10% at E = 10 10 GeV depending on the choice of nuclear PDF.
So far, we have considered calculations based on the standard integrated parton distribution functions and the collinear framework. However as discussed above, neutrino production at high energy probes the region of very small values of x of the gluon distribution, which is not very well constrained at present. The standard DGLAP evolution, which is based on the resummation of large logarithms of scale, does not provide constraints on the small x region. Therefore, it is worthwhile to explore other approaches which resum the potentially large logarithms α s ln 1/x. There are two approaches at present, the dipole model and the k T factorization. The dipole model [25-31, 34, 35, 37-45] is particularly convenient for including corrections due to parton saturation. Parton saturation in this approach is taken into account as multiple rescatterings of the dipole as it passes through the nucleus. The dynamics is encoded in the dipole cross section, which can be either parametrized or obtained from the nonlinear evolution equation. Below we shall explore improvements to the previous calculation based on the dipole model [10], which include using more modern parametrizations for the dipole scattering cross section. Another approach to evaluating the prompt neutrino flux is based on k T  [11]. The data points for the total charm cross section from pp collisions at RHIC and LHC energies are from refs. [15,[60][61][62][63][64][65][66][67][68][69], while the lower energy data are from a compilation of fixed target data in ref. [70]. Right: Energy dependence of the charm and bottom total cross section in nucleon-nucleon collision obtained in NLO pQCD approach using NLO CT14 PDFs and the EPS09 NLO nuclear modification factor R A i (solid blue curve) [51] and (M F , M R ) = (2.1, 1.6)m Q . The upper and lower curves correspond to the same variation of the factorization and renormalization scales as in the left panel.
factorization [46][47][48][49]. In this approach the dynamics of the gluon evolution is encoded in the unintegrated parton densities, which include information about the transverse momentum dependence of the gluons in addition to the longitudinal components. We shall be using the unified BFKL-DGLAP evolution approach, with nonlinear effects, to compute the evolution of the unintegrated PDFs, which should provide for a reliable dynamical extrapolation of the gluon density towards the small x regime.

Dipole model
The color dipole model [25-27, 32, 33, 36] is an alternative approach to evaluating the heavy quark pair production cross section. The advantage of this framework is that gluon saturation at small x can be included in a relatively straightforward way, as a unitarization of the dipole-proton scattering amplitude. The partonic interaction cross section of the gluon with the target can be described in the regime of high energy by a two-step process. First, a gluon fluctuation into a qq pair is accounted by a wave function squared, then this dipole interacts with the target with a dipole cross section. In this framework, the partonic cross section for qq production can be written as [25] σ gp→qqX (x, M R , Q 2 ) = dz d 2 r |Ψ q g (z, r, M R , Q 2 )| 2 σ d (x, r) , for gluon momentum squared Q 2 and renormalization scale M R . The wave function squared, for pair separation r and fractional momentum z for q = c and q = b, is in terms of the modified Bessel functions K 0 and K 1 . The dipole cross section σ d can be written in terms of the color singlet dipole σ d,em applicable to electromagnetic scattering [27,32] Using eqs. (2.3,2.4) in the expression given by eq. (2.2), the heavy quark rapidity distribution in proton-proton scattering is given by [10] dσ(pp → qqX) dy where we use Similarly, the Feynman x F distribution in the dipole model is given by [10], in terms of the qq invariant mass squared M 2 qq and center of mass energy squared s. A LO gluon PDF is used for the value of large x 1 , while the dipole cross section encodes the information about the small x dynamics of the target, including the saturation effects.
In ref. [10], the dipole cross section parametrized by Soyez [37] was used to evaluate the prompt atmospheric lepton flux. This parametrization was based on the form discussed by Iancu, Itakura and Munier [35] which approximated the solution to the nonlinear Balitsky-Kovchegov (BK) [71,72] evolution equation. In the present calculation we use updated PDFs and dipole model parametrizations. In the flux evaluations and comparisons with LHCb data, we have updated the fragmentation fractions (see Appendix A) compared to earlier work [10]. For g(x 1 , M F ), we use the CT14 LO PDF [58]. There has been significant progress in the extractions of the dipole scattering amplitudes by including the running coupling constant (rcBK), and now, more recently, full NLO corrections to the BK equation. Here, we use Albacete et al.'s AAMQS dipole cross section result that includes heavy quarks in the rcBK formalism [39] which has been fitted to the inclusive HERA data. We compare the cross section and flux calculations using this parametrization with the calculations based on the dipole cross section by Soyez.
Finally, we use a third dipole model that is phenomenologically based. Starting from a parametrization of the electromagnetic structure function F 2 (x, Q 2 ) guided by unitarity considerations by Block et al. [45], one can show that the dipole cross section for electromagnetic scattering is approximately for z 0 2.4 [43,44]. We refer to this approximate form with the parametrization of F 2 from ref. [45] as the "Block dipole." This dipole cross section does well in describing electromagnetic, weak interaction and hadronic cross sections [73], and yields a flux similar to the AAMQS and Soyez calculations. Fig. 5 shows the cross sections for charm and bottom pair production from pp interactions calculated from the various dipole models introduced above with the gluon factorization scale varied between m c and 4m c . While all the cc cross sections are comparable at E > ∼ 10 6 GeV, for the bb cross sections, there is a difference by a factor of 1.8 (1.6) at E = 10 6 (10 8 ) GeV between the Soyez (lowest) and the AAMQS (highest) results.
The dipole cross section is applicable for x 2 1. The AAMQS dipole cross section is provided for x < 0.008. The unphysical sharp increase as a function of energy for the AAMQS result near E = 10 3 GeV is an artifact of this cutoff in x. We checked that this default x max of AAMQS has no important effect for E > 10 6 GeV, relevant energies where the prompt fluxes are dominant. As in previous work [10], we fix α s at α s = 0.373 for charm and α s = 0.215 for bottom quark production. These values of α s come from taking M R ∼ m q . We take a central factorization scale equal to M F = 2m q . These choices give reasonable cross sections as fig. 5 shows. In this approach, the cross section scales linearly with α s . Within the constraints of the cross section measurements and other experimental results, e.g., LHCb, α s can be varied with different renormalization scales. Rather than make this scale variation, we keep α s fixed for all the dipole model calculations presented here.
The comparison shown in fig. 5 of the total charm and total bottom cross sections in pp collisions shows good agreement with the data at high energies. However, at low energies dipole models underestimate the cross section because of the aforementioned limitation of x max , and the fact that quark and anti-quark contributions to the cross section are not included in this model. At high energies, initial state gluons dominate, while at low energies, this is not the case.
Nuclear effects can be incorporated in the dipole model by modifying the saturation scale. This approach, as discussed by Armesto, Salgado and Wiedemann (ASW) [74], involves a relative scaling of the free proton saturation scale by an A dependent ratio (AR 2 p /R 2 A ) 1/δ where the power δ = 0.79 is a phenomenological fit to γ * A data [74] and R A = 1.12 A 1/3 − 0.86/A 1/3 fm is the nuclear radius and R p is the proton radius. This method is used in ref. [10], where the Soyez dipole cross section is described in terms of the saturation scale which depends on r and x, however, the ASW approach cannot be used if the dipole is not parametrized in terms of a saturation scale. The method used here is the Glauber-Gribov formalism, where The nuclear profile function T A (b) depends on the nuclear density ρ A and is normalized to unity: We use a Gaussian distribution for nuclear density, with a 2 = 2R 2 A /3. This agrees well with a three parameter Fermi fit [75], used in other studies [76,77].
The nuclear corrections in the dipole model are smaller than in the NLO pQCD approach with nCTEQ15-14 PDFs, however, they are similar to the EPS09 nuclear corrections. For the Block dipole model nuclear effects range from 1% at 10 3 GeV to about 88% at 10 10 GeV, while for the AAMQS at 10 10 GeV it is 93%, and for the Soyez dipole model it is approximately 90%. The nuclear corrected cross sections for charm and bottom pair production are presented in fig. 7 with the cross sections from other approaches.

k T factorization
In this subsection we discuss the calculation of the heavy-quark production cross section using the approach of k T factorization. As mentioned previously, since the kinematics of the process is such that the dominant contribution to the neutrino flux comes from forward production of the heavy quark, the values of the longitudinal momenta of partons in this process are highly asymmetric. The longitudinal momentum fraction x of the parton participating in this process from the target side (the air nucleus) is very small, and hence one needs to extrapolate the parton densities beyond the region in which they are currently constrained by experimental data. On the other hand, we know that in the regime of small x and relatively low scales, one should take into account potentially large logarithms α s ln 1/x. Such contributions are resummed in the framework of k T factorization and BFKL evolution [78][79][80][81]. The k T factorization approach to heavy quark production in hadron-hadron collisions has been formulated in [46,47]. The framework involves matrix elements for the gg → QQ process with off-shell incoming gluons. For the forward kinematics relevant here, we shall be using an approximation in which the large x parton from the incoming cosmic ray particle is on-shell and the low x parton from the target is off-shell. This is referred to as the hybrid formalism, in which on one side the integrated collinear parton density is used, and on the other side the unintegrated gluon density with explicit k T dependence is used (for a recent calculation in the color glass condensate framework of the hybrid factorization see [82]). The k T factorization formula for the single inclusive heavy quark production with one off-shell gluon reads (2.14) In the above formula, x F is the Feynman variable for the produced heavy quark, x 1 g(x 1 , M F ) is the integrated gluon density on the projectile side,σ off (z,ŝ, k T ) is the partonic cross section for the process gg * → QQ, where g * is the off-shell gluon on the target side, and f (x 2 , k 2 T ) is the unintegrated gluon density. For the unintegrated gluon density, we have used the resummed version of the BFKL evolution which includes important subleading effects due to DGLAP evolution and the kinematical constraint [83][84][85]. These terms are relevant since they correspond to the resummation of subleading terms in the small x expansion. As a result, the calculation with resummation should be more reliable than the calculation based on purely LL or NLL small x terms. We have used the latest fits, where the unintegrated parton density has been fitted to high precision experimental data on deep inelastic scattering from HERA [86]. In addition, we have considered two cases, with or without parton saturation effects included for σ cc and σ bb , shown in fig. 6. Parton saturation was included through a nonlinear term in the parton density in the evolution [84][85][86]. Both calculations of the total integrated charm cross section, as compared with the BERSS calculation, are consistent with the perturbative calculation for high energies ≥ 10 4 GeV. The calculation without parton saturation effects is higher than with saturation. At low energies, the calculation within k T factorization tends to be below the NLO perturbative calculation within the collinear framework. This is understandable as the k T factorization can be thought of as a higher order computation with respect to the LO collinear framework, but only in the region of high energies. At low energies the ln 1/x resummation is not effective anymore, and k T factorization becomes closer to the LO collinear calculation. In order to match to NLO collinear in this region one would need to include other NLO effects in the calculation or supplement the k T factorization calculation with the energy dependent K factor.
We also analyzed the impact of nuclear corrections in the k T factorization approach. The nuclear effects in this approach are encoded in the unintegrated gluon parton density through the nonlinear term in the evolution equation as described in [86]. The strength of the nonlinear term in the nuclear case is enhanced by the factor A 1/3 with respect to the proton case.

QCD predictions for charm and bottom quark total and differential cross section
In fig. 7, we show results for the energy dependence of the total charm and bottom cross sections obtained using the three different QCD models: perturbative, dipole and k T factorization (linear evolution, nonlinear evolution). For comparison we also show our calculation based on k T factorization with nuclear effects included. We find good agreement with the experimental data with all models for LHC energies. However, at lower energies only the perturbative NLO approach gives a good agreement with the data. In the calculation of the prompt neutrino flux, the higher energy cross sections are relevant. In order to evaluate the neutrino flux, one needs to convolute the differential cross section for charm production with the steeply falling (with energy) cosmic ray flux. This evaluation is sensitive to charm production in the forward region. The differential charm quark distributions are different for different QCD models, some being in better agreement with the data than others, as we discuss below. Similar to the procedure in refs. [12] and [14], we constrain our QCD parameters by comparing our results for charm production in pp collisions with LHCb measurements in the forward region at √ s = 7 TeV and √ s = 13 TeV. In refs. [15,16], the charm cross section has been measured by the LHCb experiment in the rapidity range 2.0 ≤ y ≤ 4.5 and for the charmed hadron transverse momentum p T ≤ 8 GeV. In table 2, we show Figure 6. The integrated charm cross section in pp collisions from k T factorization, using the unintegrated gluon from linear evolution from resummed BFKL (solid blue, upper curve), and non-linear evolution (dashed magenta, lower curve). Both calculations were based on the unintegrated gluon PDFs taken from [86]. Shown for comparison is the perturbative cross section from ref. [11] (black short-dashed curve) and data points as in fig. 3.
the experimental values for the total charm cross section measured by the LHCb collaboration compared with our theoretical calculations. The experimental results from LHCb [15,16]  We also show the calculation using the dipole model (DM) and k T factorization. The dipole model uncertainty band comes from the three different dipole models and the scale variation in the gluon PDF from M F = m c to M F = 4m c . For the central values, we take the average of the cross sections for all the models considered, with M F = 2m c .
. Total cc and bb cross sections as a function of the incident proton energy. The different curves correspond to: NLO perturbative (solid blue) obtained with nCTEQ15 parton distributions, dipole model calculation based on the Block parametrization (dashed-magenta), k T factorization with unintegrated PDF from linear evolution (dashed-dotted green), k T factorization with unintegrated PDF from non-linear evolution for nucleon (short-dashed violet) and k T factorization with unintegrated PDF from non-linear evolution for nitrogen (dashed orange). Comparison is made with the results from previous NLO calculation, denoted by BERSS (short-dashed black curve), ref. [11] and data points as in fig. 3.
The upper limit of the uncertainty band corresponds to the Block dipole with M F = 4m c while the lower one is the Soyez dipole with M F = 1m c . Our results which include theoretical uncertainties are in agreement with the LHCb rapidity distributions at 7 TeV and at 13 TeV.
In refs. [15,16], data are presented for transverse momentum and rapidity distributions. Imposing a cut on transverse momentum, p T < 8 GeV where possible (see below), we show dσ/dy for 2 ≤ y ≤ 4.5 evaluated using perturbative NLO, dipole model and k T factorization. We also show the transverse momentum distributions in rapidity ranges y = 2 − 2.5, y = 3 − 3.5 (scaled by 10 −2 ) and y = 4 − 4.5 (scaled by 10 −4 ) where possible. All the calculations were performed by computing the differential distribution of charm quarks, multiplied by the fragmentation fraction for c → D 0 , and finally a factor of two was included to account for antiparticles. The results are shown in figs. 8, 9, 10 respectively. The highest rapidity bin from LHCb does not include the p T to 8 GeV, but the distribution falls off rapidly. The dipole model already includes  The experimental data are from LHCb measurements [15,16].
the full p T range, but again, the steep distribution in p T means the dipole result is a good approximation.
In fig. 8 we show NLO differential distributions of charm pairs evaluated using the free proton nCTEQ15-01 PDFs. The blue shaded band shows the prediction for the range of scales proportional to m T , while the dashed magenta lines show the predictions for the scale dependence proportional to m c , the scale range taken to be the same as in fig. 3. The m T range used in our flux evaluations is consistent with the LHCb results, as is the very large range coming from m c dependent scales. Based on this comparison, in our discussion of the prompt flux we use our error band that corresponds to the range of scales proportional to m T .
Fig. 9 presents the differential cross sections for charm pair production from the dipole models for 2.0 < y < 4.5. In dipole model calculations, there is no explicit p T dependence so we do not make any p T cuts. As the LHCb data show, the differential cross section decreases with p T . For example, already between p T = 7 − 8 GeV, the differential cross section in the y = 2 − 2.5 bin is less than 3% of the differential cross section in the p T = 1 − 2 GeV bin for √ s = 7 TeV [15], so the lack of a p T cut should not introduce a large error in the rapidity comparison.
The cross section in the dipole picture is naturally written in a mixed representation where the momentum dependence in the transverse plane is Fourier transformed to coordinate space, while the longitudinal momentum dependence is kept in momentum space (corresponding to the integration variables r and z in eq. (2.2). It is possible   [58] compared with LHCb data [15,16]. The shaded blue region shows the range of scale dependence proportional to m T , while the dashed magenta outer histograms show the scale dependence proportional to m c = 1.27 GeV. The range of scales is the same as in fig. 3 to obtain a formula for the differential cross section dσ Gp→QQX /d 2 k T , which, when integrated over d 2 k T gives eq. (2.2) [33], and in principle it could be possible to obtain the p T dependence of the cross section in this way. In the context of calculating the prompt neutrino flux, this does not seem to be a fruitful approach, but it has, however, been used to demonstrate that if the dipole cross section is calculated in LO QCD, then the LO pQCD approach is exactly equivalent to the dipole approach [33].
In fig. 9 the blue band is the range from the different dipole models with the factorization scale M F = m c to M F = 4m c for the Block   The calculation using the k T factorization formalism is shown in fig. 10. It can be seen that the differential distribution in rapidity is quite sensitive to the resummation of ln 1/x terms, in particular to the parton saturation effects. As expected, the calculation with saturation effects included is substantially below the calculation without it. To illustrate the sensitivity to the small x effects we have varied the upper limit on the integral over the transverse momentum of the off-shell gluon. This illustrates how a significant contribution comes from the lack of transverse momentum ordering in this process, and illustrates the sensitivity to transverse momenta of partons larger than the typical transverse momentum of the produced charm quark. We see that the calculation is quite sensitive to the change of the upper limit in integrals over the transverse momentum of the off-shell small x gluon. The plots with wider bands (lower plots in fig. 10) were performed with the variation of the upper limit on this integral between (m T , k max ) where k max is essentially kinematical limit in the subprocess. The smaller bands (upper plots in fig. 10) correspond to the variation of (2.5m T , k max ). We see that apparently the LHCb data exclude the most extreme limits of these calculations. For  Figure 10. Rapidity distribution for pp → D 0 /D 0 X at √ s = 7 TeV (left) and at 13 TeV (right) obtained using k T factorization formalism, and data from LHCb experiment [15,16]. Blue-bands correspond to the gluon density without saturation, magenta to the calculation including saturation. Bands represent the variation in the k T integration upper limit corresponding to (m T , k max ) (lower plots) and (2.5m T , k max ) (upper plots). See text for more explanation.
the evaluation of the neutrino flux, we have thus included the linear and non-linear calculation with scale choices of 2.5m T and k max correspondingly. They correspond to the inner band (gray shaded area) in the upper plots in fig. 10.
In Fig. 11 we show transverse momentum distributions obtained within the k T factorization formalism as compared with the our pQCD NLO results presented in fig.  8. We see that both computations are in agreement with each other and the data for this distribution.
Figs. 8, 9 and 10 all show common trends in the comparison of theory with experiment. In each case, the ratio of the rapidity distributions at √ s = 13 TeV and at √ s = 7 TeV for fixed QCD parameters is smaller than the data. However, including the theoretical uncertainty due to the choice of factorization scale, one finds reasonable dσ/dp c Figure 11. Transverse momentum distribution for pp → D 0 /D 0 X at √ s = 7 TeV (left) and at 13 TeV (right) obtained using k T factorization formalism, with data from LHCb experiment [15,16]. Red-lines correspond to the calculation from k T factorization using gluon with saturation (lower lines) and without saturation (upper lines). This calculation is compared with the NLO calculation with m T scale variation shown in fig. 8. agreement with the experimental data.
In figs. 12, 13, 14, differential distributions dσ/dx c are shown for two energies E = 10 6 and E = 10 9 GeV, for perturbative NLO (x c = x E x F ), for the dipole model (x c = x F ) and for the calculation using k T factorization (x c = x F ), respectively.
Additionally, each of the calculations is compared with BERSS calculation [11]. The difference in fig. 12 between the BERSS results (black dotted curves) and the current central NLO calculation (magenta dashed curves) can be attributed to the different choice of the PDFs (CTEQ10 vs nCTEQ15). As we show below, this PDF choice reduces the flux by about 1%−6% between 10 3 −10 8 GeV relative to the CTEQ10 choice.
In fig. 12 we also show nuclear effects on the differential distribution (blue solid curve) in the perturbative NLO calculation. We observe that the nuclear corrections evaluated here are non-negligible for higher energies. LHCb data on charm production for proton-lead collisions will be able to constrain nuclear effects for heavy nuclei in the future, as noted in ref. [87]. Fig. 13 shows that the x F distributions in the dipole model are harder than the BERSS pQCD distributions. In fig. 14, which is from the k T factorization approach, one can clearly see the effect of resummation of ln 1/x effects, as well as the impact of gluon saturation. The linear BFKL-DGLAP evolution tends to give large contribution at large x F , which increases at larger energies. On the other hand the calculation with the non-linear evolution tends to give lower values due to the suppression of the gluon density at low x. We see that the effect of the nuclear corrections are also non-negligible in this approach, and further reduce the cross section for higher energies and large values of x F . Finally, in Fig. 15 we compared calculations from all approaches which include the nuclear corrections. The NLO perturbative and k T factorization seem consistent with each other, on the other hand the dipole calculation is somewhat higher than the other two.  The prompt fluxes are evaluated using the semi-analytic Z-moment method. This procedure is described in detail in, e.g., refs. [88] and [89]. This one-dimensional method consists of using spectrum weighted differential cross section for the production of hadrons, and for decays of hadrons to neutrinos, as inputs to approximate low energy and high energy solutions to the coupled cascade equations for p, N, h, ν.  Figure 13. The differential cross section dσ/dx F as a function of x F from the dipole models for cc production, evaluated with α s = 0.373 and µ F = 2m c using the CT14 LO PDF set. The charm mass is used 1.4 GeV for the Soyez dipole and 1.27 GeV for the AAMQS and the Block dipoles. The differential cross section from ref. [11] is presented for comparison. Figure 14. Left: The differential cross section dσ/dx F as a function of x F for two energies E = 10 6 GeV and E = 10 9 GeV from k T factorization, with linear evolution (solid upper blue), and non-linear evolution (lower dashed magenta). Shown for comparison is the perturbative differential cross section from ref. [11]. Right: Comparison of the k T factorization with nonlinear evolution for the proton case (dashed magenta) and the nitrogen (solid black). equations for particle j and column depth X are  Figure 15. The comparison of the differential cross section dσ/dx F as a function of x F from NLO pQCD (Blue), the dipole model (Magenta) and the k T factorization with non-linear evolution (Green) at energies of E = 10 6 GeV and E = 10 9 GeV. All calculations contain nuclear corrections.
The Z-moment method approximates the source term for k → j with interaction length λ k for φ k (E, X) = φ 0 k (E)f (X). The factorization of the X dependence in the flux is a good approximation for the Earth's atmosphere, where we approximate the target nucleon density with an exponential atmosphere where h 0 = 6.4 km and ρ 0 h 0 = 1300 g/cm 2 . The column depth is then given by X( , θ) = ∞ d ρ(h( , θ)), where h( , θ) is the height at distance from the ground and zenith angle θ. We shall be focusing on vertical fluxes, θ = 0. Using the assumption of the exponential dependence of density on height in the atmosphere, the approximate solutions can be conveniently written in terms of the interaction lengths Λ k = λ k /(1 − Z kk ), giving f (X) = exp (−X/Λ k ). For particle k decays in the relativistic limit, λ dec k = E k τ k ρ/m k . As a result, the high and low energy lepton fluxes at Earth scale with the cosmic ray flux and can be expressed as where each Z-moment and effective interaction length Λ k depends on the energy of the prompt lepton. We have here defined the critical energy k = (m k c 2 h 0 /cτ k )g(θ) for hadron k, which separates the high-and low-energy regimes. The angular dependence of the flux enters through the critical energy, with the function g(θ) 1/ cos θ for small angles close to vertical, but for angles near horizontal, it is more complicated due to the geometry of the Earth and the atmosphere. We will not need the details here, but more information can be found in e.g. [89]. For both h c and h b , the transition between low energy and high energy fluxes is at energies E ν ∼ 10 8 GeV. We evaluate the flux by taking the sum over hadron contributions: .

(3.10)
This approach is applicable to ν = ν µ , ν e and ν τ plus their antiparticles (and to µ + + µ − , which are stable nearly massless leptons to first approximation). The prompt electron neutrino flux and prompt muon flux are essentially equal to the prompt muon neutrino flux. This comes from the nearly identical kinematics in the semileptonic charmed hadron and b-hadron decays for the eν e and µν µ final states. Tau neutrinos, however, are produced at a lower rate, with the dominant source being h = D s [8,9] where D s → ν τ τ comes with a branching fraction of 5.54 ± 0.24% [55].
Tau decays are also prompt, with cτ τ = 87.03 µm (as compared to cτ Ds = 149.9 µm [55].). Tau energy loss in the atmosphere is negligible, so we include the chain decay D s → τ → ν τ as well as the direct decay D s → ν τ . Since most of the D s energy goes to the tau lepton, and approximately 1/3 of the tau energy goes to the tau neutrino, the chain decay contribution is larger than the direct contribution of D s decays to the prompt tau neutrino flux. For the chain decay D s → τ → ν τ , In the results below, we approximate all τ decays as prompt. This overestimates the tau neutrino flux above E ν ∼ 10 7 GeV [9], however, at these high energies, the tau neutrino flux is quite low anyway.

Cosmic ray flux, fragmentation and decays
A broken power law (BPL) approximation of the cosmic ray nucleon flux in terms of the energy per nucleon E is of the form nucleons cm 2 s sr GeV = 1.7 (E/GeV) −2.7 E < 5 · 10 6 GeV = 174 (E/GeV) −3 E > 5 · 10 6 GeV . (3.12) The BPL has been used to evaluate the prompt flux in many earlier references. While the all particle spectrum resembles a broken power law, recent analyses have shown that it poorly represents the nucleon spectrum, even though the composition of the cosmic rays is not completely known [53,54,90]. Recent three component models [90] with extragalactic protons (called H3p here) or an extragalactic mixed composition (H3a) come from an analysis of cosmic ray measurements paired with fits to functional forms for the spectrum by composition. The four component model by Gaisser, Stanev and Tilav (GST*) [53] is labeled here by GST4. To compare with other flux calculations, we use the BPL for reference. The H3a flux dips more precipitously at high energies than the H3p flux. We show below that the H3p and GST4 inputs to the prompt flux lead to similar predictions. Our three approaches to charm production provide the x E or x F distribution of the charmed quark, and similarly for the b quark. We use the Kniehl and Kramer fragmentation functions for charm [91], using the LO parameters with the overall normalization scaled to account for updated fragmentation fractions determined from a recent review of charm production data in ref. [92]. The original normalizations in ref. [91] had the fragmentation fractions add to 1.22 rather than 1.00. Using the updated fragmentation fractions to rescale the fragmentation functions, the sum of fractions is 0.99. For the B meson fragmentation, we use the power law form of Kniehl et al. from ref. [93], rescaled to match the fragmentation fractions of ref. [94]. The details of the fragmentation functions and parameters are shown in Appendix A.1. Decay distributions and parameters for the branching fractions and effective hadronic masses for semi-lepton heavy meson decays are listed in Appendix A.2.
The Kniehl and Kramer (KK) provide LO and NLO fragmentation functions [91], the later, in principle, being more suitable for our NLO pQCD calculation. However, the dipole model approach is a LO calculation, so we use the LO KK fragmentation, consistent with the way fragmentation was used in ERS [10]. We note that in our previous NLO pQCD work, BERSS [11], we also used the LO KK fragmentation functions. Therefore, in order to compare our new NLO pQCD results and results from different dipole models and k T factorization with the previous work in ERS and BERSS, we use the LO parametrization of the KK fragmentation function even for the NLO pQCD calculation. Since the dominant contribution is at threshold, the fragmentation functions are not evolved. We find that if we use NLO KK fragmentation functions, it results in only 2% lower fluxes than obtained with the LO KK fragmentation functions when both are rescaled to match the fragmentation fractions in [92].
The results presented below use z = E h /E Q as the fractional variable in the fragmentation functions. There are other choices possible, however, this definition can be used for all three approaches. The flux has some sensitivity to the charmed hadron fragmentation functions. We have also used the Braaten et al. [95] (BCFY) fragmentation functions for quark fragmentation to pseudoscalar and vector states, updating the discussion of Cacciari and Nason in ref. [96], for D 0 and D ± , and approximate forms for D s and Λ c . In comparison to the Kniehl-Kramer fragmentation results, the flux evaluated with the BCFY fragmentation functions may be as much as 30-50% larger depending on parameter choices.
The fragmentation process is thus somewhat uncertain. An alternative to using fragmentation functions is to use Monte Carlo event generators such as Pythia [97] to hadronize the partonic state. This typically leads to a harder spectrum of D-mesons than with fragmentation functions, see e.g. [98,99]. The reason is that the hadron that contains the charm (or anticharm) quark can pick up some momentum from the beam remnant of the interaction, and the charmed hadron can in fact have a larger longitudinal momentum than the original charm quark. In Pythia this is modeled as a drag effect from the string. This effect is not possible in e + e − collisions, and in hadronhadron collisions it is larger for more forward production. Fragmentation functions, however, are fitted from e + e − data and more central hadron-hadron data, so if such an effect is real, it may not show up in the fit.
Garzelli et al. [12] computed the charm cross section at NLO using the NLO Monte Carlo POWHEG BOX [24,[100][101][102], and used Pythia for fragmentation. The resulting prompt flux is compared with the BERSS flux [11], which is computed using fragmentation functions, and is found to be around 30% larger than BERSS. We believe that one reason for this is the different fragmentation method used. We have checked, using POWHEG BOX, that the produced spectrum of D-mesons is indeed harder when the Pythia fragmentation is used rather than the KK fragmentation functions applied to the same spectrum of charm quarks. The difference becomes larger for larger x F , but is small for small x F . The situation is therefore somewhat unclear, and ideally more forward data would be needed to resolve the question. The most forward data that we have considered is the LHCb data discussed above, but this is not forward enough to be sensitive to differences in fragmentation. The prompt flux, however, could be sensitive to such effects.
Finally, to evaluate the prompt atmospheric lepton fluxes, we need the Z-moments Z pp , Z DD , Z BB and interaction lengths. We make the same assumptions as in ref. [11] (BERSS). We use expressions for dσ/dx E that factorize into energy and x E dependent functions. The energy dependence ensures that the differential cross section grows with energy following the relevant cross section growth with energy. For pA scattering, the cross section parametrization of EPOS 1.99 is used [103], with the x E dependent function proportional to (1 − x E ) n with n = 0.51. For DA and BA scattering, we approximate the energy dependence with the cross section for kaon-nucleon scattering determined by the COMPAS group (see ref. [55]), scaled by A 0.75 for nuclear corrections, and use the same (1 − x E ) n behavior for both B and D scattering with air nuclei, with n = 1.

Prompt muon neutrino flux
We begin with the muon neutrino fluxes predicted by each of the three approaches.  [57] and used in our earlier flux evaluation [11].
The new evaluation of the prompt flux has several differences relative to the BERSS result shown in black. Here, with the updated charm fragmentation fractions [92], the flux is reduced by about 20%. The bb contribution was not included in BERSS. It contributes 5-10% of the flux at E ν ∼ 10 5 − 10 8 GeV, as shown in the left panel of fig. 17. Finally, in this work we have included nuclear effects on the prompt neutrino flux, which reduces the flux by 20%-30% for energies between E = 10 5 − 10 8 GeV for the nCTEQ15-14 PDFs, with the largest effect at the highest energy. We note that the nuclear effect is more pronounced on the prompt neutrino flux than it is for the total charm cross section, due to the fact that the flux calculation probes forward charm production (small x of the parton from the target air nucleus) where nuclear suppression is larger. For reference, the total charm cross section section suppression due to nuclear effects is between 4% and 13% at 10 5 − 10 8 GeV with the nCTEQ15-14 PDFs. The gluon PDF in fig. 1 helps illustrate this point. The cross section is dominated by the small x F region, where the parton momentum fractions are nearly equal, so probing less the shadowing region. The ratio of the flux with nuclear effects to the flux using free protons (nCTEQ15-14 PDFs compared to nCTEQ15-01 PDFs) is shown as a function of energy in the right panel of fig. 17.
The combination of all these effects results in our NLO pQCD prompt flux estimate being 30% lower than BERSS at 10 3 GeV, about 40% lower at 10 6 GeV and almost 45% lower at 10 8 GeV, when we use nCTEQ15-14 PDF as parton PDFs in the air.
When we use CT14 PDFs plus EPS09 for nuclear effects, our results are only moderately affected by nuclear corrections. In the left panel of fig. 18, we show the fluxes, and in the right panel, the ratio of the flux with nuclear effects to the one without. At very high energies, the CT14 PDFs predict a similar flux to the one obtained with nCTEQ15 PDFs, with the nuclear correction being somewhat smaller than for nCTEQ15 case. Nuclear corrections are uncertain for a larger range of x. The EPS09 suppression factors are frozen at R A i (x min , Q) for x < x min = 10 −6 , halting a decline in energy for the ratio of fluxes with nuclear corrected and free nucleon targets.
The dipole model results are shown in fig. 19, together with our ERS dipole result from ref. [10] for the broken power law. Compared to the ERS result, we have used updated PDFs (LO CT14) and included the bb contribution, and we have considered two other dipole models beyond the Soyez model used in ref. [10]. In comparing the Soyez dipole calculations, the updated fragmentation fraction reduces the overall flux by Figure 17. Left: The NLO pQCD flux prediction from bottom hadrons. The fluxes from B hadrons have a ratio of about 2 %, 7 (6)% and 9 (7) % to those from charm at 10 3 , 10 6 and 10 8 GeV, respectively, for nitrogen (Proton) PDF. Right: Nuclear effect in the prompt neutrino flux evaluated in the NLO pQCD approach. Figure 18. The central prompt neutrino flux prediction using the CT14 PDFs with EPS nuclear corrections (left), and the ratio of the fluxes with and without the nuclear corrections (right), as a function of neutrino energy. As in fig. 16, the upper and lower limits correspond to variation in the QCD scales and the uncertainty from the different PDF sets.
approximately 20%. Relative to the ERS calculation, we have updated the Z pp and Z DD moments, as discussed in detail in ref. [11], which gives a further reduction (about 35%) in the flux prediction. The AAMQS dipole and phenomenological Block dipole give very similar results and are the upper part of the uncertainty bands. Nuclear corrections to the dipole model flux predictions reduce the flux by about 10% at E ν ∼ 10 5 GeV and reduce by about 20% at E ν ∼ 10 8 GeV. Fig. 20 shows the ν µ +ν µ flux predictions in the k T formalism with linear and Figure 19. Dipole model flux predictions for ν µ +ν µ as a function of energy, compared to the ERS dipole results [10]. The error bands come from the three dipole models, and varying the factorization scale dependence from M F = m c to M F = 4m c . nonlinear evolution for the unintegrated gluon density. Nuclear corrections for nitrogen are included in this approach. The predicted high energy flux in the k T factorization formalism is consistent the other approaches, with the exception of the low energy where the flux is somewhat smaller. The low energy deficit reflects the same deficit of the cross section shown in fig. 6 since the k T factorization model applies to small x physics and therefore applies to high energies. At the high energies shown, the linear k T approach is about 7 times larger than the non-linear k T flux prediction, reflecting the range of impact that small-x effects can have on the high energy prompt flux. Figure 21. Comparison of the muon neutrino plus antineutrino fluxes using all the approaches: NLO perturbative QCD with nCTEQ15 (blue) and EPS09 (orange), dipole model (magenta), k T factorization (green) with the other calculations (black): BERSS [11], ERS [10], GMS [12] and GRRST [14].
Finally, in fig. 21, we compare the three approaches using the broken power law with the BERSS [11], ERS [10], GMS [12] and GRRST [14] results. Relative to the BERSS flux, the dipole model predicts a larger low energy flux, while the k T factorization model based on the linear evolution predicts a larger high energy flux. On the other hand the flux based on the k T factorization with nuclear corrections is consistent with the lower end of the NLO pQCD calculation. Our new perturbative result lies below the BERSS band for most of the energy range, due to a combination of the nuclear shadowing and the rescaling of the fragmentation fractions to sum to unity. The total uncertainty range of our predictions from the different approaches is compatible with the other recent evaluations of GMS and GRRST, with the latter giving a larger error band.

Prompt tau neutrino flux
To finish our discussion of prompt neutrino fluxes, for completeness we show the tau neutrino plus antineutrino flux in fig. 22 in the NLO pQCD framework. Using the central scale choice in the NLO pQCD calculation of charm pair production with nCTEQ15-14 PDFs, together with scale variations to obtain the uncertainty band, and the KK fragmentation functions, the sum of the direct and chain contributions from D s and the semileptonic B 0 and B + (and charge conjugate) decays to ν τ , the prompt tau neutrino flux is shown in fig. 22 for the broken power law and H3p cosmic ray fluxes [52]. As above, the KK fragmentation fraction of 12.3% of charm to D s has been updated to 8.1% following ref. [92].
The flux is about 10% of the prompt muon neutrino plus antineutrino flux. As noted above, the very high energy flux shown here overestimates the tau neutrino flux because we have approximated the tau decays as all prompt. The steeply falling shaded green band shows the range of the tau neutrino flux coming from ν µ → ν τ oscillations for ν µ 's from pion and kaon decays, with the upper edge coming from ν µ → ν τ conversion through the diameter of the Earth. Secondary production of tau neutrinos from interactions of atmospheric leptons in the Earth are quite small [104][105][106].
The chain decay D s → τ → ν τ dominates the flux. The tau neutrino flux from the chain decay is approximately a factor of 4 larger at low energies than from the direct decays, and becoming larger at higher energies. The fraction of the tau neutrino flux from B's is shown in the right panel of fig. 22. We have not included the B s and Λ b decays as the branching fractions to ν τ are not measured. They may contribute as much to the prompt tau neutrino flux as the B 0 + B + contributions.

Comparison with IceCube limit
In fig. 23 we show our results for the prompt neutrino flux obtained in the three different QCD approaches, compared with the conventional neutrino flux and the IceCube limit on the prompt neutrino flux using 3 years of data [107]. In the left figure, we scale the flux by E 2 and in the right figure, by E 3 . The IceCube result is an upper limit at 90% C.L. that places a bound on the normalization of a flux with the same spectral shape as the ERS model [10], and is thus not completely model independent. The upper limit corresponds to 0.54 times the ERS flux, rescaled from the broken power law CR flux used by ERS to the H3p CR flux, using the method proposed in [90]. For comparison, the red band in fig. 23 represents the 1σ error on the measured astrophysical neutrino flux. It is important to note that the evaluation of this IceCube limit is not independent of the modeling of the astrophysical neutrino flux, which in this case is taken as an unbroken power law, and the normalization of the ERS flux is taken as a free parameter in a likelihood fit to the data, yielding the displayed upper limit.
From fig. 23 we note that the IceCube limit is in tension with all dipole model predictions, and very close or at the border of the upper limit of the k T factorization approach. On the other hand both the NLO pQCD prediction which includes nuclear effects via the nuclear parton distributions, and the nonlinear k T calculation, are below the IceCube limit. We note, however, that the nuclear effects in the dipole model and with the EPS09 pQCD approach are smaller than in the nCTEQ15-14 pQCD approach. IceCube data may help distinguish between nuclear suppression models at small-x.

Discussion and conclusions 4.1 LHC and IceCube
As figs. 8,9,10 show, rapidity distributions measured at 7 and 13 TeV [15,16] seem to be somewhat in tension within all three approaches if one considers a fixed prescription for the scales independent of energy. The theoretical error bands, however, do accommodate the data as noted in ref. [14]. Figs. 8, 9, 10 compare the distributions of charm quarks with the measured D 0 distributions. In the case of the k T factorization approach the 7 TeV data seem to be more consistent with the calculation with the Figure 23. A comparison of our neutrino fluxes from different QCD models with the IceCube upper limit from ref. [107]. The results are for the H3p cosmic ray spectrum. We also show the conventional vertical neutrino flux. nonlinear gluon density, or the lower band of the calculation with linear gluon density, whereas the data at 13 TeV are more in line with the evaluation with the linear evolution. This is rather counterintuitive and perhaps could suggest that the calculation with nonlinear evolution is disfavored by the data. However, given the spread of the uncertainty of the calculation it is not possible to make decisive statement at this time and more studies are necessary. The dipole model evaluation favors the Soyez form for √ s = 7 TeV and the AAMQS or Block form for √ s = 13 TeV for our fixed value of α s . The central pQCD predictions seem to indicate that the distributions don't rise quickly enough with increasing √ s. In ref. [13], the NLO pQCD prediction of the ratio of dσ/dy in the forward region for LHCb for √ s = 13 TeV to √ s = 7 TeV was predicted to be on the order of 1.3-1.5, which we also see in fig. 8. The data show the ratio to be closer to a factor of 2. Nevertheless, for all three approaches, the LHCb data fall within the theoretical uncertainty bands.
We have calculated the rapidity and p T distribution using our theoretical QCD parameters, i.e. the range of factorization scales for a given charm mass which was determined from the energy dependence of the total charm cross section. We have found our results to be consistent with LHCb data. The range of m T dependent factorization scales in the pQCD evaluation adequately cover the range of LHCb data, while the range of m c dependent factorization scales overestimate the uncertainty. In the case of the dipole models, in which there is no explicit p T dependence, we have only made comparison with the LHCb rapidity distributions.
The IceCube Neutrino Observatory and other high energy neutrino detectors may be useful in getting a handle on forward charm production. Indeed, the high energy prompt lepton flux depends on charm production at even higher rapidity than measured by LHCb, as can be seen by the following argument. In both the high and low energy forms of the prompt lepton fluxes, the Z-moments for cosmic ray production of charm, e.g., Z pD 0 (E), depend on the lepton energy E. To evaluate the Z-moment for charm production, the energy integral over E in eq. (3.6) can be cast in the form of an integral over x E = E/E that runs from 0 → 1, account for incident cosmic rays (p) with energy E producing, in this case, D 0 with energy E. Fig. 24 shows the fraction of the Z-moment integral in eq. (3.6) for x E = 0 → x max for two different energies using NLO pQCD with the central scale choice and the H3p cosmic ray flux. For E = 10 6 GeV, about 10% of the Z-moment comes from x E < x c = 3.6 × 10 −2 , while for E = 10 7 GeV, this same percentage comes from x E < x c = 1.5 × 10 −2 . We can use the value of x E > x c that gives 90% of the Z-moment as a guide to what are the useful kinematic ranges in high energy pp collider experiments.
We approximate in terms of the hadronic center of mass rapidity, which leads to for 90% of the Z-moment evaluation. For E = 10 6 GeV, this indicates that the Zmoment is dominated by y cm > 4.9 with Finally, let us note that there could potentially be another source of charm in the proton. It has been suggested [108,109] that there could exist heavy quark pairs in the Fock state decomposition of bound hadrons. This would be an additional nonperturbative contribution and is usually referred to as intrinsic charm to distinguish it from the perturbative and radiatively generated component considered in this work. Intrinsic charm parameterized in PDFs has been explored in e.g., refs. [110,111]. There are studies that explore how to probe intrinsic charm in direct and indirect ways [112,113]. Intrinsic charm, if it exists, would be mostly concentrated at high values of x of the proton and may therefore be another contribution to the very forward production relevant to the prompt lepton flux. Its unique features were recently discussed in [114][115][116]. We note however, that the current IceCube limit on the prompt flux is already quite constraining and leaves a rather narrow window for a sizeable intrinsic charm component. Eventual IceCube observations (rather than limits) of the prompt atmospheric lepton flux, may be unique in its ability to measure or constrain the physics at high rapidities.

Summary
In this work we have presented results for prompt neutrino flux using several QCD approaches: an NLO perturbative QCD calculation including nuclear effects, three different dipole models and the k T factorization approach with the unintegrated gluon density from the unified BFKL-DGLAP framework with and without saturation effects. Numerical results are listed in tables 8, 9, 10 and 11. The energy dependence of the total charm cross section, measured from low energies (100 GeV) to LHC energies (13 TeV), is best described with NLO pQCD approach. On the other hand the dipole and k T factorization approaches are theoretically suited to describe heavy quark production at high energies, however, they need additional corrections at lower energies.
We have included theoretical uncertainties due to the choice of PDFs, choice of scales and nuclear effects, constrained by the total charm cross section measurements for energies between 100 GeV up to 13 TeV. We have shown that differential cross sections for charmed mesons obtained with these QCD parameters are in good agreement with the LHCb data. We have found that the prompt neutrino flux is higher in case of the dipole model and the k T factorization model (without saturation) than the NLO pQCD case. The former seem to be numerically consistent with the previous ERS [10] results, while NLO pQCD is smaller than BERSS [11]. For the nCTEQ15-14 evaluation, this is mostly due to the nuclear effects. In particular, we have found that the nuclear effect on the prompt neutrino flux is large in case of the pQCD approach with nCTEQ15-14 PDFs, as large as 30% at high energies, while this effect is smaller (∼ 20%) for the dipole model approach. The EPS09 nuclear corrections suppress the pQCD flux calculations with free nucleons by only ∼ 10%. The nuclear corrections are also significant in the k T factorization approach, as large as 50% at high energies, thus lowering the flux to the level comparable with that obtained using the NLO pQCD with nuclear PDFs.
Contributions from bb are on the order of 5-10% to the prompt flux of ν µ +ν µ in the energy range of interest to IceCube. For completeness, we have also evaluated the flux of ν τ +ν τ from D s and B decays.
We have also shown results for different cosmic ray primary fluxes and show how the shape of the particular choice affects the neutrino flux. As before, the updated fluxes for the primary CR give much lower results than the simple broken power law used in many previous estimates.
Finally we have compared our predictions with the IceCube limit [2]. We have found that the current IceCube limit seems to exclude some dipole models and the upper limit of the k T factorization model (without any nuclear shadowing), while our results obtained with the NLO pQCD approach with nCTEQ15-14 and the calculations based on the k T factorization with nuclear corrections included are substantially lower and thus evade this limit.
Since it is very important to determine the energy at which prompt neutrinos become dominant over the conventional neutrino flux, we expect that the calculation of the conventional flux might be improved by using the two experiments at the LHC that have detectors in the forward region, the Large Hadron Collider forward (LHCf) experiment [117] and the Total, Elastic and Diffractive Cross-section Measurement experiment (TOTEM) [118]. The LHCf experiment measures neutral particles emitted in the very forward region (8.8 < y < 10.7), where particles carry a large fraction of the collision energy, of relevance to the better understanding of the development of showers of particles produced in the atmosphere by high-energy cosmic rays. The TOTEM experiment takes precise measurements of protons as they emerge from collisions in the LHC at small angles to the beampipe, thus in the forward region. In addition to measuring the total and elastic cross section, TOTEM has measured the pseudorapidity distribution of charged particles at √ s = 8 TeV in the forward region (5.3 < |η| < 6.4).
These measurements could be used to constrain models of particle production in cosmic ray interactions with the atmosphere and potentially affect the conventional neutrino flux, which is coming from pion decays. Future IceCube measurements have a good chance of providing valuable information about the elusive physics at very small x, in the kinematic range which is beyond the reach of the present colliders. Keeping in mind the caveats involved in the current IceCube treatment of the atmospheric cascade and the incoming cosmic ray fluxes, the observation or non-observation of the prompt flux may give important insight into the QCD mechanism for heavy quark production.
First, the nuclear gluon distribution in the region x ≤ 0.01 is currently not constrained with collider or fixed target experiments. On the other hand, we expect that the upcoming 6 year IceCube data will be sensitive to our pQCD flux results, especially those obtained with the EPS framework that includes nuclear effects.
Second, from our study we find that the IceCube limit shown in fig. 23 already severely constrains the dipole model approach, even with the lowest cosmic ray flux (H3p). While it is possible that a modified dipole approach, such as a next-to-leading order calculation, would yield a lower charm cross section that is not in tension with IceCube, the dipole model calculation used here is not flexible enough to modify so that it evades the limit, i.e., this tension cannot be solved by adjusting dipole parameters, because they are constrained by the LHCb data.

Note added
After submitting our paper, the IceCube Collaboration released their upper limit on the prompt atmospheric muon neutrino flux as 1.06 times the ERS flux based on the 6 year data [119]. We find that our results presented in this work are below this new IceCube limit.

A.1 Fragmentation
The parameters for the charm fragmentation functions are shown in Table 3 for the fragmentation function of energy fraction z, of the form discussed by Kniehl and Kramer in ref. [91], (A.1) We use the LO parameters of ref. [91], with the overall normalization N rescaled to account for updated fragmentation fractions determined from a recent review of charm production data in ref. [92].  Table 3. Parameters for the charm quark fragmentation from [91], with the normalization N rescaled to match the fragmentation fractions in ref. [92].
Alternate fragmentation functions for charmed hadrons is provided by Braaten et al. [95] (BCFY) fragmentation functions for quark fragmentation to pseudoscalar and vector states, with input from Cacciari and Nason in ref. [96] for D 0 and D ± . The parameter r in their fragmentation functions  Table 4. Following ref. [96] for the BCFY fragmentation functions of the form eq. (A.2), normalized to the fragmentation fractions of ref. [92].
For the B meson fragmentation, we use the power law form of Kniehl et al. from ref. [93], rescaled to match the fragmentation fractions of ref. [94]. The functional form of the fragmentation functions for B mesons is with the parameters in  Table 5. Parameters for the b quark fragmentation from Kniehl et al. [93], with the normalization rescaled to match the fragmentation fractions in [94,121]. The Λ b fragmentation function is approximated by the meson fragmentation function, normalized to a fragmentation fraction B b = 0.236.

A.2 Decay distributions
Meson decay moments are evaluated using the decay distributions where z = E ν /E h , the fraction of the hadron energy carried by the neutrino, and B h→ν is the branching fraction. Following ref. [122], we approximate charmed meson semileptonic decay distributions as a function of neutrino energies by the three-body decay distribution. (See also ref. [89].) This fractional energy distribution comes from evaluating the pseudoscalar three-body semileptonic decay to a lighter pseudoscalar meson (e.g., D → K ν ) with form factor f + (q 2 ) f + approximately constant and f − (q 2 ) 0, keeping the mass of the meson in the final state but neglecting the lepton mass. The distribution is the same for = µ, e and ν and is where r = m 2 /m 2 h for m the mass of the hadron in the final state and m h the mass of the decaying particle. Details for the decay in rest frame are described in ref. [123], and the procedure to convert the distribution to the frame where the decaying particle has energy E h can be found in ref. [88]. We use the same form of the decay distribution for B meson semileptonic decays, and to simplify the evaluation, for Λ c and Λ b decays to ν µ and ν e .
In ref. [122], an effective hadron mass s ef f X is used for charm decays to account for the contributions of both pseudoscalar mesons, vector mesons and two mesons in the final state. The effective hadron masses are used for r = s ef f X /m 2 h . The branching fractions and effective hadron masses used in our evaluations of the prompt fluxes are listed in table 6.
The three body decays of B mesons to τ ν τ require additional mass terms proportional to m 2 τ /m 2 B . Using the same approximations as above except for keeping the tau mass, we find for the neutrino distribution, and  Table 6. Parameters for the branching fractions and effective hadronic mass for charmed and b hadron decays, following ref. [122] for charmed meson decay distributions.
for the tau distribution, where r sum = (m τ + s ef f X ) 2 /m 2 b and w(a, b, c) = [a 2 + b 2 + c 2 − 2(ab + bc + ac)] 1/2 . Numerically, for the effective hadron masses in table 6, the normalization constants are N ν = N τ = 2.72 × 10 −2 . The kinematic limits for y are for r τ = m 2 τ /m 2 Ds = 0.815 and branching ratio B Ds→τ ντ = (5.54 ± 0.24)% [55]. The energy distribution of the taus in D s decays, in terms of z τ ≡ E τ /E Ds , is of the same form as eq. A.11, however, with a different range for z τ that gives the τ a larger fraction of E Ds than the ν τ energy fraction: The energy distribution of tau neutrinos in tau decays in terms of y ≡ E ν /E τ depends on the tau polarization. It is for τ → i decays. The branching fractions B i and functions g i 0 , g i 1 are listed in Table 7. Here, the function P τ = P τ (E Ds , E τ ) is the polarization of the tau from the D s decay: in the relativistic limit. The energy distributions of ν τ andν τ are the same, a result of the opposite polarizations of the τ andτ and the opposite neutrino and antineutrino helicities. We have approximated P τ for taus coming from the semi-leptonic B decays with the polarization of taus coming from D s decays.

A.3 Flux Tables
We provide here in tables 8, 9 and 10 numerical values of the distributions in the vertical direction for E 3 φ for ν µ +ν µ from cc and bb from atmospheric production by cosmic rays. Up to E ν ∼ 10 7 GeV, these fluxes are isotropic. The prompt fluxes for = ν µ , ν e and µ are equal. Table 11 shows the vertical ν τ +ν τ flux evaluated using the nCTEQ15-14 PDFs. log 10 (E/GeV)     Table 11. The NLO pQCD vertical flux of ν τ +ν τ scaled by E 3 in units of Gev 2 /cm 2 /s/sr using the nCTEQ15-14 PDFs with low-x grids, evaluated with m T dependent renormalization and factorization scales, with the H3p cosmic ray flux.