First Monte Carlo analysis of fragmentation functions from single-inclusive e + e − annihilation
N. Sato, J.J. Ethier, W. Melnitchouk, M. Hirai, S. Kumano, A. Accardi
JJLAB-THY-16-2327, KEK-TH-1920, J-PARC-TH-0060
First Monte Carlo analysis of fragmentation functionsfrom single-inclusive e + e − annihilation N. Sato, J. J. Ethier, W. Melnitchouk, M. Hirai, S. Kumano,
4, 5 and A. Accardi
1, 6 Jefferson Lab, Newport News, Virginia 23606, USA College of William and Mary, Williamsburg, Virginia 23187, USA Nippon Institute of Technology, Saitama 345-8501, Japan High Energy Accelerator Research Organization (KEK),1-1, Oho, Tsukuba, Ibaraki 305-0801, Japan J-PARC Center, 203-1, Shirakata, Tokai, Ibaraki, 319-1106, Japan Hampton University, Hampton, Virginia 23668, USA
Jefferson Lab Angular Momentum (JAM) Collaboration (Dated: July 16, 2018)
Abstract
We perform the first iterative Monte Carlo (IMC) analysis of fragmentation functions constrainedby all available data from single-inclusive e + e − annihilation into pions and kaons. The IMCmethod eliminates potential bias in traditional analyses based on single fits introduced by fixingparameters not well contrained by the data and provides a statistically rigorous determination ofuncertainties. Our analysis reveals specific features of fragmentation functions using the new IMCmethodology and those obtained from previous analyses, especially for light quarks and for strangequark fragmentation to kaons. a r X i v : . [ h e p - ph ] S e p . INTRODUCTION Understanding the generation of hadrons from quarks and gluons (partons) remains afundamental challenge for strong interaction physics. High-energy collisions of hadrons orleptons offers the opportunity to study the formation of mesons and baryons from partonsproduced in hard collisions [1, 2]. While the hard scattering process can be computed pertur-batively from the underlying QCD theory, the hadronization of the quarks and gluons occursover long distances, and provides a unique window on nonperturbative QCD dynamics [3].Within the collinear factorization framework [4], the formation of hadrons is characterizedby universal nonperturbative fragmentation functions (FFs), which in an infinite momentumframe can be interpreted as probability distributions of specific hadrons h produced witha fraction z of the scattered parton’s longitudinal momentum or energy. As in the caseof parton distribution functions (PDFs), which describe the quark and gluon momentumdistributions inside hadrons, the nonperturbative FFs are presently not calculable from firstprinciples, and must be determined phenomenologically from QCD-based analyses of high-energy scattering data or from QCD-inspired nonperturbative models [5].In addition to providing information on the fundamental hadronization process, FFs arealso indispensable tools for extracting information on the partonic structure of the nucleonfrom certain high-energy processes, such as semi-inclusive deep-inelastic scattering (SIDIS)of leptons from nucleons. Here, assuming factorization of the scattering and hadronizationsubprocesses, the SIDIS cross section can be expressed in terms of products of PDFs andFFs summed over individual flavors. The selection of specific hadrons in the final state, suchas π ± or K ± , then allows separation of the momentum and spin PDFs for different flavors.The need for well-constrained FFs, especially for kaon production, has recently beenhighlighted [6–8] in global analyses of polarized SIDIS observables used to determine thestrange quark contribution ∆ s to the spin of the nucleon. Inclusive deep-inelastic lepton–nucleon scattering data alone are incapable of determining this without additional input fromtheory, such as the assumption of SU(3) symmetry, or other observables. Kaon production inpolarized SIDIS in principle is such an observable, involving a new combination of polarized u , d and s quark PDFs, which, when combined with the inclusive data, allow each of theflavor distributions to be determined – providing the FFs are known.As pointed out by Leader et al. [7], however, the variation between the strange-to-kaon2Fs from different analyses is significant and can lead to qualitatively different conclusionsabout the magnitude and even sign of the ∆ s distribution. In particular, analysis [7, 9] ofthe polarized SIDIS data using the DSS [10] parametrization of FFs, together with inclusiveDIS polarization asymmetries, suggests a positive ∆ s at intermediate x values, x ∼ . − . s at all x obtained from inclusive DIS data alone,assuming constraints on the weak baryon decays from SU(3) symmetry [11]. Employinginstead the HKNS [12] FF parametrization, in which the strange fragmentation to kaons isseveral times smaller in some regions of z compared with that from the DSS [10] fit, yieldsa negative ∆ s consistent with the inclusive-only analyses [8]. It is crucial, therefore, tounderstand the origin of the differences in the magnitudes and shapes of the strange, as wellas other, FFs found in the different analyses before one can draw reliable conclusions aboutthe strange quark content of the nucleon extracted from analyses including SIDIS data.Differences between FFs can come from a variety of sources, including different data setsused in the analyses (single-inclusive e + e − annihilation, SIDIS, inclusive hadron productionin pp collisions), the choice of parametrization for the FFs, assumptions about FFs that arenot well constrained by data, or even the presence of local minima in the fitting procedure.Most of the analyses to date have been performed at next-to-leading order (NLO) accuracyin the strong coupling constant [6–8, 10, 12–18], although more recent studies have exploredthe effects of incorporating next-to-next-to-leading order (NNLO) corrections [19], as wellas other theoretical developments such as threshold resummation [20–22] and hadron masseffects [22].A common feature of all existing FF analyses is that they are obtained from single fits,using either e + e − single-inclusive annihilation (SIA) data alone, or in combination withunpolarized SIDIS and inclusive hadron production in pp collisions. In order to addresssome of the questions raised by the recent ambiguities in the strange quark FFs and theirimpact on the ∆ s determination, in this paper we go beyond the standard fitting paradigmby performing the first Monte Carlo (MC) analysis of FFs. In particular, we extend themethodology of the iterative Monte Carlo (IMC) approach introduced in Ref. [11] for theanalysis of spin-dependent PDFs to the case of FFs.The virtue of the IMC approach is that it allows for a full exploration of the parameterspace when sampling initial priors for any chosen parametric form for the fitting function. Itthereby eliminates any bias introduced by fine-tuning or fixing specific parameters that are3ot well contrained by the data, a practice often employed to control single fits. Furthermore,the conventional polynomial-type parametrization choice can have multiple solutions thatlead to various local minima in the χ landscape, whereas the IMC technique statisticallysurveys all possible solutions, thereby avoiding the fit being stuck in false minima.A further important advantage of the IMC technology is in the extraction of uncertaintieson the FFs. In standard analyses the theoretical errors are typically determined using theHessian [12] or Lagrange multiplier methods [10], in which a tolerance parameter ∆ χ isintroduced to satisfy a specific confidence level (CL) of a χ probability density function with N degrees of freedom. In the IMC framework, the need for tolerance criteria is eliminatedentirely and the uncertainties are extracted through a robust statistical analysis of the MonteCarlo results.As a first IMC analysis of FFs, we confine ourselves to the case of charged pion andkaon production in e + e − SIA, using all available π ± and K ± cross section data from DESY[23–26], SLAC [27–31], CERN [32–36], and KEK [37], as well as more recent, high-precisionresults from the Belle [38, 39] and BaBar [40] Collaborations at KEK and SLAC, respec-tively. Although SIA data in principle only constrain the sum of the quark and antiquarkdistributions, we also make use of flavor-tagged data [33] which allow separation of hadronproduction from heavy and light quarks. In addition, the availability of data over a range ofkinematics, from relatively low center-of-mass energies Q ≈
10 GeV up to the Z -boson pole, Q ≈
91 GeV, allows for the separation of the up- and down-type FFs due to differencesin the quark–boson couplings in the γ and Z channels [18]. To ensure proper treatmentof data at z ∼
1, we systematically apply correct binning by integrating over each z bins,rather than taking bin averages as in previous analyses. We also studied the z cuts on thedata in different channels that need to be applied at low z values, below which the collinearframework breaks down and our analysis is not expected to be reliable.Note that our aim here is not so much the definitive determination of FFs, which wouldrequire inclusion of all possible processes that have sensitivity to FFs, but rather to explorethe application of the IMC methodology for FFs to determine the maximal information thatcan be extracted from the basic e + e − SIA process alone. The lessons learned here will beused in subsequent analyses of the entire global set of SIA and other high-energy scatteringdata to provide a more definitive determination of the individual FFs.We begin in Sec. II by reviewing the formalism for the e + e − annihilation into hadrons,4ncluding a summary of the SIA cross sections at NLO and Q evolution of the fragmentationfunctions. To improve the computational efficiency we perform the numerical calculationsin moment space, recontructing the momentum dependence of the fragmentation functionsusing inverse Mellin transforms. The methodology underpinning our global analysis is pre-sented in Sec. III, where we describe the parametrizations employed and the treatment ofuncertainties. This section also outlines the essential features of the IMC method used toperform the fits to the data, highlighting several improvements in the methodology comparedto that introduced originally in the global analysis of the JAM spin-dependent PDFs [11].The experimental data sets analyzed in this study are summarized in Sec. IV, and the resultsof our analysis presented in Sec. V. We compare the fitted cross sections with all available e + e − data, for both inclusive and flavor-tagged cross sections, finding good overall χ valuesfor both pion and kaon production. We illustrate the convergence of the iterative procedurefor the favored and unfavored FFs, the latter being partially constrained by the flavor-taggeddata. The shapes and magnitudes of the FFs from our IMC analysis are compared and con-trasted with those from previous global fits, highlighting important differences in the lightquark sector and for quark fragmentation to kaons. Finally, in Sec. VI we summarize ourfindings and preview future extensions of the present analysis. II. FORMALISMA. Cross section and fragmentation functions
The e + e − → hX cross section is typically measured as a function of the variable z = 2 p h · q/Q , where p h is the momentum of the detected hadron h and q is the momentumof the exchanged photon or Z -boson with invariant mass Q = (cid:112) Q . In the e + e − center-of-mass frame, z = 2 E h /Q can be interpreted as the momentum fraction of the parent quarkcarried by the produced hadron. For a given hadron h the experimental z distribution isusually given as F h ( z, Q ) = 1 σ tot dσ h dz ( z, Q ) , (1)5hich we shall refer to as the empirical fragmentation function for a given hadron of type h . In Eq. (1) the total inclusive e + e − cross section σ tot can be calculated at NLO as σ tot ( Q ) = (cid:88) q πα Q ˜ e q (cid:0) a s ( µ ) (cid:1) + O ( a s ) , (2)where α = e / π is the electromagnetic fine structure constant and a s ( µ R ) ≡ α s ( µ R ) / π ,with the strong coupling constant α s evaluated at the ultraviolet renormalization scale µ R .The index q runs over the active quark flavors allowed by the hard scale Q , and we introducethe shorthand notation for the charges˜ e q = e q + 2 e q g qV g eV ρ ( Q ) + (cid:0) g e A + g e V (cid:1) (cid:0) g q A + g q V (cid:1) ρ ( Q ) . (3)Here the quark vector and axial vector couplings are given by g qV = − sin θ W and g qA = + for the q = u, c flavors, while for the q = d, s, b flavors these are g qV = − + sin θ W and g qA = − . Similarly, the electron vector and axial vector couplings are given by g eV = − + 2 sin θ W and g eA = − , respectively. Because the weak mixing angle sin θ W is ≈ / ρ and ρ arise from γZ interference and Z processes,respectively, and are given by ρ ( Q ) = 14 sin θ W cos θ W Q ( M Z − Q )( M Z − Q ) + M Z Γ Z , (4a) ρ ( Q ) = 1 (cid:0) θ W cos θ W (cid:1) Q ( M Z − Q ) + M Z Γ Z , (4b)where M Z and Γ Z are the mass and width of the Z boson, respectively.Within the collinear factorization framework, the empirical fragmentation function F h ( z, Q ) can be approximately calculated in terms of quark fragmentation functions intohadrons, F h ( z, Q ) ≈ F h coll ( z, Q ) = (cid:88) i (cid:2) H i ⊗ D hi (cid:3) ( z, Q , µ , µ ) + O ( a s ) , (5)where “ ⊗ ” refers to the standard convolution integral [ H ⊗ D ]( z ) = (cid:82) z ( d ˆ z/ ˆ z ) H (ˆ z ) D ( z/ ˆ z ),and the sum runs over all parton flavors i = q, ¯ q, g . Here H i is the short-distance hard crosssection calculable in fixed-order perturbative QCD, and D hi is the partonic fragmentationfunction. As discussed below, the quark contributions H q depend on the charges ˜ e q , whilethe gluon contribution is independent of the charges.6t NLO in the MS scheme (which we use throughout in this analysis), the hard crosssection can be written H i (ˆ z, Q , µ , µ ) = H (0) i (ˆ z, Q , µ , µ ) + a s ( µ R ) H (1) i (ˆ z, Q , µ , µ ) + O ( a s ) , (6)where ˆ z is the partonic energy fraction carried by the outgoing hadron. As in Eq. (2), µ R isthe renormalization scale stemming from regularization of the ultraviolet divergences in thevirtual graphs that contribute to H (1) i , while µ FF is a factorization scale associated with theFF D hi . Note that the dependence of the convolution integral in Eq. (5) on the scales µ R and µ FF is a remnant of the fixed-order perturbative QCD approximation to F coll , which willbe cancelled by inclusion of higher order terms in the perturbative series. At leading orderin a s , the 2 → z = z , so that H (0) i is proportional to δ (ˆ z − z ). Athigher orders, additional QCD radiation effects open up the phase space for the outgoingfragmenting parton such that ˆ z varies between z and 1.The partonic FF D hi can be interpreted as the number density to find a hadron of type h in the jet originating from the parton i with momentum fraction ˆ z [41]. As for PDFs, FFsare sensitive to ultraviolet divergences, and after renormalization they acquire dependenceon the scale µ FF . (The subscript “FF” denotes the final state factorization scale, in contrastto the initial state factorization scale in PDFs.) In practice, to optimize the perturbativeexpansion of the hard cross section, we set µ R = µ FF = Q . However, for completeness weleave the dependence of µ R and µ FF in Eq. (5) and below explicit. In general, variationof the scales around Q allows one to assess the uncertainty in the perturbative expansion.For instance, in Ref. [19] a significant reduction of the scale dependence was found with theinclusion of the NNLO corrections. B. Scale dependence
In perturbative QCD the scale dependence of the FFs is described by the evolutionequations, dD hi (ˆ z, µ ) d ln( µ ) = (cid:2) P ij ⊗ D hj (cid:3) (ˆ z, µ ) , (7)where P ij are the timelike i → j splitting functions. Since the FFs cannot be calculatedfrom first principles, the ˆ z dependence is fitted to the data at some input scale µ = Q .7he latter is chosen at the lowest possible value where a perturbative QCD description canbe applied in order to minimize errors induced by backward evolution from the truncationof the perturbative series.The simplest approach to solving the evolution equations (7) is to use one of severalnumerical approximation techniques to solve the integro-differential equations directly inˆ z space [42]. Alternatively, as discussed in Ref. [11], it can be more efficient to solve theequations in Mellin moment space, where the N -th Mellin moment of a function f ( z ) isdefined as f ( N ) = (cid:90) dz z N − f ( z ) , (8)and similarly for all other moments of functions denoted in boldface. In this framework theconvolution integrals in Eqs. (6) and (7) can be rendered as ordinary products of the Mellinmoments, F h coll ( N, Q ) = (cid:88) i H i ( N, Q , µ , µ ) D hi ( N, Q , µ , µ ) + O ( a s ) , (9)and d D hi ( N, µ ) d ln( µ ) = P ij ( N, µ , µ ) D hj ( N, Q , µ , µ ) . (10)The evolution equations for D hi can be solved using the methods described in Ref. [43], andthe hadronic fragmentation function in z -space can be obtained using the inverse Mellintransform, F h coll ( z, Q ) = 12 πi (cid:90) C dN z − N F h coll ( N, Q ) . (11)The main advantage of the Mellin techniques is the improvement in speed in the evaluationof the observables and evolution equations. Another advantage is that the experimentalcross sections are typically presented as averaged values over bins of z . Such averaging,between z min and z max , can be simply done analytically, (cid:10) F h coll ( z, Q ) (cid:11) z bin = 1( z max − z min ) 12 πi (cid:90) C dN (cid:0) z − N max − z − N min (cid:1) − N F h coll ( N, Q ) , (12)without deteriorating the numerical performance. In contrast, such advantage does not existif one evaluates F h coll ( z, Q ) and solves the DGLAP evolution equations directly in z space844]. In practice, at small z the bins sizes are quite small and taking the central z valuesmight be appropriate. However, at large z the bin sizes increase and, depending on theprecision of the measured cross sections, the averaging step becomes important.For clarity, we express the Mellin moments of the hard factor in Eq. (9) in terms ofunnormalized hard factors (cid:102) H i , H q ( N, Q , µ , µ ) = ˜ e q (cid:80) q (cid:48) ˜ e q (cid:48) (cid:102) H q ( N, Q , µ , µ )(1 + 4 a s ( µ )) , (13a) H g ( N, Q , µ , µ ) = (cid:102) H g ( N, Q , µ , µ )(1 + 4 a s ( µ )) , (13b)where the charge factors for the gluon moments cancel. The perturbative expansion of (cid:102) H i is then given by (cid:102) H q ( N, Q , µ , µ ) = 1 + a s ( µ ) (cid:102) H (1) q ( N, Q , µ , µ ) + O ( a s ) , (14a) (cid:102) H g ( N, Q , µ , µ ) = a s ( µ ) (cid:102) H (1) g ( N, Q , µ , µ ) + O ( a s ) , (14b)where the gluon contribution begins at NLO. Physically, this corresponds to gluon frag-mentation into hadrons from real QCD radiation that occurs at NLO. For completeness, inAppendix A we list the formulas for (cid:102) H (1) q,g at NLO.To solve the evolution equations in Eq. (9), we follow the conventions of Ref. [43], whichwe briefly summarize here. For convenience we work in a flavor singlet and nonsinglet basis,in which we define the flavor combinations D h ± = D hu ± − D hd ± , (15a) D h ± = D hu ± + D hd ± − D hs ± , (15b) D h ± = D hu ± + D hd ± + D hs ± − D hc ± , (15c) D h ± = D hu ± + D hd ± + D hs ± + D hc ± − D hb ± , (15d) D h ± = D hu ± + D hd ± + D hs ± + D hc ± + D hb ± − D ht ± , (15e) D h ± = D hu ± + D hd ± + D hs ± + D hc ± + D hb ± + D ht ± , (15f)where D hq ± are the Mellin moments of the charge conjugation-even and -odd FFs D hq ± ( z, Q ) = D hq ( z, Q ) ± D h ¯ q ( z, Q ). Depending on the number of active flavors n f , oneneeds to consider only the equations up to D ± n f − , otherwise the system becomes degenerate.9he evolution equations in this basis can be expressed as ∂ D h ± j ∂ ln µ = P ± NS D h ± j , (16a) ∂ D h − ∂ ln µ = P − NS D h − (16b) ∂∂ ln µ D h + D hg = P qq P qg P gq P gg D h + D hg , (16c)with the splitting functions in Mellin space P ij listed in Appendix B. An important ob-servation here is that all the “+” FFs maximally couple to the gluon FFs, while the “ − ”functions decouple completely. In particular, if one consider observables that depend onlyon “+” combinations, then the “ − ” components can be ignored.In our analysis we use an independent implementation of the evolution equations in Mellinspace as described in Ref. [43], finding excellent agreement with existing evolution codes. III. METHODOLOGYA. Input scale parametrization
In choosing a functional form for the FFs, it is important to note that the SIA observablesare sensitive only to the charge conjugation-even quark distributions D hq + ( z, Q ) and thegluon FF D hg ( z, Q ). These couple maximally in the Q evolution equations, while thecharge conjugation-odd combinations D hq − ( z, Q ) decouple entirely from both D hq + ( z, Q ) and D hg ( z, Q ). In our analysis we therefore seek only to extract the D hq + and gluon distributions,and do not attempt to separate quark and antiquark FFs. This would require additionaldata, such as from semi-inclusive deep-inelastic hadron production, which can provide afilter on the quark and antiquark flavors.As a reference point, we consider a “template” function of the formT( z ; a ) = M z α (1 − z ) β (cid:82) dz z α (1 − z ) β , (17)where a = { M, α, β } is the vector of shape parameters to be fitted. The denominator ischosen so that the coefficient M corresponds to the average momentum fraction z .Using charge conjugation symmetry, one can relate D h + q + = D h − q + , D h + g = D h − g , (18)10or all partons. For pions we further use isospin symmetry to set the u + and d + functionsequal, while keeping the remaining FFs independent. Since the u + and d + distributionsmust reflect both the “valence” and “sea” content of the π + , we allow two independentshapes for these, while a single template function should be sufficient for the heavier flavorsand the gluon, D π + u + = D π + d + = T( z ; a πud ) + T( z ; a (cid:48) πud ) , (19a) D π + s + , c + , b + , g = T( z ; a πs, c, b, g ) . (19b)The additional template shape for the u + or d + increases the flexibility of the parametrizationin order to accomodate the distinction between favored (“valence”) and unfavored (“sea”)distributions, having different sets of shape parameters a πud and a (cid:48) πud .For the kaon the s + and u + FFs are parametrized independently because of the massdifference between the strange and up quarks. Since these contain both valence and seastructures, to improve the flexibility of the parametrization we use two template shapeshere, and one shape for each of the other distributions, D K + s + = T( z ; a Ks ) + T( z ; a (cid:48) Ks ) , (20a) D K + u + = T( z ; a Ku ) + T( z ; a (cid:48) Ku ) , (20b) D K + d + , c + , b + , g = T( z ; a Kd, c, b, g ) . (20c)The total number of free parameters for the kaon FFs is 24, while for the pions the numberof parameters is 18.For the heavy quarks c and b we use the zero-mass variable flavor scheme and activate theheavy quark distributions at their mass thresholds, m c = 1 .
43 GeV and m b = 4 . Q evolution we use the “truncated” solution in Ref. [43], which is more consistent withfixed-order calculations. Finally, the strong coupling is evaluated by solving numerically the β -function at two loops and using the boundary condition at the Z pole, α s ( m Z ) = 0 . B. Iterative Monte Carlo fitting
In all previous global analyses of FFs, only single χ fits have been performed. In thiscase it is common to fix by hand certain shape parameters that are difficult to constrainby data in order to obtain a reasonable fit. However, since some of the parameters and11istributions are strongly correlated, this can bias the results of the analysis. In addition,there is no way to determine a priori whether a single χ fit will become stuck in any oneof many local minima. The issues of multiple solutions can be efficiently avoided throughMC sampling of the parameter space, which allows exploration of all possible solutions.Since this study is the first MC-based analysis of FFs, we briefly review the IMC procedure,previously introduced in the JAM15 analysis of polarized PDFs [11], and highlight severalimportant new features.In the IMC methodology, for a given observable O the expectation value and varianceare defined by E[ O ] = (cid:90) d m a P ( a | data) O ( a ) , (21)V[ O ] = (cid:90) d m a P ( a | data) ( O ( a ) − E[ O ]) , (22)respectively, where a is the m -component vector representing the shape parameters of theFFs. The multivariate probability density P ( a | data) for the parameters a conditioned bythe evidence ( e.g. , the data) can be written as P ( a | data) ∝ L (data | a ) × π ( a ) , (23)where π ( a ) is the prior and L (data | a ) is the likelihood . In our analysis π ( a ) is initially setto be a flat distribution. For L (data | a ) we assume a Gaussian likelihood, L (data | a ) ∝ exp (cid:18) − χ ( a ) (cid:19) , (24)with the χ function defined as χ ( a ) = (cid:88) e (cid:88) i (cid:32) D ( e ) i N ( e ) i − T ( e ) i α ( e ) i N ( e ) i (cid:33) + (cid:88) k (cid:16) r ( e ) k (cid:17) . (25)Here D ( e ) i and T ( e ) i represent the data and theory points, respectively, and α ( e ) i are theuncorrelated systematic and statistical experimental uncertainties added in quadrature. Thenormalization uncertainties are accounted for through the factor N ( e ) i , defined as N ( e ) i = 1 − (cid:88) k r ( e ) k β ( e ) k,i D ( e ) i . (26)12ere β ( e ) k,i is the k -th source of point-to-point correlated systematic uncertainties in the i -thbin, and r ( e ) k the related weight, treated as a free parameter. In order to fit the r ( e ) k values,a penalty must be added to the definition of the χ , as in the second term of Eq. (25).Clearly the evaluation of the multidimensional integrations in Eqs. (21) and (22) is notpractical, especially when O is a continuous function such as in the case of FFs. Insteadone can construct an MC representation of P ( a | data) such that the expectation value andvariance can be evaluated as E[ O ] = 1 n n (cid:88) k =1 O ( a k ) , (27)V[ O ] = 1 n n (cid:88) k =1 ( O ( a k ) − E[ O ]) , (28)where the parameters { a k } are distributed according to P ( a | data), and n is the number ofpoints sampled from the distribution P ( a | data).Our approach to constructing the Monte Carlo ensemble { a k } is schematically illustratedin Fig. 1. The steps in the IMC procedure can be summarized in the following workflow:1. Generation of the priors
The priors are the initial parameters that are used as guess parameters for a givenleast-squares fit. The resulting parameters from the fits are called posteriors . Duringthe initial iteration, a set of priors is generated using a flat sampling in the parameterspace. The sampling region is selected for the shape parameters α > − . β > β restricts thedistributions to be strictly zero in the z → α and β are selected to cover typical ranges observed in previous analysis [10, 12, 16]. Note,however, that the posteriors can be distributed outside of the initial sampling region,if this is preferred by the data.For each subsequent iteration, the priors are generated from a multivariate Gaussiansampling using the covariance matrix and the central parameters from the priors ofthe previous iteration. The central parameters are chosen to be the median of thepriors, which is found to give better convergence compared with using the mean. Thissampling procedure further develops the JAM15 methodology [11], where the priorswere randomly selected from the previous iteration posteriors. This allows one to13onstruct priors that are distributed more uniformly in parameter space as opposedto priors that are clustered in particular regions of parameter space. The latter canpotentially bias the results if the number of priors is too small.2. Generation of pseudodata sets
Data resampling is performed by generating pseudodata sets using Gaussian smear-ing with the mean and uncertainties of the original experimental data values. Eachpseudodata point (cid:101) D i is computed as (cid:101) D i = D i + R i α i , (29)where for each experiment D i and α i are as in Eq. (25), and R i is a randomly gener-ated number from a normal distribution of unit width. A different pseudodata set isgenerated for each fit in any given iteration in the IMC procedure.3. Partition of pseudodata sets for cross-validation
To account for possible over-fitting, the cross-validation method is incorporated. Eachexperimental pseudodata set is randomly divided 50% /
50% into “training” and “val-idation” sets. However, data from any experiment with fewer than 10 points are notpartitioned and are entirely included in the training set.4. χ minimization and posterior selection The χ minimization procedure is performed with the training pseudodata set usingthe Levemberg-Marquardt lmdiff algorithm [45]. For every shift in the parametersduring the minimization procedure, the χ values for both training and validationare computed and stored along with their respective parameter values, until the bestfit for the training set is found. For each pseudodata set, the parameter vector thatminimizes the χ of the validation is then selected as a posterior.5. Convergence criterion
The iterative approach of the IMC is similar to the strategy adopted in the MC VEGASintegration [46]. There, one constructs iteratively a grid over the parameter space suchthat most of the sampling is confined to regions where the integrand contributes themost, a procedure known as importance sampling . Once the grid is prepared, a largeamount of samples is generated until statistical convergence of the integral is achieved.14n Ref. [11] the convergence of the MC ensemble { a k } was estimated using the χ distribution. While such an estimate can give some insight about the convergence ofthe posteriors, it is somewhat indirect as it does not involve the parameters explicitly.In the present analysis, we instead estimate the convergence of the eigenvalues of thecovariance matrix computed from the posterior distributions. To do this we constructa measure given by V = (cid:89) i (cid:112) W i , (30)where W i are the eigenvalues of the covariance matrix. The quantity V can be inter-preted in terms of the hypervolume in the parameter space that encloses the posteriors,and is analogous to the ensemble of the most populated grid cells in a given iterationof the VEGAS algorithm [46]. The IMC procedure is then iterated starting from step1, until the volume remains unchanged.6. Generation of the Monte Carlo FF ensemble
When the posteriors volume has reached convergence, a large number of fits is per-formed until the mean and expectation values of the FFs converge. The goodness-of-fitis then evaluated by calculating the overall single χ values per experiment accordingto χ e ) = (cid:88) i (cid:32) D ( e ) i − E [ T ( e ) i ] /E [ N ( e ) i ] α ( e ) i (cid:33) , (31)which allows a direct comparison with the original unmodified data.Finally, note that while the FF parametrization adopted here is not intrinsically more flexiblethan in other global analyses, the MC representation is significantly more versatile andadaptable in describing the FFs. Indeed, the resulting averaged central value of the FFs asa function of z is a linear combination of many functional shapes, effectively increasing theflexibility of the parametrization. 15 V. DATA SETS
In the current analysis we use all available data sets from the single-inclusive annihilationprocess e + e − → hX , for h = π ± and K ± mesons. Table I summarizes the various SIAexperiments, including the type of observable measured (inclusive or tagged), center-of-massenergy Q , number of data points, and the χ values and fitted normalization factors for eachdata set. Specifically, we include data from experiments at DESY (from the TASSO [23–25]and ARGUS [26] Collaborations), SLAC (TPC [27–29], HRS [30], SLD [31] and BaBar [40]Collaborations), CERN (OPAL [32, 33], ALEPH [34] and DELPHI [35, 36] Collaborations)and KEK (TOPAZ [37] and Belle [38, 39] Collaborations). Approximately half of the 459 π ± data points and 391 K ± data points are near the Z -boson pole, Q ≈ M Z , while the mostrecent, high-precision Belle and BaBar data from the B -factories are at Q (cid:39) . z region, and reveal clearer scaling violation effects compared with the previous higher-energymeasurements.In the TPC, OPAL, DELPHI and SLD experiments, light-quark and heavy-quark eventswere separated by considering the properties of final-state hadrons. In the SLD experiment,for example, events from the primary c and b quarks were selected by tracks near theprimary interaction point. For each secondary vertex, the total transverse momentum andinvariant mass were obtained, after which the data were separated into c - and b -tagged eventsdepending on the masses and transverse momenta. Some events without the secondaryvertex were considered as light-quark ( u, d, s )-tagged if a track did not exist with an impactparameter exceeding a certain cutoff value. Other tagged data sets used different techniquesfor selecting the quark-tagged events. In the OPAL experiment, separated probabilitiesfor u , d and s quark fragmentation were also provided, which in practice provide valuableconstraints on the flavor dependence in the light-quark FFs.For the Belle measurements [38], the data are provided in the form dσ h /dz , and caremust be taken when converting this to the hadronic FF in Eq. (1). The fragmentationenergy scale Q/ Q/ ABLE I: Single-inclusive e + e − annihilation experiments used in this analysis, inluding the typeof observable (inclusive or tagged), center-of-mass energy Q , number of data points N dat , averagefitted correlated normalization (when different from “1”), and χ values for pion and kaon produc-tion. Note that the normalization factors for the TASSO data, indicated by ( ∗ ) in the table, arein the range 0.976 – 1.184 for pions and 0.891 – 1.033 for kaons. For the BaBar pion data [40]the “prompt” data set is used in the fit discussed in this paper, with normalization and χ valuesobtained using the “conventional” data set in parentheses. experiment ref. observable Q pions kaons(GeV) N dat norm. χ N dat norm. χ ARGUS [26] inclusive 9.98 35 1.024(1.058) 51.1(55.8) 15 1.007 8.5Belle [38, 39] inclusive 10.52 78 0.900(0.919) 37.6(21.7) 78 0.988 10.9BaBaR [40] inclusive 10.54 39 0.993(0.948) 31.6(70.7) 30 0.992 4.9TASSO [23–25] inclusive 12-44 29 ( ∗ ) 37.0(38.8) 18 ( ∗ ) 14.3TPC [27–29] inclusive 29.00 18 1 36.3(57.8) 16 1 47.8 uds tag 29.00 6 1 3.7( 4.6) b tag 29.00 6 1 8.7( 8.6) c tag 29.00 6 1 3.3( 3.0)HRS [30] inclusive 29.00 2 1 4.2( 6.2) 3 1 0.3TOPAZ [37] inclusive 58.00 4 1 4.8( 6.3) 3 1 0.9OPAL [32, 33] inclusive 91.20 22 1 33.3(37.2) 10 1 6.3 u tag 91.20 5 1.203(1.203) 6.6( 8.1) 5 1.185 2.1 d tag 91.20 5 1.204(1.203) 6.1( 7.6) 5 1.075 0.6 s tag 91.20 5 1.126(1.200) 14.4(11.0) 5 1.173 1.5 c tag 91.20 5 1.174(1.323) 10.7( 6.1) 5 1.169 13.2 b tag 91.20 5 1.218(1.209) 34.2(36.6) 4 1.177 10.9ALEPH [34] inclusive 91.20 22 0.987(0.989) 15.6(20.4) 18 1.008 6.1DELPHI [35, 36] inclusive 91.20 17 1 21.0(20.2) 27 1 3.9 uds tag 91.20 17 1 13.3(13.4) 17 1 22.5 b tag 91.20 17 1 41.9(42.9) 17 1 9.1SLD [31] inclusive 91.28 29 1.002(1.004) 27.3(36.3) 29 0.994 14.3 uds tag 91.28 29 1.003(1.004) 51.7(55.6) 29 0.994 42.6 c tag 91.28 29 0.998(1.001) 30.2(40.4) 29 1.000 31.7 b tag 91.28 29 1.005(1.005) 74.6(61.9) 28 0.992 134.1 TOTAL:
459 599.3(671.2) 391 395.0 χ /N dat =1.31(1.46) χ /N dat =1.01 < . × Q/
2. For each bin the measured yields are reduced by these fractions to excludeevents with large ISR or FSR contributions. To convert the dσ h /dz data with the ISR/FSRcut to the total hadronic FF in Eq. (1) one therefore needs to correct the theoretical totalcross section σ tot by multiplying it by the ISR/FSR correction factor, which is estimated to17e 0.64616(3) [38, 39].For the BaBar experiment [40], two data sets were provided, for “prompt” events, whichcontain primary hadrons or decay products of lifetimes shorter than 10 − s, and “conven-tional” events, which include decays of lifetimes (1 − × − s. For pions the conventionalcross sections are ∼ −
15% larger than the prompt cross sections, while for kaons theseare almost indistinguishable. The prompt data are numerically close to the LEP and SLDmeasurements after taking into account Q evolution, although the conventional ones aretechnically closer to most previous measurements which included all decays. In our analysis,we consider both data sets, and assess their impact on the fits phenomenologically.Finally, our theoretical formalism is based on the fixed-order perturbation theory, anddoes not account for resummations of soft-gluon logarithms or effects beyond the collinearfactorization which may be important at small values of z . To avoid inconsistencies betweenthe theoretical formalism and the data, cuts are applied to exclude the small- z region fromthe analysis. In practice, we use a cut z > . Z -bosonmass and z > .
05 for the data at Q ≈ M Z . For kaon data, below z ≈ . z > . Q kaon data sets from ARGUS andBaBar. V. ANALYSIS RESULTS
In this section we present the main results of our IMC analysis. We first establish thestability of the IMC procedure by examining specific convergence criteria, and then illustratethe results for the fragmentation functions through comparisons with data and previous anal-yses. Programs for generating the FFs obtained in this analysis, which we dub “JAM16FF”,can be downloaded from Ref. [47].
A. IMC convergence
We examine two types of convergence tests of the IMC procedure, namely, the iterativeconvergence of the priors (the “grid”), and the convergence of the final posterior distribu-tions. As discussed in Section III B, the convergence of the priors can be tested by observing18he variation of the volume V with the number of iterations, as shown in Fig. 2. For eachiteration 200 fits are performed. During the initial ∼
10 iterations, the volume changes some9 orders of magnitude, indicating a very rapid variation of the prior distribution. After ∼ fits. In Fig. 3 we illustrate the statistical properties of the final posterior distribution byshowing averaged ratios of FFs with smaller samples (100, 200, 500 and 1000) relative to thetotal 10 samples (the averaged error bands are displayed only for the 200 and 10 samples).Using 200 posterior samples, one obtains uncertainty bands that are comparable with thosewith 10 samples. For the central values most of the FFs with 200 samples agree well withthe 10 samples. Some exceptions are the D πs + , D πg , D Kd + and D Kg FFs; however, here thedifferences are in regions where the FFs are poorly determined and the relative error bandsare large. For practical applications these effects will be irrelevant, and using a sample of200 posteriors will be sufficient to give an accurate representation of FFs. Unless otherwisestated the results presented in the following use 200 fits from the final sample.
B. SIA cross sections
In Fig. 4 the normalized yields of the final posteriors versus χ per datum for the trainingand validation sets are presented using the full sample of 10 fits. In the ideal Gaussian limit,the distributions are expected to peak around 2 [11]. In practice, inconsistencies betweendata sets shift the peak of the distribution to larger χ /N dat values. This is evident for thepion production case in Fig. 4, where the χ /N dat distribution peaks around 2.5. In contrast,for kaon production the distribution peaks around 2.1. We stress, however, that even if thepeak occurs at 2, it does not imply consistency among the data sets (or data vs. theory),since the larger experimental uncertainties in the kaon data sets compared with the pioncan induce such behavior.The ratios of experimental SIA cross sections to the fitted values are shown in Figs. 5and 6 for pions and kaons, respectively. For the pion production data, at the lower energies Q (cid:46)
30 GeV there is good overall agreement between the fitted cross sections and the data,19ith the exception of a few sets (TPC, HRS and TOPAZ) that differ by ∼ − ≈
10% normalization, which may be related to the overallnormalization correction from initial state radiation effects [39] or other corrections. Thisshould not, however, affect the z dependence of the extracted FFs.A relatively good description is also obtained of the data at higher energies, Q = M Z ,which generally have smaller uncertainties, although some discrepancies appear at higher z values. In particular, an inconsistency is apparent between the shapes of the DELPHI[35, 36] and SLD [31] spectra at z (cid:38) . uds -tagged data, with theDELPHI data lying systematically above the fitted results and SLD data lying below. Forthe heavy quark tagged results the agreement with DELPHI and SLD data is generallybetter, with only some deviations at the highest z values where the errors are largest. TheOPAL tagged data [32, 33] are the only ones that separate the individual light quark flavors u , d and s from the heavy flavors. The latter have rather large χ values for both pionand kaon data sets, particularly the b –tagged sample of the pion case. While the unfavored d –tagged kaon sample is well described in the fit, the unfavored s –tagged pion data appearless consistent with the theory. In all cases the OPAL tagged data require a normalizationof ≈ z -integrated cross sectionfrom z min to 1.The total χ /N dat for the resulting fit to all pion data sets is ≈ .
31. Using the con-ventional BaBar pion data set instead of the prompt gives a slightly worse overall fit, with χ /N dat = 1 .
46, with the difference coming mostly from the BaBar and TPC inclusive datasets. The Belle data, on the other hand, are better fitted when the conventional BaBar dataset is used. Since the conventional BaBar data lie ∼
10% higher than the prompt, whichthemselves lie slightly below the Belle data, the Belle cross sections require a normalizationshift that is closer to that needed for the conventional BaBar data.For the kaon cross sections, the overall agreement between theory and experiment isslightly better than for pions, mostly because of the relatively larger uncertainties on the K data. At low energies, as was the case for pions, the TPC data [27–29] lie ≈
10% belowthe global fit. Interestingly, though, the Belle kaon data [38, 39] do not require as large anormalization shift as was needed for the Belle pion data in Fig. 5. At energies near the20 -boson pole, Q = M Z , the deviations at large z between the theoretical and experimentalcross sections are not as prominent as for pions, with only the SLD heavy quark tagged data[31] exhibiting any significant disagreement. The OPAL flavor-tagged data [32, 33] generallyprefer an ≈ −
15% normalization for all quark flavors. The DELPHI inclusive and lightquark tagged data [35, 36], which do not include an overall normalization parameter, appearto systematically lie ≈
10% below the fitted results across most of the z range. Fits toother high energy data sets generally give good agreement, and the χ /N dat value for thecombined kaon fit is found to be 1 . C. Fragmentation functions
The fragmentation functions resulting from our IMC analysis are shown in Fig. 7 at theinput scale, which is taken to be Q = 1 GeV for the u , d , s and g flavors and at the massthresholds Q = m q for the heavy c and b quarks. The curve bundles represent randomsamples of 100 posteriors from the full set of fitted results, with the central values andvariance bands computed from Eqs. (21) and (22) using the 200 posteriors selected for thefinal JAM16FF results [47]. Generally the pion FFs have a larger magnitude than the kaonFFs, with the exception of the strange quark, where the s + to kaon distribution D K + s + islarger than that for the pion, D π + s + , over most of the z range. As expected, the u + and d + FFs to π + , which correspond to sums of favored and unfavored distributions and reflect thevalence structure of the pion, are dominant at intermediate and large values of z , z (cid:38) . u and d quarks),these are in fact identical, D π + u + = D π + d + . The s + to pion distribution, in contrast, is smaller inmagnitude, with a peak value at x ∼ . − . ≈ / D π + u + = D π + d + FFs, but one shape forall other pion distributions, Eqs. (19).For the heavy quark FFs to pions, the characteristic differences between the s + , c + and b + distributions generally reflect the different masses of the quarks, with larger masscorresponding to softer distributions. The c + and b + FFs, in particular, are large at low z values, z (cid:46) .
1, and comparable to the light-quark FFs evolved to the same scale. (Note21hat the heavy quark distributions exist only above the mass threshold,
Q > m q .) Thegluon FF D π + g is less singular, but is strongly peaked at z ≈ .
25 at the input scale. Itsuncertainties are also larger than those for the favored distributions, as their effects on theSIA cross sections are of higher order in α s .For the fragmentation to kaons, one of the most conspicuous differences with pions is thelarge magnitude of the strange FF D K + s + at intermediate and high values of z , where it iscomparable to the u + and d + FFs to pions. Reflecting the valence quark structure of K ± ,the D K + u + FF is also similar in size, but because of the mass difference between the strangeand nonstrange light quarks there is no reason for the favored u + and s + fragmentation tokaons to be equal. In fact, we find D K + s + (cid:38) D K + u + at high values of z . Unlike for pions, the d + fragmentation to kaons is unfavored, D K + d + (cid:28) D K + u + , with relatively large uncertainties,peaking at z ∼ . s + fragmentation to π + . The heavyquark FFs to kaons are also sizeable compared with the light quark functions, but peak atslightly larger z values than the corresponding pion FFs. The gluon FF to kaons, D K + g ,peaks at rather high z values, z ≈ .
85, at the input scale, consistent with the findings ofsome earlier analyses [12], and is very small in magnitude.The unusual shapes of some of the FFs, such as the gluon to π + and K + or the unfavoredlight quark FFs, lead to the natural question of whether these are robust distributions orpossibly artifacts of the fitting procedure. We can address this by observing snapshots in theIMC chain, as illustrated in Fig. 8, where the FFs from selected iteration steps are plottedat the input scale as a function of z . The first and last rows in Fig. 8 show the initialand final steps in the IMC procedure, respectively. In addition to the posterior shapes anduncertainties, we display in each row the prior distributions as individual curves. Afterperforming the initial iteration, the large spread in the prior FFs due to the flat samplingof the parameter space is reduced significantly, especially for distributions that are morestrongly constrained by the SIA data. For the FFs that are less directly constrained bythe data, more iterations are needed before convergence is reached, as illustrated by the s + to π + distribution, for example. We find that after ≈
30 iterations all of the distributionsbecome stable, which is consistent with the convergence of the volumes observed in Fig. 2.Although the peaks in some of the FFs, such as D π + ,K + g and D π + s + , are prominent at theinput scale, after Q evolution these become largely washed out. This is illustrated in Fig. 9,where the FFs are evolved to a common scale for all FFs that are above the quark threshold,22amely at Q = 1 ,
10 and 100 GeV and at the Z -boson pole, Q = M Z . Recall that thelowest Q in any of the data sets is ≈
100 GeV , so the shapes at Q = 1 and 10 GeV arenot directly compared with experimental data and are shown for illustration only.Compared with parametrizations from other global FF analyses, our fitted FFs are quali-tatively similar for the most part, but reveal important differences for specific distributions,as Fig. 10 illustrates. For pions, our u + and d + distributions are close to the HKNS [12]and DSS [10] results at large z , but are ∼ −
30% larger in magnitude at low z values, z (cid:46) .
3. The strange quark to pion FF peaks at somewhat larger z than the nonstrange,with a magnitude similar to that in previous fits. The peak in the gluon FF at z ≈ . z tail.The comparison between the various parametrizations for the kaon FFs is quite instruc-tive, especially for the light quark flavors and the gluon. The favored D K + u + and D K + s + FFs inour fit turn out to be of comparable magnitude, with the u + closer to the HKNS results and s + closer to DSS. In particular, for the u quark to kaon FF our result is ≈ −
50% largerthan HKNS, but some 2–3 times greater than DSS over the range 0 . (cid:46) z (cid:46) .
9. On theother hand, the strange to kaon FF lies between the HKNS and DSS results at intermediate z values, but coincides with the DSS at z (cid:38) .
5. Interestingly, we do not observe the largeexcess of s to K fragmentation over u to K found in the DSS analysis, which has importantphenomenological consequences for the extraction of the polarized strange quark PDF fromsemi-inclusive DIS data [7, 8].Recall that in our analysis we use two shapes for the favored D K + u + and D K + s + FFs,Eqs. (20), and one shape for all other kaon distributions. In contrast, previous analyses[10, 12] parametrized the u + and (the unfavored) ¯ u functions separately, assuming that atthe input scale D K + ¯ u = D K + d = D K + ¯ d . In contrast, with the IMC procedure in the presentanalysis we do not impose any relation between the ¯ u and ¯ d FFs, parametrizing only the q + distributions as constrained by data.For the gluon to kaon FF we find a similarly hard distribution as in earlier analyses,peaking at rather large z values, z ∼ . D K + g compared with D π + g can be understood in terms of the higher energy needed for a gluon tosplit to an s ¯ s pair than to a u ¯ u or d ¯ d pair in the pion case [12].Despite the striking shape of the gluon FF at the input scale, it is almost entirely washed23ut after Q evolution to the Z -boson scale, as Fig. 11 illustrates. Here the FFs D ( z ) (ratherthan zD ( z )) are compared for the HKNS [12], DSS [10] and the AKK [16] parametrizations.Viewed on a logarithmic scale, the qualitative features of the shapes of FFs are similar acrossall the parametrizations, especially the HKNS, DSS and the present fit. The AKK resultsgenerally lie above the other parametrizations in the low- z region, while more variation isobserved at higher z values. VI. CONCLUSION
We have performed the first Monte Carlo based QCD analysis for parton to hadronfragmentation functions within collinear factorization, using all existing single-inlusive e + e − annihilation data into pions and kaons. In particular, we include the recent high-precisionSIA data from the Belle [38, 39] and BaBar [40] Collaborations, which significantly extendthe kinematical coverage to large values of z .Our analysis is based on the iterative Monte Carlo approach, first adopted in the recentQCD analysis of polarized PDFs [11], which provides a robust determination of expectationvalues and uncertainties for the FFs. We further extended this methodology by samplingnew priors from previous iterations using a multivariate Gaussian distribution, implementinga new strategy for assessing the convergence of the IMC chain by considering the covari-ance matrix of the posterior distributions. This allowed us to sample fairly the parametereigenspace after each iteration instead of the posteriors, which can exhibit several distinctsolutions. We find that an accurate representation of the FFs can be attained with a sampleof 200 fits.We obtained a relatively good overall description of the pion and kaon SIA data atboth low and high center-of-mass energies, despite some tensions between the high-energyDELPHI and SLD pion data sets in the large- z region. For the kaon data a very good χ /N dat ∼ D π + s + and the D π + g distributions. The latter is morestrongly peaked around the maximum at small z values than either the HKNS or DSSresults, while the former has a somewhat harder z distribution. The kaon FFs, on the other24and, show greater deviation from the earlier results. Here, the favored D K + s + function issimilar in magnitude to that from the DSS parametrization [10] for 0 . (cid:46) z (cid:46)
1, but displaysimportant differences at z (cid:46) . D K + u + FF at moderate to low z values compared with the DSS fit in particular. In contrast, the gluon to kaon distribution,which peaks at very large z values, z ∼ .
85, but with a very small magnitude, is consistentwith the DSS result. The disparity between the fitted D π + g and D K + g functions is particularlystriking. At energies on the order of the Z -boson mass, the evolved distributions are muchmore similar to those of the previous analyses, with the exception of the D π + g and D π + s + FFs.The partial separation of the FFs for the various quark flavors has been possible becauseof the existence of the tagged flavor data and the Q dependence of SIA cross sections, fromlow Q ∼
10 GeV up to the Z -boson mass, selecting differently weighted combinations ofFFs in the γ and Z -exchange cross sections. To further decompose the quark and antiquarkFFs, and better constrain the gluon fragmentation, additional information will be neededfrom SIDIS and meson production in pp collisions. More immediately, it will be particularlyinteresting to examine the effect of the strange to kaon fragmentation on the extraction ofthe polarized strange quark PDF ∆ s from SIDIS data. A combined analysis of polarizedDIS and SIDIS data and SIA cross sections is currently in progress [48]. Acknowledgment
We are grateful to Hrayr Matevosyan for helpful discussions. This work was supportedby the US Department of Energy (DOE) contract No. DE-AC05-06OR23177, under whichJefferson Science Associates, LLC operates Jefferson Lab, and by the DOE contract DE-SC008791. N.S. thanks KEK and J-PARC for their hospitality during a visit where some ofthis work was performed. The work of S.K. and M.H was supported by JSPS KAKENHIGrant Number JP25105010. 25 ppendix A: Hard scattering coefficients
For completeness, in this appendix we give the hard coefficient functions in Mellin momentspace at NLO. For the quark case, the NLO coefficient is [13, 49] (cid:102) H (1) q ( N, Q , µ R = Q , µ ) = 2 C F (cid:20) S ( N ) + S ( N ) + S ( N ) (cid:18) − N ( N + 1) (cid:19) − N + 3( N + 1) −
32 1( N + 1) −
92 + 1 N + (cid:18) N ( N + 1) − S ( N ) + 32 (cid:19) ln Q µ (cid:21) , (A1)while for the gluon one has (cid:102) H (1) g ( N, Q , µ R = Q , µ ) = 4 C F (cid:20) − S ( N ) N + N + 2( N − N ( N + 1) − N − + 4 N − N + 1) + 4( N − N + N + N + 2 N ( N −
1) ln Q µ (cid:21) , (A2)where C F = 4 /
3. Here the harmonic sums S ( N ) and S ( N ) can be written in terms of theEuler-Mascheroni constant γ E , the polygamma function ψ N , and the Riemann zeta function ζ , analytically continued to complex values of N [49], S ( N ) = N (cid:88) j =1 j −→ γ E + ψ (0) N +1 , (A3) S ( N ) = N (cid:88) j =1 j −→ ζ (2) − ψ (1) N +1 , (A4)where the m -th derivative of the polygamma function ψ ( m ) N is given by ψ ( m ) N = d m ψ N dN m = d m +1 ln Γ( N ) dN m +1 . (A5)26 ppendix B: Timelike splitting functions The N -th moments of the splitting functions in the timelike region, up to O ( a s ) correc-tions, can be written for the general case when µ R (cid:54) = µ FF as [43] P ij ( N, µ , µ ) = a s ( µ ) P (0) ij ( N ) + a s ( µ ) (cid:16) P (1) ij ( N ) − β P (0)NS ( N ) ln µ µ (cid:17) , (B1)where β = 11 − n f /
3. At leading order the timelike splitting function moments are givenby the well-known expressions [50, 51] P (0)NS ± = P (0) qq = − C F (cid:20) S ( N ) − − N ( N + 1) (cid:21) , (B2a) P (0) qg = 4 n f C F N + N + 2 N ( N − N + 1) , (B2b) P (0) gq = N + N + 2 N ( N + 1)( N + 2) , (B2c) P (0) gg = − C A (cid:20) S ( N ) − − N ( N − − N + 1)( N + 2) (cid:21) − n f , (B2d)where C A = 3. Note that our notation for the off-diagonal timelike splitting functions P (0) qg and P (0) gq is opposite to that in Ref. [52]. 27t NLO accuracy, the timelike splitting function moments are given by [50, 52, 53] P (1)NS ± = − C F (cid:20) S ( N ) (2 N + 1) N ( N + 1) + 8 (cid:18) S ( N ) − N ( N + 1) (cid:19) (cid:0) S ( N ) − S (cid:48) ± ( N ) (cid:1) +12 S ( N ) + 32 (cid:101) S ± ( N ) − S (cid:48) ± ( N ) − − N + N − N ( N + 1) ∓ N + 2 N + 1) N ( N + 1) (cid:21) − C A C F (cid:20) S ( N ) − (cid:18) S ( N ) − N ( N + 1) (cid:19) (cid:0) S ( N ) − S (cid:48) ± ( N ) (cid:1) − S ( N ) − − (cid:101) S ± ( N ) + 2 S (cid:48) ± ( N ) −
29 (151 N + 236 N + 88 N + 3 N + 18) N ( N + 1) ± N + 2 N + 1) N ( N + 1) (cid:21) − n f C F (cid:20) − S ( N ) + 163 S ( N ) + 23 + 89 (11 N + 5 N − N ( N + 1) ) (cid:21) + ∆ (1)NS , (B3a) P (1) qq = P (1)NS + + n f C F (cid:20) (5 N + 32 N + 49 N + 38 N + 28 N + 8)( N − N ( N + 1) ( N + 2) (cid:21) + ∆ (1) qq , (B3b) P (1) gg = − n f C A (cid:20) − S ( N ) + 163 + 89 (38 N + 76 N + 94 N + 56 N + 12)( N − N ( N + 1) ( N + 2) (cid:21) − n f C F (cid:20) N + 4 N + N − N − N − N − N − N ( N + 1) ( N + 2) (cid:21) − C A (cid:20) S ( N ) + 32 S ( N ) (2 N + 5 N + 8 N + 7 N − N − N − N ( N + 1) ( N + 2) − S (cid:48) ( N ) ( N + N + 1)( N − N ( N + 1)( N + 2) − S ( N ) S (cid:48) ( N ) + 16 (cid:101) S + ( N ) − S (cid:48) ( N ) −
29 (457 N + 2742 N + 6040 N + 6098 N + 1567 N − N − N )( N − N ( N + 1) ( N + 2) −
29 (560 N + 1488 N + 576)( N − N ( N + 1) ( N + 2) (cid:21) + ∆ (1) gg , (B3c)28 (1) qg = 2 n f C F (cid:20)(cid:18) S ( N ) − S ( N ) − π (cid:19) ( N + N + 2)( N − N ( N + 1)+2 S ( N ) (cid:18) N − − N − N − N + 3( N + 1) − N + 1) (cid:19) − N − N + 8( N − N + 2 N + 8 N − N + 1( N + 1) − N + 1) + 92( N + 1) (cid:21) + 2 n f C F C A (cid:20)(cid:18) − S ( N ) + 5 S ( N ) − G (1) ( N ) + π (cid:19) ( N + N + 2)( N − N ( N + 1)+2 S ( N ) (cid:18) − N − + 2( N − N + 2 N − N + 1) + 1 N + 1 (cid:19) − N − + 6( N − + 179( N −
1) + 4( N − N − N − N − N + 5 N − N ( N + 1) − N + 1) − N + 1) − N + 1 − N + 2) + 449( N + 2) (cid:21) , (B3d) P (1) gq = 13 n f (cid:20) S ( N + 1) ( N + N + 2) N ( N + 1)( N + 2) + 1 N − N − N ( N + 1) − N + 1) + 43( N + 1) + 4( N + 2) − N + 2) (cid:21) + 14 C F (cid:20)(cid:16) − S ( N + 1) + 2 S ( N + 1) + 10 S ( N + 1) (cid:17) ( N + N + 2) N ( N + 1)( N + 2)+4 S ( N + 1) (cid:18) − N + 1 N + 1 N ( N + 1) + 2( N + 1) − N + 2) (cid:19) − N + 5 N − N + 4 N ( N + 1) − N ( N + 1) − N ( N + 1)+ 4( N + 1) − N + 1) + 23 N + 1 − N + 2 (cid:21) + 14 C A (cid:20)(cid:18) S ( N + 1) − S ( N + 1) − S ( N + 1) + 2 G (1) ( N + 1) − π (cid:19) × ( N + N + 2) N ( N + 1)( N + 2) − S ( N + 1) (cid:18) − N + 1 N + 1 N ( N + 1) + 4( N + 1) − N + 2) (cid:19) − N −
1) + 4 N + 83 N + 269 N − N ( N + 1) + 223 N ( N + 1) + 16( N + 1) + 683( N + 1) − N + 1) + 8( N + 1) ( N + 2) − N + 2) + 3569( N + 2) (cid:21) , (B3e)29here the terms∆ (1)NS = C F (cid:20) − S ( N ) + 3 + 2 N ( N + 1) (cid:21) (cid:20) S ( N ) − π − N + 1 N ( N + 1) (cid:21) , (B4a)∆ (1) qq = 12 n f C F (cid:20) −
809 1 N − N + 12 N − N + 8( N + 1) + 28( N + 1) − N + 1 + 323 1( N + 2) + 2249 1 N + 2 (cid:21) , (B4b)∆ (1) gg = 12 n f C F (cid:20) −
163 1( N − + 809 1 N − N − N + 12 N + 8( N + 1) − N + 1) + 4 N + 1 −
163 1( N + 2) − N + 2 (cid:21) − n f C A (cid:20) S ( N ) − N − + 1 N − N + 1) + 1( N + 2) − π (cid:21) + C A (cid:20) − S ( N ) S ( N ) + 8 S ( N ) (cid:18) N − − N + 1( N + 1) − N + 2) + π (cid:19) + (cid:18) S ( N ) − π (cid:19) (cid:18) N − − N + 1 N + 1 − N + 2 + 1112 (cid:19) − N − + 223 1( N − − N − N − N − N − N −
143 1 N − N + 1) + 143 1( N + 1) − N + 1) ( N + 2) − N + 1)( N + 2) − N + 2) −
223 1( N + 2) (cid:21) (B4c)are present specifically for the timelike functions [52]. In Eqs. (B3) the sum S (cid:48) m ± ( N ) = 2 m − N (cid:88) j =1 − j j m (B5a)has the analytic continuation S (cid:48) m + ( N ) −→ S m ( N ) , (B5b) S (cid:48) m − ( N ) −→ S m ( N − ) , (B5c)30ith S ( N ) = N (cid:88) j =1 j −→ ζ (3) + ψ (2) N +1 , (B6) (cid:101) S ± ( N ) = − ζ (3) ± (cid:20) S ( N ) N − ζ (2)2 (cid:0) ψ (0)( N +1) / − ψ (0) N/ (cid:1) + Li( N ) (cid:21) , (B7) G (1) ( N ) = ψ (1)( N +1) / − ψ (1) N/ . (B8)The last term in Eq. (B8) involves an integral over the dilogarithm function,Li( N ) ≡ (cid:90) dx x N − Li ( x )1 + x , (B9a)and can be approximated using the expansion [49]Li( N ) ≈ . N + 1 − . N + 2 + 1 . N + 3 − . N + 4 + 0 . N + 5 . (B9b)31
1] S. Albino, Rev. Mod. Phys. , 2489 (2010).[2] A. Metz and A. Vossen, arXiv:1607.02521 [hep-ex].[3] R. D. Field and R. P. Feynman, Nucl. Phys. B136 , 1 (1978).[4] J. C. Collins, D. E. Soper and G. F. Sterman, Adv. Ser. Direct. High Energy Phys. , 1 (1988).[5] H. H. Matevosyan, A. W. Thomas and W. Bentz, Phys. Rev. D , 074003 (2011); ibid. D ,114010 (2011).[6] E. Leader, A. V. Sidorov and D. B. Stamenov, Phys. Rev. D , 114018 (2010).[7] E. Leader, A. V. Sidorov and D. B. Stamenov, Phys. Rev. D , 014002 (2011).[8] E. Leader, A. V. Sidorov and D. B. Stamenov, Phys. Rev. D , 054017 (2015).[9] D. de Florian, R. Sassot, M. Stratmann and W. Vogelsang, Phys. Rev. D , 034030 (2009).[10] D. de Florian, R. Sassot and M. Stratmann, Phys. Rev. D , 114010 (2007).[11] N. Sato, W. Melnitchouk, S. E. Kuhn, J. J. Ethier and A. Accardi, Phys. Rev. D , 074005(2016).[12] M. Hirai, S. Kumano, T.-H. Nagai and K. Sudoh, Phys. Rev. D , 094009 (2007).[13] S. Kretzer, Phys. Rev. D , 054001 (2000).[14] S. Albino, B. A. Kniehl and G. Kramer, Nucl. Phys. B725 , 181 (2005).[15] S. Albino, B. A. Kniehl and G. Kramer, Nucl. Phys.
B734 , 50 (2006).[16] S. Albino, B. A. Kniehl and G. Kramer, Nucl. Phys.
B803 , 42 (2008).[17] D. de Florian, R. Sassot, M. Epele, R. J. Hernandez-Pinto and M. Stratmann, Phys. Rev. D , 014035 (2015).[18] M. Hirai, H. Kawamura, S. Kumano and K. Saito, arXiv:1608.04067 [hep-ph].[19] D. P. Anderle, M. Stratmann and F. Ringer, Phys. Rev. D , 114017 (2015).[20] S. Albino, B. A. Kniehl and G. Kramer, Phys. Rev. Lett. , 192002 (2008).[21] D. P. Anderle, F. Ringer and W. Vogelsang, Phys. Rev. D , 034014 (2013).[22] A. Accardi, D. P. Anderle and F. Ringer, Phys. Rev. D , 034008 (2015).[23] R. Brandelik et al. [TASSO Collaboration], Phys. Lett. B , 444 (1980).[24] M. Althoff et al. [TASSO Collaboration], Z. Phys. C , 5 (1983).[25] W. Braunschweig et al. [TASSO Collaboration], Z. Phys. C , 189 (1989).[26] H. Albrecht et al. [ARGUS Collaboration], Z. Phys. C , 547 (1989).
27] H. Aihara et al. [TPC Collaboration], Phys. Rev. Lett. , 577 (1984).[28] X.-Q. Lu, Ph.D. thesis, The Johns Hopkins University (1986).[29] H. Aihara et al. [TPC Collaboration], Phys. Rev. Lett. , 1263 (1988).[30] M. Derrick et al. [HRS Collaboration], Phys. Rev. D , 2639 (1987).[31] K. Abe et al. [SLD Collaboration], Phys. Rev. D , 072003 (2004).[32] R. Akers et al. [OPAL Collaboration], Z. Phys. C , 181 (1994).[33] G. Abbiendi et al. [OPAL Collaboration], Eur. Phys. J. C , 407 (2000).[34] D. Buskulic et al. [ALEPH Collaboration], Z. Phys. C , 355 (1995).[35] P. Abreu et al. [DELPHI Collaboration], Nucl. Phys. B444 , 3 (1995).[36] P. Abreu et al. [DELPHI Collaboration], Eur. Phys. J. C , 585 (1998)[37] R. Itoh et al. [TOPAZ Collaboration], Phys. Lett. B , 335 (1995).[38] M. Leitgab et al. [Belle Collaboration], Phys. Rev. Lett. , 062002 (2013).[39] M. Leitgab, Ph.D. thesis, University of Illinois at Urbana-Champaign (2013).[40] J. P. Lees et al. [BaBar Collaboration], Phys. Rev. D , 032011 (2013).[41] J. C. Collins, Foundations of perturbative QCD , Cambridge University Press (2011).[42] M. Hirai and S. Kumano, Comput. Phys. Commun. , 1002 (2012).[43] A. Vogt, Comput. Phys. Commun. , 65 (2005).[44] M. Stratmann and W. Vogelsang, Phys. Rev. D , 114007 (2001).[45] J. J. More, B. S. Garbow and K. E. Hillstrom, User Guide for Minpack-1 , ANL-80-74.[46] G. P. Lepage, J. Comput. Phys. , 192 (1978).[47] See .[48] J. J. Ethier et al. , in preparation (2016).[49] M. Gl¨uck, E. Reya and A. Vogt, Z. Phys. C , 471 (1990).[50] E. G. Floratos and C. Kounnas, Nucl. Phys. B192 , 417 (1981).[51] T. Weigl and W. Melnitchouk, Nucl. Phys.
B465 , 267 (1996).[52] M. Gl¨uck, E. Reya and A. Vogt, Phys. Rev. D , 116 (1993).[53] G. Curci, W. Furmanski and R. Petronzio, Nucl. Phys. B175 , 27 (1980). ampler priors fitfitfit posteriorsoriginal datapseudodatatrainingdatafitparameters fromminimization steps validationdatavalidationposterioras initialguessprior FIG. 1: Workflow of the iterative Monte Carlo fitting strategy. In the upper diagram (red lines)an iteration begins at the prior sampler and a given number of fits are performed generating anensemble of posteriors. After the initial iteration, with a flat sampler, the generated posteriors areused to construct a multivariate Gaussian sampler for the next iteration. The lower diagram (withblue lines) summarizes the workflow that transforms a given prior into a final posterior.
10 20 30 40 50 iterations − − − − − L og ( V / V ) πK FIG. 2: Normalized IMC volume versus number of iterations for pions (red lines) and kaons (bluelines). The approximate convergence of the volumes are indicated by the colored regions. . . . . z . . . . . R a t i o t o10 fi t s u + π . . . . z . . . . . d + . . . . z . . . . . s + . . . . z . . . . . R a t i o t o10 fi t s c +
100 fits200 fits500 fits 0 . . . . z . . . . . b + fits10 fits 0 . . . . z . . . . . g . . . . z . . . . . R a t i o t o10 fi t s u + K . . . . z . . . . . d + . . . . z . . . . . s + . . . . z . . . . . R a t i o t o10 fi t s c +
100 fits200 fits500 fits 0 . . . . z . . . . . b + fits10 fits 0 . . . . z . . . . . g FIG. 3: Fragmentation functions computed from 100 (pink), 200 (black), 500 (green), 10 (yellow)and 10 (red for pions, blue for kaons) fits, normalized to the latter. The uncertainties for the 200(black shaded) and 10 results are indicated by the bands. χ /N dat . . . . . . n o r m a li ze d y i e l d π traningvalidationcombined χ /N dat . . . . . . K FIG. 4: Normalized yield of IMC fits versus χ /N dat for the training (blue forward hashed),validation (green backward hashed), and combined (red dotted) samples for π (left panel) and K production (right panel). . . . . z . . . . . d a t a / t h e o r y ARGUS π . . . . z . . . . Belle . . . . z . . . . BaBar . . . . z . . . . TASSO . . . . z . . . d a t a / t h e o r y TPC . . . z . . . . TPC(uds) . . . z . . . . TPC(c) . . . z . . . . TPC(b) . .
12 0 . z . . d a t a / t h e o r y HRS . . . z . . TOPAZ . . . z min . . OPAL . . . z min . . OPAL(u) . . . z min . . d a t a / t h e o r y OPAL(d) . . . z min OPAL(s) . . . z min OPAL(c) . . . z min OPAL(b) . . . z . . d a t a / t h e o r y ALEPH . . . z . . . DELPHI . . . z . . . DELPHI(uds) . . . z . . . DELPHI(b) . . . z . . . d a t a / t h e o r y SLD . . . z . . . SLD(uds) . . . z . . . . SLD(c) . . . z . . . . SLD(b)
FIG. 5: Ratio of experimental single-inclusive e + e − cross sections to the fitted values versus z (or z min for OPAL data [32, 33]) for pion production. The experimental uncertainties are indicated bythe black points, with the fitted uncertainties denoted by the red bands. For the BaBar data [40]the prompt data set is used. . . . z . . . . . d a t a / t h e o r y ARGUS K . . . . z . . . . Belle . . . . z . . . . BaBar . . . . z . . . TASSO . . . . z . . . d a t a / t h e o r y TPC .
15 0 . . z . . . . TOPAZ . . . z min . . OPAL . . . z min . . OPAL(u) . . . z min . . d a t a / t h e o r y OPAL(d) . . . z min . . . . OPAL(s) . . . z min . . OPAL(c) . . . z min . . OPAL(b) . . . z . . . . d a t a / t h e o r y ALEPH . . . z . . . . DELPHI . . . z . . . . DELPHI(uds) . . . z . . . . DELPHI(b) . . . z . . d a t a / t h e o r y SLD . . . z . . SLD(uds) . . . z . . . . SLD(c) . . . z . . . . SLD(b)
FIG. 6: As in Fig. 5, but for kaon production. . . . . . . . . z D ( z ) π + K + u + . . . . . . . . π + K + d + . . . . . . . . π + K + s + . . . . z . . . . z D ( z ) π + K + c + . . . . z . . . . π + K + b + . . . . z . . . . π + K + g FIG. 7: Fragmentation functions for u + , d + , s + , c + , b + and g into π + (red bands) and K + (bluebands) mesons as a function of z at the input scale ( Q = 1 GeV for light quark flavors and gluon, Q = m q for the heavy quarks q = c and b ). A random sample of 100 posteriors (yellow curves for π + , green for K + ) is shown together with the mean and variance (red and blue bands). . . . . z . . . . i n i t i a l u + . . . . z . . . . d + . . . . z . . . . s + . . . . z . . . . c + . . . . z . . . . b + . . . . z . . . . g . . . . z . . . . i n t e r m e d i a t e π + K + . . . . z . . . . . . . . z . . . . . . . . z . . . . . . . . z . . . . . . . . z . . . . π + K + . . . . z . . . . i n t e r m e d i a t e . . . . z . . . . . . . . z . . . . . . . . z . . . . . . . . z . . . . . . . . z . . . . . . . . z . . . . fi n a l . . . . z . . . . . . . . z . . . . . . . . z . . . . . . . . z . . . . . . . . z . . . . FIG. 8: Iterative convergence of the π + (red bands) and K + (blue bands) fragmentation functionsfor the u + , d + , s + , c + , b + and g flavors (in individual columns) at the input scale. The first rowshows the initial flat priors (single yellow curves for π + and green curves for K + ) and their corre-sponding posteriors (error bands). The second and third row are selected intermediate snapshotsof the IMC chain, and the last row shows the priors and posteriors of the final IMC iteration. . . . z . . . . z D ( z ) u + π + K + . . . z . . . . d + π + K + . . . z . . . . s + π + K + Q = 1 GeV Q = 10 GeV Q = 100 GeV Q = M Z . . . z . . . . z D ( z ) c + π + K + . . . z . . . . b + π + K + . . . z . . . . g π + K + FIG. 9: Evolution of the u + , d + , s + , c + , b + and g fragmentation functions to π + (red curves)and K + (blue curves) with the scale, from the input scale Q = 1 GeV (solid) to Q = 10 GeV (dot-dashed), Q = 100 GeV (dashed) and Q = M Z (dotted). . . . . z . . . . z D ( z ) u + π + K + . . . . z . . . . d + π + K + . . . . z . . . . s + JAMHKNSDSS0 . . . . z . . . . z D ( z ) c + π + K + . . . . z . . . . b + π + K + . . . . z . . . . g π + K + FIG. 10: Comparison of the JAM fragmentation functions (solid curves) for π + (red curves) and K + (blue curves) with the HKNS [12] (dashed curves) and DSS [10] (dotted curves) parametrizationsat the input scale Q = 1 GeV for the light quark and gluon distributions, and Q = 10 and20 GeV for the c + and b + flavors, respectively. . . . z − − − − D ( z ) u + Q = m Z π + K + . . . z − − − − d + π + K + . . . z − − − − s + JAMHKNSDSSAKK . . . z − − − − D ( z ) c + π + K + . . . z − − − − b + π + K + . . . z − − − − g π + K + FIG. 11: Comparison of the JAM fragmentation functions (solid curves) for π + (red curves) and K + (blue curves) with the HKNS [12] (dashed curves), DSS [10] (dotted curves) and AKK [16](dot-dashed curves) evolved to a common scale Q = M Z . Note that the fragmentation functions D ( z ) are shown rather than zD ( z ).).