Improved constraints on parton distributions using LHCb, ALICE and HERA heavy-flavour measurements and implications for the predictions for prompt atmospheric-neutrino fluxes

The impact of measurements of heavy-flavour production in deep inelastic ep scattering and in pp collisions on parton distribution functions is studied in a QCD analysis at next-to-leading order. Recent combined results of inclusive and heavy-flavour produc- tion cross sections in deep inelastic scattering at HERA are investigated together with heavy-flavour production measurements at the LHC. Differential cross sections of charm- and beauty-hadron production measured by the LHCb collaboration at the centre-of-mass energies of 5, 7 and 13 TeV as well as the recent measurements of the ALICE experiment at the centre-of-mass energies of 5 and 7 TeV are explored. These data impose additional constraints on the gluon and the sea-quark distributions at low partonic fractions x of the proton momentum, down to x ≈ 10−6. The impact of the resulting parton distribution function in the predictions for the prompt atmospheric-neutrino fluxes is studied.


Introduction
The fundamental structure of the nucleon is described by the theory of strong interactions, quantum chromodynamics (QCD). In the collinear factorisation, the nucleon structure is expressed in terms of parton distribution functions (PDFs), defined as probability densities for partons to carry a fraction x of the nucleon momentum at a factorisation scale µ f . While the scale evolution of the PDFs is calculated in perturbative QCD (pQCD) using the DGLAP equations [1][2][3][4][5][6][7], the x dependence must be constrained from the experimental measurements. The constraining power of experimental data on particular parton distributions is to a large extent defined by the acceptance of the experiment. Measurements of neutral current (NC) and charged current (CC) cross sections in deep inelastic scattering (DIS) at HERA [8] probe the x range of 10 −4 < x < 10 −1 , impose most significant constraints on the light quark PDFs and probe the gluon distribution via scaling violations. Additional constraints on the flavour separation of the quark sea and on the gluon distribution at low and high x are obtained by using the measurements from fixed target experiments and in proton-(anti)proton collisions. Heavy-flavour production in proton-proton (pp) collisions at the LHC is dominated by gluon-gluon fusion, therefore corresponding measurements probe the gluon distribution directly [9][10][11][12]. The measurements of forward charm [13] and beauty [14] production by the LHCb experiment at the centre-of-mass energy √ s = 7 TeV were used for the first time by the PROSA collaboration [9] to improve constraints on the gluon distribution at 5 × 10 −6 < x < 10 −4 , in the region hardly covered by any other measurements to that date. The resulting PDFs (PROSA 2015) were further used to predict the prompt neutrino flux from the decays of charmed mesons produced via cosmic ray interactions in the Earth's atmosphere [15], which constitute an irreducible background JHEP04(2020)118 in searches for the extraterrestrial neutrino flux by IceCube. The first QCD analysis [9] using forward heavy-flavour production at the LHCb has triggered attention of several PDF groups, which confirmed the significant improvement of the constraints on the gluon distribution by using forward heavy-flavour production measurements obtained with the LHCb experiment. In particular, one group has updated the PDF set NNPDF3.0 [16] by including the LHCb charm measurements into the NNPDF3.0 fit using the Bayesian reweighting method and providing the NNPDF3.0+LHCb PDF set [10]. Further LHCb measurements collected at √ s = 5 TeV, 7 TeV and 13 TeV were included in the PDF studies of refs. [11,12], where the normalised distributions and cross-section ratios between measurements at different centre-of-mass energies were used, following the suggestion of ref. [17], to quantify the impact of the LHCb D-meson measurements on the NNPDF3.0 and NNPDF3.1 PDF sets. In all analyses [10][11][12], the pQCD predictions for charm production in the forward region based on the FONLL approach [18] were used.
Recent improvements in the precision of the HERA measurements [8,19], new experimental data on heavy flavour production at the LHC at different √ s [20][21][22][23], together with new developments in the theory and improvements of the phenomenological tools, offer possibilities for stronger constraints on the gluon distribution at low x. These improvements in experimental measurements and the theory are explored in the QCD analysis presented in this paper, which updates the earlier PDF result [9]. The results, referred to as PROSA 2019, are used to update the predictions for the prompt atmospheric-neutrino fluxes.

