Determination of Fragmentation Functions from Charge Asymmetries in Hadron Production

Fragmentation functions (FFs) are fundamental quantities in high-energy physics that describe the probability density in the transition of a parton into a specific hadron, differential in the light-cone momentum fraction of the parton carried by the hadron. The study of FFs has evolved from early parton models, such as the Field-Feynman model [1−3], to the modern framework of QCD collinear factorization [4, 5]. FFs are essential non-perturbative inputs for investigating the internal structure of nucleons and correlations within nuclei. This includes their application in determining polarized parton distribution functions (PDFs) [6, 7] and nuclear PDFs [8−11]. The precise determination of FFs is particularly critical for the upcoming era of high-precision nuclear physics, driven by the development of electron-ion colliders (EICs) [12, 13]. Furthermore, FFs extracted from high-energy data offer a significant opportunity to understand hadronization and confinement.

Global fits that combine perturbative QCD calculations with diverse experimental measurements have made it possible to extract FFs, especially of light charged hadrons. Notable efforts in this area include DSS [14], HKNS [15], AKK [16], NNFF [17], MAPFF [18], JAM [19], and NPC23 [20−22]. These analyses were performed at next-to-leading order (NLO) in QCD, utilizing different data sets and theoretical prescriptions. Recently, significant progress has been made toward determining FFs at next-to-next-to-leading order (NNLO). These efforts have used data solely from single-inclusive electron-positron annihilation (SIA) [23, 24], or combined SIA and semi-inclusive deep-inelastic scattering (SIDIS) data at approximate [25, 26] or full [27] NNLO. However, the relationship between FFs extracted from high-energy data and non-perturbative QCD models remains largely unexplored.

In this work, we propose a direct extraction of the non-singlet (NS) FFs, defined as the difference between quark and anti-quark FFs, from the charge asymmetry measured in pion and kaon production in SIA and SIDIS processes. The availability of world data, combined with state-of-the-art NNLO perturbative QCD calculations, enables a precise and robust determination of the NS FFs for charged pions and kaons. We perform a pion-only fit and demonstrate that a simple three-parameter model for the NS FFs accurately describes all SIDIS data. The scaling index β of the FFs at large momentum fractions is extracted, with a thorough assessment of uncertainties. Our results support the prediction of $ \beta\sim 1 $ from Nambu-Jona-Lasinio (NJL) model [28−31] and disagree with the predictions of $ \beta\sim 2 $ based on either perturbative QCD [32−35] or Dyson-Schwinger equations [31, 36, 37]. In a joint fit of both pion and kaon, we determine a strangeness suppression factor of about 0.5 and observe consistency in the NS FFs for fragmentation into pions and kaons. We note there exist previous determinations on non-singlet FFs of kaons from u and d quarks utilizing SIA data on charged and neutral kaons production with assumptions of isospin symmetry [38−40]. In comparison, our work represents a direct determination of NS FFs with all available world data and less model-dependence.

We consider charged hadron production via fragmentation in neutral current (NC) or (anti-)neutrino charged current (CC) SIDIS on target nucleus H, and SIA at the Z boson mass pole. Based on QCD collinear factorization [5], theoretical predictions on charge asymmetry of hadron production cross sections can be expressed as convolutions of PDFs $ f_{i/H} $, FFs $ D^{h^+}_{i} $ and non-singlet coefficient functions $ \mathcal A_{NS} $,

where i and h label parton flavors and hadron species, respectively. We have assumed charge conjugation symmetry for the fragmentation functions, i.e., $ D^{h^+}_{i}(z)=D^{h^-}_{\bar{i}}(z) $, and have suppressed the dependence on QCD scales in the formulas above. The kinematic variables x, z and θ are Bjorken variable, hadron energy fraction, and hadron polar angle, respectively. The non-singlet fragmentation functions are defined as $ D^{h^+}_{i^-}(z)\equiv D^{h^+}_{i}(z)-D^{h^+}_{\bar{i}}(z) $. In the case of CC SIDIS, $ i'={u,c} $ for $ i={d,s} $ when assuming diagonal CKM matrix. The non-singlet coefficient functions $ \mathcal A_{NS} $ have been calculated to NNLO for SIDIS [41−45] and SIA [46−53].

