Confronting current NLO parton fragmentation functions with inclusive charged-particle spectra at hadron colliders

The inclusive spectra of charged particles measured at high transverse momenta ($p_T\gtrsim$2GeV/c) in proton-proton and proton-antiproton collisions in the range of center-of-mass energies $\sqrt{s}=200-7000$GeV are compared with next-to-leading order perturbative QCD calculations using seven recent sets of parton-to-hadron fragmentation functions (FFs). Accounting for the uncertainties in the scale choices and in the parton distribution functions, we find that most of the theoretical predictions tend to overpredict the measured LHC and Tevatron cross sections by up to a factor of two. We identify the currently too-hard gluon-to-hadron FFs as the probable source of the problem, and justify the need to refit the FFs using the available LHC and Tevatron data in a region of transverse momenta, $p_T\gtrsim$10GeV/c, which is supposedly free from additional non-perturbative contributions and where the scale uncertainty is only modest.


Introduction
The inclusive production of large-transverse-momentum (p T 2 GeV/c) hadrons at protonproton (p-p) and proton-antiproton (p-p) colliders provides a ground for testing the factorization theorem of Quantum Chromodynamics (QCD) [1,2] that predicts the universality and evolution of the two non-perturbative elements in the theoretical cross-sections: parton distribution functions (PDFs) and parton-to-hadron fragmentation functions (FFs). While the PDFs can nowadays be tested and constrained by a multitude of different processes in deep-inelastic scattering and hadronic collisions [3], the variety and kinematic reach of data by which the FFs can be determined is more limited [4,5]. In this context, the inclusive hadron measurements at the LHC extending to unprecedentedly large values of center-of-mass energy ( √ s) and hadron p T [6], are particularly useful for studying the FFs and their universality. From an experimental perspective, the best precision and widest kinematic reach in the p T -differential cross sections are achieved when no particle identification is performed. The pre-LHC measurements for unidentified charged-hadron spectra in p-p and p-p collisions range from fixed-target experiments at √ s 60 GeV [7,8] to collider experiments covering a wide range of center-of-mass energies √ s = 200 − 1960 GeV [9,10,11,12,13,14,15,16]. However, most of these data are restricted to moderate values of p T 10 GeV/c, and the accuracy for p T ≥ 10 GeV/c is rather poor. The importance of the new LHC data to overcome such limitations is underscored by the confusion [17,18,19] triggered by the original CDF data [16] which seemed to display deviations from next-to-leading order (NLO) perturbative QCD (pQCD) calculations up to three orders of magnitude at p T 150 GeV/c, but which was later on identified as an experimental issue (see the erratum of Ref. [16]).
Interestingly, the NLO pQCD predictions presented along with the recently published CMS [20,21] and ALICE [22] inclusive charged hadron spectra, appear to overshoot the data by up to a factor of two in the kinematical region where effects such as e.g. intrinsic transverse momentum of the colliding partons (intrinsic k T ) [23,24,25], soft-gluon resummation [26,27], or small-z instabilities of the FFs [28] should not play a major role. Especially since the recent LHC measurements of the p T -differential cross sections for inclusive jets [29,30] and prompt photons [31,32] are in perfect agreement with the NLO pQCD expectations [33,34,35], the data-vs-theory discrepancies for inclusive charged-hadrons come totally unexpected. Resolving such inconsistencies is also of relevance for other QCD analyses such as those related to the suppression of high-p T hadrons in ultrarelativistic heavy-ion collisions where the p-p spectra are required as baseline measurements [36].
In this paper, we present a systematic comparison of the theoretical predictions for unidentified charged-hadron production to experimental data with a special emphasis on the latest LHC measurements. Our aim is to demonstrate that such a process in hadron colliders is predominantly sensitive to the gluon-to-hadron FFs which are presently not well determined and, consequently, large discrepancies among the modern sets of FFs exist. These differences not only translate into a significant scatter in the corresponding predictions for the cross sections, but none of the current FF sets can consistently reproduce the current LHC and Tevatron data at p T 10 GeV/c. As the data measured by different experiments at the same collision energies are in mutual agreement, it seems excluded that such discrepancies are due to an experimental issue. Instead, this hints to a severe problem in the gluon-to-hadron FFs in most of the existing sets. We conclude that the gluon FF, which is currently mildly constrained by charged-hadron spectra from hadronic collisions at RHIC and SppS energies, should be refitted by using the LHC and Tevatron hadron spectra in the region p T 10 GeV/c, where the theoretical scale uncertainties appear tolerable and which should be free from additional, non-perturbative hadron production mechanisms.

