Energy evolution of the moments of the hadron distribution in QCD jets including NNLL resummation and NLO running-coupling corrections

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 alpha_s strong coupling. The evolution equations are solved using a distorted Gaussian parametrisation, which successfully reproduces the spectrum of charged hadrons of jets measured in e+e- 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 e+e- collisions in the range sqrt(s)~2-200 GeV to the NMLLA+NLO* predictions allow one to extract a value of the QCD parameter Lambda_QCD, and associated two-loop coupling constant at the Z resonance alpha_s(m_Z^2)= 0.1195 +/- 0.0022, in excellent numerical agreement with the current world average obtained using other methods.


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 Λ QCD ≈ 0.2 GeV. The quantitative description of the distribution of hadrons of type h in a jet is encoded in a (dimensionless) fragmentation function (FF) which can be experimentally obtained, e.g. in e + e − collisions at c.m. energy √ s, via D h (ln(1/x), s) = dσ(ee → hX) σ tot d ln(1/x) , where x = 2 p h / √ s is the scaled momentum of hadron h, and σ tot the total e + e − hadronic cross section.
Its integral over x gives the average hadron multiplicity in jets. Writing the FF as a function of the (log of the) inverse of x, ξ = ln(1/x), 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" [1,3,4]. Thus, due to QCD colour coherence and interference of gluon radiation, not the softest partons but those with intermediate energies (E h ∝ E 0.3 jet ) 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 Q 0 is the collinear cut-off parameter of the perturbative expansion which can be pushed down to the value of Λ QCD (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 E jet = 10 GeV (1 TeV).
The measured fragmentation function (1) corresponds to the sum of contributions from the fragmentation D h i of different primary partons i = u, d, · · · , g: and, although one cannot compute from perturbation theory the final parton-to-hadron transition encoded in D h i , the evolution of the "intermediate" functions D bc a describing the branching of a parton of type a into partons of type b,c can be indeed theoretically predicted. The relevant kinematical variables in the parton splitting process are shown in Fig. 1 for the splitting a(k) → b(k 1 ) + c(k 2 ), such that b and c carry the energy-momentum fractions z and (1 − z) of a respectively. The Sudakov parametrisation for k 1 and k 2 , the four-momentum of partons b and c, can be written as with the light-like vector n 2 = 0, and time-like transverse momentum k 2 ⊥ > 0 such that, k ·k ⊥ = n ·k ⊥ = 0. Then, the scalar product k 1 · k 2 reads: Writing now the 4-momenta k = E, k , k 1 = zE, k 1 , k 2 = (1 − z)E, k 2 one has, | k 1 |= zE, | k 2 |= (1 − z)E for on-shell and massless partons k 2 i ≈ 0. From energy-momentum conservation: such that, replacing Eq. (4) in (3), one finally obtains: In the collinear limit, one is left with k ⊥ ≈ z(1 − z)Q, where Q = Eθ is the jet virtuality, or transverse momentum of the jet. The calculation of the evolution of D bc a inside a jet suffers from two types of singularities at each order in the strong coupling α s : collinear ln θ-singularities when the gluon emission angle is very small (θ → 0), and infrared ln(1/z)-singularities when the emitted gluon takes a very small fraction z 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 α s ln k 2 ⊥ /µ 2 n where k ⊥ 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, g → gg and q(q) → gq(q), for small values of x and θ [α s ln(1/z) ln θ] n ∼ O(1) [8,9], (iii) Single Logarithms (SL) [4,10] [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 e + e − 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 O √ α s coming from the emission of hard-collinear gluons and quarks at large x ∼ 1 and small θ i (g → gg, q(q) → gq(q) and g → qq). 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 x 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 D ≃ C(α s (t)) exp t γ(α s (t ′ ))dt , t = ln Q where C(α s (t)) = 1 + √ α s + α s . . . are the coefficient functions, and γ = 1 + √ α s + α s . . . is the so-called anomalous dimension, which in Mellin space at LLA reads, where ω is the energy of the radiated gluon and N c the number of colours. At small ω or x, the expansion of the FF expression leads to a series of half-powers of α s , γ ≃ √ α s + α s + α 3/2 s + . . ., while at larger ω or x in DGLAP, the expansion yields to a series of integer powers of α s , γ ≃ α s + α 2 s + α 3 s + . . . for FFs and PDFs. In the present work we are mostly concerned with series of half-powers of √ α s generated at small ω, which can be truncated beyond O (α s ) in the high-energy limit.
In this paper, the set of next-to-MLLA corrections of order O (α s ) for the single inclusive hadron distribution in jets, which further improve energy conservation [17,18], including in addition the running of the coupling constant α s 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 N and the ratio r = N g /N q 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 Λ QCD 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 Λ QCD . By fitting all the measured e + e − jet distributions in the range of collision energies √ s ≈ 2-200 GeV [6,[23][24][25][26][27][28][29][30][31][32][33][34][34][35][36][37] a value of Λ QCD can be extracted which agrees very well with that obtained from the NLO coupling constant evaluated at the Z resonance, α s (m 2 Z ), in the minimal subtraction (MS) 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 Λ QCD 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 O (α s ) terms from the splitting functions in Mellin space. In subsection 3.1, the new NMLLA+NLO * anomalous dimension, γ NMLLA+NLO * ω , 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 γ NMLLA+NLO * ω . 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 (Q 0 → Λ QCD ) 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 formulae 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 e + e − jet data is carried out in Sect. 4, from which a value of Λ QCD 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.

Evolution equations for the low-x parton fragmentation functions
The fragmentation function of a parton a splitting into partons b and c satisfies the following system of evolution equations [4,5] as a function of the variables defined in Fig. 1: where P ac (z) are the regularised DGLAP splitting functions [13][14][15], which at LO are given by with C F = (N 2 c −1)/2N c and N c respectively the Casimirs of the fundamental and adjoint representation of the QCD colour group SU (3) c , T R = 1/2, and n f is the number of active (anti)quark flavours. The regularisation of the splitting functions in Eq. (6) is performed through the + distribution ‡ in Eqs. (7) and (8). The α s 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 (x = 1) and ‡ The plus distribution applied to a function F (x), written [F (x)]+, is defined as (1)) for any function g(x).
hadronises into a single hadron. The equations (6) are identical to the DGLAP evolution equations but for one detail, the shift in ln z in the second argument of the fragmentation function x z D b c x z , ln z + ln Eθ , that for hard partons is set to zero, ln z ∼ 0, 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 (α s ln 2 ) n , 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 Y provides the parton-energy dependence of the fragmentation process, and the λ specifies, in units of Λ QCD , the value of the hadronisation scale Q 0 down to which the parton shower is evolved. Standard parton showers Monte Carlo codes, such as PYTHIA [45], use Q 0 values of the order of O (1 GeV) whereas in the limiting spectrum [4], that will be used here, it can be taken as low as λ → 0, i.e. Q 0 → Λ QCD . Applying the Mellin transform to the single inclusive distribution in Eq. (6) and introducingξ with k ⊥ ≈ zEθ in the soft approximation (z ≪ 1), one is left with the integro-differential system of evolution equations for the non-singlet distributions where 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 xD(x, Y ) is then obtained through the inverse Mellin transform: where the contour C lies to the right of all singularities in the ω-complex plane. In the high-energy limit (Q ≫ Q 0 ) and hard fragmentation region (Y ≫ξ or x ∼ 1), one can replace in the r.h.s. of Eq. (13) the following expansion § : Thus, replacing Eq. (17) into (13) one obtains which allows for the factorisation of α s (Y )D(ω, Y ), and leads to the equation more suitable for analytical solutions. Truncating the series at higher orders translates into including corrections O (α s ) which better account for energy conservation, particularly at large x. In Mellin space, the expansion can be made in terms of the differential operator Ω ≡ ω + ∂/∂Y such that, up to the second term in Ω, one is left with NMLLA corrections of order O (α s ) [11]. Explicitly, the inclusion of higherorder corrections from the second term of α s (Y −ξ)D(ω, Y −ξ) ≈ α s D −ξ∂(α s D)/∂Y , followed by the integration over the splitting functions (7)-(8) in x space in the r.h.s. of Eq. (13), is equivalent to the expansion P (Ω) = P (0) + P (1) Ω in Mellin space in the r.h.s. of (19), where P (0) and P (1) are constants. The expansion of the matrix elements P (Ω) in Ω can be obtained from the original expressions of the Mellin transformed splitting functions [46], as given in Eqs. (114a)-(114d) in Appendix A, which leads to the following expressions: where the finite terms for Ω → 0 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 δP qq (Ω)Ω = 0. By including those proportional to Ω, one is in addition considering the set of higherorder corrections O (α s ) 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] P ±± (Ω) = 1 2 P gg (Ω) + P qq (Ω) ± (P gg (Ω) − P qq (Ω)) 2 + 4P gq (Ω)P qg (Ω) .
Substituting Eqs. (20a)-(20d) into (21) and performing the expansion again up to terms O (Ω), 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 D ± -basis the respective solutions read: where the ratios in front of D ± are the coefficient functions that will be evaluated hereafter. Notice that in the D ± basis, the off-diagonal terms P +− (Ω) = 0 and P −+ (Ω) = 0 vanish for LO splitting functions, while this is no longer true for time-like splitting functions obtained from the MS factorisation scheme beyond LO [47], as explained in [21] for multiparticle production. Following this logic, D ± 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 *