Input data sets and used theory predictions
The main objective of the present QCD analysis is to demonstrate the constraining power of the updated measurements of heavy-flavour production in DIS and pp collisions for the determination of the PDFs of the proton. The QCD analysis is performed at next-toleading order (NLO) using the xFitter framework [24]. The updated combinations of the inclusive DIS cross sections [8] and of charm and beauty production cross sections [19] are used together with the measurements of charm and beauty hadroproduction in pp scattering at the LHC. The latter include the measurements of charm hadroproduction by the LHCb collaboration at √ s = 5 TeV [21], 7 TeV [13] and 13 TeV [20], and by ALICE at √ s = 5 TeV [23] and 7 TeV [22]. The measurements of beauty hadroproduction by LHCb at √ s = 7 TeV [14] are also used. The cross sections measured by LHCb and ALICE in each p T range are normalised in rapidity y, is the cross section in the central LHCb rapidity bin, 3 < y < 3.5. This bin is chosen for normalisation since all relevant LHCb data sets have a measurement in this rapidity range. Choosing a different bin for normalisation would not change the results. When normalising cross sections in this way, ALICE measurements at |y| < 0.5 are divided by the LHCb cross-section measurement in 3 < y < 3.5. The advantage of using the normalised cross section, demonstrated in the earlier PROSA analysis [9], is a significant reduction of the scale dependence in the theoretical prediction, while retaining the sensitivity to the PDFs, as illustrated in figure 1 Figure 1. NLO QCD predictions for D 0 meson production at LHCb at √ s = 13 TeV (top row), cross section ratios at √ s = 13 TeV and √ s = 5 TeV, denoted as R 13/5 (middle row), and normalised cross sections, denoted as R y (bottom row) with their PDF uncertainties obtained using the HERA-PDF2.0 PDF set with three active flavours (HERAPDF20 NLO FF3A EIG and HERAPDF20 NLO FF3A VAR LHAPDF sets) and different scale choices for two representative p T regions (the cross sections in the top right panel are scaled by a factor of 200 for presentation purposes). The lower panels indicate the ratio of the predictions to the nominal one. figure, this observable exhibits a smaller scale dependence than ratios of cross sections at different centre-of-mass energies which were used in refs. [11,12].
In the presented QCD analysis, bin-to-bin correlations in the input measurements are taken into account as described in the following. The treatment of correlated experimental uncertainties for the HERA data follows that of the original publications [8,19].

