Automated calculation of Jet fragmentation at NLO in QCD

Cross sections for any infrared and collinear (IRC) safe observables in a standard subtraction scheme at NLO have the following schematic form:

where $|M|_{B}^{2}$, $|M|_{V}^{2}$ and $|M|_{R}^{2}$ represent square of matrix elements at leading order (LO), one-loop level and in real corrections, respectively. $|{\cal I}|_{m+1}^{2}$ denotes the local subtraction terms constructed in $D=4-2\epsilon$ dimensions when using dimensional regularizations, and

are the integrated subtraction terms over the phase space of real radiations. We have taken a single differential cross section in observable $F$ as an example for a process with $m(m+1)$ finale state particles at LO (in real corrections). One can imagine $F$ being the transverse momentum of either colorless particles or a clustered jet produced. The measure function $\hat{F}$ for the observable applies on either the Born kinematics and flavors $\{p_{m};f_{m}\}$, or those in real corrections $\{p_{m+1};f_{m+1}\}$, and in real subtractions $\{\tilde{p}_{m};\tilde{f}_{m}\}$. For local subtraction schemes, for instance, in CS dipole subtraction Catani:1996jh ; Catani:2002hc or FKS subtraction Frixione:1995ms ; Frixione:1997np , contributions from the second line of Eq. (1) can be evaluated immediately in four dimensions due to cancellations of both infrared and collinear singularities. Note that in the infrared or collinear limits the measure function of IRC safe observable equals for configurations $\{p_{m+1};f_{m+1}\}$ and $\{\tilde{p}_{m};\tilde{f}_{m}\}$. Both the virtual corrections and integrated subtraction terms carry poles in $\epsilon$ which again cancel among each other which renders the first line of Eq. (1) being finite in four dimensions.

For fragmentation processes, the complications are due to uncancelled collinear singularities from splitting of final state partons and thus observables are not collinear safe. However, those singularities are universal and can be absorbed into definitions of bare fragmentation functions similar to the mass factorization in scattering with initial hadrons. Considering the transverse momentum distribution of a tagged hadron, a typical observable in parton fragmentations, one can attempt to use the same formalism as for IRC safe observable and calculate

The bare fragmentation functions for finding a tagged hadron ($h$) with a momentum fraction of $x$ of the mother parton ($i$) can be expressed in terms of the physical ones using mass factorization with $\overline{\rm MS}$ scheme as

where the time-like convolution kernel can be expressed as

where $\mu_{R}$ and $\mu_{D}$ are the renormalization and fragmentation scale respectively, and $\gamma_{E}$ is the Euler-Mascheroni constant. The one-loop regularized splitting functions $P^{+(0)}_{ji}$ are given in Appendix C and coincide with those of space-like splitting functions while the two-loop results are available from Ref. Stratmann:1996hn .

However, one can not take the four dimensional limit in each of the two phase space integrals and carry out numerical calculations due to the aforementioned collinear singularities. Additional subtraction terms and their integrals are needed for the NLO calculation, for instance, as given in Ref. Catani:1996jh ; Catani:2002hc . In this study we propose a minimal modification of Eq. (3) by using additional slicing of radiation phase space to single out the collinear singularities instead of including further local subtractions. We denote that as a hybrid scheme since it involves both methods of local subtractions and slicing of phase space. The master formula for the same distribution is given by

where we have inserted two $\Theta$ functions to partition into the unresolved and resolved collinear regions with a cutoff $\lambda$. The slicing variable $C$ can be chosen as either the minimum of the usual angular separations of all QCD partons $\Delta\theta=min\{\Delta\theta_{ij}\}$ in the center of mass frame or the minimum of the boost-invariant angular separation of all QCD partons $\Delta R=min\{\Delta R_{ij}\equiv\sqrt{\Delta\phi_{ij}^{2}+\Delta y_{ij}^{2}}\}$. The phase space integral of $m+1$-body above the cutoff is free of infrared and collinear singularities and can be calculated numerically in four dimensions. The integral below the cutoff can be factorized using the collinear approximations. The results contain collinear singularities and are included in terms $|\tilde{\cal J}|_{m}^{2}$ which we calculate in the following.