The pQCD framework for inclusive hadron production
The cross section for the inclusive production of a single hadron h 3 with a momentum p 3 in the collision of two hadrons h 1 and h 2 carrying momenta p 1 and p 2 respectively, can be expressed, differentially in transverse momentum p T and (pseudo)rapidity η, as [24,37] In this expression, f h k i (x k , µ 2 fact ) denote the PDFs of the colliding hadrons evaluated at parton fractional momenta x k and scale µ 2 fact . We will use the ct10nlo [38] parametrization throughout this work. The parton-to-hadron FFs are denoted by D l→h 3 (z, µ 2 frag ) where z is the fraction of the parton momentum carried out by the outgoing charged hadron. The PDFs and FFs are convoluted with the partonic coefficient functions dσ for which we use their fixed-order NLO O(α 3 s ) expressions [37,39,40], treating the partons and hadrons as massless particles. In practice, we evaluate these cross sections employing the INCNLO [37,41] program which we have modified to improve the convergence at small values of p T . The fixed-order calculations are supposed to be adequate for p T 1 GeV/c but still sufficiently away from the phase-space boundary p max T ∼ √ s/2 (at midrapidity, η ≈ 0), where soft-gluon resummations [26,27] become relevant due to large logarithmic contributions from an incomplete cancellation of the infrared divergences.
Truncating the partonic coefficient functions to O(α 3 s ) leads to the well-known scale dependence of the pQCD calculations. For inclusive hadron production, there are three independent scales: the renormalization scale µ ren , factorization scale µ fact , and the fragmentation scale µ frag . The sensitivity of the computed cross sections to the variations of these scales is typically taken as an indication of the size of the missing higher-order corrections. Our default choice is µ ren = µ fact = µ frag = p T , and we take the scale uncertainty as the envelope enclosed by the following 16 scale variations [42] (2) We omit the combinations in which µ ren and µ fact or µ ren and µ frag are pairwise scaled by a factor of two in opposite directions due to the appearance of potentially large contributions of the form log(µ 2 ren /µ 2 fact ) and log(µ 2 ren /µ 2 frag ) in the calculation. The next-to-NLO (NNLO) calculations are expected to definitely reduce the scale dependence. However, although the PDF analyses can be nowadays carried out partly at NNLO level [43,44,45,46], the time-like splitting functions needed in the NNLO evolution of the FFs are not yet fully known [47,48], nor are the NNLO coefficient functions needed in Eq. (1), although the latter could finally emerge through the work currently carried out for jets [49,50,51]. Table 1 lists the seven commonly used sets of NLO parton-to-charged-hadron FF parametrizations: Kretzer (kre) [52], kkp [53], bfgw [54], hkns [55], akk05 [56] dss [57,58], and akk08 [59], that we employ in our calculations. In a few of these analyses the charged hadron FFs are constructed as a sum of the individual FFs for pions (π ± ), kaons (K ± ), plus (anti)protons; but e.g. in dss there is still a small "residual" contribution on top of these. There are other sets of FFs with restricted particle species [60,57,61], but we do not consider them here as we focus only on the inclusive sum which should be better constrained than the FFs for individual hadron species. For reviews, see e.g. [4,5]. Table 1: Characteristics of the existing sets of parton-to-charged-hadron FFs. The hadron species included, use of different collision systems, attempts to estimate the FF errors, minimum value of z considered, and the available Q 2 -range are indicated.

