The prompt atmospheric neutrino flux in the light of LHCb

The recent observation of very high energy cosmic neutrinos by IceCube heralds the beginning of neutrino astronomy. At these energies, the dominant background to the astrophysical signal is the flux of `prompt' neutrinos, arising from the decay of charmed mesons produced by cosmic ray collisions in the atmosphere. In this work we provide predictions for the prompt atmospheric neutrino flux in the framework of perturbative QCD, using state-of-the-art Monte Carlo event generators. Our calculation includes the constraints set by charm production measurements from the LHCb experiment at 7 TeV, and has been recently validated with the corresponding 13 TeV data. Our results for the prompt flux are a factor of about 2 below the previous benchmark calculation, in general agreement with two other recent estimates, and with an improved estimate of the uncertainty. This alleviates the existing tension between the theoretical prediction and IceCube limits, and suggests that a direct direction of the prompt flux is imminent.


Introduction
The IceCube experiment [1] at the South Pole has recently made the first detection of high energy cosmic neutrinos from the Southern sky with deposited energies between 30 and 2000 TeV and arrival directions consistent with isotropy [2][3][4]. Although these are mainly ν e and ν τ charged-current and neutral-current ('cascade') neutrino interactions, the 37 events are consistent with expectations for equal fluxes of all three neutrino flavours [5]. Subsequently cosmic ν µ charged-current ('track') events have also been seen from the Northern sky [6] with comparable flux [7]. At these high energies, the 'conventional' atmospheric neutrino flux, from the decays of pions and kaons produced by the collisions of cosmic rays with nuclei in the atmosphere [8][9][10][11], is suppressed due to energy loss before the decays occur. However charmed mesons decay almost instantaneously and therefore, despite their smaller production cross-section, the 'prompt' neutrino flux from their decays (∝ E −2.7 ν ) should dominate over the conventional flux (∝ E −3.7 ν ) at high energies. Thus the prompt component is the most relevant background for the similarly hard spectrum expected for the astrophysical neutrino flux [12,13]. Indeed the statistical significance (5.7σ) with which an atmospheric origin can be rejected for the 37 IceCube events is limited by the uncertainty in the expected atmospheric prompt neutrino flux.
Many calculations of the prompt neutrino flux have been presented [14][15][16][17][18][19][20][21][22][23][24], however so far IceCube has not detected it and set only an upper limit of 1.52 times the central value of the benchmark ERS calculation [16] at 90% CL [25]. In a recent analysis this limit has been lowered further by a factor of 3 to only 0.54 times the above benchmark calculation [26]. This motivates a re-evaluation with state-of-the-art tools and inputs, providing, in particular, an improved estimate of all theoretical uncertainties. The major uncertainty is in the calculation of charm production at high-energies due to higher-order corrections and, especially, the imprecise knowledge of the gluon parton distribution function (PDF) at small-x Bjorken . A recent breakthrough has been the availability of charm hadroproduction data from the LHCb experiment [27,28] which covers the kinematical range directly relevant to the calculation of the atmospheric prompt neutrino flux.
With this motivation, we have recently validated state-of-the-art perturbative QCD calculations with the LHCb forward charm production data at 7 TeV [27,28], and included these measurements into the NNPDF3.0 global analysis [29]. We were thus able to construct a new PDF set, NNPDF3.0+LHCb, which is tailored for calculation of the prompt neutrino flux. 1 We benchmarked three other codes, FONLL [32], POWHEG [33][34][35] and aMC@NLO [36], finding good agreement both amongst themselves and with the LHCb data. Our calculations have subsequently been found to be in good agreement with the 13 TeV LHCb charm production data [37], which probe even smaller values of x.
In this work we calculate the atmospheric prompt neutrino flux using the canonical set of cascade equations implemented in the 'Z-moment' framework (see [16] and references therein). Charm cross-sections and decays are obtained using the NLO Monte Carlo generator POWHEG with the NNPDF3.0+LHCb PDF set as input. We consider several parameterisations of the cosmic ray flux, including the most recent models [38,39], and study the dependence of our result on the choice of input PDF set.
We compare our calculation with previous results, in particular the benchmark ERS calculation [16], as well as the recent BERSS [20] and GMS [24] analyses. We emphasise that our calculation is the only one which is directly validated with LHCb data. All four calculations are consistent within our theoretical uncertainty band, with ERS being at the edge of the upper limit. Our central value is similar to the BERSS result, while the GMS result is a little higher. We also compare our result to the recent IceCube limit on the prompt neutrino flux, finding that our central value is consistent. Moreover our lower limit indicates that the prompt neutrino flux will soon be detected, thus enabling reliable subtraction of any contamination of the astrophysical neutrino candidates. This paper is organised as follows. In § 2 we discuss the various inputs that enter the calculation of the prompt neutrino flux, including the parameterisations of the cosmic ray flux, the solution of the cascade questions, and the calculations of the various Z moments. The results for the prompt flux are presented in § 3, where we compare with other recent determinations as well as with the latest IceCube limit. We also study the dependence of our result on the input PDF set and on the cosmic ray flux parameterisation. Our results are summarised in § 4 and are made publicly available in the form of an interpolation code which returns the prompt neutrino flux and its uncertainty for all adopted models of the cosmic ray flux.

Calculation of the prompt neutrino flux
First we present the parameterisations of the cosmic ray flux used in this work. We review the cascade equations for the propagation of particles in the atmosphere, and their solu-tion using Z-moments. Then we discuss the calculation of the various Z-moments, with emphasis on the charm production cross-section and the associated theory uncertainties.

The incoming cosmic-ray flux
The flux of cosmic rays incident on the atmosphere is dominated by protons and has been measured by a variety of experiments (see recent reviews [38][39][40][41]). At energies 10 3 GeV relevant for calculating the prompt neutrino flux, a traditional parameterisation has been the broken Power Law (BPL) with a 'knee' at E p = 5 × 10 6 GeV, assuming all cosmic rays are protons: Recently, more elaborate parameterisations of the cosmic ray flux have been provided, with emphasis on including composition data and improving the description above the 'knee' in the spectrum. One such set [43] follows Hillas' proposal [44] for accommodating three populations of cosmic rays: one associated with acceleration by supernova remnant shocks, a second galactic component from unspecified sources, and finally an extra-galactic component which dominates at the highest energies. The assumption is that 5 groups of nuclei, i, are contained in each of these three source components, j, such that the cosmic ray spectrum for the nuclear species i can be written as where R c,j is the magnetic rigidity for the source component j, and a ij and γ ij are the corresponding normalization constants and spectral indices [43]. We construct an equivalent 'all-proton' spectrum 2 by re-weighting the various nuclei: with A i the atomic number of species i, and, to obtain the total cosmic ray flux, we sum over each of the 5 nuclear components: Thus we do not need to consider an effective nuclear attenuation length, since collective effects in nucleus-nucleus collisions can be safely ignored when calculating the nucleon interaction lengths inside the projectile. (Henceforth we drop the subscript p on E p except when essential.) We consider two types of 'all-proton' spectra, one where the third extra-galactic population contains contributions from all 5 nuclear groups, and another where only protons contribute, denoted respectively by H3A and H3P [43]. These parameterisations have been extended [38,39] to include both additional heavy nuclear species, H14a, and to include a fourth population, H14b.
All the above spectra are compared in figure 1 with the flux rescaled by E 3 so that the region above the 'knee' is a horizontal line for the BPL spectrum and the difference from the other more recent parameterisations is emphasised. The latter are similar up to ∼ 10 8 GeV, but differ significantly thereafter. Although now outdated, the results with the BPL spectrum are required for comparison with older calculations of the prompt neutrino flux.

The cascade equations and their solution
We now review the cascade formalism in the framework of the Z-moment approach [14,45] which is used to simulate the propagation of high energy particles and their decay products through the atmosphere. The aim is to solve a series of coupled differential equations dependent on the slant depth X(l, θ) measuring the atmosphere traversed by a particle: where ρ is the density of the atmosphere dependent on the distance from the ground l (along the particle trajectory) as well as on the zenith angle θ. We adopt an isothermal model of the atmosphere, appropriate for atmospheric depths 10-40 km within which the bulk of particle production occurs: The horizontal depth of the atmosphere is X 3.6 × 10 4 gm cm −2 while its vertical depth is 1.3 × 10 3 gm cm −2 . As in previous calculations, we are concerned with small angles about the vertical, θ 0, where the conventional atmospheric neutrino flux arising from the decays of charged pions and kaons is the smallest.
For a particle of species j and energy E j that has traversed a slant depth X, the cascade equation for the corresponding flux φ j (E j , X) is where λ j is the interaction length, λ dec j is the decay length, and S kj are '(re)generation functions' describing the production of particle j from particle k (where the sum includes k = j). This says that as a particle traverses the atmosphere, its flux will decrease when the particle undergoes an interaction (thus losing energy) or decays, as well as increase from the decay or interaction of other particle species. The (re)generation function is: where dn(k → j; E k , E j ) is the differential transition rate between particle species k and j.
Assuming that the particle flux factorises into components dependent respectively on the energy E and the slant depth X, it can be rewritten more simply as with the key property that the moment Z kj , is independent of the slant depth X (which cancels in the ratio of fluxes). Under this factorisation assumption, the cascade equations describing the flux of the various relevant species (protons p, mesons m, and leptons l) as they propagate through the atmosphere can be written as a set of coupled differential equations: 3 where in the last equation the sum is restricted to the charmed hadrons that contribute to the prompt flux. Here d m (E) = cβγτ m is the decay length of a particle with proper lifetime τ m .
Although the solution of these equations is in general quite involved, there are simple asymptotic solutions. The first equation for the proton flux (2.11) can be trivially integrated to give where we have defined the nucleon attenuation length as This depends in general on the nucleon's interaction length in the atmosphere λ p (E): where A = 14.5 is the average atomic number of air molecules, N 0 is Avogadro's number, and the total inelastic air-nucleon cross-section is denoted by σ pA . Concerning the total proton-air cross-section, several parameterisations are available [46][47][48][49][50]. We use the QGSJet0.1c model [47] which fits available data well through the relevant energy range, including recent measurements made at the LHC [51] and the Pierre Auger Observatory [52]. The prompt neutrino flux in fact depends very weakly on the modelling of σ pA (E) [24].
Given eq. (2.14) the cascade equation (2.12) for the meson flux can be solved by neglecting, in the low energy limit, the interaction and regeneration terms and, in the high energy limit, the decay terms. This yields two asymptotic solutions: where the dependence on the energy of the attenuation lengths Λ p and Λ m is implicit.
Because of the additional dependence on the decay length in the low energy solution, these fluxes effectively scale with the proton flux (2.14) as follows: These intermediate solutions for meson fluxes are subsequently used as inputs in the corresponding low and high energy solutions for the leptonic decay of a meson m → l, with either l = µ or ν. The final vertical flux of leptons expected at the detector can then also be described by two asymptotic solutions, valid in the low-and high-energy regimes respectively: where m is a critical energy below which the probability of a meson to decay is greater than it is to interact: The smaller the critical energy, the longer the decay length, hence the more energy the particle will lose by interactions in the atmosphere before it actually decays. In eqs. (2.21) and (2.22), φ l,m represents the flux of lepton l from the decays of the meson m.
Pions and kaons have a critical energy of O(10 2 ) GeV. However, heavy quark mesons, such as B and D, are characterised by much larger critical energies of O(10 7 ) GeV and therefore mainly decay before losing energy in interactions with the atmosphere. This is why the 'prompt' neutrino flux from their decays is expected to dominate over the 'conventional' flux from π, K at high energies. The contribution of B mesons is usually neglected (despite a larger critical energy as compared to D mesons) because the b-quark pair production cross-section is < 10% of that of c-quark pairs up to ∼ 100 PeV. However at such high energies, the contribution from D mesons would be damped as they start interacting before decaying, hence we show the prompt neutrino flux only up to 10 7.5 GeV, The final step in solving the cascade equations in the Z-moment approach is the geometrical interpolation of the low-and high-energy asymptotic solutions (2.21) and (2.22) which yields the final expression for the prompt lepton (neutrino) flux: .

(2.24)
In the sum over mesons contributing to the prompt flux we include (the leptonic decays of) D 0 ,D 0 , D ± , D ± s and Λ ± c . In fact D 0 ,D 0 , and D ± account for the bulk of charm production in the atmosphere, with the other mesons contributing only around 10%.

Computation of the Z moments
We need to compute the various Z-moments in order to evaluate the prompt flux (2.24), of which the most crucial is the nucleon to meson moment, Z pm , which depends on the charm production cross-section in pp collisions. We now discuss how to compute Z ml , Z pp , Z mm and Z pm , and the corresponding theoretical uncertainties.
The generic Z-moment (2.10) defined earlier can be written for a (re)generation moment accounting for the interaction of a proton or a meson with an air nucleus, as while the decay moments that account for the leptonic decays of mesons are given by Here the differential distributions dn(i → f ; E , E)/dE correspond to the number n of final state particles f with energies between E and E + dE produced in an interaction where the initial state particle has energy E .
We now outline how each of the moments has been computed in this work: • For the calculation of the leptonic decay moment Z ml (2.26), we use the fact that the energy distribution of leptons from charmed meson decays obeys a scaling law: where F m→l (E) is the energy spectrum of the lepton l from the decay of the meson m, computed in the rest frame of the latter. Defining the scaling variable x E = E/E , we obtain Exploiting the fact that the meson flux φ m (E m ) scales in the high-energy (2.19) and lowenergy (2.20) limits we find that for the BPL cosmic ray spectrum, the leptonic moment reduces to a relatively simple expression where the exponent is β low = {1.7, 2} in the low-energy solution, and β high = {2.7, 3} in the high-energy solution, for energies below and above the 'knee' respectively. The calculation of the leptonic energy spectrum F (x E ) from charmed meson decays is performed with the Pythia8 [53] event generator and a boost is applied to transform F (x E ) from the laboratory to the meson rest frame. For the leptonic branching fractions of charmed mesons, we use the Particle Data Group recommended values [54] for inclusive decays: B(D ± → ν l X) = 0.161, B(D 0 → ν l X) = 0.065, B(D ± s → ν l X) = 0.065, and B(Λ ± c → ν l X) = 0.028. These values are adopted for both muon and electron neutrinos. The uncertainty on the branching fractions is well below 10% for D 0 and D ± , which are the most abundantly produced hadrons due to their large fragmentation fractions. Our result for the Z ml moments using the BPL cosmic ray spectrum are quite consistent with those reported earlier [17].
In figure 2 we compare the low-energy solution for the leptonic moment Z low ml (E) using the BPL cosmic ray spectrum for the four charmed mesons that contribute to the prompt flux, and where a sum over charge conjugate states is understood. Note that the decays of D 0 and D ± contribute the bulk of the prompt leptonic flux. We also show a comparison of Z low ml (E) for D ± only, using the different parameterisations of the cosmic ray flux to illustrate the large variations.
• When calculating the regeneration moments Z pp and Z mm that account for the interactions of protons and mesons in the atmosphere yielding a final state containing the same particle species, we follow previous studies [20,24] in adopting scaling laws for the protonproton and meson-proton cross-sections, viz.  where, as before, x E = E/E is the fraction of the original energy retained by the incoming particle after interaction with an air nucleus in the atmosphere and the exponents are n 1 = 0.51 and n 2 = 1.0 Eq. (2.31) is based on the approximation that the cross-section for charmed meson scattering off nucleons can be related to the corresponding kaon-proton cross-section. The attenuation length of charmed mesons will be given under the same approximation as [15] Λ m (E) where the dependence of the kaon-proton cross-section on energy is from [54].
In figure 3 we compare the proton and meson regeneration moments Z pp (E) and Z KK (E) for the 5 cosmic ray flux parameterisations. As for the leptonic moments, differences become appreciable only at high energies above the 'knee'.
• Finally we discuss the calculation of the proton-meson moment Z pm , which is the main ingredient of the present work, as it contains the information on charm production in high energy collisions. The number distribution can be related to the differential charm production cross-section as: We assume that the charm production cross-section scales with the mean atomic number of air A as compared to the corresponding pp cross-section where D is a generic charmed meson. This approximation is justified since even for forward D production in pP b collisions, the nuclear modification of the differential D hadron crosssection results in a suppression of at most 10% [55]. Although such effects are expected to increase in strength with atomic number, it is reasonable to ignore them when air is the target. This approximation is also supported by recent B production data on pP b collisions at the LHC [56] which show no evidence for nuclear modification effects.
Since we assume that ratios of fluxes are independent of the slant depth X to first approximation, we can write down a simplified version of the Z pm moment (2.25) in terms of the charm production cross-section as follows: Our calculation of the differential charm production cross-section in pp collisions at high energies has been discussed in detail earlier [57]. As explained therein, it is based on perturbative QCD as implemented in the NLO Monte Carlo event generator POWHEG [35], benchmarked with the corresponding FONLL [32] and aMC@NLO [36] calculations. The input PDF set is NNPDF3.0+LHCb [29,57], which includes the constraints on the small-x gluon from the LHCb 7 TeV charm production cross-sections. The parton showering and fragmentation are modeled with Pythia8 [53] using the Monash 2013 tune [58]. This is consistent with the semi-analytical fragmentation implemented in FONLL, tuned to LEP data [59]. For the fragmentation probabilities, which describe the transition f (c → D) for the different types of charmed mesons, rather than using the default Pythia8 values we use the recent LHCb measurements [28]: .080, and f (c → Λ c ) = 0.094. Using this framework, we have computed the moment Z pm (E) for a wide range of energies, from 10 3 to 10 7.5 GeV. This requires the calculation of the charm production cross-section for incoming proton energies up to E p = 10 10.5 GeV in the laboratory frame. 4 The POWHEG calculation is done in the center-of-mass frame for a wide range of √ s values, then boosted to the laboratory frame. In each case we have computed all the associated theoretical uncertainties from missing higher-orders, PDFs, and from the value of the charm mass [57] as follows:  • The charm quark pole mass is varied as m c = (1.5 ± 0.2) GeV, • Renormalisation and factorisation scales are varied independently by a factor of 2 around the central scale µ 0 = p T + m 2 c , with the constraint 1/2 < µ F /µ R < 2.
• PDF uncertainties are included at 68% CL, • Finally the total theory uncertainty is obtained by addition in quadrature of these 3 components so may be considered as a crude '1σ' band.
As with the other moments, the calculation of Z pm (E) is performed for all 5 cosmic ray flux parameterisations. In figure 4 we show the central theory prediction for Z pm (E) for the BPL spectrum for the four relevant charmed mesons (left plot) and then, for the D 0 andD 0 mesons only, using all parameterisations (right plot).

Results
This section contains our main result, the updated calculation of the prompt neutrino flux. We discuss its dependence on the various inputs, in particular the adopted cosmic ray flux parameterisation and PDF set used. We compare our result with other recent calculations and also provide the spectral index of the prompt flux as a function of energy. Figure 5 shows the prompt neutrino flux up to 10 7.5 GeV using the BPL cosmic ray spectrum. Since PDF uncertainties have been substantially reduced using the LHCb data, the error band is dominated by the 'scale uncertainties' of the NLO perturbative QCD calculation which can be reduced only when the corresponding NNLO result is available [60]. However at energies above a PeV, PDF uncertainties still make an important contribution to the total error band. We also show the central value of the ERS calculation [16], which has been used as a benchmark in several IceCube analyses but is now in tension with the 90% CL upper limit labeled '0.54×ERS' [26]. The central value of our calculation is a factor of 2 smaller, and just below the IceCube limit on the prompt neutrino flux. Note  Figure 5. The prompt neutrino flux using the BPL cosmic ray spectrum as input. The error band includes all relevant sources of theoretical uncertainties: from PDFs (68% CL), missing higher orders and the charm mass, as discussed in the text. The ERS benchmark calculation [16] is shown for comparison, as is the recent 90% CL upper limit on the prompt flux from IceCube [26]. Conv.

The prompt neutrino flux
Prompt Figure 6. The prompt neutrino flux and its uncertainty using the H3A cosmic ray spectrum as input, compared to the conventional neutrino flux at IceCube [61].
that this limit should be interpreted with some care, since it depends e.g. on the specific parameterisation of the cosmic ray flux in the analysis.
In figure 6 we compare the prompt neutrino flux with the conventional neutrino flux from the decays of pions and kaons, using the same cosmic ray spectrum (H3A). We use the updated calculation [11] for the South Pole location as implemented in the code NeutrinoFlux used by the IceCube collaboration [61] . Whereas for the conventional flux the location of the experiment is important (as this determines the geomagnetic rigidity cut-off which filters incoming cosmic rays), this is irrelevant for the prompt flux which arises from the interaction of much higher energy cosmic rays). The cross-over energy where the prompt component begins to exceed the conventional one is about 4 × 10 6 GeV.
As discussed in §. 2.1, an essential component of any calculation of the prompt neutrino flux is the parameterisation of the incoming cosmic ray flux, which is rather uncertain at the relevant high energies. Since cosmic rays with energies (100 − 1000)E ν contribute to the prompt neutrino flux at a given E ν , a prediction of the prompt flux up to 10 7.5 GeV requires knowledge of the cosmic ray flux up to at least 10 10.5 GeV.  Figure 8. Comparison of our baseline calculation and its uncertainty using the NNPDF3.0L set [57], with the corresponding central results using other PDFs as input: ABM11 [62], CT14 [63], HERAPDF1.5 [64] and MMHT14 [65]. All calculations assume the BPL cosmic ray spectrum.
NNPDF3.0L), with the central prediction obtained using other PDF sets: 5 ABM11 [62], CT14 [63], HERAPDF1.5 [64] and MMHT14 [65], in all cases at NLO. For each PDF set, the POWHEG calculation has been set up to include the required scheme modification terms; for instance, when n f = 5 PDFs are used as input, the scheme transformation terms from n f = 3 to n f = 5 are included [55,57].
Results for the prompt flux using different PDF sets and the BPL cosmic ray spectrum are shown in figure 8 where the total theory uncertainty is shown for NNPDF3.0L only. All PDF sets yield results in good agreement, except for MMHT14 which yields a substantially larger flux at energies above 10 5 GeV.
Thus the choice of PDF set is (with the exception of MMHT14) not important for the central value of the calculated flux. However it should be emphasised that the theory uncertainty band, which is shown here only for NNPDF3.0L, would have been much larger had we not included the LHCb charm hadroproduction data to reduce the uncertainty in the small-x gluon [57]. Thus our estimate of the uncertainty in the prompt neutrino flux is more robust than for all other calculations to date, and accordingly we advocate its use for inferring a lower limit which can be used as a prior in analyses of experimental data.

Comparison with previous calculations
In figure 9 we compare our result with the central values from the ERS [16], BERSS [20] and GMS [24] calculations, all using the BPL cosmic ray flux as input. The relative differences would change only mildly if a different cosmic ray flux parameterisation was used as input.
The central values of these three previous calculations are contained within the total theory uncertainty band of our result. Our central value is close to BERSS, but system-  Figure 9. Comparison of our calculation (GRRST) with the central values from ERS [16], BERSS [20] and GMS [24], all calculated using the BPL cosmic ray spectrum.
atically smaller than GMS, while the benchmark ERS result is at the upper end of the theory uncertainty band. Note that the BERSS calculation is based on the CT10 NLO PDF set [66] while the GMS calculation uses the ABM11 PDF set [62], neither of which incorporate the recent LHCb charm hadroproduction data. The ERS calculation was not based on pQCD at all, but the empirical 'colour dipole model'. It is evident that there is now some stability in calculations of the prompt neutrino flux and that in particular a theoretical lower limit can be set (subject of course to the large systematic uncertainty in the parameterisation of the incoming cosmic ray flux).

Spectral index of the prompt neutrino flux
It is useful to extract the local spectral index of the prompt neutrino flux, defined as: in order to compare with the standard expectation that γ 2.7. Both are shown in figure 10 which illustrates that above 10 5 GeV the naïve scaling is not obeyed. The BPL, H3P and H3A cosmic ray fluxes all yield a a prompt neutrino spectrum which falls off more steeply, while with the H14a and H14b fluxes a harder spectrum is obtained (it is worth keeping in mind that at very high energies, above ∼ 50 PeV, charmed mesons too will begin to lose energy by interaction with air nuclei before decaying, and at this point the fall-off of the prompt neutrino flux with E ν will start to become similar to that of the conventional flux.).
This indicates that a extraction of the prompt flux from a fit to data (including both the conventional flux and a cosmic signal) requires the full calculation of φ ν (E ν ) as a prior, with the overall normalisation left free but bounded by the total uncertainty band shown in figure 5. At a minimum, the lower limit on the prompt neutrino flux should be used as a prior, rather than allowing it to be zero as in current analyses [26].

Outlook
We have presented predictions for the flux of prompt neutrinos arising from the decays of charmed mesons produced in the collisions of high energy cosmic rays in the atmosphere. Our calculation of charm production at high energy makes extensive use of NLO Monte Carlo event generators and parton distribution functions. The novelty of our approach is that it has been validated with the 7 TeV charm cross-sections measured by the LHCb experiment, and found to be consistent with the more recent 13 TeV measurements.
As input we have used the NNPDF3.0+LHCb PDF set, where the inclusion of the LHCb 7 TeV data substantially reduces the PDF uncertainties in the small-x gluon. We include theory uncertainties arising from PDFs, missing higher-orders, and the value of m c .
We have studied the dependence of our result on the choice of input cosmic ray fluxes, including the most recent parameterisations, and on the choice of input PDF set. Our predictions have been compared with other calculations, in particular with ERS [16], BERSS [20] and GMS [24]. All three calculations are within the uncertainty band of our result, though our central value is the lowest. Our result is just consistent with the current experimental upper limit, suggesting that the prompt neutrino flux will be detected soon.
Our result for the prompt neutrino flux φ ν (E ν ) and its uncertainty, evaluated for 5 different input cosmic ray flux parameterisations, is available in terms of a C++ interpolation code from https://promptnuflux.hepforge.org. The interpolation tables can be used for neutrino energies between 10 3 and 10 7.5 GeV.
Since our calculations of charm hadroproduction have been validated with LHCb data, failure to detect the prompt neutrino flux in the range indicated would imply a flaw in the input assumptions, e.g. the cosmic ray flux parameterisation or possibly the Z-moment approach itself (e.g. the scalings in eqs. (2.30)-(2.32)). The latter can only be addressed by performing a full Monte Carlo simulation as in [17] with updated QCD tools and data.