Automated calculation of jet fragmentation at NLO in QCD

We present FMNLO, a framework to combine general-purpose Monte Carlo generators and fragmentation functions (FFs). It is based on a hybrid scheme of phase-space slicing method and local subtraction method, and accurate to next-to-leading order (NLO) in QCD. The new framework has been interfaced to MG5_aMC@NLO and made publicly available in this work. We demonstrate its unique ability by giving theoretical predictions of various fragmentation measurements at the LHC, followed by comparison with the data. With the help of interpolation techniques, FMNLO allows for fast calculation of fragmentation processes for a large number of different FFs, which makes it a promising tool for future fits of FFs. As an example, we perform a NLO fit of parton fragmentation functions to unidentified charged hadrons using measurements at the LHC. We find the ATLAS data from inclusive dijet production show a strong constraining power. Notable disparities are found between our gluon FF and that of BKK, DSS and NNFF, indicating the necessities of additional constraints and data for gluon fragmentation function.


Introduction
Fragmentations of quarks and gluons into hadrons have been the central topic of QCD since only hadrons are observed experimentally.QCD factorization ensures separation of the short and long distance effects into matrix elements on production of partons and the fragmentation functions(FFs) [1][2][3].In its simplest form, fragmentation functions describe probability distribution on the fraction of momentum of the initial parton carried by the identified hadron.Due to its non-perturbative essential, fragmentation functions are usually extracted from fits to a variety of experimental data.However, dependence of the fragmentation functions on the momentum transfers or the so-called fragmentation scale follows the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation with time-like splitting kernels.For extraction of FFs at next-to-next-to-leading order (NNLO), three-loop evolution kernels at O(α 3 s ) in the strong coupling are needed, which have been calculated in Ref. [4][5][6][7][8][9].The experimental data used to extract FFs includes Single-Inclusive Annihilation (SIA) on lepton colliders, Semi-Inclusive Deep-Inelastic Scattering (SIDIS), and hadron production on hadron-hadron colliders.The corresponding parton-level cross sections for SIA have been calculated at NNLO in Ref. [10][11][12][13].For SIDIS, the next-to-leading order (NLO) corrections are given in Ref. [14][15][16][17][18][19], and approximate NNLO and N 3 LO corrections have also been obtained from expansion of the resummed expressions [20,21].For pp collisions, the corresponding NLO corrections for single-inclusive production of a hadron are given by Ref. [22][23][24][25][26][27].
Various FFs [28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45] have been extracted from SIA, SIDIS and pp collisions.However, there exist several limitations in the current tools on calculations of parton fragmentations at NLO.First, the available processes of hard scattering are limited and are usually implemented case-by-case.Furthermore, interactions in the hard processes are usually constrained to be SM interactions, thus are restrained from applications to study of various new physics beyond the SM (BSM).Besides, even direct calculations at NLO are too costly in computation time to be included into a global analysis of fragmentation functions.
In this work we provide a solution, dubbed FMNLO, by introducing a hybrid scheme of NLO calculations utilizing a phase-space slicing of collinear regions in combination with the usual local subtraction methods.Due to its simplicity we are able to realize the hybrid scheme based on the widely used program MG5 aMC@NLO [46,47].That ensures numerical calculations on partonic cross sections for arbitrary hard processes of fragmentations at NLO within SM and BSMs.We further generate the nominal and convoluted fragmentation functions using the HOPPET program [48,49] for fast convolutions with the partonic cross sections.
The rest of our paper is organized as follows.In Sec. 2, we present the FMNLO framework, which combines partonic cross section calculation and FFs at NLO QCD, in a way suitable for Monte Carlo calculation.We also show how the corresponding calculation can be boosted with interpolation techniques.In Sec. 3, our framework and its implementation are validated by comparing our results with analytic predictions for SIA at lepton colliders and the predictions of other program for inclusive hadron production at hadron colliders.Then we utilize FMNLO to study three cases of hadron production at the LHC and compare our NLO QCD predictions with the experimental measurements in Sec. 4. In Sec. 5, more pp collision measurements at the LHC are considered, and a NLO fit of parton fragmentation functions to unidentified charged hadron is performed, followed by comparisons of our fitted FFs with existing FFs.Finally our summary and conclusions are presented in Sec. 6.