JHEP04(2020)118
The correlated uncertainties in the ALICE and LHCb measurements reported in the original publications [13,14,[20][21][22][23] and listed in the respective tables as exact uncertainty values in each kinematic bin in p T and y, are treated as fully correlated, and the uncorrelated uncertainties are obtained by subtracting the correlated ones from the total uncertainties, in quadrature. Because of this treatment of systematics, most of these correlated systematic uncertainties cancel in the calculation of the normalised cross sections. In case of the LHCb cross section ratio measurements, the uncertainties cancel completely. Further systematic uncertainties, reported as error intervals, see e.g. table 2 of ref. [21], are assumed uncorrelated, since no details about their size in individual p T and y bins are provided. For different final state measurements within one experiment, the tracking and luminosity uncertainties are treated as correlated. Furthermore, all experimental uncertainties are treated as uncorrelated among measurements at different centre-of-mass energies. The uncorrelated uncertainties in the normalised cross sections are propagated as correlated uncertainties to the respective complementary rapidity bins. It is worthwhile to note, that the details of the experimental uncertainties and their correlations in each individual kinematic range is of great importance and therefore most detailed information about the systematic correlations in the experimental measurement is required.
In the presented QCD analysis, the scale evolution of partons is calculated through the DGLAP equations at NLO, as implemented in the QCDNUM programme [25]. The description of the inclusive HERA data in the PDF fits improves in the kinematic range of small x and low virtuality Q 2 , by including higher twist effects [26,27] or, alternatively, small x resummation [28,29]. These upgrades are left for future analyses, once all necessary theoretical ingredients have become available. The changes in the PDFs when varying the Q 2 min = 3.5 GeV 2 cut imposed on the HERA data, 2.5 ≤ Q 2 min ≤ 5.0 GeV 2 , are found to be small with respect to other uncertainties. Therefore we are confident that inclusion of higher-twist terms would not modify the results of this analysis in a substantial way.
The theoretical predictions for the heavy-quark and inclusive HERA data are obtained using the OPENQCDRAD [30] code in the fixed-flavour-number scheme (FFNS) with three active flavours in the proton using the MS mass scheme, following ref. [19]. Similar to the earlier PROSA analysis [9], the theoretical predictions for the fully differential heavy-quark hadroproduction in pp collisions, available at NLO in FFNS, are used. These are calculated using the MNR code [31], with the single-particle inclusive distributions computed using the pole mass scheme for the heavy quarks, and translated into the MS mass scheme expressions using the MS mass m Q (m Q ) and following ref. [32]. The MS mass scheme is then consistently used in the calculations for all used processes.
The factorisation and renormalisation scales are chosen to be Q for inclusive DIS, and µ r = µ f = Q 2 + 4m Q (m Q ) 2 for heavy quark production in DIS, respectively, with m Q (m Q ) representing the heavy-quark mass in the MS scheme. For heavy quark produc- T is assumed. The calculations for heavy quark hadroproduction are supplemented with phenomenological non-perturbative fragmentation functions to describe the transition of heavy quarks into hadrons. The fragmentation of charm quarks into D mesons is described by the -4 -

JHEP04(2020)118
Kartvelishvili function D Q (z) ∝ z α K (1 − z) with α K = 4.4 ± 1.7 as measured at HERA [33,34], and for the fragmentation of beauty quarks α K = 11 ± 3 is used as measured at LEP [35], following the previous PROSA analysis [9]. Studies of the uncertainties related to the fragmentation in ref. [36] for a determination of the charm-quark mass in the MS scheme from deep inelastic scattering at HERA data have shown that the dominant effect is captured by varying α K within its uncertainties. This treatment of charm quark fragmentation is independent of the choice of a particular renormalisation scheme for the heavy quark mass. The latter is needed in a determination of the initial condition for the perturbative heavy quark fragmentation function, which is known to NNLO [37]. The subsequent range of evolution in the case of charm quark fragmentation into D mesons from the scale of hadronisation to scales of order of the charm quark mass is very short, so that the modelling with the non-perturbative Kartvelishvili function D Q (z) is justified.
The main QCD analysis is performed in the FFNS and the sensitivity of the heavy quark measurements to the PDFs and to the masses of the charm and beauty quarks is fully explored by treating m c (m c ) and m b (m b ) as free parameters in the fit. The fit is also performed in the variable flavour number scheme (VFNS) to allow for incorporation in shower Monte Carlo event generators and applications in e.g. underlying event tuning at the LHC.

PDF parametrisation
The PDFs are parametrised at the starting evolution scale of µ 2 f 0 = 1.9 GeV 2 , similar as in ref. [8] and ref. [38], as follows: Here, xg(x), xu v (x) and xd v (x) represent the gluon, up and down valence quark distributions, respectively. The sea quark distribution is defined as xΣ(x) = xu(x) + xd(x) + xs(x), with xu(x), xd(x), and xs(x) denoting the up, down, and strange antiquark distributions, respectively. For the up-and down-type antiquark distributions, xU(x) and xD(x), relations xU(x) = xu(x) and xD(x) = xd(x) + xs(x) are assumed. The normalisation parameters A uv , A dv , and A g are determined by the QCD sum rules. The strangeness fraction f s = xs/(xd + xs) is fixed to f s = 0.4 as in the HERAPDF2.0 analysis [8]. Additional constraints B U = B D and A U = A D (1 − f s ) are imposed to ensure the same normalisation for the xu and xd distributions as x → 0. The term F g log x was proposed in [38] to provide a flexible functional form at low x and replace the 3-parameter extra term in ref. [8].
The predicted and measured cross sections together with their corresponding uncertainties are used to build a global χ 2 , minimised to determine the initial PDF parameters. The χ 2 definition follows that of eq. (32) in ref. [8]. In the minimisation, performed using the MINUIT package [39], the experimental uncertainties in the heavy-quark normalised JHEP04(2020)118 cross sections are treated as additive, and the treatment of the experimental uncertainties for the HERA DIS data follows the prescription given in ref. [8].
The parameters in eq. (3.1) are selected by first parametrising each PDF as and setting all D and E parameters to zero. Additional parameters in each resulting PDF are included in the fit one at a time. The improvement in χ 2 of the fits is monitored and the procedure is stopped when no further improvement is observed. The inclusion of the F g parameter does not lead to significant change in χ 2 , in particular, its fitted value is consistent with 0 within uncertainty, however the variation of F g significantly affects the fit uncertainties.
To ensure that the gluon PDF at low x is not over-constrained in the fit, different functional forms in the parametrisation were tested, as used in the ABMP16 [27], CT14 [40], HERAPDF2.0 [8] and Bonvini-Giuli (BG) [38] PDF fits: 1 These functional forms are characterised by 3 (HERAPDF2.0 no 'flexible' g), 4 (ABMP16) or 5 (CT14, HERAPDF2.0, BG) parameters controlling the gluon PDF, cf. 4 parameters in the presented nominal parametrisation of eq. (3.1). The resulting gluon distributions are presented in figure 2. The parametrisations of ABMP16, HERAPDF2.0 without the flexible gluon, and BG provide very similar results to that of the nominal parametrisation in eq. (3.1). Note that also the HERAPDF2.0 analysis considered the parametrisation without the flexible gluon, referred to as an 'alternative' gluon parametrisation [8], provided primarily for predictions of cross sections at very low x, such as very high-energy neutrino cross sections.
The fit using the HERAPDF2.0 and CT14 parametrisations yielded a gluon distribution with a sharp turnover to negative values at x ∼ 10 −6 , i.e. at the edge of the kinematic reach of the used measurements. Using such PDFs would lead to a negative prediction for the total charm hadroproduction cross sections at √ s 20 TeV, similar to the observation of ref. [26]. Therefore these parametrisations are discarded (despite they provide an improved χ 2 , by 22 and 7 units when using the HERAPDF2.0 and CT14 parametrisations, respectively). JHEP04(2020)118

PDF uncertainties
The PDF uncertainties are investigated according to the general approach of the HERA-PDF2.0 analysis [8], with the fit, model, and parametrisation uncertainties taken into account.
The fit uncertainties arising from the uncertainties in the measurements are estimated by using the Hessian method, adopting the tolerance criterion of ∆χ 2 = 1, and correspond to 68% confidence level.
To investigate the impact of model assumptions on the resulting PDFs, alternative fits are performed and the differences to the central result are considered as model uncertainties. The strangeness fraction is varied as 0.3 ≤ f s ≤ 0.5 and the value of Q 2 min imposed on the HERA data as 2.5 ≤ Q 2 min ≤ 5.0 GeV 2 . The FFNS strong coupling constant is assumed as 0.105 < α [42]). The variation of the fragmentation parameters α K = 4.4 ± 1.7 for charm hadrons [33,34] and α K = 11 ± 3 for beauty hadrons [35] is performed. The scales µ f and µ r for heavy quark production are varied independently and simultaneously up and down by a factor of two, excluding variations of the two scales in opposite directions. Note that for the normalised cross section predictions, the simultaneous variation of the µ f and µ r scales in the same direction results in the largest deviation in the resulting PDFs and is considered as one PDF uncertainty eigenvector.
The parametrisation uncertainty is estimated by extending the functional form of each PDF in eq. (3.1) with additional parameters D and E, see eq. (3.2), which are added or removed one at a time and do not impact the χ 2 . Furthermore, the shape of the gluon PDF is extended by adding a +G g log 2 x term [38]. This modification does not result in -7 -JHEP04(2020)118 an improvement in χ 2 and therefore is not considered in the nominal parametrisation. The variation of the starting scale, 1.6 < µ 2 f0 < 2.2 GeV 2 , is also taken into account as contribution to the parametrisation uncertainty. The parametrisation uncertainty is constructed at any given scale as an envelope built from the maximal differences between the PDFs resulting from all the parametrisation variations and the central fit at each x value.
The total PDF uncertainty is obtained by adding experimental, model, and parametrisation uncertainties in quadrature.