For NLO calculations, there only exist single unresolved collinear regions. Precisely, giving the cutoff is small enough, there is no overlap of phase space between any collinear regions of two partons $\{kl\}$. We can write the integral below the cutoff as following when neglecting power corrections that vanish when taking the limit of $\lambda$ to zero,

We have parameterized the phase space of radiations with the invariant mass square $s$ and the momentum fraction $z$ in each of the collinear regions. We further use the fact that both the square of real matrix elements and its subtractions can be written as unregularized splitting functions, $P^{(0)}_{ji}(z,\epsilon)\equiv P^{(0)}_{ji}(z)+\epsilon P^{\prime(0)}_{ji}(z)$, times the square of Born matrix elements. It is understood that the subscripts except the labels ($q,\,\bar{q},\,g$) represent flavor of that parton ($q,\,\bar{q},\,g$) when appearing in the splitting functions. We arrive at a rather compact form for $|\tilde{\cal J}|_{m}^{2}$ after carrying out integrals in $s$ and $z$, which is given by

with $\tilde{D}_{h/i}(x,\mu_{D})\equiv\sum_{j}I_{ji}\otimes D_{h/j}(x,\mu_{D})$, and the kernel of residuals $I_{ji}(z)$ can be expressed using the unregularized splitting functions,

Substituting $|\tilde{\cal J}|_{m}^{2}$ back to the master formula Eq. (6) one can find that the remaining collinear divergences or poles in $\epsilon$ cancel with those from mass factorizations. After that all remaining pieces are ready for numerical calculations performed directly in four dimensions as will be explained in the next section.

In above derivations, we have chosen the slicing variable $C=\Delta R$ as an example. It can be easily transformed into the case of $C=\Delta\theta$ by exchanging $\lambda p_{T,i}$ with $\lambda E_{i}$ in Eq. (7). We emphasize that the hybrid scheme proposed above can apply equally to processes without initial state hadrons, e.g., lepton collisions or particle decays, as well as lepton-hadron or hadron-hadron collisions. In latter cases, that implies the usual NLO subtraction terms related to initial hadrons and mass factorizations of parton distributions are included implicitly in the derivations. It is interesting to compare our scheme with the two-cutoff method for NLO calculations of fragmentations introduced in Ref. Harris:2001sx . We note that our kernels of residuals $I_{gq,qg}(z)$ coincide with similar quantities therein since the corresponding splittings are free of soft divergences. For $I_{qq,gg}(z)$, there is no simple correspondence of the two methods since the soft divergences are handled differently.

The advantage of above scheme is the capability of easy implementations into various existing programs for NLO calculations designed for IRC safe observables. We demonstrate that for calculations of differential distribution in the energy fraction $x_{h}\equiv 2E_{h}/Q$ carried by the tagged hadron. Schematically our master formula can be recast as

with

and $d\sigma^{(0)}_{m}/dx_{i}$ is the LO partonic differential cross section with respect to the energy fraction carried by the $i$-th parton. The partonic cross sections $d\tilde{\sigma}^{(1)}_{m}/dx_{i}$ and $d\tilde{\sigma}^{(1)}_{m+1}/dx_{i}$ can be identified by comparing to various contributions in Eq. (6). In summary the differential cross sections at NLO and at hadron level can be expressed as a convolution of the original fragmentation functions ($D$) and its integrals ($\bar{D}$ and $\tilde{D}$) with various partonic cross sections.

We have constructed a fast interface specialized for our calculations as following. First of all the fragmentation function and its integrals at arbitrary scales can be approximated by an interpolation on a two-dimensional grid of $x$ and $Q$,