In this study we have analyzed a comprehensive set of data on identified $ \pi^+ $, $ \pi^- $, $ K^+ $ and $ K^- $ production from SIDIS and SIA processes.The HERMES experiment measured pion and kaon production in NC SIDIS on both proton and deuteron targets at a center-of-mass energy of 7.26 GeV, with z values up to 0.8 [54].The COMPASS experiment measured pion and kaon production in NC SIDIS on isoscalar targets [55, 56] and, more recently, proton [57] targets at a center-of-mass energy of 17.33 GeV, with z values up to 0.85.The ABCMO experiment measured pion production in (anti-)neutrino CC SIDIS on proton targets at a center-of-mass energy of about 8.8 GeV, with z values up to 1 [58].Both HERMES and COMPASS measurements are differential in the momentum transfer Q, and we include only data with $ Q>2 $ GeV to ensure the validity of our perturbative calculations.The ABCMO data on hadron multiplicity are integrated over $ Q>1 $ GeV.The SLD collaboration measured charged pion and kaon multiplicities in the light quark-jet hemisphere at the Z-pole, using longitudinally polarized electron beams [59, 60].We calculated corresponding theoretical predictions at NNLO in QCD [61] using the projected-to-Born method [62].The reported multiplicity data for $ h^+ $ and $ h^- $ are converted to charge asymmetry data by taking their differences, accounting for correlated systematic uncertainties where known.This study marks the first time that COMPASS data from proton targets, SLD charge asymmetry data, and ABCMO neutrino SIDIS data have been included in a global analysis of FFs at NNLO.

We focus on charged pion and kaon production at large-z values ($ z>0.3 $). A simple three-parameter functional form for the non-singlet FFs at the initial scale $ Q_0 $ is well-motivated,

In our nominal fits we set $ Q_0 $=1.3 GeV. In alternative four-parameter form $ a_0 $ is replaced with $ a_0+a_1\sqrt z $. We conduct two independent analyses: one fitting pion FFs and another jointly fitting pion and kaon FFs. We consistently assume isospin symmetry for pion FFs, $ D^{\pi^+}_{u^-}=-D^{\pi^+}_{d^-} $ and set the non-singlet FFs for non-constituent quarks to zero at the initial scale ($ D^{\pi^+}_{(s/c/b)^-}=0 $ and $ D^{K^+}_{(d/c/b)^-}=0 $). In the joint fit, we further assume that all non-singlet FFs share the same α and β parameters, allowing only the normalization to differ due to quark masses. This results in three free parameters for our nominal pion-only fit ($ D^{\pi^+}_{u^-} $) and five for our nominal joint fit ($ D^{\pi^+}_{u^-} $, $ D^{K^+}_{u^-} $ and $ D^{K^-}_{s^-} $). The non-singlet FFs are evolved to higher scales using three-loop time-like splitting kernels [51, 63−65], implemented in a modified version of HOPPET [66], to ensure NNLO consistency. Differential cross sections are calculated at NNLO in QCD using the FMNLO program [67, 68]. Unless otherwise specified, we use the CT18 NNLO PDFs with $ \alpha_S(M_Z)=0.118 $ [69] for calculations involving initial hadrons. The renormalization, factorization, and fragmentation scales are set equal to the momentum transfer Q for both SIA and SIDIS, with scale variations incorporated into the covariance matrix for the $ \chi^2 $ calculations [27].

We first perform a fit to the non-singlet FFs of pions using data from HERMES, COMPASS and ABCMO, as the SLD charge asymmetry data for pions have excessively large uncertainties. The extracted parameters and $ \chi^2 $ values are presented in Tab. I. Our nominal fit yields a global $ \chi^2 $ of 252.4 for 249 data points, indicating an excellent description of the charge asymmetry data with NNLO calculations and a three-parameter non-perturbative FF. Detailed comparisons of the theory and data are also available (see Appendix). The uncertainties of FF are determined using the Hessian method with a tolerance of $ \Delta\chi^2=2.3 $. This tolerance value is estimated based on the agreements of individual data sets in the fit, as detailed in Ref. [27]. The scaling index at large-z, β, is determined to be 0.69 with an uncertainty of about 0.2. We also test predictions from various models, as listed in Table 1. In the first model, the parameter α is fixed to 0, while the other two parameters are fitted. This results in a best-fit $ \chi^2 $ that is 5 units higher than our nominal fit, with β close to 1. The resulting FFs and $ \chi^2 $ are similar if we fix $ \beta=1 $, as suggested by the NJL model [31]. Conversely, fixing $ \beta=2 $, as suggested by [35−37], increases the $ \chi^2 $ by about 74 units. We emphasize that it is possible these predictions on the scaling power are only valid in a region of x very close to one, e.g., $ z > 0.9 $, which are not well probed by current world data. The preference for $ \beta=1 $ over $ \beta=2 $ aligns with observations from global analyses of pion PDFs [70, 71], and could be altered by the inclusion of threshold resummation effects [36, 72, 73]. In the final two models, the non-singlet pion FFs are fully fixed, based on the Field-Feynman model [2] and the CSM FFs [74]. Note that the CSM FF was refitted using our four-parameter form, and an additional constant term was added for the FF from Field-Feynman model. Neither of these models provides a good description of the charge asymmetry data.