PROSA 2019 parton distributions
The quality of the overall fit can be judged based on the global χ 2 divided by the number of degrees of freedom, n dof . For each data set included in the fit, a partial χ 2 divided by the number of measurements (data points), n dp , is provided. The correlated part of χ 2 quantifies the influence of the correlated systematic uncertainties in the fit. The global and partial χ 2 values for each data set are listed in table 1, illustrating a satisfactory agreement among all the data sets (note that the χ 2 expression does not include theoretical uncertainties, but we account for them in the evaluation of the theory uncertainties associated to the PDFs). The central values and the uncertainties of the fitted PDF parameters are given in table 2. The fitted masses of the heavy quarks are m c (m c ) = 1.230 ± 0.031 GeV, m b (m b ) = 3.977 ± 0.100 GeV. These values are a bit lower than, but consistent with, those obtained from the HERA data only [19]. The corresponding full set of other potential systematic uncertainties was not evaluated here.
The resulting PROSA 2019 PDFs with their total uncertainties at the scale µ 2 f = 10 GeV 2 are shown in figure 3. These are compared to the result of the PROSA 2015 fit [9]. In figure 4 (left), the gluon distribution normalised to the one from the PROSA 2015 fit is shown. The two results are in a very good agreement and a significant improvement in the precision of the gluon PDF is achieved at x < 10 −4 , as compared to the PROSA 2015 fit. This improvement is attributed mainly to the presence of the new LHCb data in the fit which extends the coverage to lower values of x. The valence and sea quark PDFs are in good agreement with the result of the HERAPDF2.0 analysis [8] and the observed differences in these distributions to the PROSA 2015 analysis are attributed to the update of the DIS measurements [43] used in ref. [9] to the final combination [8] of the HERA data. It turns out that once the LHCb data are included in the fit, the inclusion of ALICE data provide no significant additional constraints on the PDFs, but since the ALICE data cover the central range of y, the consistent description of these data serves as a non-trivial self-consistency check of the fit using normalised cross sections.
The relative total, fit, model and parametrisation uncertainties for the gluon PDF are shown in figure 4 (right). The total uncertainties are dominated by the model uncertainties, with the largest contributions arising from the scale variations in predictions for heavy-quark hadroproduction. Reduction of these uncertainties would require theoretical calculations at higher order.