where we choose the interpolation variables $y(x)=x^{0.3}$, $w(Q)=\ln(\ln(Q/0.3\,{\rm GeV}))$ and the interpolation order $n=4$. $D^{k,l}$ is the value on the $k$-th node in $x$ and $l$-th node in $Q$. The spacing $\delta y$($\delta w$) has been chosen so as to give $N_{x}=50$ ($N_{Q}=16$) grid points evenly distributed for the typical kinematic regions considered. We use an $n$-th order polynomial interpolating function $I_{i(j)}^{(n)}$ and the starting grid point $k(l)$ is determined such that $x(Q)$ is located in between the $k(l)+1$-th and $k(l)+2$-th grid points. Substituting the interpolated functions to Eq. (2.2) we arrive at

In practice we calculate partonic cross sections with MG5_aMC@NLO Alwall:2014hca ; Frederix:2018nkq and extract matrix of coefficients $G$, $\bar{G}$, and $\tilde{G}$ for a series of $x_{h}$ values. The matrices need to be calculated once and stored using histograms in MG5_aMC@NLO. For arbitrary choices of fragmentation functions we use HOPPET Salam:2008qg ; Salam:2008sz to carry out DGLAP evolution and convolution of fragmentation functions. Thus the final hadronic cross sections can be obtained via matrix multiplications efficiently without repeating the calculations of NLO partonic cross sections which are time consuming. Experimental measurements provide bin-averaged cross sections rather than differential cross sections at a single value of $x_{h}$. They can be constructed again using interpolations from differential cross sections on a dense grid of $x_{h}$ which we choose to be the same as the one used for $x$ interpolation of fragmentation functions. We have verified that the prescribed interpolations on both fragmentation functions and hadronic cross sections give a precision better than a few per mille in general. We emphasize that the above fast interface can also work for any hadronic differential cross sections related to longitudinal momentum in fragmentations with minimal modifications. A driver for running MG5_aMC@NLO to generate the NLO coefficient tables for a variety of distributions and associated fast interface have been made available as explained in Appendix A.

We consider two scenarios of lepton collisions for benchmark purpose. We focus on the NLO predictions for distribution of the energy fraction carried by the tagged hadron, namely $d\sigma/dx_{h}$. In the first case, the LO hard process involves the annihilation of an electron-positron pair into quarks through virtual photons. These quarks subsequently undergo fragmentation to produce the tagged hadrons. In the second case, the LO hard process involves the annihilation of a muon-anti-muon pair into two gluons via the coupling with the SM Higgs boson. These gluons then undergo fragmentation to produce the tagged hadrons. For above processes the NLO predictions can be calculated analytically with results collected in Appendix C. For simplicity we use a toy model of the fragmentation functions in the calculations

with $N_{i}=1$ and $9/4$ for (anti-)quarks and gluons respectively. We choose a center of mass energy $Q=200$ GeV and set the renormalization and fragmentation scales to $\mu_{R}=\mu_{D}=Q$.

We show comparisons of our numerical results, denoted as FMNLO, and those using NLO analytical formulas in Fig. 1 for the di-quark and in Fig. 2 for di-gluon production respectively. We include two groups of results from FMNLO using a cutoff parameter of $\lambda=0.01$ and $0.04$ to check the consistency of our hybrid scheme. The upper panel shows the NLO predictions on distributions normalized to the LO total cross sections. The middle and lower panel show the ratios of the three NLO predictions to the analytical results at NLO and the ratios of the NLO results to the LO ones respectively. We find very good agreement between our predictions and the analytical results for both channel. For instance, the NLO predictions with $\lambda=0.04$ differ with the analytical ones by at most two per mille in the range of $x_{h}$ from $0.01$ to $1$. We have checked that these differences are indeed due to the interpolations used and can be reduced if a denser $x$-grid is used. The differences between FMNLO predictions with $\lambda=0.01$ and $0.04$ for the virtual photon case are mostly due to fluctuations of Monte Carlo (MC) calculations. We further compare predictions with a variety of $\lambda$ choice ranging from 0.001 to 0.08 and conclude that the choice of $\lambda=0.04$ is sufficiently small to ensure convergence and stability of MC integration. It is worth noting that the numeric effects of the size of a few per mille are much smaller than the typical experimental uncertainties or scale variations of NLO predictions.

