Energy evolution of the moments of the hadron distribution in QCD jets

including NNLL resummation and NLO running-coupling corrections

Redamy Pérez-Ramos***e-mail:  and David d’Enterriae-mail:

Department of Physics, University of Jyväskylä, P.O. Box 35 (YFL), F-40014 Jyväskylä, Finland

CERN, PH Department, CH-1211 Geneva 23, Switzerland


The moments of the single inclusive momentum distribution of hadrons in QCD jets, are studied in the next-to-modified-leading-log approximation (NMLLA) including next-to-leading-order (NLO) corrections to the strong coupling. The evolution equations are solved using a distorted Gaussian parametrisation, which successfully reproduces the spectrum of charged hadrons of jets measured in collisions. The energy dependencies of the maximum peak, multiplicity, width, kurtosis and skewness of the jet hadron distribution are computed analytically. Comparisons of all the existing jet data measured in collisions in the range  2–200 GeV to the NMLLANLO predictions allow one to extract a value of the QCD parameter , and associated two-loop coupling constant at the Z resonance  = 0.1195  0.0022, in excellent numerical agreement with the current world average obtained using other methods.

1 Introduction

One of the most ubiquitous manifestations of the fundamental degrees of freedom of Quantum Chromodynamics (QCD), quark and gluons, are the collimated bunches of hadrons (“jets”) produced in high-energy particle collisions. The evolution of a parton into a final distribution of hadrons is driven by perturbative dynamics dominated by soft and collinear gluon bremsstrahlung [1, 2] followed by the final conversion of the radiated partons into hadrons at non-perturbative scales approaching  0.2 GeV. The quantitative description of the distribution of hadrons of type in a jet is encoded in a (dimensionless) fragmentation function (FF) which can be experimentally obtained, e.g. in collisions at c.m. energy , via

where is the scaled momentum of hadron , and the total hadronic cross section. Its integral over gives the average hadron multiplicity in jets. Writing the FF as a function of the (log of the) inverse of , , emphasises the region of relatively low momenta that dominates the spectrum of hadrons inside a jet. Indeed, the emission of successive gluons inside a jet follows a parton cascade where the emission angles decrease as the jet evolves towards the hadronisation stage, the so-called “angular ordering” [3, 1, 4]. Thus, due to QCD colour coherence and interference of gluon radiation, not the softest partons but those with intermediate energies () multiply most effectively in QCD cascades [4]. As a result, the energy spectrum of hadrons as a function of takes a typical “hump-backed plateau” (HBP) shape [4, 5], confirmed by jet measurements at LEP [6] and Tevatron [7] colliders, that can be written to first approximation in a Gaussian form of peak and width :


where is the collinear cut-off parameter of the perturbative expansion which can be pushed down to the value of (the so-called “limiting spectrum”). Both the HBP peak and width evolve approximately logarithmically with the energy of the jet: the hadron distribution peaks at  2 (5) GeV with a dispersion of  0.7 (1.4) GeV, for a parton with  = 10 GeV (1 TeV).

The measured fragmentation function (1) corresponds to the sum of contributions from the fragmentation of different primary partons :

and, although one cannot compute from perturbation theory the final parton-to-hadron transition encoded in , the evolution of the “intermediate” functions describing the branching of a parton of type into partons of type , can be indeed theoretically predicted. The relevant kinematical variables in the parton splitting process are shown in Fig. 1 for the splitting , such that and carry the energy-momentum fractions and of respectively. The Sudakov parametrisation for and , the four-momentum of partons and , can be written as


with the light-like vector , and time-like transverse momentum such that, . Then, the scalar product reads:


Writing now the 4-momenta , , one has, , for on-shell and massless partons . From energy-momentum conservation:


such that, replacing Eq. (4) in (3), one finally obtains:


In the collinear limit, one is left with , where is the jet virtuality, or transverse momentum of the jet.

Figure 1: Relevant kinematical variables in the parton splitting process : is the energy of the leading quark or gluon of virtuality ; and () are the energy fractions of the intermediate offsprings and which finally fragment (at virtualities ) into hadrons carrying a fraction of the parent parton momentum.