Fit in VFNS
The fit in the VFNS is performed using the APFEL library [45] interfaced to xFitter. The theoretical predictions for the HERA data are computed using the FONLL-B scheme [46] with the pole charm and beauty quark masses set to m pole c = 1.4 GeV and m pole c = 4.5 GeV respectively. However, no VFNS calculation for heavy-quark pp hadroproduction is interfaced to public QCD analysis tools like xFitter. To use the MNR calculations with the VFNS, the functionality of the APFEL library is exploited, allowing to choose arbitrary heavy-quark matching thresholds [47]. These thresholds are set as: The kinematic requirements p T < 5 GeV and p T < 16 GeV are imposed on the LHC charm and beauty data, respectively, to ensure that not more than 3 (4) flavours are considered when calculating predictions for charm (beauty) data when choosing µ r = µ f = Q 2 + 4m Q (m Q ) 2 . The strong coupling constant is set to α n f =5 s (M Z ) = 0.118 [42], while all other settings are the same as in the FFNS fit. The specific matching thresholds in eq. (5.1) are chosen to ensure that a sufficient amount of the LHC charm and beauty data is still included in the fit. The choice of the matching thresholds is arbitrary and coincides with the renormalisation scheme choice of ref. [47]. The results are proven to remain stable under variations of 3.1 ≤ µ Q /m pole Q ≤ 6, whereby the p T cuts in the charm and beauty cross-section measurements of the LHC are modified accordingly.
In the VFNS variant of the PDF fit, χ 2 = 2114 is obtained for n dof = 1714, indicating a similar quality of data description as compared to the fit in the FFNS (consistent with the observations in refs. [8,19]). The resulting PDFs are available in the LHAPDF format at the PROSA web-site [44]. No PDF uncertainties are provided with this set.
The performance of the PROSA 2019 VFNS PDFs is tested by computing predictions for the inclusive and multi-jet production in DIS [48][49][50][51][52] and jet [53] and top quarkantiquark production [54,55] in pp collisions. The results collected at the PROSA website [44] are found to be similar to those using HERAPDF2.0 PDF.
In figure 5 we compare the gluon distribution obtained in our fits in FFNS and VFNS with the NNPDF3.1 PDFs [56] and with the PDFs from ref. [12] obtained using NLO calculations which also exploit the LHCb measurements of D-meson cross-sections at 5, 7 and 13 TeV. The different fits are in good agreement.  -12 -