FF set
Species Fitted data Error estimates zmin Q 2 (GeV 2 ) Kretzer (kre) [52] π The main constraints in global fits of FFs come from the inclusive hadron production in e + e − collider experiments, from which the data are abundant [62]. These experiments are, however, mainly sensitive to the quark FFs leaving the gluon-to-hadron FFs largely unconstrained. Similarly, the data for semi-inclusive hadron production in deeply inelastic scattering [63], used in the dss fit, predominantly originates from the quark fragmentation. This leaves the gluon FFs prone to parametrization-bias and other theoretical assumptions. To improve on this, dss and akk08 include various datasets (in different combinations) from hadronic collisions at RHIC [10,11,64], SppS [12,13,14] and Tevatron [15,65]. The bulk of these measurements is, however, concentrated at rather small values of transverse momentum p T 5 GeV/c (dictating the hard scales of the process) where the sensitivity to the gluon FFs is certainly present, but where the NLO pQCD calculations are not well under control due to the large scale uncertainty (see later). Therefore, the extraction of the gluons FFs based on these data cannot be considered completely safe, and, therefore, we believe that the "older" sets of FFs with only e + e − data should not be blindly discarded. Although the uncertainties in FFs due to the experimental errors (or lack of data) have been addressed in some FF extractions (see also Ref. [66]), only hkns provides the necessary information to fully propagate these uncertainties to further observables. The bfgw package provides three alternative sets with somewhat different gluon FFs, but the mutual variation between these sets translates only to a few-percent difference in the observables we discuss here. The error analysis of dss was performed via the method of Lagrange multipliers [67] which does not allow the end user to estimate the propagation of their FF errors. For these reasons, in what follows, we will present the central results from all the parametrizations including also the hkns error bands. In Fig. 1, we present a comparison of the up-quark and gluon FFs into charged hadrons for all the available FFs at a common scale Q = 20 GeV. The spread among the gluon FFs is significantly larger than in the case of quarks -a clear indication of lack of definitive constraints -in particular for moderate and large z > 0.3 values. Indeed, even the hkns error band is not broad enough to cover all different sets above z ∼ 0.5. We believe this is mainly due to the small amount of fit parameters that could be left free in the absence of strict gluon constraints from the e + e − data alone, as discussed above. To understand at which z values the hadron production at LHC energies probes the FFs, and to what extent this depends on the hardness of the gluon FFs, we examine the shape of the differential z distributions. A couple of typical examples corresponding to the LHC kinematics are shown in Fig. 2. The old Kretzer and the modern dss FFs are considered here as the hardness of their gluon FFs is quite different. The z distributions appear to be rather broad and to depend significantly on the specific set of FFs used. The spectra in the range p T = 5-20 GeV/c probe average hadron fractional momenta z ≈ 0.4-0.6 at √ s = 900 GeV, decreasing to z ≈ 0.3-0.5 at √ s = 7 TeV. The behaviour of these z distributions can be understood approximating thê p 3T dependence of the convolution between the PDFs and the partonic cross sections dσ by a power law, as follows [68]: where C is a p T -independent constant, so that one can write Due to the factor z n−1 , the contributions from small values of z are efficiently suppressed, and the average value of z becomes much larger than the kinematic lower limit z min ≈ 2p T / √ s at midrapidity. However, towards larger √ s, the p T distributions are flatter (the power n is smaller), and the z distributions become on average shifted towards smaller values of z. In any case, the contributions from the small-z region, z 0.05 − 0.1, which is more difficult to treat in DGLAPbased approaches [28,6], remain very small and the use of the standard FF framework is well justified. The relative contributions from the quark and gluon fragmentations are plotted in Fig. 3 for √ s = 900 GeV (left) and √ s = 7000 GeV (right). At small p T , the gluon fragmentation clearly dominates but towards large values of p T the quark fragmentation becomes eventually predominant. In any case, the gluon contribution is always significant and therefore the LHC promises to be a good "laboratory" to determine the gluon FFs in the region of p T where the NLO pQCD calculations can be considered to be well under control.

