SUSY-like relation of the splitting functions in evolution of gluon and quark jet multiplicities

. We show the new relationship [1] between the anomalous di-mensions, resummed through next-to-next-to-leading-logarithmic order, in the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations for the ﬁrst Mellin moments D q ,g ( µ 2 ) of the quark and gluon fragmentation functions, which correspond to the average hadron multiplicities in jets initiated by quarks and gluons, respectively. So far, such relationships have only been known from supersymmetric (SUSY) QCD. Exploiting available next-to-next-to-next-to-leading-order (NNNLO) information on the ratio D + g ( µ 2 ) / D + q ( µ 2 ) of the dominant plus components, the ﬁt of the world data of D q ,g ( µ 2 ) for charged hadrons measured in e + e − annihilation leads to α (5) s ( M Z ) = 0 . 1205 + 0 . 0016 − 0 . 0020 .

In QCD [2], the inclusive production of single hadrons involves the notion of fragmentation functions D a (x, µ 2 ), where µ is the factorization scale. Owing to the factorization theorem, the D a (x, µ 2 ) functions are universal in the sense that they do not depend on the process by which parton a is produced. By local parton-hadron duality [3], there should be a local correspondence between parton and hadron distributions in hard-scattering processes. So, D a (x, µ 2 ) are genuinely nonperturbative and need to be determined by fitting experimental data. However, once D a (x, µ 2 0 ) are assumed to be known, their µ 2 dependences are governed by the timelike DGLAP evolution equations [4,5], whose splitting functions P ba (x) are known at next-to-next-to-leading order [6]. The scaling violations, i.e., the µ 2 dependences, of D a (x, µ 2 ) may be exploited in global data fits to extract the strong-coupling constant α s = g 2 s /(4π), leading to very competitive results [7] as for the world average [8]. The DGLAP equations are conveniently solved in Mellin space, where D a (N, µ 2 ) = dx x N−1 D a (x, µ 2 ) with N = 1, 2, . . . and similarly for P ba (x): where , with n f being the number of active quark flavors, is the quark singlet component. The quark non-singlet component is irrelevant for the following. After solving the DGLAP equations in Mellin space, one returns to x space via the inverse Mellin transform, analytically continuing N to complex values.
The first Mellin moment D a (µ 2 ) ≡ D a (1, µ 2 ) is of particular interest in its own right because it corresponds to the average hadron multiplicity n h a of jets initiated by parton a. There exists a wealth of experimental data on n h q , n h g , and their ratio r = n h g / n h q for charged hadrons h taken in e + e − annihilation at various center-of-mass energies √ s, ranging from 10 to 209 GeV (for a comprehensive compilation of experimental publications, see Ref. [9]), which allows for a high-precision determination of α s [9,10]. This provides a strong motivation for us to deepen our theoretical understanding of D a within the QCD formalism as much as possible, which is actually limiting the error in the value of α s thus extracted.
The study of D a is a topic of old vintage; the LO value of r, C −1 = C A /C F with color factors C F = 4/3 and C A = 3, was found four decades ago [11]. Subsequent analyses [10,12] were performed using the generating-functional approach in the modified leading-logarithmic approximation (MLLA) [13]. The description of the µ 2 dependences of D a at fixed order in perturbation theory are spoiled by the fact that P ba ≡ P ba (1) are ill defined and require resummation, which was performed for the leading logarithms (LL) [14], the next-to-leading logarithms (NLL) [15], and the next-to-next-to-leading logarithms (NNLL) [16]. In Refs. [9,17], Eq. (1) is first diagonalized for arbitrary value of N at LO, and then the NNLL resummation is incorporated. Unfortunately, this two-step procedure, which has been standard practice in the literature so far [18,19], fails to fully exploit the available knowledge on the higher-order corrections and yields an approximation, the uncertainty of which is difficult to estimate reliably.
In Ref. [1] (see also [20]), we exposed a relationship between the NNLL-resummed expressions for P ba , which has gone unnoticed so far. Its existence in QCD is quite remarkable and interesting in its own right, because a similar relationship is familiar from SUSY QCD, where C = 1 [5,13,16,21].
(2) is written in a form that allows us to glean a novel relationship: which is independent of n f . Eq. (4) generalizes the case of SUSY QCD [5,13,16,21] from C = 1 to arbitrary C values. The corresponding relation in N = 1 SUSY [5] is known to be violated beyond LO [6]. It will be interesting to see if Eq. (4) also holds beyond O(γ 3 0 ), et least in the case of the schemes, which preserve supersymmetry properties, such as the dimensional reduction. The choice of a scheme in above consideration is not so important because a difference in the results of various schemes is exactly canceled in Eq. (4).
Perhaps, the result (4) may be relate with Lipatov observation [22] on an integrability in the high-energy limit of QCD. Of course, the Lipatov observation is based on the resummation of the large ln(1/x) terms in the space-like kinematics. Here we have a similar resummation in the time-like kinematics and the possible relation is not so obvious and should need strong investigations.
We solve Eq. (1) for N = 1 exactly by exploiting Eq. (4). To this end, we diagonalize the NNLL DGLAP evolution kernel as (6) Eq. (1) for N = 1 thus assumes the form where the second term contained within the square brackets stems from the commutator of µ 2 d/dµ 2 and U, and Owing to Eq. (4), the square root in Eq. (6) is exactly canceled, and we have simple expressions for P ±± Inserting Eq. (9) in Eq. (5), we have Using the QCD β function, after a small algebra we may cast Eq. (1) in its final form, The initial conditions are given by Eq. (8) for µ = µ 0 in terms of the three constants α s (µ 2 0 ), D s (µ 2 0 ), and D g (µ 2 0 ). The solution of Eq. (12) is greatly facilitated by the fact that one entry of the matrix on its right-hand side is zero. We may thus obtain D − as the general solution of a homogeneous differential equation, where The correction ∝ γ 0 in Eq. (14) originates from the extra term in Eq. (7) and represents a novel feature of our approach. In Refs. [9,17] and analogous analyses for parton distribution functions [23,24], the minus components do not participate in the resummation.
We are then left with an inhomogeneous differential equation for D + . The general solutioñ D + of its homogeneous part reads where K (1) with d + = −4C A K (1) + /β 0 and b 1 = β 1 /(2C A β 0 ). Adding toD + a special solution of the inhomogeneous differential equation for D + , we find its general solution to be The final expressions for D − and D + in Eqs. (13) and (15), respectively, are fully renormalization group improved because all µ dependence resides in γ 0 . Using Eqs. (5) and (8), we now return to the parton basis, where it is useful to decompose D a = D + a + D − a into the large and small components D ± a proportional to D ± , respectively. Defining r ± = D ± g /D ± s and using Eqs. (2), (3), and (9), we then have D ± s = ∓D ± and . (19) Recalling that n h q = D s and n h g = D g , we thus have Eq. (19) differs from Eqs. (53) and (54) in Ref. [9], On the other hand,r + in Eq. (21) agrees with the result for r obtained in Ref. [25] in the approximation of putting D − a = 0 and extended to through O(γ 3 0 ) in Refs. [26,27], which is in line with the reasoning in Chapter 7 of Ref. [13]. By the same token, we may accommodate the higher-order corrections [26,27] by including within the curly brackets in Eq. (21) the termsc 2 γ 2 0 +c 3 γ 3 0 . The exact results forc 2 andc 3 can be found in [1]. For n f = 5, The difference between r ± andr ± is an artifact of the different diagonalization procedures adopted here and in Ref. [9]. In fact, taking the limit N → 1 in D a (N, µ 2 ) and diagonalizing the DGLAP equations are noncommuting operations. Consequently, our components D a differ from those in Ref. [9], D a withr ± = D ± g /D ± s , by terms of O(γ 0 ). Specifically, we have In fact, this transformation convertsr ± into r ± and allows us to extend our result for r + through O(γ 3 0 ); the counterpart of Eq. (22) reads We denote the approximation of using Eq. (24) on top of Eqs. (19) and (20) as NNNLO approx +NNLL.
Power-like corrections were found to be indispensable for a realistic description of the experimental data of n h q , n h g , and r [27,28]. Following Refs. [27,28], we include them by multiplying r + in Eq. (24) with the factor (to obtain it we collect all terms ∼ γ 0 /µ in the Appendix 1 of [27]) where µ cr is a critical scale parameter to be fitted. In the MLLA approach, µ cr = K cr Λ QCD usually serves as the initial point of the evolution, which is implemented with the basic variables Y = ln(µ/µ 0 ) and λ = ln K cr . The most frequent choice, λ = 0, corresponds to the limitingspectrum approximation [3]. Other recent choices include λ = 1.4 and λ = 2.0 [10]. Since logarithmic and powerlike corrections become comparable at small values of µ 2 , a judicious choice of µ is important to prevent strong correlations. Motivated by Refs. [11,29,30], we choose µ 2 = R 2 Q 2 + 4M 2 eff , where R is the jet radius, Q 2 = √ s, and M eff is the effective gluon mass. We adopt R = 0.3 as a typical value from Ref. [29] and M eff (Q 2 ) = m 2 /[1 + (Q 2 /M 2 ) γ ] with m = 0.375 GeV, M = 0.557 GeV, and γ = 1.06 from Ref. [30].
We are now in a position to perform a global fit to the available measurements of n h q and n h g for changed hadrons h in e + e − annihilation, which were carefully compiled in Ref. [9]. As in Ref. [9], we choose the reference scale to be Q 0 = 50 GeV, which roughly corresponds to the geometric mean of the smallest and largest of the occurring √ s values, and put n f = 5 throughout our analysis. As may be seen in Fig. 1, our NNNLO approx + NNLL fit yields an excellent description of the experimental data included in it, with a χ 2 per degree of freedom of χ 2 dof = 1.32. The fit parameters are determined to be n h (Q 2 0 ) q = 16.38 ± 0.05, n h (Q 2 0 ) g = 23.87 ± 0.07, K cr = 7.09 +1.71 −1.21 , and α (5) s (M 2 Z ) = 0.1205 +0.0026 −0.0037 , where the errors refer to the 90% confidence level (CL) and are evaluated as explained in Ref. [9]. At 68% CL, we have which nicely agrees with the present world average, α (5) s (M 2 Z ) = 0.1181 ± 0.0011 [8]. Our fit results turn out to be very insensitive to the precise choice of Q 0 . The power corrections turn out to be sizeable, with λ = 1.96 +0. 21 −0.19 , in agreement with Ref. [10]. In Fig. 2, we compare our NNNLO approx + NNLL prediction for r with the experimental data compiled in Ref. [9], which did not enter our fit. Someone can see that our results lie something below of the most of points that shows about a disagreements between the experimental data for multiplicities and their ratio r.  Figure 2. Comparison of our NNNLO approx +NNLL prediction of r(µ 2 ) with experimental data excluded from the fit.
In summary, we shown an unexpected, SUSY-like relationship [1] between the NNLLresummed first Mellin moments of the timelike DGLAP splitting functions in real QCD, Eq. (4), which is n f independent. Also incorporating the appropriately transformed O(γ 2 0 ) and O(γ 3 0 ) corrections to r + as well as power-like corrections, we shown also (performed in [1]) a global fit to the world data of charged-hadron multiplicities in quark and gluon jets produced by e + e − annihilation and so extracted the competitive new value of α (5) s (M 2 Z ) in Eq. (26), which nicely agrees with the present world average. Notice that there is only a 1% difference between the result (26) and one in [9], that shows an independence of the results