JHEP04(2020)118 6 Predictions for prompt atmospheric-neutrino fluxes
Various applications in high-energy astroparticle physics could benefit from accurate PDFs in the low-x region. One of the most interesting cases is the evaluation of the prompt flux of atmospheric neutrinos, originating from the semileptonic decays of heavy-flavoured hadrons produced in the interactions of cosmic rays (CR) with nuclei in the atmosphere. The prompt atmospheric-neutrino flux represents a relevant background for searches of highly energetic cosmic neutrinos, which are supposed to be produced in the vicinity of far astrophysical sources and in the Galactic Plane [57]. Such searches are conducted at Very Large Volume Neutrino Telescopes such as ANTARES [58], IceCube [59] and KM3NeT [60], which register and analyse the features of the track and cascade events induced by the charged-current and neutral-current weak interactions of the impinging neutrinos with the water/ice nuclei. To date, no direct measurement of the prompt atmospheric-neutrino flux is available. Therefore, the most precise theoretical predictions for these fluxes are needed for the reliable interpretation of the experimental data in order to disentangle the cosmic neutrino component from the atmospheric background [61]. In this paper, the predictions for the prompt atmospheric-neutrino fluxes are calculated, in general following the method detailed in ref. [15]. It is assumed, that pA and AA interactions leading to charm production can be described in terms of pp interactions (superposition model) in pQCD. For the proton structure description, the PROSA 2019 PDF fit is used among other PDFs. Production and decay of the D ± , D 0 ,D 0 , D ± s , Λ ± c in the atmosphere is considered dominant, since the contribution of other charmed hadrons, as well as b-flavoured hadrons, amounts to 5-15% of the dominant one [62]. In the computation of charmed-hadron production cross sections, the renormalisation and factorisation scales are chosen as µ R = µ F = µ 0 = p 2 T + 4m 2 c , consistent with the scale choice adopted in the theory predictions of D-and B-meson production at LHCb and ALICE used in the PDF fit. Note that this scale choice differs from the one of ref. [15] (PROSA 2015), where µ R = µ F = p 2 T + m 2 c was used, consistent with [9]. While the difference between the two scale choices reduces with increasing p T , at low p T the present scale choice is motivated by faster convergence of the perturbative series to NNLO for the total pp → cc + X cross section at the LHC energies, as reported in ref. [63].
In the present work, the central value of the pole mass of the charm quark, m pole c = 1.43 GeV is used, corresponding to m c (m c ) = 1.23 GeV in the PDF fit (see table 2), as obtained using 1-loop conversion. It is worthwhile to note that this value is somewhat larger then the one 2 used in the PROSA 2015 computation. The uncertainty due to the choice of the charm quark mass is evaluated by varying the pole mass by ± 0.15 GeV around the central value.