We compare our calculations with the INCNLO program INCNLO for the case of unidentified charged hadron production via QCD at hadron collisions. We consider the scenario of $pp$ collisions at LHC $7$ TeV and predictions for the transverse momentum distribution of the charged hadrons $d\sigma/dp_{T,h}$. The charged hadrons are required to have rapidity $|y|<2.0$. On various theoretical inputs we use CTEQ6M NLO parton distributions Pumplin:2002vw and the BKK NLO fragmentation functions bkk1 ; bkk2 ; bkk3 . Furthermore, for simplicity, we fix both the renormalization and factorization scales to $100$ GeV, and set the fragmentation scale to $p_{T,h}$, namely transverse momentum of the charged hadron.

In Fig. 3 we demonstrate independence of our NLO predictions on the choice of the cutoff parameter. We show predictions of the double differential cross sections for $p_{T,h}$ in three kinematic bins from $60$ to $220$ GeV, for several choices of $\lambda$ from $0.002$ to $0.08$. The variations are within $1\%$ in general mostly due to uncertainties of MC integration. In practice, we recommend using a value of $0.02\sim 0.04$ for numerical stability. In Fig. 4 we present comparisons of our NLO predictions with those from INCNLO1.4 INCNLO for a finer binning. The agreements of the three predictions with $\lambda=0.02$, $0.04$, and $0.08$ are similar to that in Fig. 2. From the middle panel we can find our predictions with $\lambda=0.04$ agree with INCNLO predictions at a few per mille in general. In the lower panel ratios of the NLO to LO predictions for various conditions is presented. The NLO corrections can reach to 70% in lower $p_{T,h}$ regions which are much larger than the discrepancies mentioned.

In Ref. 1801.04895 the CMS collaboration measured parton fragmentation based on hard scattering events in $pp$ collisions ($\sqrt{s}=5.02$ TeV) consisting of an isolated photon in association with jets. The photon is required to have a transverse momentum $p_{T,\gamma}>$ 60 GeV and a pseudo-rapidity $|\eta_{\gamma}|<$ 1.44. Jets are clustered with anti-$k_{T}$ algorithm Cacciari:2008gp with $R=0.3$ and are required to have $p_{T,j}>$ 30 GeV and $|\eta_{j}|<$ 1.6. They select jets that have an azimuthal separation to the photon $\Delta\phi_{j\gamma}>7\pi/8$ and analyze the charged-particle tracks inside the jet with transverse momentum $\vec{p}_{T,h}$ in Ref. 1801.04895 . The charged tracks are required to have $p_{T,h}>1$ GeV. The transverse momentum of the photon $\vec{p}_{T,\gamma}$ serves as a good reference of the initial transverse momentum of the fragmented parton. Thus $\xi_{T}^{\gamma}\equiv\ln[-p_{T,\gamma}^{2}/(\vec{p}_{T,\gamma}\cdot\vec{p}_{T,h})]$ is a good probe of the momentum fraction carried by the charged hadron. The results are presented in a form of $1/N_{j}dN_{trk}/d\xi_{T}^{\gamma}$ which is simply a linear combination of the quark and gluon fragmentation functions at the LO evaluated at a momentum fraction of $e^{-\xi_{T}^{\gamma}}$. Fig. 5 comprises three panels that present a comparison between NLO predictions obtained from different fragmentation function sets and experimental data. The first panel displays the results derived from the BKK and NNFF1.1 sets. Upon examination, it becomes evident that the results from the BKK set closely resemble the experimental data and exhibit a good agreement in the lower $\xi_{T}^{\gamma}$ region, ranging from $0.5$ to $2.5$. As $\xi_{T}^{\gamma}$ increases, the discrepancy enlarges, which can be attributed to the lack of fitted data in the small $x=p_{T,h}/p_{T,\gamma}$ ($<0.01$) regions bkk1 ; bkk2 ; bkk3 . This observation is further supported by the second panel, where the results are normalized to the experimental data. The NNFF1.1 results match the experimental data in the lower and higher regions, while a significant deviation is observed in the middle region. Moreover, the error band indicates that the theoretical uncertainties increase with higher values of $\xi_{T}^{\gamma}$. Finally, in the last panel, we give the LO and NLO predictions based on NNFF1.1, normalized to NNFF1.1 results at NLO with nominal scale choice. The ratios indicate that the NLO corrections contribute insignificantly, less than 20%, in the region $\xi_{T}^{\gamma}<3.5$. However, for $\xi_{T}^{\gamma}>3.5$, the contribution switches from negative to positive, and the ratio rises rapidly to nearly 2.