Alternative fits under varying conditions are summarized in Tab. I to investigate various systematic effects. First, we use NNPDF4.0 NNLO PDFs [75] instead. This yields a slightly lower $ \chi^2 $, but the extracted parameters are fully consistent with our nominal fit using CT18 PDFs. In the other two fits, we raise the lower cut on z to 0.5 or include only data from proton targets. In both cases, the scaling index β only changes slightly. Furthermore, we attempt a four-parameter fit or lower the initial scale to $ Q_0= $ 1 GeV. The $ \chi^2 $ remains almost unchanged, while β increases or decreases by about 0.13. In the final scenario, we perform a NLO fit, which favors a β value smaller by 0.16. The NLO fit shows a slightly smaller $ \chi^2 $ than our nominal NNLO fit because the scale variations, accounted for in the $ \chi^2 $ calculation, are larger at NLO. Additionally, we perform a dedicated scan of β and plot the change in the global $ \chi^2 $ relative to our nominal best-fit as a function of β in Fig. 1. The horizontal dashed line indicates the global $ \chi^2 $ change corresponding to the 68% C.L. interval. The various markers represent the Field-Feynman model, the best-fit with α fixed to 0, and best-fits with alternative NNLO PDFs, including MSHT20 [76] and NNPDF4.0.

Figure 1.  (color online) The global and individual $ \chi^2 $ variations as a function of the β parameter. The best-fit $ \chi^2 $ values from alternative fits and several models are also shown.

In Fig. 1 we also plot the $ \chi^2 $ changes for individual data sets as a function of β during the scan. All three data sets show consistency, with preferred β values falling within the uncertainty range of the global determination. The COMPASS measurements on proton and isoscalar targets, with a total of 209 data points, provide the strongest constraints, followed by the ABCMO (28 data points) and HERMES (12 data points) measurements.

Our non-singlet pion FFs, $ zD_{u^-}^{\pi^+} $, are shown in Fig. 2 as a function of z at a fragmentation scale of 1.3 GeV. This includes results from the nominal NNLO fit with Hessian uncertainties and from alternative fits under various conditions. Our results are compared to previous determinations, including the Field-Feynman model, CSM FFs, FFs from the BDSSV22 [25], MAP10 [26], and NPC23 [27] global analyses at NNLO, and the JAM20 [19] global analysis at NLO. We also show predictions based on MC simulations with PYTHIA8 [77]. For the simulation, we select electron-positron collisions at a center of mass energy of 11 GeV; the results are insensitive to this choice. Our NLO fit favors a slightly harder distribution compared to the NNLO fit. The CSM results include an uncertainty that reflects possible variation in the value of the hadron scale, and exhibit a very different shape with respect to our nominal results. The FFs from the Field-Feynman model deviate significantly in the region $ z>0.8 $ due to the constant term not vanishing at $ z=1 $. The BDSSV22 and NPC23 FFs are in good agreement with our nominal results, considering their respective uncertainties. The PYTHIA8 predictions show differences in the shape of the FFs, and the MAP10 (JAM20) FFs are significantly below (above) all other predictions.

Figure 2.  (color online) The NS FF $ D_{u^-}^{\pi^+} $ at a scale of 1.3 GeV, as a function of z from our nominal and alternative fits, compared with predictions from the Field-Feynman model, MC generators, and previous FF determination.