JHEP04(2020)118
The predictions for the prompt (ν µ +ν µ ) fluxes using PROSA 2019 PDFs are presented in figure 6. Those are obtained by using different hypotheses for the primary CR all-nucleon flux [64,65], which are derived from the measured CR all-particle spectrum [66], under specific assumptions for the CR composition. In particular, these assumptions concern the proton and nuclear groups included in the derivation of the spectra, their spectral indices, their rigidity, the number of populations of galactic origin and their sources, the presence or not of an additional proton population of extragalactic origin, as detailed in the aforementioned references.
The QCD uncertainties in the resulting prompt (ν µ +ν µ ) fluxes encompass the uncertainties in the charm quark mass, PDF and those related to the scale choice, with the latter being the dominant uncertainty. The quoted scale uncertainty bands are obtained at fixed PDF, i.e. they do not include the contribution related to scale variation in the PDF fit. The effect of varying the (µ R , µ F ) scales when comparing theoretical predictions with experimental data in the PDF fit process, according to the method detailed in section 4, is instead accounted for in the PDF uncertainty band. As expected, the uncertainty on prompt neutrino fluxes due to the variation of the charm quark mass value in the constant interval around its central constant value, decreases with energy: at small E ν, lab this uncertainty dominates over the PDF uncertainty, whereas at E ν, lab ∼ 10 7 -10 8 GeV, both uncertainty contributions become similar. At high E ν, lab , the PDF uncertainties are reduced with respect to those of the PROSA 2015 computation.
The different contributions of the PROSA 2019 PDF uncertainty in the flux prediction are shown in figure 7. All the uncertainties increase with increasing E ν, lab , which corresponds to the decreasing x of the target parton, probed.
Prompt neutrinos with energy E ν, lab are mostly produced by air collisions of CR protons with laboratory energies (10 -100) times larger. Therefore, neutrinos with energies of some PeV, i.e. the most energetic neutrinos seen so far by IceCube, are mostly obtained by collisions up to the LHC centre-of-mass energies. On the other hand, neutrinos with higher energies can be the result of collisions at energies not yet probed at accelerators. It is worthwhile to note that the PDF uncertainties for 10 6 < E ν, lab < 10 8 GeV are calculated assuming the PDFs can be extrapolated to x-values lower than the kinematic reach of the data used in the PDF fit, x ≈ 10 −6 . For an x-range lower than 10 −6 , there are no measurements probing this region to date. 3 However, the computation of the prompt neutrino flux at the highest E ν, lab energies involves a non-negligible contribution from initial state partons with x lower than 10 −6 , as shown in figure 4 of [67]. The agreement of the results based on the PROSA 2019 and the PROSA 2015 PDF sets can be considered as a consistency test of the extrapolation procedure, assuming no New Physics contribution in the probed x-range. Furthermore, at the neutrino energies of E ν, lab 10 5 GeV, the assumptions on the CR composition become very important (see figure 6, bottom righthand plot), having an impact on both the shape and the normalisation for the prompt atmospheric-neutrino flux. In particular, at the highest E ν, lab , corresponding to the low- 3 We recall that the kinematic formula relating the projectile/target parton x with the pT and y of a produced heavy-quark with mass mQ in pp → QQ collisions at a laboratory energy Ep is x =  est x values, the spread between central predictions obtained using as input different CR primary all-nucleon spectra amounts to a factor of O(5-10), that is much larger than the extrapolated PDF uncertainty.
In figure 8, predictions for prompt atmospheric-neutrino fluxes using different descriptions of the proton structure, are compared among each other. The predictions using FFNS PROSA 2019, PROSA 2015 and ABM11 PDFs with corresponding α S (M Z ), have been obtained using as a basis matrix elements for cc hadroproduction at NLO in FFNS (N f =3), matched, according to the POWHEG formalism [68,69], to the PYTHIA8 parton shower and hadronisation algorithms 4 [70]. Each pp collision can produce more than two D-hadrons, because charm quarks can be produced both in the hard-scattering and during the parton shower processes. The predictions using PROSA 2019 at high energies are somewhat lower than those using PROSA 2015 and ABM11 PDFs. In the same figure, the predictions obtained in the general-mass VFNS framework of ref. [71] (GM-VFNS), using as input VFNS PDFs (CT14nlo and the PROSA 2019 VFNS) are shown. The NLO QCD corrections are included in the partonic cross section, whereas the transition from partons to hadrons is described by fragmentation functions evolving with the factorisation scale [72], a procedure which resums logarithms of p T /m c at next-to-leading-logarithmic accuracy. Both central predictions using the GM-VFNS shown in the plot are compatible among each other, but show shape differences with respect to the FFNS ones. Part of these differences are related to the different treatment of the transition of partons into hadrons (parton shower + hadronisation on the one hand, vs. fragmentation functions on the other hand). Also, a different factorisation scale is used in the GM-VFNS predictions. 5 For comparison, the upper limit on the prompt neutrino flux obtained in the IceCube analysis [73] of up-going muons from the northern hemisphere is also shown and is well described by the predictions.
The various predictions shown in figure 8 were all obtained under identical assumptions for all inputs used in the solution of the cascade equations for the evaluation of prompt neutrino fluxes, except for the explained differences in the evaluation of D-hadron production. On the other hand, in figure 9 the presented flux prediction using the PROSA 2019 PDFs is compared to those obtained by other groups. The BPL primary CR all-nucleon spectrum [57] is used as input for this comparison because of its very simple form which has allowed an easy incorporation of this spectrum in the computation of many different authors. Although the general methodology for the calculation of prompt neutrino fluxes is the same, the calculation by different authors are obtained in a completely independent way and, thus, at least in principle, might differ in many respects, not only related to the methodology for the computation of charm hadroproduction, but even for other assumptions in the solution of cascade equations (e.g. the details of the atmospheric model, of the p-Air total inelastic cross-section and of the proton and hadron regeneration processes [63]).