The calculation of the evolution of inside a jet suffers from two types of singularities at each order in the strong coupling : collinear -singularities when the gluon emission angle is very small (), and infrared -singularities when the emitted gluon takes a very small fraction of the energy of the parent parton. Various perturbative resummation schemes have been developed to deal with such singularities: (i) the Leading Logarithmic Approximation (LLA) resums single logs of the type where is the transverse momentum of the emitted gluon with respect to the parent parton, (ii) the Double Logarithmic Approximation (DLA) resums soft-collinear and infrared gluons, and , for small values of and  [8, 9], (iii) Single Logarithms (SL) [4, 10] account for the emission of hard collinear gluons (), , and (iv) the Modified Leading Logarithmic Approximation (MLLA) provides a SL correction to the DLA, resumming terms of order  [4]. While the DLA resummation scheme [10] is known to overestimate the cascading process, as it neglects the recoil of the parent parton with respect to its offspring after radiation [9], the MLLA approximation reproduces very well the data, although Tevatron jet results require further (next-to-MLLA, or NMLLA) refinements [11, 12]. The MLLA [4], partially restores the energy-momentum balance by including SL corrections of order coming from the emission of hard-collinear gluons and quarks at large and small (, and ). Such corrections are included in the standard Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) [13, 14, 15] splitting functions which describe the parton evolution at intermediate and large in the (time-like) FFs and (space-like) parton distribution functions (PDFs). The first comparison of the MLLA analytical results to the inclusive particle spectra in jets, determining the energy evolution of the HBP peak position was performed in [16].

The solution of the evolution equations for the gluon and quark jets is usually obtained writing the FF in the form

where are the coefficient functions, and is the so-called anomalous dimension, which in Mellin space at LLA reads,

where is the energy of the radiated gluon and the number of colours. At small or , the expansion of the FF expression leads to a series of half-powers of , , while at larger or in DGLAP, the expansion yields to a series of integer powers of , for FFs and PDFs. In the present work we are mostly concerned with series of half-powers of generated at small , which can be truncated beyond in the high-energy limit.

In this paper, the set of next-to-MLLA corrections of order for the single inclusive hadron distribution in jets, which further improve energy conservation [17, 18], including in addition the running of the coupling constant at two-loop or next-to-leading order (NLO) [19], are computed for the first time. Corrections beyond MLLA were considered first in [20], and more recently in [21], for the calculation of the jet mean multiplicity and the ratio in gluon and quark jets. We will follow the resummation scheme presented in [20] and apply it not just to the jet multiplicities but extend it to the full properties of parton fragmentation functions using the distorted Gaussian (DG) parametrisation [22] for the HBP which was only used so far to compute the evolution of FFs at MLLA. The approach followed consists in writing the exponential of Eq. (1) as a DG with mean peak and width , including higher moments (skewness and kurtosis) that provide an improved shape of the quasi-Gaussian behaviour of the final distribution of hadrons, and compute the energy evolution of all its (normalised) moments at NMLLA+NLO accuracy, which just depend on as a single free parameter.

Since the evolution of each moment is independent of the ansatz for the initial conditions assumed for the jet hadron spectrum, and since each moment evolves independently of one another, we can obtain five different constraints on . By fitting all the measured jet distributions in the range of collision energies  2–200 GeV [23, 24, 25, 26, 27, 28, 29, 6, 30, 31, 32, 33, 34, 34, 35, 36, 37] a value of can be extracted which agrees very well with that obtained from the NLO coupling constant evaluated at the resonance, , in the minimal subtraction () factorisation scheme [38, 39, 40]. Similar studies –at (N)MLLA+LO accuracy under different approximations, and with a more reduced experimental data-set– were done previously for various parametrizations of the input fragmentation function [41, 42, 43, 44] but only with a relatively modest data-theory agreement, and an extracted LO value of with large uncertainties.

The paper is organised as follows. In Sect. 2 we write the evolution equations and provide the generic solution including the set of terms from the splitting functions in Mellin space. In subsection 3.1, the new NMLLANLO anomalous dimension, , is obtained from the evolution equations in Mellin space, being the main theoretical result of this paper. In subsection 3.2 the Fong and Webber DG parametrisation [22] for the single-inclusive hadron distribution is used and the energy evolution of its moments (mean multiplicity, peak position, width, skewness and kurtosis) is computed making use of . In subsection 3.3, the results of our approach are compared for the quark and gluon multiplicities, recovering the NMLLA multiplicity ratio first obtained in [17]. The energy-evolution for all the moments in the limiting spectrum case () are derived in subsection 3.4, and the role of higher-order corrections contributing to the resummed components of the DG which improve the overall behaviour of the perturbative series, are discussed in subsection 3.5, and the final analytical formulæ are provided. Subsection 3.6 discusses our treatment of finite-mass effects and heavy-quark thresholds, as well as other subleading corrections. The phenomenological comparison of our analytical results to the world jet data is carried out in Sect. 4, from which a value of can be extracted from the fits. Our results are summarised in Sect. 5 and the appendices provide more details on various ingredients used in the calculations.

2 Evolution equations for the low- parton fragmentation functions

The fragmentation function of a parton splitting into partons and satisfies the following system of evolution equations [4, 5] as a function of the variables defined in Fig. 1:


where are the regularised DGLAP splitting functions [13, 14, 15], which at LO are given by