We now turn to the relevant calculations for the production process of Z boson in association with jets at LHC. In Ref. 2103.04377 the CMS collaboration measured parton fragmentation based on hard scattering events of the above process, where $\sqrt{s}=5.02$ TeV. The $Z$ boson is required to have a transverse momentum $p_{T,Z}>$ 30 GeV, and no jet reconstructions are performed. They also analyzed all charged-particle tracks with an azimuthal separation to the $Z$ boson $\Delta\phi_{trk,Z}>7\pi/8$ in Ref. 2103.04377 . The charged tracks are required to have $p_{T,h}>1$ GeV and $|\eta_{h}|<$ 2.4. Different from the production process of an isolated photon in association with jets, we use a distribution of $1/N_{Z}dN_{trk}/dp_{T,h}$ to show our results.

In Fig. 6, the BKK and NNFF1.1 results are depicted as previously mentioned. It is apparent from the first panel that the BKK data exhibits better agreement with the experimental data in the whole kinematic region. In the second panel, we find that, in most regions, the experimental data lies within the error band of the BKK results, with a maximum deviation of approximately 20%. Meanwhile, the NNFF1.1 results show a greater discrepancy, particularly in the middle region. In the third panel, it can be seen that, in most regions, the NLO corrections are negative, but they diminish as $p_{T}$ increases. The maximum corrections at NLO is approximately 20%.

In this subsection, we present the third example of the calculations mentioned above. In Ref. 1906.09254 the ATLAS collaboration measured parton fragmentation based on hard scattering events in $pp$ collisions ($\sqrt{s}=13$ TeV) consisting of two or more jets. Jets are clustered with anti-$k_{T}$ algorithm with $R=0.4$ and are required to have $p_{T,j}>$ 60 GeV and $|\eta_{j}|<$ 2.1. The two leading jets are required to satisfy a balance condition $p_{T,j1}/p_{T,j2}<$ 1.5, where $p_{T,j1(2)}$ are the transverse momentum of the (sub-)leading jet. They also analyzed charged-particle tracks inside the jet classified according to its transverse momentum and pseudo-rapidity (forward or central) in Ref. 1906.09254 . The charged tracks are required to have $p_{T,h}>0.5$ GeV and $|\eta_{h}|<$ 2.5. The results are presented in a differential cross section of $1/N_{j}dN_{trk}/d\zeta$ with $\zeta\equiv p_{T,h}/p_{T,j}$ and $p_{T,j}$ being the transverse momentum of the jet probed¹¹1We note that the distributions presented in the experimental publication have been multiplied by the bin width of each data points..

We present our NLO predictions and compare them to the ATLAS measurement using the central jet of the two leading jets and with $p_{T,j}\in[200,\,300]$ GeV in Fig. 7. The data are displayed as mentioned before. From the first two panels, we find both the NNFF1.1 and BKK results fit well in the high $\zeta$ region. However, the BKK data aligns more closely with the experimental data. In the lower $\zeta$ region, it can be seen that the first three bins of the NNFF1.1 data exhibit a closer resemblance. And the error band of the BKK results in these regions is considerable. In the third panel, it is apparent that the NLO correction is more significant compared to the previous two experiments, and the ratio can reach nearly 2. The NLO correction transitions from negative in the lower region of $\zeta$ to positive in the higher region.