JHEP04(2020)118
Notwithstanding these possible further sources of discrepancies, we observe that our predictions turn out to be consistent with those by other authors, within uncertainties. Due to this similarity, the experimental collaborations, in many of their works, limit themselves to consider as input only very few (if not only one) of the theoretical predictions available. The result of ref. [74] shows the largest differences with the presented result due to using the charm cross-section calculation at LO only. The ERS dipole model prediction [75], that is mostly used by the experimental collaborations in their data analysis, is also consistent with the prediction of this paper, within uncertainties. The uncertainties in the ERS prediction are smaller compared to the QCD-based prediction, however the way of the uncertainty estimate in both calculations can not be directly compared. Indeed the dipole approach is expected to effectively resum logarithmic contributions of the form α S ln(1/x) in the PDF evolution, a property which could lead to a reduction of the PDF uncertainties associated to the target parton at low x, whereas the resummation of these logarithms is not included in the DGLAP evolution. The uncertainties associated to the projectile parton distributions in ref. [75] were estimated by comparing two different central PDF sets, without considering the PDF uncertainty associated to each of those. Furthermore, the factorisation scale variation in ref. [75] is performed in a limited range of µ F = m c , µ F =2 m c .
In figure 10, the prediction for the prompt neutrino flux based on the superposition model for both the projectile CR and the target nucleon of the air, obtained using the PROSA 2019 proton PDF set, is compared to the calculation of ref. [62] which uses nuclear PDFs to describe the target nucleon (nitrogen) and the proton PDFs for the projectile CR. The H3p CR all-nucleon spectrum is adopted as an input in our calculation, to be consistent with the choice of ref. [62]. In general, nuclear PDF fits are at a less advanced stage of development with respect to the proton PDF fits, due to the fact that less experimental data on collisions involving at least one nucleus are available with respect to the proton case, and that the theory for describing these collisions is also less advanced, with persistent difficulties in disentangling the different possible sources of cold nuclear matter effects. Additionally, the study of p-A collisions at the Large Hadron Collider has been performed by mostly using Pb beams, while the atmosphere involves much lighter nuclei (N, O), which necessarily requires an important extrapolation. However, at present stage, it is remarkable to observe that predictions using proton PDFs and the superposition model turn out to agree with those using nuclear PDFs, at least within present uncertainty. This might point to the conclusion that the approximation of using proton PDFs and the superposition model, instead of nuclear PDFs (which, in principle, would be more appropriate because the air is made by nuclei instead of being made by unbound protons and neutrons) can still be considered as well justified, at least considering the present status of uncertainties.

Summary
In this paper, improved constraints on the parton distributions are presented, as obtained in a QCD analysis at NLO using DIS and pp collision data. In particular, the recent measurements of the LHCb and ALICE experiments of hadroproduction of charm and beauty-flavoured hadrons in different kinematic ranges (forward and central) provide ad-   Figure 6. Predictions for the prompt atmospheric-neutrino fluxes and their uncertainties related to scale variation, charm mass and PDF uncertainties. Each of the first five panels refers to a different CR primary all-nucleon spectrum (GST-3, GST-4, H3a, H3p [64,65], and Broken-Power-Law (BPL) [57], respectively), chosen among those most widely used in literature. In the bottom panel on the right-hand side the central predictions using all these primary all-nucleon spectra are compared among each other, and with those obtained with the more recently introduced Nijmegen [76] and Global Spline Fit (GSF) [77,78] all-nucleon spectra. At the highest E ν lab , the largest predictions (GST-4 and BPL) are seven time larger than the smallest one (H3a).   Figure 8. Predictions for prompt neutrino fluxes from this paper as compared to other predictions previously obtained by some members of our group. Predictions obtained with the PROSA PDF fit of ref. [9], using the same µ R = µ F scale, but a slightly different charm mass value (m c = 1.4 GeV) are shown by dotted (pink) lines; the predictions of ref. [63], using the ABM11 PDF fit [79], are shown by double-dotted-dashed (blue) lines. Finally the predictions using the VFNS version of the PDF fit described in this paper, in association with a GM-VFNS calculation of charm hadroproduction [71], are also reported and compared to those of ref. [71] itself that made use of the CT14nlo PDF fit [40]. The experimental upper limit on prompt (ν µ +ν µ ) flux, extracted from IceCube in the analysis of ref. [73], is also reported. The solid red line is the limit inferred by the available IceCube data, whereas the dashed line is its extrapolation to different E ν , still computed by the IceCube collaboration.   Figure 10. Predictions for prompt atmospheric-neutrino fluxes obtained in the presented analysis, compared to those by other authors [62]. The all-nucleon flux H3p CR is used for consistency with ref. [62]. For the description of the Nitrogen (air) target, the nuclear PDFs EPS09 [83] and nCTEQ15 [84] are used.