Comparison of NLO pQCD to high-p T charged-hadron collider data
In this section, we compare the data from various experiments with the NLO calculations using the seven FF sets listed in Table 1. Our main attention will be on the latest data from CMS [20,21] and ALICE [22] for p-p collisions at the LHC, as well as the CDF measurements [16] in p-p collisions at Tevatron. We do not include the similar ATLAS measurements [69] here as their results have not been given as invariant cross sections, but only in terms of the absolute yields. However, we have checked that the shapes of the ATLAS p T -spectra are in agreement with those measured by the other LHC experiments. In order to compare with some earlier hadron-collider data, included in the akk08 and dss global fits, we consider also the p-p results from UA1 [12,13,14] as well as the p-p spectra measured by STAR [10]. As an example of the p T -differential cross sections, Fig. 4 presents the CMS measurements for inclusive charged hadrons at √ s = 0.9, 2.76, 7 TeV at midrapidity. The data, spanning five orders of magnitude in perturbatively-accessible values of p T , are compared to the NLO calculations using two sets of FFs, dss and Kretzer. While dss clearly overshoots the data (by up to a factor of 2), the Kretzer FFs do a much better job in describing the spectra both in shape and absolute normalization. The apparent difference between the two parametrizations derives from the fact that the large-z Kretzer gluon FFs are much softer than those of dss (Fig. 1).
The central result of our paper is shown in Fig. 5 which presents a comprehensive comparison of the charged-hadron world-data in the TeV-range to the NLO predictions using the seven most recent sets of FFs. The panels show the ratio of the data to the predictions obtained with the Kretzer FF (data points), as well as the ratios of cross sections obtained with various FFs over those from Kretzer (curves). The light-blue band denotes the scale uncertainty envelope (cf. Eq. (2)), and the dark-blue band the uncertainty derived from the ct10nlo PDF (90% confidence level) error sets. The hkns error band is shown in light-brown color. The scale uncertainty below p T ≈ 10 GeV/c is prohibitively large making it difficult to draw any strict conclusion regarding the level of agreement between the data and the calculations there (although the central NLO prediction obtained with the Kretzer FFs agrees very well with the data). The scale dependence, however, stabilizes below ±20% beyond p T ≈ 10 GeV/c and it is rather this region where the NLO calculations are to be fully trusted. In all cases, the PDF errors are almost negligibly small in comparison to the scale uncertainty.   [20,21], ALICE [22] (diamonds), CDF (squares) [16], and UA1 (triangles) [13] at √ s = 900-7000 GeV, over the corresponding NLO calculations using the Kretzer FFs. The curves show the NLO predictions obtained with other FF sets: kkp (pink scarcely dashed), dss (green dashed), bfgw (brown long-dashed), hkns (purple dashed-dotted), akk08 (yellow dotted-dashed), and akk05 (red long-dashed short-dashed) relative to Kretzer FFs. The point-to-point systematic and statistical errors are indicated by colored rectangles (gray for CMS and CDF, brown for ALICE, green for UA1) and error bars. The boxes at the beginning of the p T -axis mark the overall normalization uncertainty. The light-blue bands correspond to the scale uncertainty envelopes while the dark-blue ones indicate the variations derived from the ct10nlo PDF error sets. The hkns uncertainties are shown by the light-brown bands.
All the LHC data are in mutual agreement within their systematic and statistical uncertainties, although especially the CMS data at 2.76 TeV seem to show larger fluctuations than those inferred from the typical size of the quoted point-to-point experimental uncertainty. The UA1 spectrum at √ s = 0.9 TeV is mostly above the CMS and ALICE results and is only barely compatible with them. However, as the shape of the p T distribution is not incompatible with the rest, this disagreement could well be an issue of the experimental determination of the overall normalization in the oldest measurement.