In this study, we analyzed several recent publications on fragmentation function measurements at the LHC over the past five years. Relevant information including the kinematic coverage are summarized in Table. 1. We focus our analysis solely on data obtained from $pp$ collisions. These studies were conducted at a center of mass energy of 5.02 TeV, with the exception of the ATLAS inclusive dijet analysis which used a higher energy of 13 TeV. The measurements can be separated into three categories including using an isolated photon or a $Z$ boson recoiling against the fragmented parton, or using the clustered jet as a reference of the fragmented parton.

In the case of the isolated-photon-tagged jets, the CMS 2018 analysis 1801.04895 measured the normalized distribution $1/N_{j}dN_{trk}/d\xi_{T}^{\gamma}$ that has been explained in previous sections. They probe a region of momentum fractions of the parton carried by hadrons from $0.01\sim 0.6$ based on the definition of $\xi_{T}^{\gamma}$. The ATLAS 2019 analysis 1902.10007 has different setups as the CMS analysis. We highlighted a few of them as below. Firstly the photons and the jet are required to have a transverse momentum in [80, 126] GeV and [63, 144] GeV respectively. The pseudo-rapidity of the photons and the jets have been extended to 2.37 and 2.1 compared to the CMS measurement. In addition they measured the normalized distribution $1/N_{j}dN_{trk}/dp_{T,h}$ in the region $p_{T,h}\in[1,100]$ GeV, that is used in our fit.

For measurements involving $Z$-tagged jets, the CMS 2021 analysis 2103.04377 measured the normalized distribution $1/N_{Z}dN_{trk}/dp_{T,h}$ with setups detailed in previous sections. In the ATLAS 2020 analysis 2008.09811 , the same distribution was measured for three transverse momentum regions of the $Z$ boson, namely [15, 30] GeV, [30, 60] GeV, and beyond 60 GeV. Besides, the requirement on azimuthal separation between the charged track and the $Z$ boson is $\Delta\phi_{trk,Z}>3\pi/4$ instead. The covered $p_{T,h}$ region is [1, 30] GeV and [1, 60] GeV for CMS and ATLAS respectively. Lastly in the ATLAS 2019 analysis of inclusive dijets at high energy and with high luminosity 1906.09254 , they measured the normalized distribution $1/N_{j}dN_{trk}/d\zeta$ detailed in previous sections. That covered momentum fractions of the parton carried by hadrons from 0.002 to 0.67, as well as a wide range of the transverse momentum of the partons by utilizing jets in finned bins of $p_{T,j}$ from 100 to 2500 GeV. Furthermore, the distributions are measured independently for the central and forward jet of the two leading jets to increase further the discrimination on fragmentation of gluon and quarks.

Despite of the wide coverage on momentum fraction or transverse momentum of the hadrons from above measurements, on the theoretical predictions it requires a careful evaluation on validity of the factorization framework and on stability of the perturbation expansions. There have been previous studies d_Enterria_2014 ; thennpdfcollaboration2019charged showing difficulties on fitting to experimental data in certain kinematic regions indicating large higher-order corrections or even violations of collinear factorization. In this study we take a conservative approach by selecting only those data points corresponding to momentum fractions $x>0.01$ at LO and data points with transverse momenta of the hadrons $p_{T,h}>4$ GeV. Furthermore, we exclude the jet transverse momentum region of [100, 200] GeV for the inclusive dijet measurements since that corresponds to a low transverse momentum of the hadrons in general. Similarly, we exclude the $\xi_{T}^{\gamma}$ regions greater than 3 for the CMS isolated-photon-tagged measurement. These kinematic selections reduce our total number of data points from 569 to 318 as can be seen from Table. 1. In principle one can perform a scan on above kinematic selections and study the stability of the fitted fragmentation functions which we leave for future investigation.