Anomalous dimension at NMLLA +NLO *
Our NMLLA+NLO * scheme involves adding further corrections O (α s ) 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 D ± result in the following expressions for D + and D − : The leading contribution to D − after setting b 2 = 0 in Eq. (27) reads: The exponent b 1 /(4N c β 0 ) = O 10 −2 √ α s induces a very small (non-Gaussian) correction, which can be neglected asymptotically, for Y + λ ≫ λ. Thus, the (+) trajectory (22a) provides the main contribution to the single inclusive distribution D(ξ, Y ) = xD(x, Y ) at small x ≪ 1, after applying the inverse Mellin transform (16). Hard corrections proportional to a 1 and a 2 account for the energy balance in the hard fragmentation region and are of relative order O √ α s and O (α s ) respectively with respect to the O (1) DLA contribution. The NLO expression (9) results in corrections ∝ β 0 at MLLA, and ∝ β 0 , β 1 at NMLLA which provide a more accurate consideration of running coupling effects at small x ≪ 1 [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 ∝ β 1 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 Q to the lowest possible energy scale Q 0 , at which the parton-to-hadron transition occurs. In Eq. (30), γ(ω, α s (y)) is the anomalous dimension that mixes g → gg and g → qq splittings and is mainly dominated by soft gluon bremsstrahlung (g → gg). Introducing the shorthand notation γ ω = γ(ω, α s (Y )), the MLLA anomalous dimension has been determined in the past [4,22], setting a 2 = 0 and β 1 = 0 in Eq. (26), and is given by where γ 2 0 is the DLA anomalous dimension amounting to The first term of Eq. (31) is the DLA main contribution, of order O( √ α s ), which physically accounts for soft gluon multiplication, the second and third terms are SL corrections O(α s ) accounting for the energy balance (∝ a 1 ) and running coupling effects (∝ β 0 ). 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 , and the third one, proportional to a 2 , is O α 2 s such that after factorising the whole equation by O (α s ) one is left with the relative orders of magnitude in √ α s . 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 ∝ β 1 as well. The equation can be solved iteratively (perturbatively) by setting the MLLA anomalous dimension written in Eq. (31) in the main and subleading contributions of Eq. (33), to find: which is the main theoretical result of this paper. Terms proportional to a 2 1 , a 1 β 0 and β 2 0 are of order O(α 3/2 s ), and were previously calculated in the (N)MLLA+LO scheme described in [42]. Those proportional to β 1 and a 2 are computed for the first time in our NMLLA+NLO* framework. Indeed, the single correction ∝ β 1 is obtained replacing Eq. (9) in the l.h.s. of (33), which leads to the equation,

Distorted Gaussian (DG) parametrisation for the fragmentation function
The distorted Gaussian (DG) parametrisation of the single inclusive distribution of hadrons in jets at small x (or ω → 0) was introduced by Fong and Webber in 1991 [22], and in x-space it reads: where, δ = (ξ −ξ)/σ, N is the asymptotic average multiplicity inside a jet, andξ, σ, s, and k are respectively the mean peak position, the dispersion, the skewness, and kurtosis of the distribution. The distribution should be displayed in the interval 0 ≤ ξ ≤ Y which depends on the jet energy, and the values of Q 0 and Λ QCD . The three scales of the process are organised in the form Q ≫ Q 0 ≥ Λ QCD . The formula (35) reduces to a Gaussian for s = k = 0 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 non-zero skewness and kurtosis, we compare in Fig. 2 the shapes of four different single-inclusive hadron distributions of width σ = 1.4 and mean position atξ = 3.5 in the interval 0 ≤ ξ 7 typical of jets at LEP-1 energies: (i) an exact Gaussian, (ii) a skewed Gaussian with s = −0.5, k = 0, (iii) a kurtic Gaussian with s = 0, k = −0.5, and (iv) a DG including both "distorting" s, k 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 NMLLA+NLO * anomalous dimension, Eq. (34), computed here. Defining K n as the n-th moment of the single inclusive distribution: the different components (normalised moments) of the DG are given by ¶ : such that after plugging Eq. (30) into (29) and what results from it into (36), one is left with which is more suitable for analytical calculations since it directly involves the anomalous dimension expression (34).
Multiplicity. The multiplicity is obtained from the zeroth moment, i.e. the integral, of the single-particle distribution. Setting ω = 0 in Eq. (34), one obtains We list also k5 which is needed to obtain the maximum peak position ξmax fromξ, as discussed below.
from which the mean multiplicity N (Y, λ) can be straightforwardly derived by integrating over y: where As expected, the mean multiplicity (40) including the two-loop α s exactly coincides with the expression obtained in [20]. This cross-check supports the validity of our "master" NMLLA+NLO * formula (34) for the anomalous dimension at small ω, which is not surprising as the gluon jet evolution equation solved in [20] for the mean multiplicity coincides with Eq. (26) after setting ω = 0 and N (Y, λ) = D + (0, Y, λ). The first term in Eq. (41) is the DLA rate of multiparticle production, the second and third terms provide negative corrections that account for energy conservation and decrease the multiplicity. However, the third term, proportional to β 1 , is positive and can be large since it accounts for NLO coupling corrections. Though, due to energy conservation, one may expect the multiplicity to decrease in the present scheme running coupling effects take over and can drastically increase the multiplicity as well as single inclusive cross-sections at the energy scales probed so-far at e + e − colliders. Only at asymptotically high-energy scales, that is for Q 0 ≫ Λ QCD , the energy conservation becomes dominant over running coupling effects, thus inverting these trends. The ratio of multiplicities in quark and gluon jets are discussed in Sect. 3.3 and compared with the calculations of [20]. Performing the numerical evaluation for n f = 5 quark flavours we obtain the final expression for the multiplicity: Peak position. The energy evolution dependence of the mean peak position is obtained plugging Eq. (34) into (30), and the latter into Eq. (29) in order to get the K n moments of the distribution from Eq. (36). Thus, for n = 1 one obtains The smallness of the constant in front of the NMLLA correction proportional to (ln(Y + λ)− ln λ) should not drastically modify the MLLA peak position and should only affect it at small energy scales.
The position of the mean peak is related to the corresponding maximum and median values of the DG distribution by the expressions [48]: As will be seen below the dependence in n f is very weak and will not affect the final normalisation of the distribution.
for which we need the fifth moment of the DG, k 5 , which reads: The final numerical expressions for the mean and maximum peak positions, evaluated for n f = 5 quark flavours, read: Width. The DG distribution dispersion σ follows from its definition in Eq. (38) for n = 2. The full expression for the second moment K 2 (Y, λ) can be found in Appendix B, Eq. (118), from which taking the squared root, followed by the Taylor expansion in (1/ √ y + λ or √ α s ) and keeping trace of all terms in (1/(y + λ) or α s ), the NMLLA+NLO * expression for the width is obtained: where the functions f i are also defined in Appendix B. The new correction term, proportional to (1/(Y + λ)), is of order O (α s ) and decreases the width of the distribution and so does λ for the truncated cascade with Q 0 > Λ QCD . The numerical expression for the width (for n f = 5 quark flavours) reads: Skewness. The NMLLA term of the third DG moment, K 3 , turns out to vanish like the leading order one [48]. According to the definition in Eq. (38), the skewness s = K 3 σ −3 presents an extra subleading term which in this resummation scheme comes from the expansion of the second contribution to σ −3 , proportional to 1/ (Y + λ), as written in Eq. (122) of Appendix B, such that In [22], only the first term of this expression was provided, the subleading contribution given here is thus new. This NMLLA+NLO * correction to Eq. (50) increases the skewness of the distribution, while for increasing λ it should decrease again, thus revealing two competing effects. The net result is a displacement of the tails of the HBP distribution downwards to the left and upwards to the right from the peak position and depending on the sign given by both effects (Fig. 2). The final numerical expression for the skewness (for n f = 5 quark flavours) reads: Kurtosis. The evolution of the kurtosis follows from the expressions for the fourth DG moment, given in Eqs. (120) and (123) of Appendix B. As shown in the same appendix, the proper Taylor expansion in powers of (1/ √ Y + λ) which keeps trace of higher-order corrections and leads to: where the functions f i can be again found in Appendix B. The new NMLLA+NLO * correction for the kurtosis affects the distribution by making it smoother in the tails and wider in the hump region. The final numerical expression for the kurtosis (for n f = 5 quark flavours) reads: for Q 0 ≫ Λ QCD . Physically, for higher values of the shower energy cut-off Q 0 , the strength of the coupling constant decreases and the probability for the emission of soft gluon bremsstrahlung decreases accordingly, making the multiplicity distribution and the peak position smaller. The difference between the MLLA and NMLLA+NLO * resummed distributions is, as mentioned above, mainly due to runningcoupling effects, proportional to β 1 , at large ξ (small x) which is not unexpected because in this region they are more pronounced due to the ln(xEθ) dependence in the denominator of the strong coupling. On the other hand, energy conservation plays a more important role in the hard fragmentation region x ∼ 1 (ξ ∼ 0), where the NMLLA+NLO * DG is somewhat suppressed compared with the MLLA DG.

Multiplicities for the single inclusive D g and D q distributions
In this section we determine the coefficient function involved in Eq. (25a) that provide higher-order corrections to the quark/gluon multiplicity ratio. As shown through Eq. (28), the D − (ω, λ) component is negligible and thus the solutions for the gluon and quark single inclusive distributions can be directly obtained from D + in the form Making use of the expressions (20a)-(20d) and (22a)-(22b), and expanding in ω results in where the numerical values of the constants, for n f = 5 quark flavours, read The c i numerical constants in Eq. (55) were obtained in [4]. Performing the inverse Mellin-transform back to the x-space, or making the equivalent replacement Ω → ∂ ∂ξ + ∂ ∂Y , one has which in a more compact form can be rewritten as A clear difference is observed in the quark and gluon jet initiated distributions given by the colour factor C F /N c = 4/9 and the role of higher-order corrections which prove more sizable for the NMLLA+NLO * scheme over the whole phase space 0 ≤ ξ ≤ Y , as observed in the right panel of Fig. 4. In [4] however, the role of O ( √ α s ) corrections, proportional to c A , which allowed for a direct comparison between the MLLA D + (ξ, Y ) and the hadronic energy-momentum spectrum (for a complete review see [10]). Asymptotically (Q → ∞), the solution of the original Eq. (61) has a Gaussian shape near its maximum: where, as a result of the expansion, Notice that up to the order O (α s ), the multiplicity ratio does not involve corrections proportional to β 1 , which only appear beyond this level of accuracy [20]. Up to the NMLLA order in O (α s ), Eq. (63) coincides with the expression found in [49], which gives further support to the calculations carried out in our work. A more updated evaluation of the mean multiplicity ratio, including two-loop splitting functions, was given recently in [21].

Limiting spectrum for the DG parametrisation
The so-called limiting spectrum, λ → 0, implies pushing the validity of the partonic evolution equations down to (non-perturbative) hadronisation scales, Q 0 ≈ Λ QCD [1]. Such an approach provides a minimal (and successful) approach with predictive power for the measured experimental distributions. We derive here the evolution of the distorted Gaussian moments for this limit which involves formulae depending only on Λ QCD as a single parameter.

Multiplicity.
Among the various moments of the DG parametrisation, only its integral (representing the total hadron multiplicity) needs an extra free parameter to fit the data. The "local parton hadron duality" (LPHD) hypothesis is a powerful assumption which states that the distribution of partons in inclusive processes is identical to that of the final hadrons, up to an overall normalization factor, i.e. that the mean multiplicity of the measured charged hadrons is proportional to the partonic one through a constant K ch , Thus, in the limiting spectrum the mean multiplicity reads which is in agreement with the mean multiplicity first found in [20], supported by the improved solution of the evolution equations accounting for the same set of corrections.
Peak position. For the limiting spectrum, the mean peak position Eq. (43) can be approximated as follows:ξ thanks to the fortuitous smallness O(10 −3 ) of the NMLLA correction toξ at high-energy where Y + λ ≫ λ. Notice that, as shown in [22], the MLLA version of Eq. (67) up to the second order is finite. The origin of the third ∝ ln Y correction in this resummation framework comes from the truncated expansion of the anomalous dimension Eq. (34) in O(α s ), which is proportional to 1/Y by making (−∂γ ω /∂ω) at ω = 0, and hence yields the ∝ ln Y term after integrating over Y . Therefore, we assume that Eq. (67) is valid for Q ≫ Q 0 ≈ Λ QCD .
The maximum of the peak position for the limiting spectrum DG can be obtained via Eq. (44) which involves the mean peak position as well as the other higher-order moments. In a generic form, the moments of the distorted Gaussian associated with the dispersion (48), skewness (50), kurtosis (52), and k 5 (45), are finite for n ≥ 2 for the limiting spectrum and can be written as where the constants K n and the functions f i (λ → 0) → 1 are written in Appendix B. In other words, the second λ-dependent part of K n in Eq. (68) can be dropped as λ → 0 for sufficiently high energy scales, Y + λ ≫ λ, where α s (Y + λ) ≪ α s (λ) in the r.h.s. of Eq. (68). Performing the same approximation in Eq. (68) as λ → 0, the expressions for the rest of moments of the fragmentation functions in the limiting spectrum are derived below. Thus inserting Eqs. (70a), (70b), (70c) and (70d) into (44), we obtain : Width. The width of the DG distribution in the limiting spectrum is obtained from Eq. (48): Skewness. The skewness of the DG distribution in the limiting spectrum reads, from Eq. (50), Kurtosis. The kurtosis can be derived from Eq. (52): Accordingly, we give the last component, k 5 , following from Eq. (45): Final DG (limiting spectrum) expression. In order to get the DG in the limiting spectrum, one should replace Eqs. (66)-(70c) into Eq. (35). We note that in our NMLLA+NLO * framework, the K ch from the DG can be smaller than that found in [20] since it should fix the right normalisation enhanced by second-loop coupling constant effects. Notice also that setting subleading corrections to zero, we recover the results from [22] as expected. In Fig. 5, the MLLA and NMLLA+NLO * distorted Gaussians are displayed in the limiting spectrum approximation for a jet virtuality Q = 350 GeV in the interval 0 ≤ ξ ≤ Y , for Y = 7.5.
We can see a sizable difference between the MLLA D + (ξ, Y ) and the NMLLA+NLO * D + (ξ, Y ) evolutions, which is mainly driven by the two-loop ∝ β 1 correction in the mean multiplicity and other moments of the DG, as mentioned above. The account of energy conservation can be observed at low ξ, i.e. for harder partons. Similar effects have been discussed in [50] where an exact numerical solution of the MLLA evolution equations was provided with one-loop coupling constant. Numerical solutions of exact MLLA equations provide a perfect account of energy conservation at every splitting vertex of the branching process in the shower. For this reason, accounting for higher-order corrections O(α n/2 s ) to the truncated series of the single inclusive spectrum of hadrons should follow similar features and trends to that provided by the numerical solutions of [50] (see also [51]), although our NMLLA+NLO * solution incorporates in addition the two-loop coupling constant.
In Fig. 6 we display the same set of curves as in the Fig. 4 with the right normalisation given by the coefficient functions for quark and gluon jets. The overall corrections provided by the coefficient functions slightly decrease the normalisation of the spectrum in a gluon jet as well as its width σ. In the quark jet, upon normalisation by the colour factor C F /N c , the normalisation is decreased while the width is slightly enlarged. In order to better visualise the less trivial enlargement for the width, we can for instance consider e + e − -annihilation into hadrons at the LEP-2 centre of mass energy √ s = 196 GeV for a quark jet of virtuality Q = √ s/2 = 98 GeV with Y = ln( √ s/(2Λ QCD )) ≈ 6.0 for Λ QCD = 0.25 GeV. If the resulting distribution D q (ξ, Y ) is refitted to a DG and compared with the D + (ξ, Y ), the enlargement of the width compared with that given by D + (70a) can reach 10%. This latter effect is mainly due to the positive O (α s ) correction to the coefficient function C g q given by the larger numerical coefficient c (0) q = 0.487. Similar effects have been discussed in [50]. In conclusion, we will directly fit the D + (ξ, Y ) distribution to the data of final state hadrons in the limiting spectrum approximation.

Higher-order corrections for the DG limiting spectrum
The exact solution of the MLLA evolution equations with one-loop coupling constant entangles corrections which go beyond O ( √ α s ), though the equations are originally obtained in this approximation [5].
The exact solution resums fast convergent Bessel series in the limiting spectrum λ → 0. Using the DG parametrisation it is possible to match the exact solution in the vicinity of the peak position δ ≪ 1 after determining the DG moments: ξ 1 =ξ, ξ 2 = ξ 2 , ξ 3 = ξ 3 , ξ 4 = ξ 4 , related to the dispersion, skewness and kurtosis through [52]: where ξ n is determined via discussed in more detail in Appendix C. Similarly, these extra corrections, which better account for energy conservation and provide an improved description of the shape of the inclusive hadron distribution in jets, will be computed and added hereafter to all the NMLLA+NLO * DG moments, as it was done in [4] for the particular case of the mean peak position,ξ, but extended here also to all other components: Eqs. (67), (70a), (70b) and (70c).
Multiplicity. The extra "hidden" corrections discussed in Appendix C result in one extra term for the multiplicity in the DG limiting spectrum, which is inversely proportional to Y and amounts to: However, we can use directly the full-NLO result obtained in [20] for the multiplicity. In this case the extra correction amounts to: , for n f = 3, and (76) although the terms ∝ 1 √ Y and ∝ 1 Y are almost constant and practically compensate to each other at the currently accessible energies. Peak position. The mean peak value of the DG distribution,ξ, truncated as done in Eq. (43) can be improved as discussed in [4]. The NMLLA correction proportional to ln Y is of relative order O( √ α s ) and is very small O(10 −3 ln Y ) compared to the second term. There is one extra correction (numerical constant) toξ coming from the exact solution of Eq. (26) with a 2 = 0, written in terms of Bessel series in Appendix C. Indeed, substituting Eq. (135) into (133) (see Appendix C for a complete derivation), one obtains the extra NMLLA term toξ: from the expansion of the Bessel series through the Eq. (133) that should be added to Eq. (43). Therefore, the full resummed expression of the mean peak position reads in its complete NMLLA+NLO * form. The corresponding position of the maximum is related to the mean peak value by the expression [48]: such that Width. Similar extra corrections can be found for the dispersion by calculating ξ 2 through this recursive procedure. By making use of Eq. (74) and the full derivation presented in Appendix C, it was found in [52]: such that, with σ 2 = ξ 2 −ξ 2 given by Eq. (71), one finds the extra correction (for n f = 5) which should be accordingly added to the r.h.s. of Eq. (70a).

Skewness.
In the case of the skewness, the expression for ξ 3 reads ξ 3 such that, if one makes use of the expression (72), the extra correction reads (for n f = 5) to be added to the r.h.s. of Eq. (70b). Notice that Eq. (84) was given in [52] without accounting for terms O z −4 and O z −7 . Such terms cannot be neglected when dealing with MLLA and NMLLA corrections.
Kurtosis. Finally, for the kurtosis, we obtain the formula for ξ 4 : which can be cast into Eq. (73) to obtain the corresponding correction which reads (for n f = 5): to be also added to Eq. (70c).
For n f = 4 quark flavours, relevant for jet analysis above the charm mass threshold (m c ≈ 1.3 GeV) but below the bottom mass, one finds and for n f = 5 quark flavours relevant for jet analysis above the bottom mass threshold (m b ≈ 4.2 GeV): The MLLA expressions first computed in [22] can be naturally recovered from our results by keeping all terms up to 1/ √ Y . For n f = 5 quark flavours, they read: which clearly highlight, by comparing to the corresponding full expressions above, the new NMLLA+NLO * terms computed in this work for the first time.

Other corrections: finite mass, number of active flavours, power terms, and Λ QCD rescaling
Mass effects: In the approach discussed so far, the partons have been assumed massless and so their scaled energy and momentum spectra are identical. Experimentally, the scaled momentum distribution ξ p = ln( √ s/(2 p h )) is measured and, since the final-state hadrons are massive, the equivalence of the theoretical and experimental spectra no longer exactly holds. One can relate the measured ξ p spectrum to the expected DG distribution (which depends on ξ ≡ ξ E ) by performing the following change of variables [53]: where the energy of a hadron with measured momentum p h = ( √ s/2) · exp −ξ p is E h = p 2 h + m eff 2 , and m eff is an effective mass of O(Λ QCD ) accounting for the typical mixture of pion, kaon and protons in a jet. In Fig. (7) we compare the DG distribution in the limiting-spectrum for the typical HBP of LEP-1 jets with and without mass corrections, using Eq. (112) with m eff = 0 and m eff = Λ QCD ≈ 0.23 GeV. As expected, the net effect of the non-null mass of the measured jet particles affects the tail of the distribution at high ξ (i.e. at very low momenta) but leaves otherwise relatively unaffected the rest of the distribution. In the analysis of experimental jet data in the next Section, the rescaling given by Eq. (112) will be applied to the theoretical DG distribution for values of m eff = 0-0.35 GeV to gauge the sensitivity of our results to finite-mass effects. Since experimentally there are not many measurements in the large ξ tail (i.e. very low particle momenta) and here the distribution has larger uncertainties than in other ranges of the spectrum, the fits to the data turn out to be rather insensitive to m eff .

Number of active flavours n f :
The available experimental e + e − data covers a range of jet energies E jet ≈ 1-100 GeV which, in its lowest range, crosses the charm (m c ≈ 1.3 GeV) and bottom (m b ≈ 4.2 GeV) thresholds in the counting of the number of active quark flavours n f present in the formulae for the energy-dependence of the DG moments. Although the differences are small, rather than trying to interpolate the expressions for different values of n f in the heavy-quark crossing regions, in what follows we will use the formulaae for n f = 5 for the evolution of all moments and rescale the obtained moments of the four lower-√ s datasets from the BES experiment [23] to account for their lower effective value of n f . The actual numerical differences between the evolutions of the DG moments for n f = 4 and n f = 5 quark flavours -given by Eqs. (94)-(99) and (100)-(105) respectively -when evaluated for energies below the bottom-quark threshold are quite small: 0-10% for N (Y ), below 1% for ξ max (Y ), around 5% for the width σ(Y ), and 5-10% for the skewness s(Y ) and kurtosis k(Y ). In this respect, the most "robust" (n f -insensitive) observable is the peak position of the distribution.
Power-suppressed terms: Power corrections of order O(Q n 0 /Q n ) appear if one sets more accurate integration bounds of the integro-differential evolution equations over z, such as Q 0 Q ≤ z ≤ 1− Q 0 Q instead of 0 ≤ z ≤ 1, which actually leads to Eq. (26) after Mellin transformation with Q 0 ∼ m h , where m h is the hadron mass (for more details see review [54,55]). For the mean multiplicity, this type of corrections was considered in [17]. They were proved to be powered-suppressed and to provide small corrections at high-energy scales. Furthermore, they become even more suppressed in the limiting spectrum case where Q 0 can be extended down to Λ QCD for infrared-safe observables like the hump-backed plateau. The MLLA computation of power corrections for differential observables is a numerical cumbersome task which, for the hump-backed plateau, would add minor improvements in the very small x domain ln(1/x) → ln(Q/Λ QCD ) away from the hump region of our interest, and thus they would not introduce any significant shift to the main moments of the hadron distributions (in particular its peak position ξ max , and width σ).
Rescaling of the Λ QCD parameter: Technically, the Λ QCD parameter is a scheme-dependent integration constant of the QCD β-function. Rescaling the QCD parameter by a constant, Λ QCD → CΛ QCD , would give an equally acceptable definition. In our formalism, such a variation would translate into a ln C-shift of the constant term of the HBP peak, Eq. (81) [4], which corresponds to higher-order contributions to the solution of the evolution equations. The approach adopted here is to connect Λ QCD to α s in the MS factorisation scheme through the two-loop Eq. (9) and, at this level of NLO accuracy, there is no ambiguity when comparing our extracted α s results to other values obtained using the same definition.
where D + (ξ, Y ) is given by Eq. (112) corrected to take into account the finite-mass effects of the hadrons (for values of m eff = 0-0. 35 GeV, see below) with Y = ln[ √ s/(2 Λ QCD )]. Each fit has five free parameters for the DG: maximum peak position, total multiplicity, width, skewness and kurtosis. In total, we analyse 32 data-sets from the following experiments: BES at √ s = 2-5 GeV [23]; TASSO at √ s = 14-44 GeV [24,25]; TPC at √ s = 29 GeV [26]; TOPAZ at √ s = 58 GeV [27]; ALEPH [28], L3 [29] and OPAL [6,30] at √ s = 91.2 GeV; ALEPH [31,34], DELPHI [32] and OPAL [33] at √ s = 133 GeV; and ALEPH [34] and OPAL [35][36][37] in the range √ s = 161-202 GeV. The total number of points is 1019 and the systematic and statistical uncertainties of the spectra are added in quadrature.  In order to assess the effect of finite-mass corrections discussed in the previous Section, we carry out the DG fits of the data to Eq. (112) for many values of m eff in the range 0-320 MeV. The lower value assumes that hadron and parton spectra are identical, the upper choice corresponds to an average of the pion, kaon and (anti)proton masses weighted by their corresponding abundances (65%, 30% and 5% approximately) in e + e − collisions. Representative fits of all the single-inclusive hadron distributions for m eff = 0, 140, and 320 MeV are shown in Figures 8-10 respectively, with the norm, peak, width, skewness, and kurtosis as free parameters. In all cases the individual data-model agreement is very good, with goodness-of-fit per degree-of-freedom χ 2 /ndf ≈ 0.5-2.0, as indicated in the data/fit ratios around unity in the bottom panels. The fits to all datasets with energies above √ s = 50 GeV turn out to be completely insensitive to the choice of m eff , i.e. the moments of the DG obtained are "invariant" with respect to the value of m eff , whereas those at lower energies are more sensitive to it. The value of the effective mass that provides an overall best agreement to the whole set of experimental distributions is m eff ≈ 140 MeV, which is consistent with a dominant pion composition of the inclusive charged hadron spectra.  The general trends of the DG moments are already visible in these plots: as √ s increases, the peak of the distribution shifts to larger values of ξ (i.e. smaller relative values of the charged-hadron momenta) and the spectrum broadens (i.e. its width σ increases). In the range of the current measurements, the peak moves from ξ max ≈ 1 to ξ max ≈ 4, and the width increases from σ ≈ 0.5 to 1.2. The expected logarithmic-like energy dependence of the peak of the ξ distribution, given by Eq. (102), due to soft gluon coherence (angular ordering), correctly reproduces the suppression of hadron production at small x seen in the data to the right of the distorted Gaussian peak. Although a decrease at large ξ (very small x) is expected based on purely kinematic arguments, the peak position would vary twice as rapidly with the energy in such a case in contradiction with the calculations and data. The integral of the ξ distribution gives the total charged-hadron multiplicity N ch which increases exponentially as per Eq. (100).    its quoted uncertainties. Varying m eff from zero to 0.32 GeV yields differences in the extracted Λ QCD parameter below ±0.5% for the ξ max fits and below ±2% for the other components, which indicate the robustness of our NMLLA+NLO * calculations for the limiting-spectrum DG with respect to finite-mass effects if a wide enough range of charged-hadron and parent-parton (jet) energies are considered in the evolution fit. The point-to-point uncertainties of the different moments, originally coming from the DG fit procedure alone, have been enlarged so that their minimum values are at least 3% for the peak position, and 5% for the multiplicity and width. Such minimum uncertainties are consistent with the spread of the DG moments obtained for different experiments at the same collision-energies, and guarantee an acceptable global goodness-of-fit χ 2 /ndf ≈ 1 for their √ s-dependence. We note that not all measurements originally corrected for feed-down contributions from weak decays of primary particles. This affects, in particular, the multiplicities measured for the TASSO [24,25], TPC [26] and OPAL [6] datasets which include charged particles from K s 0 and Λ decays. The effect on the peak position (and higher HBP moments) of including secondary particles from decays is negligible (<0.5%), but increases the total charged particles yields by 8% according to experimental data and Monte Carlo simulations [45]. For these three data-sets, we have thus reduced accordingly the value of N ch . The DG skewness and kurtosis are less well constrained by the individual fits to the measured fragmentation functions and have much larger uncertainties than the rest of moments. As a matter of fact, in the case of the kurtosis our NMLLA+NLO * prediction for its energy-evolution Eq. (105), fails to provide a proper description of the data and seems to be above the data by a constant offset (Fig. 15). Whether this fact is due to missing higher-order contributions in our calculations or to other effects is not yet clear at this point. Apart from the kurtosis, the QCD coupling value extracted from all the other moments has values around α s (m 2 Z ) = 0.118, in striking agreement with the current world-average obtained by other methods [56,57]. Table 1: Values of Λ QCD and associated α s (m 2 Z ) at NLO (MS scheme, n f = 5 quark flavours) obtained from the fits of the √ s-dependence of the moments of the charged hadron distribution of jets in e + e − collisions obtained from their NMLLA+NLO * evolution. The last column provides the weighted-average of the individual measurements with its total propagated uncertainty.  Table 1 lists each value of the Λ QCD parameter individually extracted from the energy evolutions of the four DG components that are well described by our NMLLA+NLO * approach, and their associated values of α s (m 2 Z ) obtained using the two-loop Eq. (9) for n f = 5 quark flavours. Whereas the errors quoted for the different Λ QCD values include only uncertainties from the fit procedure, the propagated α s (m 2 Z ) uncertainties have been enlarged by a common factor such that their final weighted average has a χ 2 /ndf close to unity. Such a "χ 2 averaging" method [57] takes into account in a well defined manner any correlations between the four extractions of α s , as well as underestimated systematic uncertainties. The relative uncertainty of the α s (m 2 Z ) determination from the DG moments evolution is about ±1.5% for the maximum peak position, ±3.5% for the width, ±7% for the total multiplicity, and about ±11% for the skewness. The last column of Table 1 lists the final values of Λ QCD and α s (m 2 Z ) determined by taking the weightedaverage of the four individual measurements. We obtain a final value α s (m 2 Z ) = 0.1195 ± 0.0022 which is in excellent agreement with the current world-average of the strong coupling at the Z mass [56,57]. Our extraction of the QCD strong coupling has an uncertainty (±2%) that is commensurate with that from other e + e − observables such as jet-shapes (±1%) and 3-jets rates (±2%) [56,57]. In a forthcoming work, we extend the extraction of the strong coupling via the NMLLA+NLO * evolution of the moments of the hadron distribution in jet world-data measured not only in e + e − but also including deep-inelastic e ± p collisions [58].

Conclusions and outlook
We have computed analytically the energy evolution of the moments of the single-inclusive distribution of hadrons inside QCD jets in the next-to-modified-leading-log approximation (NMLLA) including next-toleading-order (NLO) corrections to the α s strong coupling. Using a distorted Gaussian parametrization, we provide in a closed-form the numerical expressions for the energy-dependence of the maximum peak position, total multiplicity, peak width, kurtosis and skewness of the limiting spectra where the hadron distributions are evolved down to the Λ QCD scale. Comparisons of all the existing jet data measured in e + e − collisions in the range √ s ≈ 2-200 GeV to the NMLLA+NLO * predictions for the moments of the hadron distributions allow one to extract a value of the QCD parameter Λ QCD and associated two-loop coupling constant at the Z resonance, α s (m 2 Z ) = 0.1195 ± 0.0022, in excellent agreement with the current world average obtained with other methods. The NMLLA+NLO * approach presented here can be further extended to full NMLLA+NLO through the inclusion of the two-loop splitting functions. Also, in a forthcoming phenomenological study we plan to compare our approach not only to the world e + e − jet data but also to jet measurements in (the current hemisphere of the Breit-frame of) deep-inelastic e ± p collisions. The application of our approach to the hadron distribution of TeV-jets produced in protonproton collisions at LHC energies would further allow one to extract α s from parton-to-hadron FFs over a very wide kinematic range. The methodology presented here provides a new independent approach for the determination of the QCD coupling constant complementary to other existing jet-based methods -relyiong on jet shapes, and/or on ratios of N-jet production cross sections-with a totally different set of experimental and theoretical uncertainties.

B NMLLA+NLO * moments K n of the distorted Gaussian
We compute here the generic for the moments of the distorted Gaussian (DG) for λ = 0 according to Eq. (38) by introducing the following functions: where the function L n was written in the form of the series, L n (B + 1, B + 2; z) = P