Theoretical framework 2.1 A hybrid scheme
Cross sections for any infrared and collinear (IRC) safe observables in a standard subtraction scheme at NLO have the following schematic form: where R represent square of matrix elements at leading order (LO), oneloop level and in real corrections, respectively.|I| 2 m+1 denotes the local subtraction terms constructed in D = 4 − 2ϵ 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 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 {p m ; fm }.For local subtraction schemes, for instance, in CS dipole subtraction [50,51] or FKS subtraction [52,53], contributions from the second line of Eq. (2.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 {p m ; fm }.Both the virtual corrections and integrated subtraction terms carry poles in ϵ which again cancel among each other which renders the first line of Eq. (2.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 3) 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 MS scheme as where the time-like convolution kernel can be expressed as where µ R and µ D are the renormalization and fragmentation scale respectively, and γ 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. [54].
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. [50,51].In this study we propose a minimal modification of Eq. (2.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 Θ functions to partition into the unresolved and resolved collinear regions with a cutoff λ.The slicing variable C can be chosen as either the minimum of the usual angular separations of all QCD partons ∆θ = min{∆θ ij } in the center of mass frame or the minimum of the boost-invariant angular separation of all QCD partons ∆R = min{∆R ij ≡ ∆ϕ 2 ij + ∆y 2 ij }.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 | J | 2 m 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 λ 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 ji (z, ϵ) ≡ P (0) ji (z), times the square of Born matrix elements.It is understood that the subscripts except the labels (q, q, g) represent flavor of that parton (q, q, g) when appearing in the splitting functions.We arrive at a rather compact form for | J | 2 m after carrying out integrals in s and z, which is given by with Dh/i (x, µ D ) ≡ j I ji ⊗ D h/j (x, µ D ), and the kernel of residuals I ji (z) can be expressed using the unregularized splitting functions, m back to the master formula Eq. (2.6) one can find that the remaining collinear divergences or poles in ϵ 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 = ∆R as an example.It can be easily transformed into the case of C = ∆θ by exchanging λp T,i with λE i in Eq. (2.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. [55].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.

Implementation
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 ≡ 2E h /Q carried by the tagged hadron.Schematically our master formula can be recast as with and dσ 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σ (1) m /dx i and dσ (1) m+1 /dx i can be identified by comparing to various contributions in Eq. (2.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 ( D and 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.3GeV)) 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 δy(δ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 (n) i(j) 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.9) we arrive at In practice we calculate partonic cross sections with MG5 aMC@NLO [46,47] and extract matrix of coefficients G, Ḡ, and 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 [48,49] 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.

Validation
In this section we demonstrate the validity of our calculation scheme and its implementation for several scenarios in both lepton collisions and hadron collisions.We note that in MG5 aMC@NLO the NLO mode of lepton-hadron collisions is not publicly available yet.
Our calculation scheme can be easily implemented once the version including lepton-hadron collisions is released.

Lepton collisions
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σ/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 µ R = µ 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 λ = 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 λ = 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 λ = 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 λ choice ranging from 0.001 to 0.08 and conclude that the choice of λ = 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.Comparison of the NLO predictions on distribution of the hadron energy fraction from FMNLO and from analytical calculations for e + e − → γ * → q q at a center of mass energy of 200 GeV.

Hadron collisions
We compare our calculations with the INCNLO program [56] 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σ/dp T,h .The charged hadrons are required to have rapidity |y| < 2.0.On various theoretical inputs we use CTEQ6M NLO parton distributions [57] and the BKK NLO fragmentation functions [28][29][30].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 λ 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 ∼ 0.04 for numerical stability.In Fig. 4 we present comparisons of our NLO predictions with those from INCNLO1.4 [56] for a finer binning.The agreements of the three predictions with λ = 0.02, 0.04, and 0.08 are similar to that in Fig. 2. From the middle panel we can find our predictions with λ = 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.
4 Applications at the LHC Prescribed calculation scheme and its numerical implementation are especially desirable for predictions of various measurements carried out at the LHC.In typical fragmentation measurements at the LHC, the requirement is often imposed that the tagged hadron is produced either within a reconstructed jet or in association with an isolated photon or a Z boson.Meanwhile, jet algorithms and various selection cuts are applied in the analyses which can be implemented easily in a MC event generator such as MG5 aMC@NLO.In the following, we show three examples of such calculations that we adapted to the corresponding LHC measurements.We focus on the measurements of spectrum of unidentified charged hadrons.Predictions on tagged hadrons with a specified flavor can be obtained easily using the same NLO grids multiplied with the corresponding fragmentation functions.
In the following calculations we use the CT14 NLO parton distribution functions [58], the BKK fragmentation functions [28][29][30] and the NNFF1.1 fragmentation functions [42,43] for unidentified charged hadrons.We set central values of the factorization and renormalization scales (µ F,0 and µ R,0 ) to the default dynamic scale used in MG5 aMC@NLO, namely the sum of the transverse mass of all final state particles divided by 2. For the fragmentation scale, we set its central value (µ D,0 ) to the maximum of the transverse momentum of all final state particles.The above central values equal in the case of only two massless particles in the final states.The scale variations are obtained by taking the envelope of theory predictions of the 9 scale combinations of µ F /µ F,0 = µ R /µ R,0 = {1/2, 1, 2} and µ D /µ D,0 = {1/2, 1, 2}.We note alternative choices on the fragmentation scale of using the transverse momentum of the jet multiplied by the jet cone size when calculating hadron fragmentation inside the jet [26].For typical jet cone sizes of ∼ 0.5 used in the LHC measurements, the choice is close to our nominal choice of the fragmentation scale.

Isolated-photon-tagged jets
In Ref. [59] the CMS collaboration measured parton fragmentation based on hard scattering events in pp collisions ( √ s = 5.02 TeV) consisting of an isolated photon in association with  jets.The photon is required to have a transverse momentum p T,γ > 60 GeV and a pseudorapidity |η γ | < 1.44.Jets are clustered with anti-k T algorithm [60] with R = 0.3 and are required to have p T,j > 30 GeV and |η j | < 1.6.They select jets that have an azimuthal separation to the photon ∆ϕ jγ > 7π/8 and analyze the charged-particle tracks inside the jet with transverse momentum ⃗ p T,h in Ref. [59].The charged tracks are required to have p T,h > 1 GeV.The transverse momentum of the photon ⃗ p T,γ serves as a good reference of the initial transverse momentum of the fragmented parton.Thus )] 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ξ γ T which is simply a linear combination of the quark and gluon fragmentation functions at the LO evaluated at a momentum fraction of e −ξ γ T .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 ξ γ T region, ranging from 0.5 to 2.5.As ξ γ T increases, the discrepancy enlarges, which can be attributed to the lack of fitted data in the small x = p T,h /p T,γ (< 0.01) regions [28][29][30].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 ξ γ T .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 ξ γ T < 3.5.However, for ξ γ T > 3.5, the contribution switches from negative to positive, and the ratio rises rapidly to nearly 2.

Z boson tagged jets
We now turn to the relevant calculations for the production process of Z boson in association with jets at LHC.In Ref. [61] the CMS collaboration measured parton fragmentation based on hard scattering events of the above process, where √ 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 ∆ϕ trk,Z > 7π/8 in Ref. [61].The charged tracks are required to have p T,h > 1 GeV and |η 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%.

QCD inclusive dijets
In this subsection, we present the third example of the calculations mentioned above.In Ref. [62] the ATLAS collaboration measured parton fragmentation based on hard scattering events in pp collisions ( √ 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 |η 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. [62].The charged tracks are required to have p T,h > 0.5 GeV and |η h | < 2.5.The results are presented in a differential cross section of 1/N j dN trk /dζ with ζ ≡ p T,h /p T,j and p T,j being the transverse momentum of the jet probed1 .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 ∈ [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 ζ region.However, the BKK data aligns more closely with the experimental data.In the lower ζ 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 ζ to positive in the higher region.

Analysis of fragmentation functions
In this section we perform a NLO fit of the parton fragmentation functions to unidentified charged hadrons using a variety of experimental data from pp collisions at the LHC.Those include processes on production of charged hadrons from inclusive dijets, in association with an isolated photon and in association with a Z boson.They can probe fragmentation of both gluon and quarks in a wide kinematic region due to different production mechanisms involved.We demonstrate that such a fit at NLO accuracy with a few hundreds of experimental data points can be accomplished easily with the help of the FMNLO framework.In the following we first briefly introduce our selection of experimental data sets and the fitting framework, and then show our best-fit and the estimated uncertainties of the fragmentation functions.

Experimental data sets
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 [59] measured the normalized distribution 1/N j dN trk /dξ γ T that has been explained in previous sections.They probe a region of momentum fractions of the parton carried by hadrons from 0.01 ∼ 0.6 based on the definition of ξ γ T .The ATLAS 2019 analysis [63] 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 pseudorapidity 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 ∈ [1, 100] GeV, that is used in our fit.For measurements involving Z-tagged jets, the CMS 2021 analysis [61] measured the normalized distribution 1/N Z dN trk /dp T,h with setups detailed in previous sections.In the ATLAS 2020 analysis [64], 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 ∆ϕ trk,Z > 3π/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 [62], they measured the normalized distribution 1/N j dN trk /dζ 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.

Experiments
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 [65,66] 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 ξ γ T 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.

Framework of the fit
The parameterization form of fragmentation functions to unidentified charged hadrons used at the initial scale Q 0 is where {α, β, 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.[54] and are implemented in HOPPET [48,49] 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 loglikelihood function (χ 2 ), which is given by [67] N 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 {α, β, a n }. σ k,µ are the correlated errors from source µ (N λ in total).We assume that the nuisance parameters λ µ follow a standard normal distribution.
By minimizing χ 2 ({α, β, a n }, {λ}) with respect to the nuisance parameters, we get the profiled χ 2 function where 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. ( 5.4) 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 σ j,µ 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 χ 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 [68] program.We use the text-book criterion of ∆χ 2 = 1 on determination of parameter uncertainties.It should be noted that tolerance conditions are usually applied for fits involving multiple data sets [67] and will lead to conservative estimation of uncertainties.In addition, we adopt the iterative Hessian approach [69] to generate error sets of fragmentation functions that can be used for propagation of parameter uncertainties to physical observable.

Results and discussions
The overall agreement between NLO predictions from our nominal fit and the experimental data can be seen from Table .2. The total χ 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 χ 2 /N pt well below 1. Description of the isolated-photon measurements is reasonable with χ 2 /N pt ∼ 2. The agreement to CMS Z-boson measurement is good while it is much worse for the ATLAS measurement with χ 2 /N pt ∼ 5.The discrepancies to data are mostly due to the low-p T,h kinematic bins (∼ 4 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 χ 2 of more than 3000 indicating the necessity of inclusion of NLO corrections.
The values of all 10 parameters of the fragmentation functions from  3. The best-fit parameters for quark and gluon from our nominal NLO fit and their uncertainties (68% C.L.) estimated using profile scans or Hessian method.The last column is the first moment of the fragmentation functions at the initial scale as calculated using the fitted functional forms.
are collected in Table .3. In addition we also calculated the first moment of the quark and gluon fragmentation functions ⟨x⟩ 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 [67].We find good agreements between uncertainties from the two methods in general.The relative uncertainties of parameters for gluon are 2 ∼ 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 DSS [36][37][38][39][40] 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 ∆χ 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.

Conclusions
In this work, we propose a new prescription for combining general-purpose Monte-Carlo generators with fragmentation functions (FFs) at NLO in QCD.This new framework, dubbed FMNLO, is based on a hybrid scheme of NLO calculations utilizing a phase-space slicing of collinear regions in combination with the usual local subtraction methods, and organizes various ingredients for fragmentation predictions in a way suitable for Monte Carlo calculations.As a proof of concept, we realize FMNLO with MG5 aMC@NLO.The corresponding code is publicly available and is introduced in Appendix A. Our scheme and its implementation are validated for several scenarios in both lepton collisions and hadron collisions.
The combination of general-purpose MC generators and FFs allows for the study of singlehadron production for various hard process at NLO in QCD with general selection cuts or jet reconstruction.As examples, we compare the predictions of FMNLO with experimental measurements of jet production with a tagged isolated photon, jet production with a tagged Z boson, and inclusive dijet production.Also, we boost FMNLO with interpolation techniques, such that for a given measurement, the time-consuming calculation of NLO partonic cross section can be reused when the fragmentation functions are changed.The combination of these two features endows it with the unique ability of making theoretical predictions for a wide range of measurements within a reasonable time consumption, which is essential for a global fit of FFs.We demonstrate this ability by performing a NLO fit of parton FFs to unidentified charged hadrons, using hadron production measurements at the LHC.Our nominal fit shows very good agreements with the LHC data.We find that the high-precision fragmentation measurements from ATLAS inclusive dijet production especially show a strong constraining power on the FFs.Our unidentified charged-hadron FFs are then compared with those from BKK, DSS and NNFF1.1.Notable disparity in gluon FF is found, indicating the necessities of additional constraints and data in gluon fit.We emphasize that our framework also works for FFs with specific flavors.
Besides its ability in extraction of FFs, the proposed scheme and its implementation open the opportunity of studying BSM effects with single-hadron production.FMNLO is also desirable for calculations of NLO hard functions needed for various predictions of QCD resummation [70][71][72][73].Furthermore, it can be generalized to calculate distributions of observable related to transverse dependent fragmentation functions which have been widely used in studies of jet substructures [8,74,75].We leave those for future investigations and updates of the program.in MG5 aMC@NLO, one has to modify madgraph/core/helas objects.py in the native MG directory to disable group of identical processes.That can be achieved by simply replace True with False at the 3486-th line of the file.We mentioned earlier that the FKS subtractions have not been implemented for NLO calculations of DIS processes in MG5 aMC@NLO.However, we are working toward an alternative solution that will come with the next release of FMNLO.
To calculate the parton fragmentation in a typical hard scattering process, two subsequent steps are followed.First, inside the directory mgen one invokes MG5 aMC@NLO with a customized analysis routine (a module) to generate the interpolation tables storing matrices of coefficient functions in Eq. 2.12.We have released modules for all processes included in this study.New modules can be added easily following existing examples.In mgen, the subdirectory common contains the common ingredients needed for all modules.Each module has a separate directory including init.sh for MG5 aMC@NLO command, pre cuts.ffor the selection of relevant phase space, and analysis HwU pp.f containing the main analysis routine.
Various input parameters are specified in the file proc.run.Each line contains a record for one input variable: a character tag with the name of the variable, followed by the variable's value.We take proc.runused for the calculation of CMS isolated-photon-tagged jet measurement as an example.
# main input for generation of NLO fragmentation grid file by MG5 process A180104895 # subgrids with name tags grid pp obs 4 cut 0.02 pta1 60.0 pta2 10000.0ptj1 30.0 ptj2 10000.0# in MG5 format set lpp1 1 set lpp2 1 set ebeam1 2510.0 set ebeam2 2510.0 set lhaid 13100 set iseed 11 set muR_over_ref 1.0 set muF_over_ref 1.0 end • process specifies the name of the directory that contains the module to be loaded.
• grid is a string indicating the name of the running job.
• obs specifies different distributions to be calculated: 1 for distribution in ζ, 2 for distribution in p T,h , 3 and 4 for distributions in ξ j T or ξ γ(Z) T .
• cut gives the slicing parameter λ and a value of 0.02 is recommended.
• pta1 and pta2 specify the lower and upper limit of the kinematic range of the transverse momentum of the photon.
• ptj1 and ptj2 specify the lower and upper limit of the kinematic range of the transverse momentum of the jet.
• Other possible inputs in this block include hrap and pth specifying the upper limit on the absolute pseudorapidity and the lower limit on the transverse momentum of the hadrons.isca determines the choice on central value of the fragmentation scale, 1 for our nominal choice of max{p T,j } (Q) for pp (e + e − ) collisions, and 2 for using p T,h (E h ).Default values of all above variables are assigned via the script file mgen.sh for individual modules if not specified in proc.run.
• The remaining inputs follow the same syntax as the normal MG5 aMC@NLO command, for instance, lpp1 and ebeam1 are type and energy of collision particle 1, lhaid specifies parton distribution of proton used for the calculation, etc.
The generation of fragmentation grid can be launched by the command ./mgen.sh proc.run.
Note that for the same process, generation of multiple grids can be grouped into a single input file by simply repeating the two blocks after the process line.Once the generation of grid is finished, it will be stored in an upper-level directory grid, for instance with a name A180104895 pp.fmg for above example.
After generation of the grid, the calculation of physical distributions can be done within seconds by running ./fmnlo in the directory data.Input parameters at this stage are specified in the file input.card.for LO and NLO respectively, followed by 0/1 for using the native evolution provided by the input fragmentation functions or evolving with HOPPET package from the initial scale Q 0 .
• The 4th line indicates the choice of fragmentation functions.A value of 0 indicates the usage of fragmentation functions from LHAPDF6 library and other integer values correspond to fragmentation functions implemented in FMNLO1.0,e.g., 1 for the NLO nominal fragmentation functions presented in this study (note one should use the NLO evolution with HOPPET concurrently), 2 for the BKK functions of unidentified charged hadrons, and 3(4) for the KKP [31](DSS) functions.The following two inputs specify the name and set number of the fragmentation function in the case of using LHAPDF6.Note that the value of the QCD coupling used in the evolution of fragmentation functions is set consistently.It can either be imported from LHAPDF6 or set by HOPPET with α S (M Z ) fixed at 0.118.Other possibilities on the choice of fragmentation function can be implemented by modifying the source file internal.f.
• The 8th line specifies choice of normalization: 0 for absolute distributions, 1 and 2 for normalized distributions to the total cross sections of corresponding order or LO respectively.The followed are the name of the pre-generated grid file for the calculation and the name of the file containing binning of the distribution.Multiple entries similar to the 8th line can be added to calculate several distributions at once.The binning is set via two-line inputs.For example, in 1801-04895.Bin, 8 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 the first line specifies the total number of kinematic bins and the second line contains all nodes of the bins in sequence.
Once ./fmnlo is executed, the format of the printout can be understood easily

B Comparison of the theory to data
In this appendix we include more details on our nominal NLO fit to the LHC data.We first show the total χ 2 profile from scans on individual parameters of the fragmentation functions in Fig. 10 and 11 for quark and gluon respectively.In each scan all other parameters are set free and are fitted to minimize the constrained χ 2 .The fragmentation function of quarks is better constrained as mentioned above.The χ 2 profile shows a parabolic shape around the minimum.The fragmentation function of gluon is less constrained especially for parameters a 0/1/2 .For instance, the increase of χ 2 can barely reach 2 units even we scan over a wide range of a 0 and a 1 .The large variations of the parameters also lead to a non-parabolic shape of the χ 2 profile in the scan region.
Our NLO predictions based on the best-fit of the fragmentation functions and their comparison to data are presented in Figs. 12 to 15.In the lower panel of each figure the predictions are normalized to the central value of the experimental data.The colored bands indicate estimated theoretical uncertainties from scale variations as explained in Sec.4.1.We show the comparison to isolated-photon production in Fig. 12 for both the CMS and ATLAS measurements.The theoretical predictions locate within 10% of the data in general with the exception of the first(last) bin in ξ γ T (p T,h ) which corresponds to very large x region.For the Z boson production shown in Fig. 13, we again find good agreements to both the CMS and ATLAS measurements.Furthermore, there is consistency observed between the two independent measurements, providing additional confidence in the accuracy of our calculations.The large χ 2 observed for the ATLAS data is mostly driven by the first p T,h bin of the highest p T,Z region (> 60 GeV), which has a high precision of a few percents.Notably it has a p T,h around 4 GeV and receives contributions from region of momentum fractions x ≲ 0.01, and may be affected by additional uncertainties from both theoretical and experimental sides.Lastly in Figs.14-15 we present comparisons to the ATLAS dijet measurements for both the central and the forward jet and for the selected bins of the transverse momentum of the jet.We find excellent agreements between our theoretical predictions and the data in the full kinematic region even without considering the theoretical uncertainties.The scale variations are almost an order of magnitude larger than the experimental uncertainties.

Figure 1 .
Figure1.Comparison of the NLO predictions on distribution of the hadron energy fraction from FMNLO and from analytical calculations for e + e − → γ * → q q at a center of mass energy of 200 GeV.

Figure 2 .
Figure 2. Comparison of the NLO predictions on distribution of the hadron energy fraction from FMNLO and from analytical calculations for µ + µ − → H * → gg at a center of mass energy of 200 GeV.

Figure 3 .
Figure3.Ratio of the NLO predictions on distribution of the hadron transverse momentum from FMNLO with different choices of cut-off parameter relative to λ = 0.04 for inclusive jet production in pp collisions with a center of mass energy of 7 TeV.Three representative bins on the transverse momentum have been selected, and a rapidity cut of |y h | < 2 of hadrons has been applied.

Figure 4 .
Figure 4. Comparison of the NLO predictions on distribution of the hadron transverse momentum from FMNLO and from INCNLO for inclusive jet production in pp collisions with a center of mass energy of 7 TeV.A rapidity cut of |y h | < 2 of hadrons has been applied.

Figure 5 .
Figure 5.Comparison of NLO predictions to CMS measurement on normalized distribution of ξ γ T for isolated photon production in pp collisions with a center of mass energy of 5.02 TeV.The two colored bands represent predictions including scale variations, based on NNFF1.1 and BKK fragmentation functions respectively.The error bars indicate the total experimental uncertainties.Theoretical predictions have been normalized to the central value of data in the middle panel.In the lower panel the two bands correspond to LO and NLO predictions based on NNFF1.1, normalized to the NLO prediction with nominal scale choice.

Figure 6 .
Figure 6.Similar to Fig. 5 but with CMS measurement on normalized distribution of p T,h for Z boson production in pp collisions with a center of mass energy of 5.02 TeV.

Figure 7 .
Figure 7. Similar to Fig. 5 but with ATLAS measurement on normalized distribution of ζ for dijet production in pp collisions with a center of mass energy of 13 TeV.

Figure 8 .
Figure 8.The u-quark fragmentation function at Q 0 = 5 GeV from our nominal NLO fit as a function of the momentum fraction x, and its comparison to the NNFF1.1,BKK and DSS results.The colored bands indicate the uncertainties as estimated with the Hessian (MC) method for our (NNFF1.1)fit.

Figure 9 .
Figure 9.The gluon fragmentation function at Q 0 = 5 GeV from our nominal NLO fit as a function of the momentum fraction x, and its comparison to the NNFF1.1,BKK and DSS results.The colored bands indicate the uncertainties as estimated with the Hessian (MC) method for our (NNFF1.1)fit.

Figure 10 .
Figure 10.Profile of total χ 2 change from scans on individual parameters of the quark fragmentation functions under nominal fit with other parameters freely varying.

Figure 11 .Figure 12 .
Figure 11.Profile of total χ 2 change from scans on individual parameters of the gluon fragmentation functions under nominal fit with other parameters freely varying.

Figure 13 .
Figure 13.Similar to Fig. 12 but for the Z boson production and in different bins of transverse momentum of the Z boson.

C. 1 )Figure 14 .
Figure 14.Similar to Fig. 12 but for the dijet production and in different bins of transverse momentum of the central jet.

Figure 15 .
Figure 15.Similar to Fig. 12 but for the dijet production and in different bins of transverse momentum of the forward jet.

Table 1 .
Summary on experimental data sets used in this analysis, including the observable measured, the number of data points before and after data selection, and the kinematic range covered.

Table 2 .
The χ 2 of individual data sets and their sum from our nominal NLO fit and alternative fits at NLO and LO without theoretical uncertainties.Numbers in parenthesis correspond to χ 2 divided by the number of data points.