The parameterization form of fragmentation functions to unidentified charged hadrons used at the initial scale $Q_{0}$ is

where $\{\alpha,\beta,a_{n}\}$ are free parameters in the fit. We choose $Q_{0}=5$ GeV and use a zero-mass scheme for heavy quarks with $n_{f}=5$. We assume fragmentation functions equal for all quarks and anti-quarks since the data sets we selected show weak sensitivity on quark flavors of the fragmented partons. The degree of polynomials is set to $p=2$ since improvements of fit by introducing higher-order terms are marginal. Thus the total number of free parameters is 10. The fragmentation functions are evolved to higher scales using two-loop time-like splitting kernels to be consistent with the NLO analysis. The splitting functions was calculated in Refs. Stratmann:1996hn and are implemented in HOPPET Salam:2008qg ; Salam:2008sz which we use in the analysis.

The quality of the agreement between experimental measurements and the corresponding theoretical predictions for a given set of fragmentation parameters is quantified by the log-likelihood function ($\chi^{2}$), which is given by 1709.04922

$N_{\operatorname{pt}}$ is the number of data points, $s^{2}_{k}$ is the total uncorrelated uncertainties by adding statistical and uncorrelated systematic uncertainties in quadrature, $D_{k}$ is the central value of the experimental measurements, and $T_{k}$ is the corresponding theoretical prediction which depends on $\{\alpha,\beta,a_{n}\}$. $\sigma_{k,\mu}$ are the correlated errors from source $\mu$ ($N_{\lambda}$ in total). We assume that the nuisance parameters $\lambda_{\mu}$ follow a standard normal distribution.

By minimizing $\chi^{2}(\{\alpha,\beta,a_{n}\},\{\lambda\})$ with respect to the nuisance parameters, we get the profiled $\chi^{2}$ function

where $\operatorname{cov}^{-1}$ is the inverse of the covariance matrix

We neglect correlations of experimental uncertainties between different data points since they are not available. However, we include theoretical uncertainties into the covariance matrix of Eq. (17) by default, assuming these to be fully correlated among points in each subset of the data shown in Table. 1. Those are data points within the same bin of the transverse momentum of either the photon, $Z$ boson or jets in the measurement. The theoretical uncertainties $\sigma_{j,\mu}$ are estimated by the half width of the scale variations from the prescription mentioned in Sec. 4.1.

The best-fit of fragmentation parameters are found via minimization of the $\chi^{2}$ and further validated through a series of profile scans on each of those parameters. The scans of the parameter space are carried out with MINUIT minuit program. We use the text-book criterion of $\Delta\chi^{2}=1$ on determination of parameter uncertainties. It should be noted that tolerance conditions are usually applied for fits involving multiple data sets 1709.04922 and will lead to conservative estimation of uncertainties. In addition, we adopt the iterative Hessian approach hessian to generate error sets of fragmentation functions that can be used for propagation of parameter uncertainties to physical observable.

The overall agreement between NLO predictions from our nominal fit and the experimental data can be seen from Table. 2. The total $\chi^{2}$ is $267.4$ for a total number of data points of 318. For the ATLAS dijet measurements which contain the majority of the data points, the agreement is quite good with $\chi^{2}/N_{pt}$ well below 1. Description of the isolated-photon measurements is reasonable with $\chi^{2}/N_{pt}\sim 2$. The agreement to CMS $Z$-boson measurement is good while it is much worse for the ATLAS measurement with $\chi^{2}/N_{pt}\sim 5$. The discrepancies to data are mostly due to the low-$p_{T,h}$ kinematic bins ($\sim 4\,{\rm GeV}$) as shown in Appendix B. For comparison we also include results from alternative fits with either excluding to the theoretical uncertainties or using LO matrix elements and LO evolution of the fragmentation functions. Impact of the theoretical uncertainties is mostly pronounced for the CMS isolated-photon and $Z$-boson measurements. The LO fit shows a total $\chi^{2}$ of more than 3000 indicating the necessity of inclusion of NLO corrections.

