Energy evolution of the moments of the hadron distribution in QCD jets
including NNLL resummation and NLO runningcoupling corrections
Redamy PérezRamos^{*}^{*}*email: redamy.perez@uv.es and David d’Enterria^{†}^{†}†email: dde@cern.ch
Department of Physics, University of Jyväskylä, P.O. Box 35 (YFL), F40014 Jyväskylä, Finland
CERN, PH Department, CH1211 Geneva 23, Switzerland
Abstract
The moments of the single inclusive momentum distribution of hadrons in QCD jets, are studied in the nexttomodifiedleadinglog approximation (NMLLA) including nexttoleadingorder (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 twoloop coupling constant at the Z resonance = 0.1195 0.0022, in excellent numerical agreement with the current world average obtained using other methods.
Contents
 1 Introduction
 2 Evolution equations for the low parton fragmentation functions

3 Evolution of the parton fragmentation functions at NMLLA +NLO
 3.1 Anomalous dimension at NMLLA +NLO
 3.2 Distorted Gaussian (DG) parametrisation for the fragmentation function
 3.3 Multiplicities for the single inclusive and distributions
 3.4 Limiting spectrum for the DG parametrisation
 3.5 Higherorder corrections for the DG limiting spectrum
 3.6 Other corrections: finite mass, number of active flavours, power terms, and rescaling
 4 Extraction of from the evolution of the distribution of hadrons in jets in collisions
 5 Conclusions and outlook
 A Mellintransformed splitting functions
 B NMLLANLO moments of the distorted Gaussian
 C Higherorder corrections to the moments of the distorted Gaussian
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 highenergy 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 nonperturbative 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 socalled “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 “humpbacked 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 :
(1) 
where is the collinear cutoff parameter of the perturbative expansion which can be pushed
down to the value of (the socalled “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 partontohadron 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 energymomentum fractions and of respectively. The Sudakov parametrisation for and , the fourmomentum of partons and , can be written as
(2) 
with the lightlike vector , and timelike transverse momentum such that, . Then, the scalar product reads:
(3) 
Writing now the 4momenta , , one has, , for onshell and massless partons . From energymomentum conservation:
(4) 
such that, replacing Eq. (4) in (3), one finally obtains:
(5) 
In the collinear limit, one is left with , where is the jet virtuality, or transverse momentum of the jet.
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 softcollinear 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 (nexttoMLLA, or NMLLA) refinements [11, 12].
The MLLA [4], partially restores
the energymomentum balance by including SL corrections of order
coming from the emission of hardcollinear gluons and quarks at large and small
(, and ).
Such corrections are included in the standard DokshitzerGribovLipatovAltarelliParisi
(DGLAP) [13, 14, 15]
splitting functions which describe the parton evolution at intermediate and large in the
(timelike) FFs and (spacelike) 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 socalled 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 halfpowers 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 halfpowers of generated at small , which can be
truncated beyond in the highenergy limit.
In this paper, the set of nexttoMLLA 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 twoloop or nexttoleading 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 quasiGaussian 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
dataset– were done previously for various parametrizations of the input fragmentation
function [41, 42, 43, 44] but only with
a relatively modest datatheory 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 singleinclusive 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 energyevolution for all the moments in the limiting spectrum case () are derived in subsection 3.4, and the role of higherorder 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 finitemass effects and heavyquark 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:
(6) 
where are the regularised DGLAP splitting functions [13, 14, 15], which at LO are given by
(7)  
(8) 
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 distribution^{‡}^{‡}‡The 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 twoloop level reads [19]
(9) 
with
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 socalled scaling violation of DGLAP
FFs in timelike evolution, and that of spacelike evolution of PDFs in
in DIS. In our framework, however, this term is responsible for the double softcollinear 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:
(10) 
where provides the partonenergy 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)
(11) 
and introducing
(12) 
with in the soft approximation (), one is left with the integrodifferential system of evolution equations for the nonsinglet distributions
(13) 
where
(14) 
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
(15) 
The parton density is then obtained through the inverse Mellin transform:
(16) 
where the contour lies to the right of all singularities in the complex plane. In the highenergy 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 .:
(17) 
Thus, replacing Eq. (17) into (13) one obtains
(18) 
which allows for the factorisation of , and leads to the equation
(19) 
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 higherorder 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:
(20a)  
(20b)  
(20c)  
(20d) 
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 higherorder 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]
(21) 
Substituting Eqs. (20a)–(20d) into (21) and performing the expansion again up to terms , yields:
(22a)  
(22b) 
where the terms proportional to are new in this framework. The set of constants involved in Eqs. (22a) and (22b) reads:
(23a)  
(23b)  
(23c)  
(23d) 
Therefore, the diagonalisation of Eq. (19) leads to two equations:
(24) 
such that in the new basis the respective solutions read:
(25a)  
(25b) 
where the ratios in front of are the coefficient functions that will be evaluated hereafter. Notice that in the basis, the offdiagonal terms and vanish for LO splitting functions, while this is no longer true for timelike 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 twoloop 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 twoloop 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 :
(26) 
(27) 
The leading contribution to after setting in Eq. (27) reads:
(28) 
The exponent induces a very small (nonGaussian) 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:
(29) 
with the evolution “Hamiltonian”:
(30) 
that describes the parton jet evolution from its initial virtuality to the lowest possible energy scale , at which the partontohadron 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
(31)  
(32) 
where is the DLA anomalous dimension amounting to
(33) 
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
(34) 
which will be solved after inserting the twoloop 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:
(35)  
(36)  
(37) 
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:
(38) 
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.
As an example of the effects of nonzero skewness and kurtosis, we compare in Fig. 2 the shapes of
four different singleinclusive hadron distributions of width and mean position at in the
interval typical of jets at LEP1 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:
(39) 
the different components (normalised moments) of the DG are given by^{¶}^{¶}¶We list also which is needed to obtain the maximum peak position from , as discussed below.:
(40) 
such that after plugging Eq. (30) into (29) and what results from it into (39), one is left with
(41) 
which is more suitable for analytical calculations since it directly involves the anomalous dimension expression (37).
Multiplicity.
The multiplicity is obtained from the zeroth moment, i.e. the integral, of the singleparticle distribution. Setting in Eq. (37), one obtains