In this section, we report the results of the joint fit to non-singlet FFs of pions and kaons using all aforementioned charge asymmetry data. In this joint fit, the SLD kaon data provide direct and dominant constraints on the FFs from strange quarks, specifically $ D_{s^-}^{K^-} $. We have raised the z cut to 0.4 in the nominal joint fit to account for the larger mass of kaons. The best-fit yields a global $ \chi^2 $ of 436.0 for 378 data points with 5 free parameters. The scaling index increases slightly to 0.79 for the joint fit.

We investigate the effects of quark mass on the FFs of light mesons, such as the strangeness suppression factor. In Fig. 3, we plot the ratios of the three NS FFs at 1.3 GeV: $ D_{u^-}^{K^+}/D_{u^-}^{\pi^+} $ and $ D_{s^-}^{K^-}/D_{u^-}^{\pi^+} $. Our best-fit predicts these ratios to be 0.49 and 1.52, respectively. The two ellipses represent the 68% and 95% C.L. region determined with the Hessian method. A tolerance of $ \Delta \chi^2=5.4 $ is estimated for our nominal joint fit due to the slightly worse agreement with the COMPASS proton data. We show predictions from alternative fits with NNPDF4.0 PDFs, a z cut of 0.5, and $ Q_0=1 $ GeV, and a NLO fit. They all fall within the 68% C.L. error ellipse of the NNLO nominal fit. The Field-Feynman model predicts the two ratios to be 0.50 and 1.0, which is close to the boundary of the 68% C.L. error ellipse. We also calculate the ratios using either NNLO and NLO FFs from global analyses or MC simulations following similar procedure as explained in the pion-only fit. In the calculations we have integrated the FFs over z from 0.4 to 1 before taking the ratios. The central prediction from NPC23 (JAM20) FFs is close to the boundary of the 95% (68%) C.L. error ellipse while that from MAP10 shows much larger values of the two ratios. Both PYTHIA8 and JETSET [78] predict a stronger strangeness suppression with $ D_{u^-}^{K^+}/D_{u^-}^{\pi^+} $ of about 0.34.

Figure 3.  (color online) Ratios of the NS FFs of kaons and pions at a scale of 1.3 GeV, from our nominal and alternative fits, compared with predictions from the Field-Feynman model, MC generators, and previous FF determinations.

In Fig. 4, we highlight the impact of the SLD measurements of charge asymmetry in the light-quark tagged hemisphere, specifically $ R^q_{K^-}-R^q_{K^+} $, where R denotes the measured hadron multiplicities. We compare our NNLO predictions from the joint fit, including Hessian uncertainties, to the SLD data, which are shown with error bars. They agree well within the large uncertainties on both sides. Predictions from the JETSET simulation also agree well with the SLD data, due to its generally lower hadron multiplicities. Predictions from the PYTHIA8 simulation, as well as NNLO predictions using NPC23 and MAP10 FFs, are well above the SLD measurements.

Figure 4.  (color online) Comparison of SLD measurements of kaon charge asymmetry in the light-quark tagged hemisphere with theoretical predictions at NNLO based on FFs from this fit and from previous determinations, as well as predictions from MC event generators.

We have presented a new approach to extracting the non-singlet fragmentation functions of charged pions and kaons at NNLO in QCD from charge asymmetries measured in both SIDIS and SIA processes. This is the first time that charge asymmetry data from SIA and neutrino SIDIS have been included in a global analysis alongside theoretical calculations at NNLO accuracy. This methodology has led to a robust and precise determination of the NS FFs of light mesons, yielding several key findings. These include a large-z scaling index close to 1, a strangeness suppression factor of approximately 0.5, and the universality of pion and kaon FFs. The extracted NS FFs were compared to those from previous global analyses and predictions from MC event generators, revealing notable differences. Our results serve as an important benchmark for testing non-perturbative QCD models and as crucial input for future electron-ion colliders.

We thank M. Lochner and L. Bonino for their communications regarding the NNLO calculations of neutrino SIDIS, C. Roberts for providing numerical results of CSM FFs and comments on the draft, T. Sjostrand for communications on the JETSET predictions, and other members of the NPC collaboration for valuable discussions.