8
The results of Fig. 5 exhibit clear systematic trends as a function of √ s and p T : As √ s increases from 0.9 TeV to 7 TeV, the experimental spectra gradually sink more and more below the theoretical predictions. However, the shape of the data-to-theory ratio remains qualitatively similar regardless of the collision energy. Especially, relative to the calculation with the Kretzer FFs, the data first show a "bump" at p T ≈ 4 GeV/c but then straighten up beyond p T ≈ 10 GeV/c. Indeed, the flatness of the data/theory ratio is worth noticing, although the absolute spectra span many orders of magnitude. This suggests that the underlying pQCD dynamics of the hadron production is indeed correctly understood and that the data-theory disagreement lies rather in the current sets of FFs. On average, the Kretzer FFs used as reference for the data/NLO ratios shown in Fig. 5, seem to do the best job in describing the data, the results from all other FFs being practically enclosed by the hkns error bands. We quantify the data-to-theory correspondence by computing the χ 2 values for each FF set (in the case of hkns only the central predictions are considered) defined by where D i corresponds to the data point with total error δ tot i (correlated and uncorrelated point-topoint uncertainties added in quadrature), and the theory values T i are specific for each set of FFs. The sum runs over all the data shown in Fig. 5 using three different cuts for the hadron transverse momenta: p min T = 1.3, 5, 10 GeV/c. Such thresholds are chosen so as to reduce the weight of the lower-p T data which would otherwise dominate the χ 2 due to their larger cross sections and associated smaller statistical uncertainties. The calculations are run for three scale choices, µ ≡ (µ ren /p T , µ fact /p T , µ frag /p T ) = 1 2 , 1 2 , 1 2 , (1, 1, 1), (2, 2, 2) , which, above p T ≈ 5 GeV/c, practically cover the larger selection of scale variations used earlier, Eq. (2). The results from this exercise are shown in Table 2, and numerically confirm what is seen in Fig. 5: The lowest χ 2 value is almost always obtained with the Kretzer FFs and the highest one with akk05. The preferred choice of scale is specific for each set of FFs and depends on the p min T cut imposed. However, on average, the choice µ = (2, 2, 2) is preferred as this set of values tends to reduce the computed cross sections and thereby moderate the data overshooting. The same conclusion is reached if the χ 2 is computed accounting separately for the correlated systematic errors. In this particular case, the only known systematic parameter is the overall normalization (given by all but UA1) and the χ 2 can be expressed [67] as where δ uncor i is the uncorrelated error, and δ norm i the quoted normalization error. The parameter f k is found by minimizing the χ 2 . The corresponding results are listed in Table 3. In comparison to the uncorrelated-uncertainties case, the χ 2 values become generally somewhat lower, since, the calculated values are usually quite above the data, and the improvement attained in the first term in Eq. (7) exceeds the growth of the latter. This procedure, however, often also leads to unnaturally large values of |1 − f k | as the disagreement between the calculations and the data tends to be much beyond the normalization uncertainty quoted by the experiments. In any case, the values in Table 3 lead to the same conclusions as those extracted from Table 2 and the final outcome is the same regardless of the way the χ 2 -function is defined, namely that NLO predictions generally overpredict the experimental charged-hadron spectra by a factor of two. Similar discrepancies have been found for high-p T neutral pion and η meson production at LHC energies [70] implying that the problem is not limited to the total (g → h + + h − ) fragmentation function but affects the identified (π ± , K ± , p, p) gluon FFs individually. Finally, we take a look at collider data at lower c.m. energies, √ s = 200-630 GeV, where the situation is different than that found at the LHC energies. Fig. 6 presents a comparison of NLO calculations to the UA1 [12,13,14] and STAR [10] charged-hadron spectra. For these datasets, the calculation obtained with the Kretzer FFs -the preferred set at the Tevatron and LHC energies -gives a bad description of the spectra. On the other hand, the dss set describes now these measurements reasonably well. This is, however, not too surprising as these datasets were included in their actual fit. In this case, also the hkns error band is wide enough to enclose these data. Looking back to Fig. 1, one can see that the lower-energy collisions prefer much harder gluons at large z than the LHC data. That is, any set of FFs that can, more or less, reproduce the lower-√ s data (preferring hard gluon FFs), will disastrously overshoot the LHC measurements (preferring softer gluon FFs). As the variation in the probed range in z is only mild as a function of √ s and p T (see Fig. 2) such a result hints that it may be difficult to tensionlessly include all these existing data into a global FF fit with a p T cut as low as e.g. p T ≥ 2 GeV/c. Indeed, as the very large scale-uncertainty indicates, the fixed-order NLO calculations are not stable below p T ≈ 10 GeV/c, and it is questionable whether these lower-p T data points should be even considered in such a fit in the absence of full NNLO corrections. It should be also recalled that within the RHIC kinematics reach, the threshold logarithms can still play a role by increasing the NLO cross sections for increasingly-high p T values, although such effects should die out towards larger √ s (at fixed p T ).
In any case, the fact that such effects cannot be seen in Fig. 6, signals that threshold resummations are likely not the main cause for the different √ s-dependence of the NLO calculations and the data.
In addition, the good agreement between fixed-order NLO calculations [33,34] and the jet data at LHC [29,30], Tevatron [71,72] and RHIC [73], and the fact that the single high-p T charged particle spectrum is dominated by leading hadrons carrying out a large fraction, z ≈ 0.5, of the parent parton (jet) energy (Fig. 2), strongly reinforces the idea that the origin of the data-theory disagreement lies in an imperfect knowledge of the final gluon-to-hadron fragmentation functions.
As a matter of fact, at low values of p T the whole picture of independent parton-to-hadron fragmentation may not be adequate, especially in the case of production of heavier baryons. Highertwist effects, where the hadron is produced directly (i.e. more exclusively) in the hard subprocess rather than by gluon or quark jet fragmentation, may contribute to the cross sections at RHIC energies in the range of transverse momenta experimentally studied [74] and hence "contaminate" the extraction of FFs in global fits that use such data. Even at LHC energies, the proton-to-pion ratio below p T ≈ 6 GeV/c (see e.g. [75]) appears to behave qualitatively very differently than the kaon-to-pion ratio or the pQCD expectations. While the kaon-to-pion ratio increases smoothly towards larger p T , the proton-to-pion ratio contains a clear "bump" around p T ≈ 3 GeV/c. To reproduce such a behaviour, additional effects outside the pQCD toolbox are called for. Assuming that the behaviour of the perturbative proton-to-pion ratio is qualitatively similar to the the kaon-to-pion ratio, one could crudely estimate that there is a roughly 5% "non-fragmentation" enhancement around p T ≈ 3 GeV/c which then diminishes towards higher p T . Subtracting such a non-perturbative contribution from the LHC data would make the disagreement between the data and e.g. dss even worse. On the contrary, in comparison to the description with the Kretzer FFs, the data-to-theory ratio would be flatter and thereby improve the compatibility. For lower √ s the surplus of (anti)protons could be even more pronounced and remain important up to higher values of p T . This could partly explain the strong √ s-dependence of the data-to-theory ratio.
In Ref. [22], it was observed that the ratios of the ALICE cross sections between different but nearby √ s become rather well reproduced -also at low p T -by the dss FFs. However, it is important to note that as the √ s dependence of the z distributions (see Fig. 2) is only mild, part of the FF dependence is bound to cancel in such ratios. Despite such cancellations, it appears that some √ s dependence still remains at low p T , supporting our conjecture regarding the presence of a √ s-dependent non-perturbative component.