with and respectively the Casimirs of the fundamental and adjoint representation of the QCD colour group , , and is the number of active (anti)quark flavours. The regularisation of the splitting functions in Eq. (6) is performed through the distributionThe plus distribution applied to a function , written , is defined as for any function . in Eqs. (7) and (8). The is the strong coupling which at the two-loop level reads [19]



being the first two coefficients involved in the perturbative expansion of the -function through the renormalisation group equation:

The initial condition for the system of evolution equations (6) is given by a delta function

running “backwards” from the end of the parton branching process, with a clear physical interpretation: when the transverse momentum of the leading parton is low enough, it can not fragment () and hadronises into a single hadron. The equations (6) are identical to the DGLAP evolution equations but for one detail, the shift in in the second argument of the fragmentation function , that for hard partons is set to zero, , in the LLA. It corresponds to the so-called scaling violation of DGLAP FFs in time-like evolution, and that of space-like evolution of PDFs in in DIS. In our framework, however, this term is responsible for the double soft-collinear contributions that are resummed at all orders as , justifying the fact that the approach is said to be modified (MLLA) with respect to the LLA.

The evolution equations are commonly expressed as a function of two variables:


where provides the parton-energy dependence of the fragmentation process, and the specifies, in units of , the value of the hadronisation scale down to which the parton shower is evolved. Standard parton showers Monte Carlo codes, such as pythia [45], use values of the order of whereas in the limiting spectrum [4], that will be used here, it can be taken as low as , i.e. . Applying the Mellin transform to the single inclusive distribution in Eq. (6)


and introducing


with in the soft approximation (), one is left with the integro-differential system of evolution equations for the non-singlet distributions




and the lower and upper indices have been omitted for the sake of simplicity. The NLO strong coupling (9) can be rewritten as a function of the new variables (12), such that


The parton density is then obtained through the inverse Mellin transform:


where the contour lies to the right of all singularities in the -complex plane. In the high-energy limit () and hard fragmentation region ( or ), one can replace in the r.h.s. of Eq. (13) the following expansion§§§Note that the MLLA solution [4] to the evolution equations corresponds to the replacement accounting for the single logarithmic corrections of relative order .:


Thus, replacing Eq. (17) into (13) one obtains


which allows for the factorisation of , and leads to the equation


more suitable for analytical solutions. Truncating the series at higher orders translates into including corrections which better account for energy conservation, particularly at large . In Mellin space, the expansion can be made in terms of the differential operator such that, up to the second term in , one is left with NMLLA corrections of order  [11]. Explicitly, the inclusion of higher-order corrections from the second term of , followed by the integration over the splitting functions (7)–(8) in space in the r.h.s. of Eq. (13), is equivalent to the expansion in Mellin space in the r.h.s. of (19), where and are constants. The expansion of the matrix elements in can be obtained from the original expressions of the Mellin transformed splitting functions [46], as given in Eqs. (139a)–(139d) in Appendix A, which leads to the following expressions:


where the finite terms for constitute the new subset to be computed for the first time in this work. The solution of the evolution equations in the MLLA were considered in [4] up to the regular terms with . By including those proportional to , one is in addition considering the set of higher-order corrections known as NMLLA that improve energy conservation [20]. The diagonalisation of the matrix (14) in order to solve (19) results into two trajectories (eigenvalues), which can be written as [4, 46]


Substituting Eqs. (20a)–(20d) into (21) and performing the expansion again up to terms , yields:


where the terms proportional to are new in this framework. The set of constants involved in Eqs. (22a) and (22b) reads:


Therefore, the diagonalisation of Eq. (19) leads to two equations:


such that in the new -basis the respective solutions read:


where the ratios in front of are the coefficient functions that will be evaluated hereafter. Notice that in the basis, the off-diagonal terms and vanish for LO splitting functions, while this is no longer true for time-like splitting functions obtained from the factorisation scheme beyond LO [47], as explained in [21] for multiparticle production. Following this logic, should first be determined in order to obtain the gluon and quark jets single inclusive distributions.

3 Evolution of the parton fragmentation functions at NMLLA +NLO

3.1 Anomalous dimension at NMLLA +NLO

Our NMLLANLO scheme involves adding further corrections from contributions proportional to in the Mellin representation of the expanded splitting functions, and considering the two-loop strong coupling, Eq. (15). We label our approach as NLO to indicate that the full set of NLO corrections are only approximately included, as the two-loop splitting functions (discussed e.g. in [21]) are not incorporated. After diagonalisation of the original evolution equations (6), the Eqs. (24) for result in the following expressions for and :


The leading contribution to after setting in Eq. (27) reads:


The exponent induces a very small (non-Gaussian) correction, which can be neglected asymptotically, for . Thus, the (+) trajectory (22a) provides the main contribution to the single inclusive distribution at small , after applying the inverse Mellin transform (16). Hard corrections proportional to and account for the energy balance in the hard fragmentation region and are of relative order and respectively with respect to the DLA contribution. The NLO expression (9) results in corrections at MLLA, and at NMLLA which provide a more accurate consideration of running coupling effects at small  [20]. In Ref. [20], the mean multiplicities, multiplicity correlators in gluon and quark jets, and the ratio of gluon and quark jet multiplicities were also studied at NMLLA, where corrections were accordingly included. Here, we extend the NMLLA analysis to all moments of the fragmentation function.

The solution of Eq. (26) can be written in the compact form:


with the evolution “Hamiltonian”:


that describes the parton jet evolution from its initial virtuality to the lowest possible energy scale , at which the parton-to-hadron transition occurs. In Eq. (30), is the anomalous dimension that mixes and splittings and is mainly dominated by soft gluon bremsstrahlung (). Introducing the shorthand notation , the MLLA anomalous dimension has been determined in the past [22, 4], setting and in Eq. (26), and is given by


where is the DLA anomalous dimension amounting to


The first term of Eq. (32) is the DLA main contribution, of order , which physically accounts for soft gluon multiplication, the second and third terms are SL corrections accounting for the energy balance () and running coupling effects (). It is important to make the difference between orders and relative orders mentioned above. Indeed, if one looks at the l.h.s. of the evolution equation (26) for , , the first term in the r.h.s is , the second one proportional to is , and the third one, proportional to , is such that after factorising the whole equation by one is left with the relative orders of magnitude in . Setting Eq. (29) in (26) leads to the perturbative differential equation


which will be solved after inserting the two-loop coupling (9) in order to include corrections as well. The equation can be solved iteratively (perturbatively) by setting the MLLA anomalous dimension written in Eq. (32) in the main and subleading contributions of Eq. (34), to find:


which is the main theoretical result of this paper. Terms proportional to , and are of order , and were previously calculated in the (N)MLLA+LO scheme described in [42]. Those proportional to and are computed for the first time in our NMLLA+NLO* framework. Indeed, the single correction is obtained replacing Eq. (9) in the l.h.s. of (34), which leads to the equation,

with . Since and , and following power counting, this correction has naturally the same order of magnitude as the other terms and should not be neglected. The other new correction adds those NMLLA contributions arising from the terms in the LO splitting functions (20a)–(20d), known to better account for energy conservation. Since this correction is multiplied by a term , the overall result is and, thus, of the same order of magnitude as the previous terms such that, the full resummed result is .

3.2 Distorted Gaussian (DG) parametrisation for the fragmentation function

The distorted Gaussian (DG) parametrisation of the single inclusive distribution of hadrons in jets at small (or ) was introduced by Fong and Webber in 1991 [22], and in -space it reads:


where, , is the asymptotic average multiplicity inside a jet, and , , , and are respectively the mean peak position, the dispersion, the skewness, and kurtosis of the distribution. The distribution should be displayed in the interval which depends on the jet energy, and the values of and . The three scales of the process are organised in the form . The formula (38) reduces to a Gaussian for and its generic expression does not depend on the approach or level of accuracy used for the computation of its evolution.

Figure 2: Comparison of various Gaussian-like hadron distributions in jets sharing the same mean position and width ( and ) but with different third and fourth moments: (i) symmetric Gaussian, (ii) skewness , (iii) negative kurtosis , and (iv) full distorted Gaussian with .

As an example of the effects of non-zero skewness and kurtosis, we compare in Fig. 2 the shapes of four different single-inclusive hadron distributions of width and mean position at in the interval typical of jets at LEP-1 energies: (i) an exact Gaussian, (ii) a skewed Gaussian with , , (iii) a kurtic Gaussian with , , and (iv) a DG including both “distorting” components above. As can be seen, the shape of the DG differs from that of the pure Gaussian, mainly away from the hump region. A negative skewness displaces the peak of the Gaussian to higher values while adding a longer tail to low , and a negative kurtosis tends to make “fatter” its width.

In order to obtain the evolution of the different DG components, we will proceed by following the same steps as in [22] but making use instead of the expanded NMLLANLO anomalous dimension, Eq. (37), computed here. Defining as the -th moment of the single inclusive distribution:


the different components (normalised moments) of the DG are given byWe list also which is needed to obtain the maximum peak position from , as discussed below.:


such that after plugging Eq. (30) into (29) and what results from it into (39), one is left with


which is more suitable for analytical calculations since it directly involves the anomalous dimension expression (37).


The multiplicity is obtained from the zeroth moment, i.e. the integral, of the single-particle distribution. Setting in Eq. (37), one obtains