This analysis incorporates charged hadron multiplicity data from several key experiments. The COMPASS collaboration provides measurements on multiplicities of charged pions and kaons in kinematic bins of Bjorken-x from 0.004 to 0.4 and inelasticity y from 0.1 to 0.7. This corresponds to a range in the squared momentum transfer, $ Q^2 $, from 0.36 to $ 60\; {\rm{GeV}}^2 $ [55, 56]. The measurements cover a region of the hadron energy fraction z from 0.2 to 0.85 with both isoscalar and proton targets. The HERMES collaboration has published charged pions and kaons multiplicities covering x from 0.023 to 0.6, z from 0.2 to 0.8, hadron transverse momentum $ P_{h,\perp} $ from $ 0 $ to $ 1.2\; {\rm{GeV}} $, and $ Q^2 $ from $ 1 $ to $ 15\; {\rm{GeV}}^2 $ [54]. Additionally, we include charged pion production data from the ABCMO collaboration, which measured (anti-)neutrino charged-current semi-inclusive deep inelastic scattering on proton target [58]. These measurements were performed at a center-of-mass energy of approximately $ 8.8\; {\rm{GeV}} $, with kinematic cuts of $ x>0.1 $ and the hadronic invariant mass $ W>3 $ GeV, and z values extending up to $ 1 $.

In Figs. A1-A4 we present a comparison between the predictions of our nominal fit and the experimental data for each of the data sets used in our analysis. The experimental value and nominal fit prediction are shown in green and red lines respectively, while the theoretical uncertainty is shown by the shaded bands. Scale variations are estimated by varying the renormalization ($ \mu_R $), factorization ($ \mu_F $), and fragmentation ($ \mu_D $) scales simultaneously $ \mu_R/\mu_{R,0}=\mu_F/\mu_{F,0}=\mu_D/\mu_{D,0}=\{1/2, 1, 2\} $ and taking the envelope. Hessian uncertainties are determined using the Hessian method. Experimental uncertainties shown are quadratic sum of statistical and uncorrelated systematic uncertainties. We also display the predictions from an alternative fit fixing $ \beta = 2 $ as a reference, which results in worse description of the global data in general.

Figure A1.  (color online) Comparisons of experimental data and predictions from nominal fit and fit fixing $ \beta=2 $ for COMPASS measurements of pion multiplicity difference with iso-scalar target. Scale uncertainties are shown in shaded bands and estimated by varying $ \mu_R/{\mu_{R,0}}=\mu_F/{\mu_{F,0}}=\mu_D/\mu_{D,0}= $$ \{1/2,1,2\} $ and taking the envelope. Hessian uncertainties are evaluated by the Hessian method. Experimental uncertainties shown are quadratic sum of statistical and uncorrelated systematic uncertainties. Each panel represents different kinematic regions in x and y.

Figure A4.  (color online) Same as Fig. A1 but for ABCMO measurements with ν and $ \bar\nu $ beams.

Figs. A1 and A2 show the comparison to COMPASS data from both isoscalar and proton targets across various kinematic bins. Fig. A3 presents the comparison with HERMES data from deuteron and proton targets in different $ Q^2 $ ranges. Fig. A4 displays the comparison to ABCMO (anti-)neutrino-proton scattering data. For the $ \bar{\nu}p $ scattering, the charged pion asymmetry is negative; we have plotted its absolute value for clarity.

Figure A2.  (color online) Same as Fig. A1 but for COMPASS measurements with proton target.

Figure A3.  (color online) Same as Fig. A1 but for HERMES measurements with deuteron and proton targets respectively.

Tabs. A1 and A2 present the statistical measures of the fit quality. Tab. A1 corresponds to the pion-only fit, while Tab. A2 shows results for the joint pion and kaon fit. In the data sets of joint fit, the SLD 2004 measurement includes charged pion and kaon multiplicities in the light quark-jet hemisphere at the Z-pole, using longitudinally polarized electron beams [59], while the 1999 measurement only includes charged kaon multiplicities [60]. We note that both of the SLD measurements presented unfolded results on hadron multiplicities corresponding to a pure light quark, according to the calculated purity of the light-quark jet. The reported hadron multiplicities have been converted back to the level of a light-quark jet, in order to compare with our theoretical predictions at higher orders [61]. For each data set, we report the number of data points $ N_{pt} $, $ \chi^2 $, $ \chi^2/N_{pt} $, and the effective Gaussian variable $ S_E $. Note in estimation of the Hessian uncertainties, we use a $ \Delta\chi^2 $ tolerance being square of the maximal effective Gaussian variables of all individual data sets [27]. For the pion-only fit, we also compare to the alternative fit fixing $ \beta = 2 $.