{'fraction_non_alphanumeric': 0.03595618063870472, 'fraction_numerical': 0.013546115767464303, 'mean_word_length': 4.341032738379307, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Parton-to-Pion Fragmentation Reloaded

We present a new, comprehensive global analysis of parton-to-pion fragmentation functions at next-to-leading order accuracy in QCD. The obtained results are based on the latest experimental information on single-inclusive pion production in electron-positron annihilation, lepton-nucleon deep-inelastic scattering, and proton-proton collisions. An excellent description of all data sets is achieved, and the remaining uncertainties in parton-to-pion fragmentation functions are estimated based on the Hessian method. Extensive comparisons to the results from our previous global analysis are performed.


I. INTRODUCTION AND MOTIVATION
The quantitative description of hard scattering processes involving identified light hadrons in the final-state requires a precise knowledge of how quarks and gluons hadronize. In the framework of perturbative QCD (pQCD), which we pursue in the following, this vital information is encoded in parton-to-hadron fragmentation functions (FFs) [1]. To match the increasing amount and precision of experimental results, the availability of reliable sets of FFs for a large variety of hadrons, in particular, for neutral and charged pions and kaons, as well as accurate estimates of their uncertainties is of the utmost relevance and the subject of this study.
Like parton distribution functions (PDFs), FFs are required in a pQCD calculation to consistently absorb certain classes of collinear parton-parton configurations related to long-distance physics, i.e., interactions happening a long time after the actual hard scattering process. As such, FFs are non-perturbative quantities, and any information about them needs to be gathered from data, preferably in a global QCD analysis combining results obtained in a large variety of processes. These fits of FFs are facilitated by assuming factorization [2], which allows one to compute the relevant short-distance hard scattering matrix elements perturbatively, and the fact that pQCD predicts the scale evolution of FFs very much in the same way as for PDFs. Fragmentation functions depend on the parton of flavor i which hadronizes, the fraction z of its four-momentum taken by the observed hadron H, and the scale Q at which they are probed in a hardscattering process. In what follows, they will be denoted as D H i (z, Q 2 ). All relevant ingredients for a global QCD analysis of FFs are fully known up to next-to-leading order (NLO) accuracy in the strong coupling α s comprising the kernels governing the time-like scale evolution [3,4], single-inclusive hadron production in electronpositron annihilation (SIA) [5][6][7] and proton-proton (pp) collisions [8], and hadron multiplicities in semi-inclusive deep-inelastic lepton-nucleon scattering (SIDIS) [7,9].
A first global QCD analysis of data collected for all these hadron production processes in terms of pion (and kaon) FFs has been performed quite some time ago in Ref. [10], commonly known as "DSS analysis", followed by similar fits for eta mesons, protons, and unidentified charged hadrons [11,12]. Compared to earlier studies based on SIA data [13,14], the DSS fit [10] fully exploited the synergy of the complementary probes of FFs (SIA, SIDIS, and pp) to arrive at a more detailed and data-driven separation of the individual FFs for different quark flavors than before. The use of charge separated, i.e., H = π ± , K ± , SIDIS multiplicities [15] was instrumental in this respect as they weight quite differently contributions of the various quark flavors i in the hadronization process. While SIA data are more straightforward to analyze and free of PDF uncertainties, the LEP [16][17][18] and SLAC [19] data used in the DSS fit only provided a handle on flavor-separated FFs when supplemented by corresponding flavor-tagged results [17,19,20], which have no unambiguous theoretical description [10,13] and heavily rely on Monte Carlo simulations to extract them experimentally [20]. One peculiar finding of the DSS global analysis [10] was an unexpectedly large charge symmetry violation between the total u-and d-quark FFs for pions, within sizable uncertainties though. In addition, single-inclusive pion data from pp collisions at BNL-RHIC [21] provided a first constraint on the gluon-to-pion FFs, which, at that time, was impossible to determine otherwise as precise enough SIA data were only available from the LEP and SLAC experiments, i.e., at a fixed scale M Z , the mass of the Z boson.
To implement the lengthy, exact NLO expressions for hadron production processes in SIDIS and pp collisions without any approximations into the theoretical framework of a global fit, DSS adopted the Mellin technique [22,23]. The gist of this method is to pre-calculate all time-consuming NLO expressions for SIDIS and pp processes once, before the actual fit is performed, and to store the required information on look-up tables; for details, see [10,23,24]. The use of Mellin moments is also most appropriate in solving the QCD evolution equations.
Uncertainties of the extracted FFs were estimated based on the robust Lagrange multiplier (LM) technique [25], but only for a specific moment of the D H i (z, Q 2 ) contributing to the momentum sum rule. More detailed studies were performed in Ref. [26], where, in addition, the applicability of the standard iterative Hessian (IH) approach [27] was explored. Comparisons with the results obtained with the LM technique revealed, however, some limitations of the IH method, mainly due to the lack of sufficiently precise experimental information at that time to warrant the assumption that any deviations from the optimum fit are quadratic in all the parameters specifying the FFs.
In the present paper we build upon the theoretical and conceptual framework developed for the DSS analysis [10] but make use of a wealth of newly available data sets, which will enable us to relax and scrutinize some of the constraints imposed on the parameter space in the DSS fit. The key assets of the new analysis are the recently published precise SIA data from BaBar [28] and Belle [29], which, in principle, should provide a novel handle on the gluon FF through QCD scaling violations of the SIA structure functions between the scale Q = M Z , rel-evant for the LEP and SLAC experiments, and the scale corresponding to the center-of-mass system (c.m.s.) energy of BaBar and Belle, Q = √ S ≃ 10.5 GeV. In addition, since the electroweak couplings of up-type and down-type quarks to the Z boson become almost equal at Q ≈ M Z , LEP and SLAC data are mainly sensitive to the total quark singlet FF for any observed hadron H. At the lower √ S of BaBar and Belle, the quarkantiquark pairs in SIA are weighted according to their electrical charge, which in our global fit should allow for some partial flavor separation of FFs.
Another important and new ingredient to the current analysis are the final SIDIS data released by the Hermes Collaboration [30] which supersede the preliminary and much less precise data utilized in the DSS fit [15]. New, still preliminary data for pion multiplicities in SIDIS are also available, for the first time, from the Compass experiment at CERN [31], which are very precise despite exhibiting a fine binning in the relevant kinematic variables. Finally, first results on single-inclusive pion spectra at high transverse momenta p T have become available from the LHC at c.m.s. energies of up to 7 TeV [32], which supplement the data from BNL-RHIC taken at √ S = 200 GeV that have been already used in the original DSS analysis [21]. We also include several recent results from the Star Collaboration for both neutral and charged pion production at √ S = 200 GeV [33][34][35][36]. We note that at variance with the original DSS analysis, we now determine the optimum normalization shifts for each data set in the fit analytically (see, e.g. Ref. [25] for a discussion of normalization shifts in PDF fits), which greatly facilitates the global fitting procedure.
The main goal of our new analysis is to extract an updated set of parton-to-pion FFs and to determine their uncertainties reliably based on the IH method [27] in light of all the newly available, precise experimental results in SIA, SIDIS, and pp collisions. This will allow us to scrutinize the consistency of the information on FFs extracted across the different hard scattering processes, i.e., to validate the fundamental notion of universality, which is at the heart of any pQCD calculation based on the factorization of short-and long-distance physics [2] sketched above.
For the time being, we have to limit ourselves to pion FFs as a similar level of improvements on the available experimental information is still lacking for kaons, most noticeable for the SIDIS process, which is crucial in determining flavor-separated FFs. Nevertheless, we strongly believe that our updated global analysis of parton-to-pion FFs is very timely for the reasons mentioned above and the fact that precise FFs are in high demand as input for global analyses of helicity PDFs [24,37] and transverse momentum dependent PDFs [38], both of which heavily draw on data with identified pions in the final-state. Other applications involve to quantify and understand possible modifications of hadron production yields in the presence of a nuclear medium, as studied in heavy ion collisions both at RHIC and the LHC [39].
Since extractions of leading order (LO) FFs have yielded a much less satisfactory description of the available pion production data in the DSS analysis [10], we only perform our global QCD fit at NLO accuracy. In any case, the need for LO FFs (and PDFs) has greatly diminished in recent years with the advent of novel theoretical tools that allow one to compute NLO cross sections largely automatically. The obtained optimum NLO parton-to-pion FFs, including the Hessian eigenvector sets, are available upon request and enable one to straightforwardly propagate our obtained uncertainties to any observable of interest.
It should be noted that the necessary time-like evolution kernels for FFs are available even at next-to-NLO (NNLO) accuracy now, with the exception of one, presumably minor, detail for phenomenological applications [40]. However, the corresponding partonic hard scattering processes have been only computed for SIA so far [41]. Nevertheless, it might be an interesting future endeavor to perform a NNLO analysis of SIA data alone and, perhaps, to investigate the impact of also available all-order resummations of potentially large logarithmic corrections near the partonic threshold [42]. This is, however, well beyond the scope of the current analysis.
The remainder of the paper is organized as follows: in the next Section we briefly summarize the main aspects of our updated global analysis, including the choice for the functional form used to parametrize the FFs at the initial scale for the QCD evolution, the selection of data sets and cuts imposed on them, and the treatment of experimental normalization uncertainties. The outcome of the new fit is discussed in depth in Sec. III. The obtained partonto-pion fragmentation functions and their uncertainties are shown and compared to the results of our previous global analysis. Detailed comparisons to the individual data sets are given to demonstrate the quality of the fit. Potential open issues and tensions among the different data sets will be discussed. We briefly summarize the main results in Sec. IV.

II. FRAMEWORK
In this Section we lay out the framework and key ingredients for our global QCD analysis of parton-to-pion FFs. We mainly focus on those aspects that differ from the original DSS analysis [10].

A. Functional Form and Fit Parameters
The functional form adopted in the DSS global analysis [10] is flexible enough to accommodate also the wealth of new experimental information included in the present fit. Therefore, we continue to parametrize the hadronization of a parton of flavor i into a positively charged pion at an initial scale of Q 0 = 1 GeV as .
(1) Here, B[a, b] denotes the Euler Beta-function, and the N i in (1) are chosen in such a way that they represent the contribution of zD π + i to the momentum sum rule. Compared to our previous analysis, the improved experimental information now allows us to impose less constraints on the parameter space spanned by the input function in Eq. (1). More specifically, as before we still have to assume isospin symmetry for the unfavored FFs of light sea quarks, i.e., and we need to relate the total u-quark and d-quark FFs by a global, z-independent factor N d+d , which quantifies any charge symmetry violation found in the fit. The fragmentation of a strange quark into a pion is now related to the unfavored FFs in Eq. (2) by rather than just using a constant as in the DSS analysis.
The charm-and bottom-to-pion FFs no longer assume γ c = γ b = 0 in Eq. (1) but can now exploit the full flexibility of the ansatz. This is not due to new flavor-tagged data but helps the global fit to accommodate the recent, very precise results from BaBar [28] and Belle [29] in SIA and from Compass [31] and Hermes [30] in SIDIS, which now constrain both the total quark fragmentation, i.e., summed over all flavors, and the individual flavorseparated, light quark FFs much better than before. As in the DSS and all other analyses [10][11][12][13][14], we include heavy flavor FFs discontinuously as massless partons in the QCD scale evolution above their MS "thresholds", Q = m c,b , with m c and m b denoting the mass of the charm and bottom quark, respectively. Conceptually, due to confinement, there has to be a heavy quark FFs present as soon as the heavy quark can be produced in the final-state of a hard scattering process. We leave it to dedicated future studies to explore and incorporate an improved theoretical framework for heavy quark-tolight hadron fragmentation functions into the global fitting procedure, following the rather elaborate schemes that have been developed for heavy quark parton densities [43] to properly include mass effects near threshold and to resum potentially large logarithms ∼ ln m 2 c,b /Q 2 for Q 2 ≫ m 2 c,b . We note that a dynamical, parameterfree generation of the heavy flavor component to light meson fragmentation functions has been developed, for instance, in Ref. [44].
In total we now have 28 free fit parameters describing our updated FFs for quarks, antiquarks, and gluons into positively charged pions, which are determined from data by a standard χ 2 minimization to be described below. The corresponding FFs for negatively charged pions are obtained by charge conjugation and those for neutral pions by assuming D π 0 i = [D π + i + D π − i ]/2. We note that none of the constraints imposed on the fit through Eqs. (2)-(4) has any impact on its overall quality.

B. Data Selection
We make use of all the currently available experimental information on single-inclusive charged and neutral pion production in SIA, SIDIS, and hadron-hadron collisions to determine the free fit parameters defined in Sec. II A.
Compared to the data sets already used in the DSS global analysis [10], we include the new results from BaBar [28] and Belle [29] in SIA at a c.m.s. energy of √ S ≃ 10.5 GeV. Both sets are very precise, with relative uncertainties of about 2 − 3%, and reach all the way up to pion momentum fractions z close to one, well beyond of what has been measured so far. We analyze both sets with n f = 4 active, massless flavors using the standard expression for the NLO SIA cross section [5]. As customary, we limit ourselves to data with z ≥ 0.1 to avoid any potential impact from kinematical regions where finite, but neglected, hadron mass corrections, proportional to M π /(Sz 2 ), might become of any importance [10,13,14]. For SIA data taken at higher √ S we use n f = 5 and z > 0.05, following the original DSS analysis. Any incompatibility of the two new precise sets of data at √ S ≃ 10.5 GeV with each other or with the old LEP and SLAC data at √ S ≃ 91.2 GeV [16][17][18][19] has the potential to seriously spoil the quality of the global fit.
In case of SIDIS, we replace the preliminary multiplicity data from Hermes [15] by their recently released final results [30]. More specifically, we use the data for charged pion multiplicities as a function of momentum transfer Q 2 in four bins of z taken on both a proton and a deuteron target. The range of average values of Q 2 covered by the data is from about 1.1 GeV 2 to 7.4 GeV 2 and 0.2 ≤ z ≤ 0.8. In addition, we include the still preliminary multiplicity data for π ± from the Compass Collaboration [31], which are given as a function of z in bins of Q 2 and the initial-state momentum fraction x. The coverage in z is the same as for the Hermes data, but due to the higher √ S of the Compass experiment the reach in x and Q 2 is wider. Experimental information is available for 0.004 ≤ x ≤ 0.7 and 1.2 ≤ Q 2 ≤ 22.4 GeV 2 . We do not have to impose any cuts on both data sets to accommodate them in the global analysis. As for SIA, having now available two precise sets of multiplicity data in SIDIS, covering slightly different but partially overlapping kinematics, makes it very important to validate their consistency in a global fit.
Finally, we add a couple of new sets of data for inclusive high-p T pion production in pp collisions to the results from the Phenix experiment [21] already included in the DSS analysis. Most noteworthy are the first results for neutral pions from the Alice Collaboration at CERN-LHC [32], covering unprecedented c.m.s. energies of up to 7 TeV. In addition, we add Star data taken at √ S = 200 GeV in various rapidity intervals for both neutral and charged pion production and for the π − /π + ratio [33][34][35][36]. As we will demonstrate and discuss in more detail in Sec. III D below, it turns out that a good global fit of RHIC and LHC pp data, along with all the other world data, can only be achieved if one imposes a cut on the minimum p T of the produced pion of about 5 GeV. Such a cut eliminates some of the pp data points included in the previous DSS analysis from the fit, in particular, all the Brahms [45] and Star [46] data at forward pseudorapidities and, hence, too small values of p T .

C. Fit Procedure and Uncertainty Estimates
The 28 free parameters describing the updated partonto-pion FFs in Eq. (1) at the chosen input scale of 1 GeV are again determined from a standard χ 2 minimization where for i = 1, . . . , N data sets, each contributing with N i data points. E j is the measured value of a given observable, δE j the error associated with this measurement, and T j is the corresponding theoretical estimate for a given set of parameters in Eq. (1). Since the full error correlation matrices are not available for some of the data sets used in the fit, statistical and systematical errors are simply added in quadrature in δE j as in all previous fits [10,11,13,14]. At variance with the original DSS fit [10], where we have introduced several extra fit parameters to account for experimental normalization uncertainties N i in (5), we now derive the optimum normalization shifts for each data set analytically from the condition ∂χ 2 /∂N i = 0, which yields Here, δN i denotes the quoted experimental normalization uncertainty for data set i. In Sec. III A we will list the so obtained normalizations N i along with the individual χ 2 values for each data set included in the fit.
In the DSS analysis [10] we assessed uncertainties in the extraction of fragmentation functions with the help of the LM technique [25] by mapping out the maximum allowed range of variation in the fit of the truncated second moments of the fragmentation functions for x min = 0.2 and Q = 5 GeV. While this method is very robust, even when some of the fit parameters are only loosely constrained by data, it has the disadvantage that uncertainties cannot be easily propagated to other observables of interest. In Ref. [26] we have therefore explored the applicability of the iterative Hessian approach [27] based on the original DSS choice of data sets, cuts, and parameters by comparing its outcome to uncertainty estimates obtained with the LM method. The main idea of the IH method is to assume a quadratic behavior of the χ 2 hyper-surface of parameter displacements and to express the χ 2 increment from its minimum value in terms of combinations of fit parameters that maximize the variation. Such an eigenvector representation of the Hessian matrix proves to be extremely suitable to compute the propagation of uncertainties to arbitrary observables in terms of a limited number of pre-calculated sets of FF functions (in fact, twice the amount of fit parameters). These sets correspond to fixed displacements along the eigenvector directions of the Hessian matrix.
With the much increased availability of precise data for the current analysis, the sole use of the computationally less demanding IH method to quantify uncertainties of FFs becomes viable and will be pursued in the following. The obtained eigenvector sets of FFs will be made available upon request from the authors along with a parametrization of the optimum fit. To define the eigenvector sets one has to choose a tolerance parameter ∆χ 2 for the increment in χ 2 which is still acceptable in the global fit. Here we proceed as follows: the tolerances for the eigenvector sets corresponding to 68% and 90% confidence level (C.L.) intervals are determined from the Gaussian probability density function for a χ 2 distribution with k degrees of freedom (d.o.f.): The ∆χ 2 related to the 68 th and 90 th percentiles are then obtained by solving χ 2 +∆χ 2 0 dχ 2 P k (χ 2 ) = 0.68 and 0.90, respectively.
Finally, we choose the NLO set of PDFs from the MSTW group [47] and the corresponding uncertainty estimates in computations of the SIDIS and pp cross sections. For consistency, we also fix the strong coupling α s to the values obtained in the MSTW fit. We note that in the x and Q 2 region relevant for our global analysis of FFs, the needed combinations of PDFs are relatively well constrained. A choice of PDFs other than the MSTW set would not alter the outcome of our fit in any significant way. We will illustrate in Sec. III D below that theoretical scale ambiguities are considerably larger than PDF uncertainties.

III. RESULTS
In this Section we present and discuss in depth the results of our global analysis of parton-to-pion FFs. First, we present the obtained fit parameters, normalization shifts, and individual χ 2 values. Next, the obtained D π + i (z, Q 2 ) and their uncertainties are shown and compared to the results of the DSS fit. The quality of the fit to SIA, SIDIS, and pp data and potential open issues and tensions among the different sets of data are illustrated and discussed in Sec. III B, III C, and III D, respectively. In Table I we list the obtained set of parameters specifying our updated, optimum parton-to-pion fragmentation functions at NLO accuracy at the input scale Q 0 = 1 GeV for the light quark flavors and the gluon, and at their respective thresholds Q 0 = m c,b for the charm and bottom quarks. Table I reveals already a notable difference to one of the findings of the DSS analysis which preferred an unexpectedly sizable breaking of the charge symmetry between u +ū and d +d FFs of about 10% [10], within large uncertainties though. This was mainly driven by the preliminary π ± multiplicities from Hermes [15] used in the fit at that time. Now, with much improved experimental information on charged pion multiplicities both from Hermes [30] and Compass [31] and new data on the ratio π − /π + in pp collisions from Star [35], the parameter N d+d in Eq. (3) prefers to stay very close to unity, i.e., very little or no breaking.

A. Parton-To-Pion Fragmentation Functions
As has been mentioned above, in case of the unfavored FFs, data now allow us to introduce some non-trivial z dependence, see Eq. (4), to parametrize a potential SU(3) breaking between theū = d and s =s FFs. This little extra freedom not only helps to accommodate all the different data sets used in the global analysis in a much better way, but it also leads to more realistic uncertainty estimates for both D π + s and D π + u . In particular, it now turns out that the uncertainties for D π + s are much bigger than for D π + u as can be inferred from Fig. 1, where we present the individual parton-to-pion FFs D π + i (z, Q 2 ) at Q 2 = 10 GeV 2 . The four leftmost panels show the opti- mum zD π + i at NLO accuracy for i = u +ū, d +d,ū = d, s =s, c =c, and the gluon g (solid lines) along with our uncertainty estimates at 68% C.L. (inner bands) and 90% C.L. (outer bands), obtained as described in Sec. II C. For better visibility, the rightmost panels give the relative uncertainties for the same set of zD π + i . The results of the previous NLO DSS fit are shown as dashed lines.
As can be inferred from Fig. 1, for the light quark flavors the old DSS results are either close to the updated fit or within its 90% C.L. uncertainty band. The best determined pion FFs is D π + u+ū , where the relative uncertainties are below 10% at 90% C.L. throughout most of the relevant z range. Only for z 0.8 the errors rapidly increase because of the lack of experimental constraints in this region. The corresponding uncertainties for D π + d+d turn out to be slightly larger as they also include possible violations of SU(2) charge symmetry through Eq. (3). We stress again, that at variance with the DSS analysis [10], the new fit does not favor any SU(2) breaking. For the unfavored FFs, D π + u = D π + d are determined well in a much more limited range of z, and uncertainties start to increase already for z 0.5. The corresponding ambiguities on D π + s = D π + s are about a factor of two larger and amount to at least 25% at 90% C.L. for z ≃ 0.3.
Bigger deviations from the DSS analysis are found for both the gluon and the charm FFs. In the latter case, this is driven by the greater flexibility of the functional form, five fit parameters rather than three, which helps with the overall quality of the global fit and cannot be pin-pointed to a particular data set. In fact, there had been no new charm (or bottom) tagged data since the LEP and SLAC era. The significantly reduced D π + g as compared to the DSS fit is a result of the new Alice pp data [32], which have a strong preference for less pions from gluon fragmentation for basically all values of z. We will discuss this finding, and possible tensions arising with the pp data from RHIC, in more detail in Sec. III D. The relative uncertainties on D π + g at Q 2 = 10 GeV 2 are about 20% at 90% C.L. up to z ≃ 0.5 and quickly increase towards larger z values.
We refrain from performing a detailed comparison to the uncertainty estimates based on the data sets available for the original DSS analysis [10,26] as they can be viewed at best as a rough approximation. Only with the quality and variety of data sets available for the current global analysis one can arrive at a first meaningful determination of uncertainties for parton-to-pion FFs, which therefore constitutes as one of the main results of this study.
We note that the new very precise SIA data from BaBar [28] and Belle [29] help to reliably constrain light quark FFs to much higher values of z than before, in particular, D π + u+ū . In combination with the LEP and SLAC data, which, at Q 2 = M 2 Z , mainly constrain the total quark singlet fragmentation function, the new precise data at √ S ≃ 10.5 GeV also help to provide some partial  Fig. 1 but now for Q 2 = M 2 Z . Note that here also a bottom-to-π + fragmentation function is shown.
flavor separation from SIA data alone, as they are sensitive to the electrical charge weighted sum of quark FFs. Multiplicities in SIDIS for identified charged pions provide further invaluable experimental input to address this question, see Sec. III C below. In case of D π + g , for the first time, some constraint can be derived from QCD scaling violations in SIA thanks to having now available two precise sets of data at different energy scales Q 2 ≈ 110 GeV 2 and Q 2 = M 2 Z . However, scaling violations for FFs in the relevant medium-to-large z range are fairly mild and also the coverage of the LEP and SLAC data is much more sparse towards high z, which, to some extent, is reflected in the still relatively large uncertainties obtained for D π + g in Fig. 1.
To demonstrate the scale evolution of FFs, we show in Fig. 2 the same zD π + i as in Fig. 1 but now at Q 2 = M 2 Z . Since we are above the bottom threshold Q = m b , we now include also our results for zD π + b = zD π + b in the middle panel of the lower row. To facilitate the comparison of the FFs computed at the two different scales, the dot-dashed lines in Fig. 2 repeat the results for the new, optimum fit at Q 2 = 10 GeV 2 shown in Fig. 1. The FFs of the DSS fit are again denoted as dashed lines. As can be seen, evolution to larger Q 2 reduces the FFs for essentially all relevant z values above about z ≃ 0.15. This trend is reminiscent of the Q 2 evolution of PDFs at not to small values of x, which is not surprising as the LO evolution kernels are essentially the same for the time-like and space-like case. The increase of the FFs at small z is phenomenologically not relevant as their range of applicability is anyhow restricted to z 0.05. The relative uncertainties, again given in the rightmost panels of Fig. 2, are largely similar to those obtained at Q 2 = 10 GeV 2 . Some of the sizable ambiguities at large z are pushed towards smaller z by evolution, most noticeable for D π + g . The overall quality of the fit is summarized in Tab. II where we list all data sets included in our global analysis, as discussed in Sec. II B, along with their individual χ 2 values and the analytically determined normalization shifts according to Eq. (6). We note that the quoted χ 2 values are based only on fitted data points, i.e., after applying the cuts mentioned in Sec. II B, and include the χ 2 penalty from the N i , i.e., the first term in Eq. (5).
Firstly, it is worth mentioning that there is a more than twofold increase in the number of available data points as compared to the original DSS analysis [10]. Secondly, the quality of the global fit has improved dramatically from χ 2 /d.o.f. ≃ 2.2 for DSS, see Tab. II in Ref. [10], to χ 2 /d.o.f. ≃ 1.2 for the current fit. A more detailed comparison reveals that the individual χ 2 values for the SIA data [16-19, 48, 49], which were already included in the DSS fit, have, by and large, not changed significantly. The description of the fully flavor separated data from Opal [20] in the fit favors a rather large normalization shift but has nevertheless deteriorated. Given that this set has only 25 data points, it is the biggest contributor to the total χ 2 . However, in general, flavor-tagged results should not be taken too literally as they lack a proper interpretation and theoretical framework beyond the lowest order as was already pointed out, e.g., in Refs. [10,13]. The biggest improvement concerns the SIDIS multiplicities from Hermes which, in their recently published version [30], are now described very well by the updated fit. Also, the preliminary charged pion multiplicities from Compass [31] and the new SIA data from BaBar [28] and Belle [29] integrate nicely into the global analysis of parton-to-pion FFs.
Finally, and as we will illustrate in detail in Sec. III D below, there is some tension among the pp data sets from RHIC and the LHC, which forced us to introduce a cut p T > 5 GeV on the pion's transverse momentum in the current fit to accommodate both of them. The obtained individual χ 2 values are all reasonable, as can be inferred from Tab. II, with the new Alice data [32] being on the high side, which largely stems from the penalty for the still sizable normalization shift. This large shift reflects the preference of the new Alice data for a smaller gluonto-pion FF than extracted by the original DSS fit based on RHIC Phenix data [21] alone. As a result of the p T cut, the number of pp data in the fit for RHIC has decreased as compared to the DSS analysis. Both the Brahms [45] and Star [46] results at forward pseudorapidities do not pass the p T cut anymore, and, hence, are excluded from the updated fit. Likewise, we do not consider the Alice data taken in pp collisions at √ s = 900 GeV [32], where only a single point would survive the cut in p T .

B. Electron-Position Annihilation Data
In Figs. 3 and 4 we present a detailed comparison of the results of our fit and its uncertainties at both 68% and 90% C.L. with the SIA data already included and newly added to the original DSS analysis [10], respectively. In general, the agreement of the fit with SIA data is excellent in the entire energy and z range covered by the experiments. For Q 2 = M 2 Z , the Delphi data [17] exhibit some mild tension with other sets at the same c.m.s. energy in the largest z bins, resulting in a somewhat higher individual χ 2 value, as can be gathered from Tab. II.
The Belle data [29], shown in the right-hand panel of Fig. 4, provide not only the finest binning in z but also reach the highest z values measured so far. Above z 0.8 one observes an increasing trend for the new fit to overshoot the data, but still within the estimated and growing theoretical uncertainties though. In this kinematic regime one expects large logarithmic corrections, which appear in each order of perturbation theory, to become more and more relevant. It is known how to resum such terms to all orders in the strong coupling [42], and it might be worthwhile to explore their relevance in a future dedicated analysis and whether they could further improve the agreement with data. Resummations also provide a window to non-perturbative contributions to the perturbative series so far little explored. The binning of BaBar data [28] is more sparse towards large z, and a similar trend as for the Belle data is not visible here.
For all the sets shown in Fig. 3, the new fit is able to follow the trend of the data even below the z values included in the analysis (the region indicated by the hatched area). Agreement with BaBar data below the cut z = 0.1 quickly deteriorates though. In this region, the data start to drop while the NLO SIA cross section continues to rise as can be seen in the left-hand panel of Fig. 4. Since the BaBar data are taken at the lowest c.m.s. energy, such an effect is not unexpected and signifies the onset of neglected hadron mass effects in the theoretical framework. In fact, this was the reason for us to choose a somewhat higher cut in z, z > 0.1, than for the other SIA data obtained at higher c.m.s. energies.  [10] (dashed lines) with data sets for inclusive pion production in SIA used in both fits, see Tab. II. The inner and outer shaded bands correspond to new uncertainty estimates at 68% and 90% C.L., respectively. Right-hand side: "(data-theory)/theory" for each of the data sets w.r.t. our new fit (symbols) and the DSS analysis (dashed lines).  [28]; also shown is the result obtained with the DSS fit [10] (dashed line). Right-hand side: same, but now for the Belle data [29]. The lower panels show "(data-theory)/theory" for each of the data sets w.r.t. our new fit (symbols) and the DSS analysis (dashed lines). The inner and outer shaded bands correspond to the new uncertainty estimates at 68% and 90% C.L., respectively.
The Belle experiment did not publish any data below z = 0.2 [29].
HERMES multiplicities π + π -FIG. 5: Comparison of our NLO results for charged pion multiplicities in SIDIS off proton (left panels) and deuteron (right panels) targets with data from the Hermes Collaboration [30]. The inner and outer shaded bands correspond to uncertainty estimates at 68% and 90% C.L., respectively. Also shown are the results obtained with the DSS FFs (dashed lines).
without any refitting or adjusting normalization shifts. The agreement with SIA data is in general very good, except for some small deviations from the recent B factory data, most noticeable in the comparison to BaBar. Contrary to the new analysis, the original DSS fit undershoots both the Belle and BaBar data at high z.
Our estimated uncertainty bands, also shown in Figs. 3 and 4, reflect the accuracy and kinematical coverage of the fitted data. They increase towards both small and large z, similar to the pattern observed for the individual D π + i in Figs. 1 and 2. One should keep in mind that the obtained bands are constrained by the fit to the global set of SIA, SIDIS, and pp data and do not necessarily have to follow the accuracy of each individual set of data.
As was already mentioned in Sec. III A, the SIA data from the LEP and SLAC experiments constrain mainly the total quark singlet fragmentation to pions as up-type and down-type quark couplings to the exchanged Z gauge boson are roughly equal at Q ≃ M Z . The new BaBar and Belle data are dominated by photon exchange and, hence, prefer up-type quark flavors. When combined, this leads to some partial flavor separation. QCD scale evolution between Q 2 ≃ 110 GeV 2 and Q 2 = M 2 Z provides some additional constraints, in particular, also for the gluon FF. The flavor-tagged LEP and SLAC data, listed in Tab. II, are still the best "direct" source of information on the charm-and bottom-to-pion FFs.
Finally, we wish to remark that despite the excellent agreement with all SIA data there are still some issues which require further scrutiny and, perhaps, more detailed comparisons among the different experimental groups. One concern is the question to what extent "feeddown" pions from weak decays contribute to the individual data sets. Different treatments of QED radiative corrections, whose main effect is to lower the "true" c.m.s. energy √ S of the collisions, might be another source of potential tension. For instance, the Belle Collaboration [29] provides only a measurement of the cross section dσ/dz, while all other experiments in SIA scale their quoted results by the total cross section σ tot for e + e − → hadrons. Since Belle cuts on radiative photon events if their energy exceeds a certain threshold, rather than attempting to unfold the radiative QED effects, one has to take this into account when normalizing the Belle data to the conventional 1/σ tot dσ/dz in a global fit.

C. Semi-Inclusive DIS Multiplicities
The most powerful constraint of flavor-separated FFs comes from charged pion multiplicities in SIDIS. Contrary to SIA, which produces π + and π − at equal rates, multiplicities are sensitive to the produced hadron's charge through the choice of the target hadron in DIS.   6: Comparison of our NLO results for π + multiplicities in SIDIS of muons off a deuteron target with preliminary data from the Compass experiment [31] for various bins in x and Q 2 . The inner and outer shaded bands correspond to our uncertainty estimates at 68% and 90% C.L., respectively. Also shown are the results obtained with the DSS FFs (dashed lines). As indicated in the plot, different constant factors are added to the multiplicities M π + µd to distinguish results for different values of Q 2 in the same x bin.
For instance, data taken on a proton target will produce more π + than π − , since u-quarks are more abundant in a proton than d-quarks, and they are also preferred in their coupling to the probing virtual photon due to their larger electrical charge.
Compared to the DSS analysis, where we only had some preliminary set of pion multiplicities on a deuteron target from the Hermes Collaboration at our disposal [15], we can now use their recently published, final set of data for both proton and deuteron targets [30]. In Fig. 5 we illustrate the quality of the new fit with respect to the Hermes data. Shown are the charged pion multiplicities M π ± e,p(d) , which are defined as the ratio of the inclusive pion yield and the total DIS cross section at the same x and Q 2 values (bins) in electron-proton (ep) or electrondeuteron (ed) scattering: The extraction of the FFs requires knowledge of the PDFs of the proton (deuteron) target for which we use the NLO parametrization of the MSTW Collaboration [47] as was already mentioned above. In the fit we consider the projection of the three-dimensional multiplicity data onto the Q 2 dependence for four different bins of the pion's momentum fraction z, which is most sensitive to the quantities we are interested in, the parton-to-pion FFs. The x integrated ratio (9) is also least sensitive to the actual choice of PDFs. We use the standard Mellin technique [23] to pre-calculate look-up tables for each data point at NLO accuracy to speed up the fitting procedure and to facilitate the uncertainty analyses significantly. We recall that at NLO, the relevant hard scattering coefficient functions for SIDIS [7,9] depend in a non-trivial way on both x and z, such that an often used naive approximation, where the x and the z dependence in Eq. (9) is assumed to completely factorize, is bound to fail. Even at LO accuracy such an assumption cannot work as soon as different quark flavors fragment differently into the observed hadron which they do for charged pions.
The agreement between the Hermes data [30] and the updated fit is remarkably better than with the preliminary results [15] used in the DSS analysis, see Fig. 4 in Ref. [10]. This is largely due to the much improved preci-   Fig. 6 but now for π − multiplicities.
sion of the final data, which also exhibit considerably less fluctuations from bin to bin, in particular, for π − . This is also reflected in the total χ 2 for the Hermes data set, which reduces from 188.2 for 64 data points on a proton target in the DSS fit [10] to 175 for 128 data points in the current analysis. In Fig. 5 we compare again also to the result of a calculation based on the DSS FFs, without any re-fitting or adjusting normalizations. As can be seen, the agreement with data is not optimal, and the theory predictions fall short of the data in all bins. Most noticeable is the disagreement at 0.2 ≤ z ≤ 0.3 and 0.6 ≤ z ≤ 0.8 for both π + and π − data. Here, the DSS result is well outside our current uncertainty estimates shown, as before, as shaded bands in Fig. 5. The use of the Hermes multiplicity data as a means of providing a reliable flavor and charge separation for pion FFs in the DSS fit was often questioned in the past because of the smallish Q 2 values of some of the data points. New, still preliminary data from the Compass Collaboration [31], taken at a higher c.m.s. energy, will shed some light on the validity of using a standard, leading-twist pQCD framework at NLO accuracy [7,9] to describe the Hermes multiplicity data for charged pions.
In the present fit we can use charged pion results from Compass obtained on a deuteron target [31]. More specifically, the data are presented as a function of z in 8 bins of x, each subdivided into various bins in Q 2 . In total 199 data points pass our cuts for both π + and π − .
The comparison of the Compass data to the results of our fit is presented in Figs. 6 -8. A very satisfactory agreement is achieved in almost all bins across the entire kinematic regime covered by data, as can be best inferred from Fig. 8, where we show "(data-theory)/theory". The obtained χ 2 /d.o.f. for both π + and π − multiplicities is close to unity, see Tab. II, demonstrating that the low energy Hermes [30] and the Compass [31] data can be described simultaneously and without spoiling the agreement with SIA results. For comparison we show again theoretical results obtained with the DSS FFs (dashed lines), which also agree well with Compass data except for some of the bins corresponding to the lowest Q 2 values. This is in line with the observations for the Hermes data above, where the deviations with DSS were found to be largest at the smallest Q 2 .
As was already mentioned in Sec. III A, the new SIDIS data now favor almost identical u +ū and d +d FFs, i.e., very little or no charge symmetry breaking. This is also preferred by data on the π − /π + ratio in pp collisions which we discuss next.

D. RHIC and LHC Data
The last of the three pillars of our global analysis of parton-to-pion FFs is the wealth of experimental information coming from hadron-hadron collisions, more   specifically, single-inclusive high-p T pion production in pp collisions at BNL-RHIC and CERN-LHC. Compared to the original DSS analysis [10], which mainly made use of the Phenix data for π 0 production at mid rapidity [21], we now have, in addition, results from the Star Collaboration for neutral and charged pions [33][34][35][36] as well as first data from the LHC [32]. Due to the complexity of the underlying hardscattering processes at NLO accuracy [8], the use of a fast, grid-based method such as the Mellin technique, to implement the relevant expressions efficiently and without the need of any approximations is indispensable here. As in various previous analyses [10][11][12]24], and for the implementation of the SIDIS multiplicities in NLO, we adopt the well-tested method based on Mellin moments as described in Ref. [23]. Since inclusive particle spectra at not too large values of p T are dominated by gluoninduced processes in pp collisions [50], the RHIC and LHC data will provide invaluable information on the otherwise only weakly constrained gluon FF D π + g . One of the main results of our updated fits is to reveal a tension between the p T spectra of neutral pions mea-sured at the RHIC experiments and by the Alice Collaboration. In some sense this was already anticipated by comparisons to expectations obtained with the previous DSS FFs, which are known to describe the RHIC data nicely down to p T ≃ 1.5 GeV [10,21] but were found to grossly overshoot recent Alice results at √ s = 7 TeV for essentially all p T values [32,51]. We have tried to accommodate both sets of pp data together by introducing additional freedom to our standard functional form in Eq. (1) but to no avail. In particular, at smallish p T values, below about 5 GeV, the two sets of data appear to be mutually exclusive in a global fit.
Since we do not want to remove either of the data sets from the analysis and, in any case, have no means of judging whether there is a potential experimental inconsistency among the different pp sets, we decided to introduce a cut on the p T of the observed pion. Including only pp data with p T ≥ 5 GeV largely resolves the observed tension between RHIC and LHC data. This is illustrated in Figs. 9 and 10 where we compare to the Phenix [21] and Alice [32] data, respectively. We note that the calculated normalization shift (6) for the Alice 7 TeV data results in a down-shift, outside the experimentally quoted normalization uncertainty, which contributes significantly to the quoted χ 2 value in Tab. II. As we have already hinted at in our discussion of the SIA data in Sec. III B, a different treatment of decay pions by the RHIC and LHC experiments might play some role for the tension observed at p T 5 GeV.
As can be seen, both data sets are well described by the global fit above the introduced p T cut which is indicated by the hatched area in both figures. One also notices the still sizable theoretical scale ambiguity at NLO accuracy, which is indicated in the lower panels of Figs. 9 and 10 and within which the data are consistent with the fit even below the imposed p T cut. The PDF uncertainties, computed with the 90% CL NLO sets from MSTW [47] and also illustrated in the same panels, are much less significant than the scale ambiguities, in particular, for the Alice data.
The resulting global fit is, as always, a compromise of all the data sets included in the analysis and, in particular, mediates between RHIC pp data preferring a larger  Fig. 9 but now for π 0 production in pp collisions at √ S = 7 TeV as measured by the Alice experiment [32].
gluon-to-pion FF and LHC Alice data favoring a smaller D π + g . The net effect is a significantly reduced D π + g as compared to the DSS fit [10], as was already discussed in Sec. III A and illustrated in Fig. 1. Because of the remaining small tension, the estimated uncertainties on D π + g are sizable, despite the available amount of rather precise experimental data from pp collisions. If both RHIC and LHC data would point to a more similar D π + g , the resulting uncertainties would likely to be somewhat smaller, however, the large theoretical scale ambiguities illustrated above still remain. We note that pp data at mid rapidity dominantly probe the gluon FF at mediumto-large z values as was, for instance, demonstrated in Ref. [50] The last two figures give a similar comparison to the Star data [33][34][35][36] for which we adopt, of course, the same p T cut as for the other pp sets. None of these results was included in the DSS analysis. In Fig. 11 we focus on the π − /π + ratio at mid rapidity [35], which is now much better described by the fit than with the DSS FFs. Scale ambiguities partially cancel in the ratio and are much less dramatic than for the individual cross sections, cf. Fig. 9. As was already mentioned, the ratio is sensitive to a potential charge asymmetry or SU(2) breaking, as parametrized by Eq. (3) in out fit. Like for the SIDIS multiplicities, the fit prefers little or no breaking, i.e.,  [35]. The inner and outer shaded bands correspond to uncertainty estimates at 68% and 90% C.L., respectively. Also shown are the results obtained with the DSS FFs (dashed line). Scale and PDF uncertainties are indicated at the base of the plot. N d+d in (3) close to unity. Figure 12 gives on overview of the three other sets of single-inclusive pion data from the Star Collaboration used in the fit [33,34,36], which span different rapidity intervals. Since we fit to the π − /π + ratio shown in Fig. 11 and π + data, we exclude results on the π − cross section to avoid double-counting. The description of the data is very good, even below the p T -cut of 5 GeV, indicating that there is a little bit less of a tension with Alice results than for the Phenix experiment. Calculations based on the DSS FFs (dashed lines) also provide a good description of data.