Summary and outlook
We have examined the LHC, Tevatron, SppS and RHIC data for inclusive unidentified chargedhadron production at √ s = 0.2 − 7 TeV against NLO pQCD calculations with seven different sets of FFs, quantifying also the systematics associated with the scale and PDF uncertainties. The spread among the predictions with different FFs is large and can be traced back to sizable mutual differences in the gluon-to-hadron FFs. None of the existing FF sets can reproduce the experimental results optimally, on average, clearly overshooting the data by up to a factor of two in the large-p T region, p T ≥ 10 GeV/c, where the fixed-order NLO pQCD calculations should be trustworthy as their scale dependence is modest. The best overall agreement with the data is obtained using the relatively old Kretzer FFs in which the gluon-to-hadron FFs are the softest of all. Below p T ≈ 10 GeV/c, the √ s-dependence of the data appears too strong to be reproduced by calculations based solely on collinear factorization (especially if also lower values of √ s are included in the comparison), although the NLO scale uncertainties become there very large preventing one from drawing definitive conclusions. However, this may not be a problem of the pQCD calculation itself as below p T ≈ 10 GeV/c there are increasing indications of additional non-perturbative effects in the case of (anti)proton production even at the LHC energies, which may have a nonnegligible impact on the total hadron yield. These observations indicate that only the region above p T ≈ 10 GeV/c of these charged-hadron data, with theoretical scale uncertainties below ±20%, should be included in forthcoming global fits of parton-to-hadron fragmentation functions.