QCD analysis of pion fragmentation functions in the xFitter framework

We present the first open-source analysis of fragmentation functions (FFs) of charged pions (entitled IPM-xFitter) computed at next-to-leading order (NLO) and next-to-next-to-leading order (NNLO) accuracy in perturbative QCD using the xFitter framework. This study incorporates a comprehensive and up-to-date set of pion production data from single-inclusive annihilation (SIA) processes, as well as the most recent measurements of inclusive cross-sections of single pion by the BELLE collaboration. The determination of pion FFs along with their theoretical uncertainties is performed in the Zero-Mass Variable-Flavor Number Scheme (ZM-VFNS). We also present comparisons of our FFs set with recent fits from the literature. The resulting NLO and NNLO pion FFs provide valuable insights for applications in present and future high-energy analysis of pion final state processes.

In the context of high-energy particle physics, Quantum Chromodynamics (QCD) is a remarkably successful theory which governs the interactions of strongly interacting particles in hard scattering processes. Accurate predictions from the QCD theory are essential to provide a detailed description of the hadron structure. This information allows for incisive tests of the Standard Model (SM), and facilitates searches for new physics.
Although the non-perturbative character of QCD imposes formidable calculational challenges, perturbative QCD circumvents this difficulty by factorizing the physics into universal parton distributions functions (PDFs) and fragmentation functions (FFs) which encode the dynamical structure of the hadrons.
The PDFs and FFs must be determined phenomenologically by a global QCD analysis of a diverse collection of hard scattering measurements involving initial-and final-state hadrons. 1 The universal feature of the PDFs and FFs allows us to extract them in a global fit if we have sufficient data to constrain these distributions. Fortunately, we have extensive and precise measurements from a variety of facilities including PETRA, SLC, LEP, BELLE, BaBar, Tevatron, HERA, and the Large Hadron Collider (LHC).
Recent determinations of the pion FFs include results by several groups such as HKNS [15], NNFF [4], JAM [6,13], SGKS [7], SGK [16], and other studies available in the literature including KKP [17,18], BKK [19], and KLC [20] These analyses differ in their choices of input functional form for the pion FFs, the fitting methodology, the method of error estimation, the input experimental data sets, and the order of the QCD perturbative expansion in the theory predictions.
The primary objective of the present study is to perform an extraction of the collinear FFs for the pion (in a vacuum) using an exhaustive set of SIA data along with the very recent pion production measurements by the BELLE collaboration (BELLE20) [21]. To help benchmark the impact of this new data set we will be especially interested to compare the BELLE20 data with the previous 2013 BELLE measurement (BELLE13) [22] as well as with the BaBar [23] data. Additionally, we will compare our analysis at NLO and NNLO accuracy.
xFitter is an open-source package that provides a framework for the determination of the PDFs, and now FFs, for many different kinds of analyses in QCD. The current xFitter version 2.1 offers an expanded set of tools and options, and incorporates experimental data from a wide range of experiments [24,25]. xFitter can analyze this data up to next-to-next-to-leadingorder (NNLO) in perturbation theory with a variety of theoretical calculations. While primarily based on the collinear factorization foundation, xFitter also provides facilities for fits of dipole models and transversemomentum dependent (TMD) PDFs.
xFitter has been used for a wide range of analyses involving different processes; examples include: pion PDFs [26]; nuclear PDFs [27]; and as well as various other studies [24,25,[28][29][30]. This work further extends the utility and capabilities of xFitter to the case of unpolarized collinear FFs of hadrons analysis. The pion FFs analysis presented in this paper represents the first step of a broader program to demonstrate this opensource xFitter tool which can serve as the basis for future analyses of both PDFs and FFs.
The structure of this paper is as follows. Sec. II reviews the theoretical formalism for e + e − annihilation into hadrons. The analysis framework and methodology for the global QCD analysis is presented in Sec. III. Here we also describe our methodology including our parametrization and the treatment of the FFs uncertainties. The SIA experimental data sets analyzed in this study are summarized in Sec. IV. Sec. V presents the details of our results, as well as comparisons other FFs available in the literature. Additionally, we discuss the impact of the new BELLE20 data on the extracted pion FFs and compare with results from the literature. Finally, we summarize our conclusions in Sec.VI and consider future applications that may prove fruitful for the xFitter framework.
To facilitate further investigations of the pion FFs, both the FF grids and the code are publicly available on the xFitter webpage (xfitter.org).
This project is a collaborative effort between the Institute for Research in Fundamental Sciences (IPM) and xFitter; hence, we shall identify our fit results as IPM-xFitter, or just IPMx for short.

II. THEORETICAL FRAMEWORK
The cross section for charged hadron production in an electron-positron annihilation process (e + e − → h ± X) is typically measured as a function of the scaling variable z = 2E h / √ s at a given center of mass energy of √ s = Q with hadron energy E h . The differential cross section is related to the hadronic fragmentation function F h ± 2 as 2 where σ 0 =4πα 2 /Q 2 and α is the electromagnetic coupling.
The factorization theorem allows us to write the nonperturbative hadronic fragmentation function F h ± 2 as a convolution of a perturbative coefficient function C i and a non-perturbative partonic fragmentation function D h ± i given by: [31][32][33] where we sum over parton flavors i = q,q, g. The coefficient functions C i have been calculated up to the NNLO accuracy in the MS scheme [34,35]. The non-perturbative partonic FFs D h ± i are universal and represent the number density for a hadron of type h ± from parton i with momentum fraction z at scale Q. It is the universal property of the FFs which will allow us to extract these quantities by parameterizing their functional form and fitting to experimental data.
To simplify the expansion of the hard scattering crosssection, we choose the renormalization scale µ R and the factorization scale µ F equal to the center-of-mass energy; 3 thus, we have µ R = µ F = √ s ≡ Q. The scale dependence of the partonic FFs is described by the DGLAP evolution equations, [36][37][38][39] where P ij are the perturbative time-like splitting functions, and the convolution integral ⊗ is We solve the integro-differential DGLAP equations directly in z space using the APFEL package [40] which provides NLO and NNLO accuracy. Now that we have outlined the key elements of the SIA cross section calculation, we next examine the framework for our analysis, including the parameterization of the non-perturbative FFs. 2 Here, we will follow the notation of Ref. [4], and the subscript on F h ± 2 suggest an analogy with the F 2 DIS hadronic structure function. 3 To be more precise, in Eq. (2) the α S (µ R ) depends on the renormalization scale µ R , and the partonic fragmentation function D h ± i (z, µ F ) depends on the factorizaton scale µ F .

III. ANALYSIS FRAMEWORK
We will obtain the Fragmentation Functions (FFs) by parameterizing their functional form in z and then performing a fit by minimizing a χ 2 function in comparison with experimental data.
In the following, we detail the analysis framework including the parameterization form, the fitting procedure, and the uncertainty analysis.

A. FFs parametrization and assumptions
We parametrize the z dependence of the FFs at an initial scale Q 0 = 5 GeV which keeps us above m b , and use the DGLAP equations to evolve to arbitrary Q scale. The flexible parametric form we use is: which has (maximally) five free parameters is the Euler beta function. For the charged pion FFs, we fit the flavor combinations i = u + , d + , s + , c + , b + and g. The beta functions in the denominator of Eq. 5 simply ensures 1 0 dz zD i = N i . There are a number of constraints we can impose to reduce the number of free parameters of the fit. From the energy sum rule, we have the relation: where h sums over all possible produced hadrons. For the pion FFs (h = π ± ) this relation provides only an upper bound, but if we expect the lighter pions carry most of the parton momentum, then we have where i = g, q,q. Note, in Table I we report N u + where u + = u +ū, hence the limit on this quantity is N u + < 2. Thus, we will use four shape parameters {α i , β i , γ i , δ i } together with the normalization parameter N i to fit our FFs.
For the π + FFs, we assume isospin symmetry for the favored (u,d) and unfavored (ū, d) components [4,14,15]: We can also use charge conjugation to relate the above π + FFs to the π − FFs: Additionally, for the s + , c + and g flavors we impose the condition γ i =δ i =0. We thus have a total of 19 free parameters to be determined by a standard χ 2 minimization strategy. We now discuss the details of the fit.

B. Fitting procedure
By comparing the theoretical and experimental measurements from our selected data sets, we determine the optimal values of the 19 independent fit parameters by minimizing the χ 2 function as implemented in the xFitter framework: [24,41,42] where the sum on i is over the number of data points. Here t i is the theory prediction, d i is the data measurement, b α are the nuisance parameters, and δ i,stat and δ i,unc are the statistical and uncorrelated systematic uncertainties.
In general, for this analysis we will add the correlated and uncorrelated uncertainties in quadrature; while this allows a generous band of uncertainties on our FFs, we still find it challenging to describe all the data as we detail in Sec. V. In the case where an individual experiment provides multiple data sets (e.g., TASSO, DELPHI, and SLD c.f., Table II), the overall normalization of these data points is correlated; we take this into account as a systematic error. The effect of this correlated systematic normalization will be evident in Fig. 8 which displays both the theory prediction as well as the "theory+shifts" showing the impact of the correlated normalization.

C. Uncertainty analysis
In addition to a central value of the FFs and PDFs, it is important to also evaluate the uncertainty band of these quantities. For the present study, we will apply the Hessian method to estimate the FFs uncertainty band [24,43,44]. The Hessian approach is based on a quadratic expansion of the χ 2 function about its global minimum and can be expressed as: where i, j run over the 19 free parameters of the fit. Here, χ 2 0 is the global minimum, ∆χ 2 is the displacement from the minimum, and y i are the parameter shifts about the minimum.
The Hessian matrix is constructed from the second derivatives of χ 2 and defined as H ij = 1 2 (∂ 2 χ 2 )/(∂y j ∂y j ). The Hessian is a symmetric n×n matrix where n=19 is the number of free parameters, and it will have n orthogonal eigenvectors v i with eigenvalues i . We find it convenient to transform from the y i parameter basis to the v i eigenvector basis as the Hessian is diagonal in this basis. Additionally, we scale the eigenvectors by the square-root of their eigenvalue to introduce a scaled eigenvector z i = √ i v i . Thus, we can express the χ 2 shift as: Finally, we can relate the original fitting parameters y i to the scaled eigenvectors z i via the equation: where V ij is the n×n matrix composed of the eigenvectors v i which diagonalizes the Hessian as in An important consideration when constructing the FFs uncertainty range is to decide upon the acceptable tolerance range T = ∆χ 2 over which we will allow the fit parameters to vary.
For a simple case of a single parameter, the confidence level is 68% for ∆χ 2 =1.
The present case of the pion FFs involves combining large data sets across many different measurements with uncertainties from both experimental and theoretical sources that can only be approximated. Thus, determination of a ∆χ 2 tolerance criteria is not straightforward. For guidance on the choice of tolerance, we reviewed the analyses used in Refs. [15,41,[45][46][47]. For the present analysis we will use ∆χ 2 =20 which was chosen based on a number of factors. We shall present the details in the following sections, but the ∆χ 2 =20 tolerance criteria gives us a reasonable overlap of the bands for our comparison of NLO and NNLO fits, our fits comparing BELLE13 and BELLE20, and also our fits using different scale choices. With a tolerance of ∆χ 2 =1, the bands for the above comparisons would not overlap, and would clearly underestimate the true uncertainties.
These data sets contain inclusive and flavor-tagged measurements of pion scaled energy (z) corresponding to the sum of light quarks (u, d, s) as well as individual charm (c) and bottom (b) quarks.
We list all the data sets used in our fits in Table. II, and indicate the experiment and the observable. The bulk of our data points are measured by the LEP and SLD Collaborations ( √ s=M Z ) and the B-factory experiments ( √ s 10). The remaining measurements range across intermediate scales of energy.
In a manner similar to Ref. [4], we will impose z cuts on the data sets to avoid non-perturbative effects not described by our calculations; such effects could arise from contributions beyond the order of our NNLO calculation, or from non-factorizable (e.g., higher twist) contributions. Specifically, for the experiments with √ s=M Z we take data in the range z ∈ [0.02, 0.9], and for all other sets we take z ∈ [0.075, 0.9].
For all our fits we used α s (M 2 Z )=0.1185 which is compatible with the world average data [47]. We also use M Z =91.1876 GeV, sin 2 θ W =0.231, Q 0 =5 GeV, and pole masses for the heavy quarks of m c =1.43 GeV and m b =4.5 GeV.
We now describe the specific features of the B-factory measurements from BaBar and BELLE which are taken at the lower scale √ s 10 GeV. As this energy is below the B-meson pair production threshold, the effects of bottom quark production are not included for this data. These data sets are important for the fit as they have high precision and extend the momentum fraction z close to one.

A. Updated BELLE Data
The BELLE Collaboration has recently released a new measurement (BELLE20) [21] which we will compare to the previous data set (BELLE13) [22].
The new data set (BELLE20) has removed the backgrounds from Υ decays, τ production and twophoton processes, and also applied an updated initial-state radiative (ISR) correction. Additionally, the updated BELLE20 data contains all systematic uncertainties separated into correlated and uncorrelated contributions.
For the present analysis, these uncertainties are added in quadrature; even with this generous error budget, we will find this data challenging to describe with our theoretical calculations. Furthermore, the number of equidistant bins in the range 0.075<z<0.9 has decreased from 70 in BELLE13 measurements to 32 in BELLE20.
In the case of the BELLE13 experiment, the measurements have been provided in the form dσ H ± /dz. To properly treat this data set consistently with the other measurements included in our analysis, we multiply the BELLE13 measurements [22] by a correction factor as explained in Ref. [4] to compensate for the kinematic cut on radiative photon events which was applied to the BELLE13 data. In contrast, no such correction is needed for the BELLE20 data.

V. RESULTS
We now present the results of our fit for the charged pion FFs. Additionally, we will also compare our resulting FFs with recent results in the literature.

IPM-xFitter Fits
As we studied the influence of the various data sets, we found it especially useful to compare the effects of B-factory data sets from BELLE and BaBar. Thus, we perform five sets of fits which we summarize in the following. Note, we include all the data in Table II except what is indicated below.
Fit A: This fit focuses on the impact of the BELLE13 data set. Thus, we exclude the BELLE20 data for a total of 386 points.
Fit B: This fit focuses on the impact of the BELLE20 data set. Thus, we exclude the BELLE13 data for a total of 348 points.
Fit C: This fit focuses on the impact of the BELLE20 data without the BaBar set. Thus, we exclude the BaBar and BELLE13 data for a total of 308 points.
Fit D: This fit focuses on the impact of cutting the low-z BELLE20 data. Thus, we exclude the BaBar and BELLE13 data, and impose a z>0.2 cut on the BELLE20 data, for a total of 304 points.
Fit E: This fit focuses on the impact of the BELLE20 and BaBar sets with cuts imposed to remove low-z data. Thus, we exclude the BELLE13 data, impose a z>0.2 cut on the BELLE20 data and impose a z>0.1 cut on the BaBar data for a total of 341 points.
Thus, in the above, Fit A will reflect the influence of the updated BELLE13 and Fit B will reflect BELLE20. We can discern the impact of the BaBar data by comparing Fits C and D with Fit E. And finally, we can see the impact of the low z cuts by comparing Fits D and E.

Fit Parameters
In Table I we present the values of the 19 fitted parameters, Eq. (5), and their uncertainties for our fits as computed by Minuit [41].
Examining the relatively small uncertainties on the shape parameters (with the possible exception of α g and α s + ), this suggests that these parameters appear to be relatively well-constrained by the SIA data.
Similarly, Table II displays the the χ 2 values for each individual experiment and the number of degrees of freedom, as well as the totals. The normalizations are treated as a correlated systematic and are also reported [24].

NLO and NNLO Comparison
Since xFitter can perform both NLO and NNLO calculations, we will compare these for the case of Fit A. The quality of our fits can be surmised by examining the total χ 2 and χ 2 per degree of freedom (χ 2 /dof ) as displayed in Table II. For Fit A we observe that the NNLO fit yields a substantially improved χ 2 /dof of 427/386 for NNLO vs. 480/386 for NLO representing an improvement of ∼11%.
If we look in more detail, Table. II also displays the χ 2 values for each individual experiment. While most of the data sets show improvement of the NNLO compared to the NLO results, there are some observables where the individual χ 2 does slightly increase such as the ones measured by the TASSO, TPC and DELPHI experiments; these increases are small compared to the overall improvement of the fit from χ 2 /dof ∼1.24 to 1.11. A similar behavior is observed in other analyses such as Ref. [4].
Finally, if we examine the uncertainty bands for the FFs of Fig. 2, the a posteriori fact the NLO and NNLO bands overlap can be taken as an indication that perturbative uncertainties are under control and this additionally suggests that our choice of tolerance T is reasonable.
While we do not explicitly show details of the other NLO fits, the results are similar. Hence, the inclusion of the higher-order QCD corrections noticeably improves the quality of our fits, and this result is in agreement with other analyses in the literature [4,16,56].
Since the NNLO yields improved results over the NLO calculation, we will focus only on the NNLO fits for the following analyses.

BELLE Data
To examine how well our fits describe the experimental data, we start by comparing Fit A which includes the BELLE13 data and Fit B which includes the recent BELLE20 data. This information is contained both in Table. II and also graphically in Fig.1.
Examining In Fig. 3 we plot the FFs for both Fits A and B. Note that in our framework, the u + and d + distributions will be equivalent; hence, we display u + as well as the combination d + +s + . We also display the uncertainty bands for the FFs computed with our tolerance of T = 20. We see that the u + =d + , charm and bottom distributions show only moderate difference between the two fits, and Fit B yields a slightly softer distribution for u + and d + . In contrast, the strange and gluon show larger difference which reflects, in part, their larger uncertainties. Compared to Fit A, Fit B moderately increases s + , c + , b + while decreasing g.

BELLE and BaBar Comparison
Comparing Fits A and B above, we observed that the inclusion of the updated BELLE20 data caused a substantial increase in the BaBar χ 2 . To explore if this is primarily a tension between the BELLE20 and BaBar data sets, we introduce Fit C which fits BELLE20 but excludes the BaBar data.
Examining the results of Fit C in Table. II we see the quality of the BELLE20 data χ 2 improves from 82/32 (Fit B) to 32/32 (Fit C), and the overall χ 2 improves from 518/348∼1.49 to 404/308∼1.31.
Additionally, in Fig. 4 we display the comparisons of our A, B, and C fits with the data for BELLE13, BELLE20 and BaBar. For BELLE13 (Fit A), the theoretical predictions are entirely consistent with the experimental data, but this is partly due to the larger uncertainties; note the vertical scale of the "theory/data" plot.
Comparing BELLE20 and BaBar, we see that the fits yield a good description of the data with the exception of the low-z region. Here, the BELLE20 and BaBar data sets appear to pull the fit in opposite directions; BELLE20 is above the theory, and BaBar is below the theory.
If we compare Fit B (with BaBar) to Fit C (without BaBar) we notice Fit C shifts toward the data throughout the z range, including the lower z region, thus reducing the χ 2 for BELLE20 from 82/32 (Fit B) to 32/32 (Fit C).

Low-z Cuts
The deviations between theory and data in the lowz region, as observed in Fig. 4, warrants a closer examination. Thus, we introduce Fit D which imposes a z>0.2 cut on the BELLE20 data again does not include the BaBar data. As summarized in Table II, the BELLE20 data set is reduced from 32 to 28 points with a significantly improved χ 2 of 9.2/28, and the total χ 2 for Fit D improves to 357/304 = 1.17.
The encouraging results of imposing the low-z cut leads us to try Fit E which includes both the BELLE20 and BaBar data, but with low-z cuts imposed. Specifically, for BELLE we impose z>0.2 which removes 4 points, and for BaBar we impose z>0.1 which removes 3 points.
(Note, we could impose a more stringent cut of z>0.2 on both data sets, but we shall find this more conservative choice already yields a notable improvement. ) We note that other analyses in the literature have also imposed kinematic cuts on the low z data. For example, the JAM19 focus was on SIDIS in the region of larger z; hence, they imposed a cut of z 0.2 on their data sets. Similarly, NNFF1 used a lower kinematic cut of z min =0.02 for Q=M Z and 0.075 for Q < M Z .
The result of these low-z cuts is dramatic. We obtain an improved 17/28 for BELLE20, 33/37 for BaBar, and a good χ 2 = 410/341∼1.20 for the total data set of Fit E.
To see the impact of these fits in more detail, Fig. 5 displays the comparisons of our fits with the data for BELLE20 and BaBar. Examining the BaBar results, we see that by cutting out the lowest 3 points in z, the theoretical predictions are generally contained within the band of the experimental uncertainty.
Similarly, for BELLE20 we see that the predictions also are generally contained within the band of the experimental uncertainty. We also see that including the BaBar data (Fit E) increases the theory curve in the lowz region, consistent with the effects observed in Fig. 4 To see the impact of the low-z cuts on the FFs, in Fig. 6 we display the results of Fit B (no z cuts) and Fit E (with z cuts). The difference is dramatic and (with the exception of the heavy charm and bottom quarks) the uncertainty bands do not overlap throughout much of the z range.
While our tolerance criteria of T = 20 gave reasonable results for the comparisons of the NLO vs. NNLO results of Fig. 2 and the comparison of BELLE13 and BELLE20 of Fig. 3, Fig. 6 underscores the fact that this does not represent the uncertainty due to the choice of data in the low z region; hence, this must be separately accounted for when computing the maximal FF uncertainty.
Regarding the impact of the z cuts on the shape of FFs, we see that this generally increases the FFs for up, down, charm and bottom while decreasing the strange and gluon FFs (which have comparably larger uncertainties).

Scale Variation
Since the results presented thus far have imposed the scale choice µ F = µ R = Q, we also want to study the effect of varying the perturbative scales, in part, to see if this could have any significant effect on the partonic FF zD π ± (z, Q) in the low-z region. In Fig. 7, we display the resulting FFs for Fit E with their uncertainties for the variation of the scale choice as a function of z. We start with a factorization scale of µ = Q = 10 GeV, and vary this up and down by a factor of 2 to display FFs for µ = Q = {5, 10, 20} GeV. Essentially, this reflects the influence of the DGLAP evolution in Q.
The variation of the uncertainty bands for the different scales is rather modest for all the flavors with the exception the strange FFs which shows a observable shift in the intermediate z region, and the gluon FFs which shows a shift in the low z region. This variation can provide an indication of higher order contributions which may not be included in our fixed-order perturbative calculation. For increasing scales, the FFs generally decrease in the small z region (z∼0.01), and this is most evident in the gluon FF.

Comparison with SIA Data
To provide further details on the quality of the QCD fits, we now compare our theoretical predictions directly with the inclusive pion production data.
In Fig. 8 we display the theoretical predictions for Fit A, Fit B and Fit E at NNLO together with the inclusive pion production data and their uncertainties. The comparisons are given as a function of z for the various energies, and the measured z region for each experiment can be read from the plots. To provide detail, we show both the absolute distributions (upper panel) and the data over theory ratios (lower panel).
The results of TASSO, DELPHI, and SLD are split across multiple data sets (c.f., Table II), but the experiments provide the correlated systematic uncertainty for the overall normalization of the data. For these data sets, the difference between the "theory" and "theory+shifts" is due to this systematic normalization. For the other data sets of Fig. 8, there is no systematic normalizaton, so the "theory" and "theory+shifts" curves will coincide.
The theoretical predictions are generally in good agreement with the experimental measurements in the central z-range. We observe some deviations in the large-z region, but as the data uncertainties are typically larger here this does not significantly impact the χ 2 of the fits. In contrast, in the low-z region where the experimental uncertainties are typically smaller we do see some deviations of the theoretical predictions which, as discussed above, can substantially impact both the FFs (e.g., Fig. 6) and the χ 2 values (e.g., Tab. II). For the BaBar and BELLE20 plots, we note Fits A and B show deviations in the low z region, and this is in contrast to Fit E where these deviations are largely absent.

Theoretical Effects at Low-z
What could be the source of the large χ 2 contributions coming from the low-z region? One concern is that we could be missing important theoretical corrections such as non-factorizable higher-twist contributions, or terms beyond the order or our perturbation theory, some of which can be estimated using a resummation procedure. The original set of z cuts discussed in Sec. IV and outlined in Ref. [4] aim to minimize the non-perturbative contribution not described by our calculations.
Additionally, Ref. [57] performs an all-order resummation of logarithmically enhanced contributions at small momentum fraction z at NNLO logarithmic order, including the dependence on the factorization and renormalization scales. Specifically, they examine terms of the form (α n s ln k z) and resum these contributions in Mellins space to reduce the theoretical uncertainty in the low z region.
This analysis suggests that although the lowest few z data points of BELLE and BaBar may have modest corrections from resummation contributions, these alone are not sufficient to address the differences we observe with the current fits. Additionally, the fact that the BELLE and BaBar data seem to pull the theory in opposite directions further complicates the resolution of this situation and may point to non-perturbative effects not incorporated in our calculations.
In conclusion, while we generally obtain a good quality for our fits across most of the z range, the description of the data in the low-z region remains an unresolved puzzle, and we must leave this for future investigation.

C. Comparison to results from the literature
We now compare our results with recent analyses from the literature. While Fit E (IPM-xFitter) is our preferred final fit, we also display Fit B to highlight the impact of the low z cuts. We compare with the results from the NNFF1 [4,58,59], JAM19 [13] and DSEHS14 [5] collaborations which are computed at NLO (JAM19 and DSEHS14) and at NNLO (NNFF1). These results are displayed in Figs. 9 and 10 at Q 2 = 100 GeV 2 and Q 2 = M 2 Z , respectively. Discretion is necessary when interpreting the very low z region as the extrapolation of the FF grids may extend beyond the region fitted in the individual analyses of Refs. [4,5,13]. For example, the JAM19 focus was on SIDIS in the region z 0.2, NNFF1 used a lower kinematic cut of z min =0.02 for Q=M Z and 0.075 for Q < M Z , and DSEHS14 also used various z cuts. Additionally, differences can arise due to the choice of the tolerance criteria as well as the method for computing the error bands such as the Hessian or Monte Carlo approach.
The DSEHS14 uses a combination of data from SIA, SIDIS, and hadron-hadron collisions in an NLO framework; additionally, our parametric form of Eq. 5 matches their initial FFs. The NNFF1 analysis is based on electron-positron SIA cross-sections for the sum of charged pion, charged kaon, and proton/antiproton production. The JAM analysis simultaneously fits both the PDFs and FFs using DIS, SIDIS, Drell-Yan and SIA data.
Comparing our up and down FFs with NNFF1 and DSEHS14, we see they are generally compatible at larger z, but differ in the low-z region; this is more pronounced for Fit E which imposes cuts on some of the low-z data. Similar conclusions apply to the case of the c + and b + heavy quarks. Comparing our gluon FF with NNFF1 and DSEHS14, again we see our FFs are generally compatible as our curves lie within the NNFF1 uncertainty band, and have overlap with DSEHS14 in the larger z region. For the strange quark there is more of a spread in results suggesting an overall increased uncertainty (note the vertical scale). This reflects, in part, the fact that our chosen data set has minimal constraints on the strange FF.
In contrast, the above FFs generally have a different behavior as compared with the JAM19 analysis. 4 The JAM19 FFs have a much steeper slope at small z for the quark flavors as compared to both our Fit B and Fit E, while the gluon generally lies above our curves for intermediate to larger z values. In addition, much smaller error bands are obtained in JAM19 analysis in respect to our results and the NNFF1 analysis. As noted above, the focus of the JAM19 study was on SIDIS in the region z 0.2, and in this region there is closer agreement among the FFs. The JAM19 study also performed a separately analysis without the SIDIS data (not shown here) and this resulted in a substantive shift of FFs. Specifically, this tends to increase the u + and decrease the g FFs in the larger z region, z 0.2. As the JAM19 results are obtain by a combined fit to both PDFs and FFs, it would be interesting to investigate further to see what difference might be separately attributable to this combined analysis and the choice of data sets.

VI. SUMMARY AND CONCLUSIONS
In this paper we have presented the first open source analysis of unpolarized pion FFs at NNLO accuracy in pQCD. Our analysis is based on a comprehensive data set from SIA experiments, including the most recent measurement of single pion production by the BELLE collaboration [21].
This analysis includes several novel features including the use of the open source xFitter framework to perform the analysis, as well as to compute the Hessian uncertainties for the FFs. We have compared both NLO and NNLO results, and also studied the impact of the B-factory data.
The SIA experimental data are reasonably well described by our QCD fits, with the exception of the low-z region. This remains as an important unresolved issue and will clearly require additional study.
The comparisons of our pion FFs to the results in the literature are also informative. While some of the results agree within uncertainties, there are a number of areas where the FFs differ, and thus warrants further examination. This paper facilitates future investigations as these results are available within the open-source xFitter program, and also as pion FF grids formatted for the LHAPDF6 library. Therefore, this analysis serves as an important step toward a broader program to improve the determination of the nuclear PDFs and FFs of hadrons as we strive to better comprehend the full character of the QCD theory.
Note added: As we completed this study, a related investigation was reported in Ref. [60] using a complementary neural-network approach with different data sets; this opens new avenues for future study.   Table II. For reference, (green) guide lines are shown for a χ 2 /dof of 1 and 3. As the results for up and down are identical, we choose to also display the d + +s + combination. Note that we are plotting the sum π + +π − .

Theory/Data
FIG. 8. The inclusive pion production as a function of z for fixed Q 2 values (displayed) for Fit A, Fit B and Fit E at NNLO. We compare with experimental data from ALEPH [50], OPAL [51], TASSO [53][54][55], TPC [48], DELPHI [52], SLD [49], BaBar [23] BELLE13 [22] and BELLE20 [21]. We show both the absolute distributions (upper panel) and the data over theory ratios (lower panel). The theory shifts arise from the correlated systematic overall normalization of the data points in the case where an individual experiment provides multiple data sets, such as TASSO, DELPHI, and SLD (c.f., Table II). FIG. 9. A comparison of our preferred Fit E [IPMx] as well as Fit B for charged pion FFs (π + +π − ) at NNLO with results from the literature at Q 2 = 100 GeV 2 . We display NNFF1 [4] at NNLO, JAM19 [13] at NLO, DSEHS [5] at NLO, with their uncertainties at Q 2 = 100 GeV 2 . Note, discretion is necessary when interpreting the very low z region as the extrapolation of the FF grids extends beyond the region fitted in the individual analyses. For example, the JAM19 focus was on SIDIS in the region z 0.2, and NNFF1 used a lower kinematic cut of zmin=0.02 for Q=MZ and 0.075 for Q < MZ . While Fit E is our preferred fit, we also display Fit B to highlight the impact of the low z cuts.