The values of all 10 parameters of the fragmentation functions from our nominal best-fit are collected in Table. 3. In addition we also calculated the first moment of the quark and gluon fragmentation functions $\langle x\rangle$ which corresponds to the total momentum fraction carried by charged hadrons at the initial scale. The values are 58.6% and 51.0% for quark and gluon respectively. We also show the estimated uncertainties of the fitted parameters in Table. 3 as from both the profile scans and the Hessian calculation. In the latter case, two error sets are generated for each of the 10 orthogonal Hessian directions, and the full uncertainties are obtained by adding uncertainties from individual directions in quadrature 1709.04922 . We find good agreements between uncertainties from the two methods in general. The relative uncertainties of parameters for gluon are $2\sim 5$ times larger than those of the quark, and also are more asymmetric, indicating a larger fraction of quark jets than gluon jet in those measurements.

We compare our fragmentation functions fitted at NLO to those from NNFF1.1, BKK and DSSdss1 ; dss2 ; dss3 ; dss4 ; dss5 as a function of the momentum fraction $x$ and at the scale $Q_{0}=5$ GeV, which is presented in Fig. 8 for $u$-quark and in Fig. 9 for gluon respectively. We should emphasize that in this comparison of our fit we use a restricted parametrization form, namely setting equal of all quark fragmentation functions, which are allowed to be different in fits of NNFF1.1, BKK and DSS, and the error criterion chosen by us is $\Delta\chi^{2}=1$. It should be noted that the small uncertainties in our work, which will be shown below, are partly attributed to the specific assumption we have made in our parametrization.

The upper panel in both figures illustrates the value of the fragmentation function multiplied by the momentum fraction $x$ as a function of $x$ and in the lower one all results are normalized to the central value of our findings. Finally, we need to mention that the colored bands represent the estimated uncertainties from the corresponding fits when they are available.

For the $u$ quark, we observe a good agreement between NNFF1.1 results and our work in the region $0.1<x<0.3$ in Fig. 8. However, significant deviations, especially in the small $x$ region, are observed. These deviations can reach up to nearly 80%. On the other hand, our results remain within the error band of NNFF1.1 results in lower regions. It is important to note that the error band of NNFF1.1 data is large in the small $x$ region (where $x<0.1$), indicating a lack of well-fitted data in that range. Comparatively, the BKK results exhibit significant deviations from our findings throughout the entire region. Particularly in the small $x$ region, the deviation can be as large as 160%. Unfortunately, the error data for BKK is not available. DSS data only fits from 0.05 to 1, and errors are not provided. In the region of $0.05<x<0.3$, a close resemblance with our work is also observed with the maximum deviation about 20%. Furthermore, both the BKK and NNFF1.1 data show a decreasing trend as $x$ increases, DSS data, in its available region, a decreasing trend is also observed, while our results demonstrate an increase as $x$ approaches approximately 0.1, followed by a decrease.

In the case of the gluon, Fig. 9 reveals significant discrepancies among the four results. Firstly, it is evident that these results do not agree except for the BKK and DSS fits, indicating a tension between different fits of the gluon fragmentation function. Secondly, the error band associated with NNFF1.1 data is notably large, suggesting a higher level of uncertainties in that dataset. Further examination reveals that in certain regions (specifically for $x<0.5$), the ratio with respect to our work is less than 4. However, beyond this range, the ratio experiences a sudden and pronounced increase. The upper panel also highlights the contrasting trends exhibited by the BKK and DSS results. Specifically, the BKK results consistently decrease as $x$ increases. But DSS results are not available at $0.01<x<0.05$, therefore, its property is not known in these regions. In contrast, our fit and the NNFF1.1 data initially display an increase, followed by a subsequent decrease. Our results reach their maximum value at approximately $x=0.05$, while the NNFF1.1 data reach their peak at around $x=0.1$. The notable disparities among different datasets emphasize the need for additional constraints and further data to improve the accuracy of the gluon fit.