IV. SUMMARY AND OUTLOOK
We have presented a new, comprehensive global QCD analysis of parton-to-pion fragmentation functions at next-to-leading order accuracy including the latest experimental information. The analyzed data for inclusive pion production in semi-inclusive electron-positron annihilation, deep-inelastic scattering, and proton-proton collisions span energy scales ranging from about 1 GeV up to the mass of the Z boson. The achieved, very satisfactory and simultaneous description of all data sets strongly supports the validity of the underlying theoretical framework based on pQCD and, in particular, the notion of factorization and universality for parton-to-pion fragmentation functions.
Compared to our previous analysis, which was based FIG. 12: Comparison of our NLO results for single-inclusive π 0 and π + production in different rapidity ranges in pp collisions at √ S = 200 GeV with the corresponding Star data [33,34,36]. The inner and outer shaded bands correspond to uncertainty estimates at 68% and 90% C.L., respectively. Also shown are the results obtained with the DSS FFs (dashed lines). Note that some of the results are rescaled for clarity.
on much less precise experimental input and to which we have made extensive comparisons throughout this work, we now obtained a significantly better fit, as measured in terms its the global χ 2 , using the same functional form with only a few additional fit parameters. While most of the favored and unfavored quark-to-pion fragmentation functions are by and large similar to our previous results, the reduced amount of pions stemming for the hadronization of gluons is a noteworthy outcome of the new analysis. This finding was driven by first data from the CERN-LHC experiments, which, surprisingly, turned out to be mutually incompatible with previously available data obtained in lower center-of-mass system energy collisions at BNL-RHIC. To remedy this tension in our fit, we were forced to introduce a lower cut on the transverse momentum of the produced pions in proton-proton collisions. We have argued that it should be worthwhile for the experiments to compare in detail their procedures to determine pion yields as, for instance, different cuts for secondary pions from decays of other, heavier mesons perhaps have some numerical impact. We believe that such a contamination from feed-down pions might show up most prominently at small transverse momenta, where we currently observe the tension between RHIC and LHC data. We wish to mention that in the quark sector, the new data do not favor any charge symmetry violation be-tween the total up-and down-quark fragmentation functions, contrary to our previous fit.
We have also performed a, what we believe, first reliable and trustworthy estimate of uncertainties for partonto-pion fragmentation functions based on the standard iterative Hessian method. This was made possible by the wealth of new data included in our updated global analysis. The obtained uncertainties are still sizable and range at best from about ten percent to twenty-five percent for the total u-quark and gluon fragmentation function, respectively, in the kinematic regions covered by data and they quickly deteriorate beyond. A new asset of the current analysis is the analytic procedure to determine the optimum normalization shift for each data set in the fit, which greatly facilitated the global fitting procedure The newly obtained pion fragmentation functions and their uncertainty estimates will be a crucial ingredient in future global analyses of both helicity and transversemomentum dependent parton densities, which heavily draw on data with identified pions in the final-state. Our results will also serve as the baseline in heavy ion and proton-heavy ion collisions, where one of the main objectives is to quantify and understand possible modifications of hadron production yields by the nuclear medium. Also, the current analysis framework will be adopted for updates of the parton-to-kaon fragmentation functions, which we will pursue once all the promised sets of new data eventually become available. Since pions and kaons constitute by far the largest fraction in frequently measured yields of unidentified charged hadrons, a precise determination of their respective, optimum sets of fragmentation functions, including reliable uncertainty estimates, is critical to determine the room left for other hadrons, such as protons, in a future global analysis of charged hadron data.
Further improvements of parton-to-pion fragmentation functions from the theory side should include an improved treatment of heavy quark-to-pion fragmentation functions, likely along similar lines as for heavy flavor parton densities. Also, the impact of higher order corrections beyond the next-to-leading order accuracy should be explored. On the one hand, it is already possible to perform a next-to-next-to-leading order analysis of electron-positron annihilation data, and, on the other hand, the theoretical framework for all-order resummations of potentially large logarithmic corrections is available. On the experimental side, RHIC and the LHC will continue to provide new data on identified hadron spectra but it should be also worthwhile to explore the potential of future accelerator projects, such as an Electron-Ion Collider currently pursued in the U.S. [52], to further our knowledge of fragmentation functions and the physics behind hadronization.