NLO QCD corrections to five-jet production at LEP and the extraction of alpha_s(M_Z)

The highest exclusive jet multiplicity studied at LEP experiments is five. In this paper we compute the next-to-leading order QCD corrections to e+e- annihilation to five jets, essentially closing the (pure) perturbative QCD studies of exclusive jetty final states at LEP. We compare fixed-order perturbative results with ALEPH data. We estimate hadronization corrections to five-jet observables using the event generator SHERPA, which employs the CKKW procedure to combine a reliable perturbative treatment of high-multiplicity jet final states with parton showers. We show that a competitive value of the strong coupling constant alpha_s(M_Z) = 0.1156 +0.0041 -0.0034 can be extracted from the distribution of the five-jet resolution parameter and the five-jet rate at LEP1 and LEP2.

Infra-red safe observables, traditionally studied in e + e − annihilation, are dominated by short-distance physics; they are computable in perturbative QCD up to corrections suppressed by inverse powers of some large energy scale related to e.g. the average relative transverse momentum of jets in a given process. Those power (or hadronization) corrections are typically estimated using event generators such as PYTHIA [39], HERWIG [40], and ARIADNE [41], under the assumption of a complete factorization of non-perturbative and perturbative physics. This factorization implies that the hadronization corrections to an infrared-safe observable O can be estimated as where O i hadr and O i part are the values of the observable O computed at the hadron and at the parton level with the event generator i. Because hadronization corrections are assumed to be factorizable, they can be used to "improve" the perturbative prediction O pt for the observable under study. Hence, one defines the quantity O impr = H i [O] O pt , and compares it to experimental data 1 .
While the use of the procedure that we just described is widespread, it is clear that it can not be fully valid. Indeed, one can imagine that, for a particular observable O, the event generator i happens to reproduce its measured value, O i hadr = O Data . As a result, O impr reads (1.2) Clearly, O impr can only be equal to O data if O pt = O i part , but this equality can not hold true for a variety of reasons, including different approximations in parton/dipole showers and fixed-order perturbative calculations, different dependencies of O pt and O i part on the renormalization scale µ, etc. Differences at the perturbative level are particularly worrisome when considering high-multiplicity final states, since standard event generators routinely used in e + e − studies are based solely on e + e − → qq and e + e − → qqg matrix elements. This implies that, since event generators are tuned to data, hadronization corrections as defined by eq. (1.1) contain both non-perturbative and perturbative effects; the latter are present to compensate for deficiencies of the partonic part of a particular event generator 2 . To the extent that both perturbative and hadronization corrections are small, the inconsistency of the whole procedure may not be very apparent, but it becomes evident if those corrections are large. It is very likely that these issues are important for five-jet production at LEP. Indeed, because e + e − → 5 jets involves a high power of α s , NLO QCD corrections are expected to be large. In addition, a correct description of five hard, wellseparated partons is difficult for traditional event generators, so that all the problems of a conventional approach to estimating hadronization corrections can be exposed by studying five-jet observables. 1 There are alternative approaches to estimate hadronization corrections that address both theoretical [42] and experimental [16,43] aspects of this procedure. 2 See ref. [44] for a related discussion.
By extracting hadronization corrections with traditional event generators -PYTHIA, HERWIG and ARIADNE -we find these corrections to be large and generator-dependent. This is unfortunate since it implies a large spread of "improved" predictions for five-jet observables when perturbative and hadronization effects are combined. For this reason, we believe it is important to obtain hadronization corrections from an event generator whose perturbative part is up to the task of describing the production of five hard, wellseparated partons. Within the context of event generation, the implementation of such a description requires a consistent matching between high-multiplicity matrix elements, and parton/dipole showers. In this paper, we employ the SHERPA event generator [45,46], which implements the CKKW [47] matching prescription 3 . By calculating hadronization corrections with SHERPA/CKKW, we find that they are relatively small (see the right pane of fig. 2), in particular in the kinematic region where perturbative QCD is reliable. This is what we expect since, with CKKW matching, the perturbative description of five-jet production provides a good approximation to the actual physical process. As a consequence, when we use SHERPA/CKKW to extract non-perturbative corrections according to eq. (1.1), the results are less contaminated by perturbative contributions, compared to the case when traditional event generators are used for this purpose. Since these considerations apply to hadron collisions as well, our results have obvious implications for jet physics at the LHC.
The goal of this paper is to investigate five-jet production at LEP, using theoretical predictions that are accurate at NLO in QCD. We discuss hadronization corrections and show that they depend significantly upon the event generator that is used to estimate them. We extract the strong coupling constant by fitting the NLO QCD predictions for the distributions of the five-jet resolution parameter y 45 and the five-jet rate R 5 . The remainder of the paper is organized as follows. In Section 2 we discuss the technical details pertinent to the computation of the NLO QCD corrections to e + e − → 5 jets. In Section 3 a phenomenological analysis of five-jet production at LEP is reported. In Section 4 the value of the strong coupling constant is extracted. In Section 5 we present our conclusions. Details of the fit procedure that we use in our analysis of the strong coupling constant are described in the Appendix.

Technical details
The computation of NLO QCD corrections to any process or observable in the context of a subtraction formalism involves the evaluation of the following contributions: the oneloop virtual corrections, the real-emission corrections and their subtraction terms, and the finite remainders obtained from the analytical integration of the subtraction terms over the degrees of freedom of the unresolved parton(s). The one-loop virtual amplitudes for e + e − → 5 jet required in this paper are computed within the generalized D-dimensional unitarity framework, as described in Refs. [48][49][50]; a few technicalities relevant to the five-jet case are given in sect. 2.1. The remaining contributions are calculated using MadFKS [51], which is also employed to perform the integration over the phase space of the final shortdistance cross sections. Details of MadFKS relevant for this computation are reviewed in sect. 2.2.

Calculation of one-loop amplitudes
Within the context of generalized D-dimensional unitarity, we compute the so-called primitive amplitudes [52], which are gauge-invariant subsets of color-ordered amplitudes. The color decomposition of amplitudes that we need in this paper coincides with the color decomposition of QCD amplitudes without any colorless vector boson γ * /Z. We use the color decomposition introduced in ref. [53]. The relation between primitive and color-ordered amplitudes that we need in this paper can be found in ref. [50].
In general, the amplitudes needed for the NLO calculation of e + e − → γ * /Z → 5 partons are very similar to the amplitudes with a W -boson and five partons. The main differences between amplitudes that involve charged and neutral currents originate from the different couplings of W -bosons and γ * /Z bosons to fermions. In particular, in the Standard Model, the coupling of the W -boson to fermions is only left-handed, while for γ * /Z one has to sum over left-and right-handed states. Indeed, the Z-boson couples to fermions via a vector and a vector-axial coupling g V γ µ + g A γ µ γ 5 . We can rewrite this coupling through left-and right-handed projection operators P L,R = (1 ∓ γ 5 )/2 as where g L,R = g V ∓g A . Therefore, all the one-loop amplitudes relevant for this paper can be obtained by considering quarks coupled to a vector current only; the only subtlety is that the vector current should couple to left-and right-handed quarks with different strength.
As we already mentioned, many of the non-trivial amplitudes that we need are identical to the amplitudes calculated for 0 → W + 5 partons case [50]. We note, however, that new amplitudes appear if γ * /Z couples directly to a loop of virtual fermions. There are two reasons for this. First, such amplitudes are not present in the W -boson case studied in ref. [50] because of charge (or flavor) conservation. Second, certain parts of those amplitudes are related to an axial anomaly and, therefore, violate the symmetry between the vector current and the axial-vector current. It is interesting to remark that within the context of generalized D-dimensional unitarity, a correct computation of the axial anomaly entails a literal implementation of the 't Hooft-Veltman prescription for dealing with γ 5 in dimensional regularization. We have calculated the required amplitudes where γ * /Z couples to the fermion loop but we neglect them in this paper since previous experience with those amplitudes shows that they are very small [22,23,54,55], especially when compared to the residual theoretical uncertainty of the five-jet rate.
The cross section involves Born and virtual amplitudes with two or four quarks in the final state. We relate amplitudes involving a Z-boson exchange to amplitudes that only involve the photon (vector) exchange: In eq. (2.2) all particles are outgoing, the dots denote all gluon momenta, and the propagator factor P(s) is given by where M Z and Γ Z are the mass and the width of the Z boson respectively. The left-and right-handed coupling of the electrons (f = e) and quarks (f = q) to the Z boson are explicitly given by where θ W is the Weinberg angle, T e 3 = T d 3 = −1/2, T u 3 = 1/2 are the values of weak isospin for quarks and leptons and Q e = −1, Q u = 2/3, Q d = −1/3 are the respective electric charges.
The amplitudes involving four quarks and a Z-boson can be written in a similar way. The amplitude that involves different quark flavors can always be written as the sum of two amplitudes, where the Z-boson couples to a particular quark line Note that when amplitudesÃ are computed, the Z-boson is allowed to couple to the quark flavor indicated by first and second argument ofÃ. The expression forÃ Z amplitude in terms ofÃ * γ is identical to eq. (2.2), so we do not repeat it here. Finally, if the flavor of the two quark lines is the same, we include the symmetry factor 1/4 and anti-symmetrise with respect to the exchange of the quarks or anti-quarks. We point out that the full one-loop matrix elements squared that we use in this paper were checked against a similar computation performed by the BlackHat collaboration [56][57][58], and complete agreement was found.

Real emission corrections with MadFKS
MadFKS is based on the (FKS) subtraction formalism of ref. [59]. The implementation of the FKS procedure is fully automated in MadFKS. In essence, MadFKS goes through the following steps. First, it determines the partonic processes that contribute to a given physical reaction, and their singularity structures. Then, it constructs the real-emission matrix Process # of FKS pairs e + e − → qqgggg 3 e + e − → qqq ′q′ gg 7 e + e − → qqqqgg 4 e + e − → qqq ′q′ q ′′q′′ 3 e + e − → qqqqq ′q′ 2 e + e − → qqqqqq 1 elements, their subtraction terms, the finite remainders, and the Born matrix elements. Finally, it proceeds to the actual computation, by sampling (possibly with multi-channeling techniques) the phase space, by evaluating the short-distance cross sections, and by returning weighted parton-level kinematic configurations, that can be used to construct as many observables as one likes. All matrix elements, except those of virtual origin, are obtained by calls to MadGraph [60] routines. The virtual matrix elements are on the other hand computed as described in sect. 2.1. We point out that MadFKS gives, for each phase-space point, a four-momentum configuration as input to the code of sect. 2.1, which returns three numbers -corresponding to the double and single (IR) pole residues, and to the finite part; the talk-to between the two codes uses the Binoth-Les Houches interface [61]. As clarified in ref. [51], MadFKS integrates quantities that are locally finite in the phase space, and in four dimensions. Therefore, the pole residues provided by the virtual amplitudes are used only to check that they match those that are known analytically from the subtraction procedure, in this way ensuring that KLN cancellation does indeed take place. The only output of the code of sect. 2.1 used in the integration of the short-distance cross sections is thus the finite part, defined according to the conventions given in the Appendix B of ref. [51].
We have performed the calculation of the five-jet cross section in a straightforward manner. We included all partonic processes, and performed explicit sums over colors and helicities, except in the case of virtual amplitudes. For the latter, we have computed separately the leading-and subleading-color contributions, and performed the sum over helicities using Monte Carlo methods. We have used five massless quark flavors. We remind the reader that the FKS formalism is particularly efficient in keeping the number of subtraction terms to a minimum. Furthermore, the real-emission matrix element minus the subtraction terms is re-organized into a sum that gets as many contributions as the subtraction terms themselves, and that are separately finite, which implies that they can be (and are) integrated independently from each other. Physically, these contributions corresponds to pairs of particles (called FKS pairs) that can give one soft and/or one collinear singularity at most. We report in table 1 the number of FKS pairs relevant to the various real-emission processes that contribute to the five-jet cross section.
With five massless flavors, the number of independent partonic subprocesses that contribute to the five-jet cross section is 25. Using the entries of table 1, this implies 81 FKS pairs in total, i.e. 81 independent integrations. On the other hand, the complexity of the kinematics is such that even in the context of an adaptive integration it may be very difficult to map correctly all the peaks of the Feynman diagrams, and thus to have a stable numerical behavior. We have therefore preferred to adopt a multi-channelling integration strategy, that in MadFKS follows the same procedure as in MadGraph [62]. In doing so, the numbers of integration channels we deal with at the real-emission and virtual level are equal to 3620 and 2×1408 respectively (the factor of two in the virtual amplitudes being due to the independent integration of the leading-and subleading-color contributions). These numbers are much larger than the 81 FKS pairs we started with; however, the Feynmandiagrams peaks can now be mapped accurately by the integration routines, and relatively small statistics is sufficient in each channel to obtain numerical stability. We conclude by stressing that the MadFKS integration channels are fully independent. Furthermore, they are not determined dynamically (e.g. by performing a preliminary integration of the cross section), but are defined a priori, by considering the topologies of the Feynman diagrams that contribute to the relevant partonic processes. The whole organization of the calculation is therefore inherently parallel.

Phenomenology of five-jet production
In this Section we present the results of our calculation. We consider e + e − → jets and define jets using the Durham jet algorithm [28] with resolution parameter y cut . The following distance between each pair of particles is used in the Durham jet algorithm where s is the center-of-mass energy of the collision squared, E i is the energy of the parton i, and θ ij is the relative angle between the partons i and j, in the e + e − center-of-mass reference frame. The pair of particles with the smallest distance is clustered together by adding their four-momenta, as long as y ij < y cut , and the procedure is then iterated. When all distances y ij are larger than y cut , the recombination stops and the number of jets in the event is defined to be equal to the number of (pseudo)-particles left at that stage. In this paper we consider two observables which we define with the Durham jet algorithm. The first observable is the differential distribution with respect to the five-jet resolution parameter y 45 , normalized to the total cross section for e + e − → hadrons, σ tot (which we compute at the NLO, i.e. at O(α s )). The resolution parameter y 45 is the maximal value of y cut such that a given event is classified as a five-jet event by the Durham jet algorithm. We note that where σ 5−jet incl is the inclusive five-jet production cross section in e + e − annihilation. The second observable that we study is the five-jet rate R 5 (y cut ). It is defined as follows where σ 5−jet excl (y cut ) is the exclusive five-jet production cross section. It is calculated by applying the Durham jet algorithm to the given event, and by requiring that exactly five jets are reconstructed.
When we compute σ −1 tot dσ/d ln y −1 45 and R 5 in perturbative QCD, we obtain a power series in the strong coupling constant where µ is the renormalization scale, b 0 = (33 − 2n f )/3 and n f = 5 is the number of quark flavors that we treat as massless. The top quark is considered to be infinitely heavy and is completely neglected in our computation. It is important to emphasize that the coefficients A 45,5 and B 45,5 depend on y 45 and y cut , respectively, but not, say, on the total center-of-mass energy squared. This feature is a consequence of the following approximations employed in our computation: 1) all particles, except the Z boson, are treated as massless; 2) the observables that we are interested in are sufficiently inclusive so that the vector and the axial currents do not interfere; 3) we neglect triangle fermion diagrams that lead to the axial anomaly so that (for equal couplings) vector and axial current contributions to the final result are equal. These three points are sufficient to ensure that A 45,5 and B 45,5 are independent of the electroweak parameters and the center-of-mass energy squared. Experimentally [16], five-jet observables are computed using the reconstructed momenta and energies of charged and neutral particles. Measurements are corrected for detector effects, so that final distributions correspond to stable hadrons and leptons, and for initial-and final-state photon radiation, which is a sizable correction for LEP2 data. Above the Z peak, relevant backgrounds are subtracted; the most important among them is W -pair production. The experimental uncertainties are estimated by varying event-and particle-selection cuts. They are below 1% at LEP1 and slightly larger at LEP2. Further details of the experimental analysis can be found in ref. [16].
In fig. 1, we compare ALEPH LEP1 data [16] for 1/σ tot dσ/d ln y −1 45 with the hadronlevel predictions of three event generators -PYTHIA, HERWIG and ARIADNE. We observe that these event generators describe experimental data fairly well; differences between data and theoretical predictions are below twenty five percent in the central 4.5 < ln y −1 45 < 9 region of the distribution, where the statistical accuracy of the data is good. This is an impressive accomplishment since σ −1 tot dσ/d ln y range of ln y −1 45 . In addition, it follows from fig. 1 that the difference between hadronization corrections, as calculated using different event generators, can be as large as 20-30%.
We attribute these features to the inability of PYTHIA, HERWIG and ARIADNE to describe hard perturbative radiation correctly. Indeed, these programs generate highmultiplicity final states starting from hard low-multiplicity processes; they produce additional jets by means of parton/dipole showers. Since these showers describe hard large-angle emissions only approximately, the so-called hadronization corrections attempt to correct for this (perturbative) deficiency. While this problem is unavoidable if traditional event generators are used to describe high-multiplicity final states, techniques exist to match parton showers and high-multiplicity matrix elements in a consistent manner, thereby improving the pure-perturbative part of event generators. One such technique is the CKKW matching procedure [47], which is implemented as default in the SHERPA event generator. The comparison of ALEPH LEP1 data with SHERPA predictions, as well as the hadronization corrections derived from SHERPA, are shown in fig. 2. Two hadronization models -Lund string [63] and cluster [64] -are employed. In the central part of the Hadr. Corr.
-ln(y 45 ) Lund Cluster Figure 2: ALEPH data for the y 45 distribution at LEP1, compared to SHERPA results. Two hadronization models -Lund string [63] and cluster [64] -are employed. The lower pane in the left plot shows the relative difference between Sherpa predictions with the two hadronization models, and ALEPH data. In the right plot, the hadronization corrections for the two models are shown.
distribution, SHERPA results agree with ALEPH data to 20 − 25%, similar to traditional event generators. Moreover, in the region of moderately small values of ln y −1 45 , where fixedorder perturbative description is reliable, the hadronization corrections are below twenty percent, in sharp contrast with estimates of hadronization corrections based on PYTHIA, HERWIG and ARIADNE. It is important to emphasize that, although in that region of ln y −1 45 traditional event generators provide slightly better description of data compared to SHERPA, this does not mean that hadronization corrections extracted with the former codes are more reliable. Indeed, traditional event generators achieve agreement with data at the price of very large hadronization corrections. This feature precludes a clear separation between long-and short-distance phenomena, which is crucial for the procedure outlined below eq. (1.1) to be meaningful.
The ALEPH data exhibit a characteristic turnover shape. This turnover means that for small values of y 45 , the result is dominated by exclusive five-jet production with very small resolution parameter, where fixed order perturbation theory fails and a resummation is required to achieve meaningful results. A resummation of α n s L 2n and α n s L 2n−1 terms, where L = ln y −1 cut , was performed for R 5 in [28], while no resummation is currently available for the five-jet resolution parameter distribution. However, there seems to be no region in L where this resummation can be valid since two conditions L ≫ 1 and α s L ≪ 1 should be satisfied simultaneously. Taking α s ∼ 0.15 as a typical value of the strong coupling constant 4 , we find that L should be smaller than 6. On the other hand, practical experience with resummations suggests that L ≫ 5 is what can be considered as a large logarithm. Clearly, 5 ≪ L < 6 leaves very little room for the validity of this approach. It should be possible to improve on the resummation by including sub-leading logarithms and matching to NLO QCD computations. However, since we do not perform any resummation in this paper, we require ln y −1 45 , ln y −1 cut < ∼ 6 for the comparison of the NLO QCD computation with data. Interestingly, a similar upper bound on ln y −1 45 appears because we neglect the mass of b-quarks in our computation. This implies that the resolution parameter times the center of mass energy should be larger than the b-quark mass, i.e. sy 45 > m 2 b , which translates into ln(y −1 45 ) < ln(s/m 2 b ) < ∼ 6, for s = M 2 Z . When fixed-order perturbative QCD calculations are compared to experimental data, the choice of the renormalization scale becomes an important issue. Traditionally, multijet observables in e + e − annihilations are computed in perturbative QCD by evaluating the strong coupling constant at the center-of-mass energy. However, for large numbers of jets this choice should be reconsidered, since the hardness of each jet decreases with their number. Dynamical renormalization scales used in event generators account for this effect by relating the choice of the renormalization scale to the event kinematics. Our choice of the renormalization scale is also motivated by dynamical considerations. To this end, we consider the clustering history of five-and six-parton configurations that results from using the Durham jet algorithm. We compute the average value of √ y 23 , where y 23 is the three-jet resolution parameter, using only phase-space weights. We find this average to be approximately equal to 0.3. Since √ y 23 s is, roughly, the relative transverse momentum of production in the range 3 < ln y −1 cut < 7.
the hardest branching in the clustering history, we select µ = 0.3 √ s as the default choice for the renormalization scale of α s that we use to describe the five-jet production. With this choice of the renormalization scale, we compare in fig. 3 ALEPH LEP1 data for σ −1 had dσ 45 /d ln y −1 45 and R 5 with our leading and next-to-leading results. We use the value of strong coupling constant α s (M Z ) = 0.130 for leading order computations and α s (M Z ) = 0.118 for next-to-leading order computations. While it is not customary to change the value of the strong coupling constant from one order in perturbation theory to the other in applications of QCD to e + e − physics, it is done routinely in the context of hadron collider physics. Our choice of the leading-order value for α s is motivated by fits of parton distribution functions and of the strong coupling constant described in ref. [65]; our next-to-leading order value for α s is close to the world average [66,67]. We do not apply hadronization corrections at this stage. In order to assess the perturbative uncertainty, a scale variation by a factor of two around the default scale µ = 0.3 M Z is performed. A close inspection of the two plots shown in fig. 3 reveals that the most important effect of the NLO QCD corrections is the reduction in the uncertainty related to the renormalization scale dependence. The renormalization scale uncertainty is reduced from [−30%, +45%] at leading order to [−20%, +25%] at next-to-leading order. When leading order and next-to-leading order results are compared for µ = 0.3 M Z , the QCD corrections increase the leading order predictions by 10 − 20% 5 . The agreement between NLO QCD results and ALEPH data is very good for both observables considered. However, for the σ −1 had dσ 45 /d ln y −1 45 distribution systematic differences appear for ln y −1 45 > 5.2, whereas the R 5 data can be described by fixed-order QCD without hadronization corrections all the way up to ln y −1 cut = 6.5.

The strong coupling constant from five-jet observables
It follows from the discussion in the previous Section that the use of the world-average value of the strong coupling constant results in good agreement between parton-level NLO predictions, and LEP1 data for the five-jet resolution parameter and the five-jet rate. Therefore, we can turn this consideration around, and extract the value of the strong coupling constant from these two five-jet observables at LEP1 and LEP2. To combine the values of the strong coupling constant extracted from different observables and at different energies, we use the procedure advocated by the LEP QCD working group [68]. Technical details of the fit procedure are reported in the Appendix. We consider ALEPH data [16] for R 5 and σ −1 tot dσ/dy 45 [64], which we therefore adopt as our default as well.
Incidentally, we observe that it gives smaller hadronization corrections in the kinematic range of interest, than the Lund string hadronization model [63], which is also available in  SHERPA. We use the Lund string model to estimate systematic uncertainties related to hadronization effects.
Since we use fixed-order perturbative results and do not perform any resummation, it is not possible to describe the data in the full kinematic ranges studied by experiments. This feature makes the choice of the kinematic range used in the fit an important but, unfortunately, somewhat a subjective issue. In general, we attempt to take the fit range as large as possible, with the condition that our computations are reliable and that the data quality is good. In the determination of the central value of α s at LEP1, we consider 3.8 ≤ − ln y 45 ≤ 5.2 (7 data points) for the five-jet resolution parameter distribution, and 4.0 ≤ − ln y cut ≤ 5.6 (8 data points) for R 5 . In order to estimate the error on α s related to our choice of the fit range, we extract the value of α s by performing a second fit, with larger ranges 3.4 ≤ − ln y 45 ≤ 5.6 (11 data points) for the five-jet resolution parameter, and 3.4 ≤ − ln y cut ≤ 6.0 for R 5 (13 data points). The difference between the values of α s obtained in the two fits is called the "fit range" error; it is supposed to quantify the uncertainty on α s due to the choice of the data points included in the fits.
At LEP2 the situation is different. Firstly, data are given with a coarser binning and, secondly, large fluctuations are present in experimental results at small values of ln y −1 45 and ln y −1 cut (for example, for some center-of-mass energies the corresponding observables are not even monotonic). Because of this, we decided to exclude those data points from our fits, effectively reducing the fit ranges. We note that those data would have had a modest impact on the final result anyhow, because they are affected by fairly large errors.  We use 4.8 ≤ − ln y 45 ≤ 6.4 for the five-jet resolution parameter (2 data points per √ s), and 2.1 ≤ − log 10 y cut ≤ 2.9 for R 5 (4 data points per √ s), to find the central values of α s .
In order to estimate the fit-range error, we employ 4.8 ≤ − ln y 45 ≤ 5.6 (1 data point per √ s), and 2.1 ≤ − log 10 y cut ≤ 2.5 (2 data points per √ s), since the choice of these ranges leads to the largest changes in the values of the strong coupling constant compared to the α s values obtained from fitting with the default ranges.
The results of our fits to LEP1 data are shown in table 2. The agreement between the two values of α s (M Z ) extracted with and without hadronization corrections is impressive; the difference is completely negligible compared to the overall uncertainties. This result could have been anticipated by inspecting fig. 2, which shows that, in the fit region, hadronization corrections are small, in particular when the default SHERPA choice, the cluster model, is used. We note that if we use the hadronization corrections as given by conventional HERWIG, PYTHIA, or ARIADNE, without matching them to high-multiplicity matrix elements, the picture changes drastically and the values of α s (M Z ) extracted with or without hadronization corrections are quite different from each other. We also note that the overall errors of the two results given in table 2 are slightly smaller when including hadronization corrections. This is due to a marginally better description of the data in the central region of the fit range -which leads to smaller value of α s (M Z ) and thus to smaller perturbative errors. However, the error reduction is partially compensated by the degradation of the fit quality when including larger values of ln y −1 45 , ln y −1 cut , where hadronization corrections increase. This feature leads to larger fit-range error compared to the no-hadronization case. Note also that if we extract the values of α s (M Z ) by fitting 6 the five-jet resolution parameter distribution and R 5 separately, we obtain 0.1168 +0.0076 −0.0060 and 0.1151 +0.0071 −0.0056 respectively. These values are consistent with the result of the combined fit shown in table 2 but have slightly larger errors. From table 2, it is clear that the sensitivity of the five-jet observables to α s is very high, as illustrated by the tiny statistical errors. This sensitivity is ultimately related to the high power of α s that enters the five-jet observables. In spite of this, the overall error is not particularly small, since the perturbative uncertainty is still quite sizable at this order in perturbative QCD.
Compared to LEP1, there are important differences when we extract α s by fitting to the LEP2 data. Firstly, because hadronization corrections are negligible at LEP1, and because these corrections decrease with energy, we do not consider them for LEP2. Secondly, for the reasons explained above, we do not consider the data points at small values of ln y −1 45 and ln y −1 cut . This fact, combined with coarser binning of data, pushes us to the region of y 45 that may be affected by large logarithms of the resolution parameter. As a result, we find larger fit-range errors at LEP2 than at LEP1. The statistical errors are also much larger at LEP2 than at LEP1, as one expects given the luminosities collected. On the other hand, since the effective strong coupling is smaller at LEP2, the perturbative uncertainty affecting five-jet observables decreases, making the α s extraction at LEP2 competitive with that done at LEP1. This effect is particularly strong for R 5 , for which the extracted central value of α s is slightly smaller than for the five-jet resolution parameter. In table 3 we present the α s values obtained by fitting separately the five-jet resolution parameter and R 5 at LEP2, since they differ from each other by a larger amount than at LEP1. Still, both values are within one standard deviation from the strong coupling constant that we obtain by performing a simultaneous fit to the two observables. We take the latter value, given in the third column of table 3, as our best determination of α s from LEP2 data.
We obtain our final estimate of the strong coupling constant by combining the values of α s (M Z ) extracted from LEP1 and LEP2 data. We assume that the statistical and systematic errors of the two results are not correlated (an assumption which is strictly correct for the former, and a very good approximation for the latter), while the perturbative errors are considered to be fully correlated. The correlation of the perturbative uncertainties is due to the fact that we estimated them by varying the renormalization scale, which results in changes of the cross sections whose pattern is independent of the center-of-mass energy. It is quite likely that a more sophisticated approach to estimating perturbative errors (see e.g. ref. [69]) will result in a smaller uncertainty on α s . Hence, the procedure that we employ in this paper is rather conservative. Using the results of tables 2 and 3, we finally obtain We note that if we perform the fit to both LEP1 and LEP2 data simultaneously, we obtain α s (M Z ) = 0.1156 +0.0045 −0.0041 , in perfect agreement with eq. (4.1). The value of α s (M Z ) that we extract from five-jet observables at LEP can be compared with other recent determinations of this quantity, shown in table 4. We see that both the central value of α s and its error, obtained from fitting five-jet observables, compare well with other determinations. On the other hand, it is interesting that α S (M Z ) in eq. (4.1) world average 0.1184 ± 0.0007 [67] e + e − → five jets 0.1156 ± 0.0038 this paper is lower than the world average. It is peculiar that a number of recent determinations of α s arrived at a similar conclusion.

Conclusions
In this paper we study the production of five jets in e + e − annihilation at LEP1 and LEP2. We improve the perturbative QCD predictions for five-jet observables, 1/σ tot dσ/dy 45 and R 5 , by computing the NLO QCD corrections. For suitably chosen renormalization scales, such corrections are between ten and twenty percent 7 . They reduce the scale uncertainty by about a factor of two with respect to the LO predictions, and lead to a better agreement between theoretical predictions and experimental data. We point out that hadronization corrections computed with event generators whose showers are not matched to high-multiplicity matrix elements (such as out-of-the-box HER-WIG, PYTHIA, and ARIADNE) are large and uncertain. For this reason, we believe it is important to describe five-jet observables in a way that incorporates high-multiplicity treelevel matrix elements. This is provided by the event generator SHERPA, which implements the CKKW procedure for matching tree-level matrix elements to parton showers. In this way, an improved description of five hard, well-separated partons is obtained, which in turn results in fairly small hadronization corrections in the range where fixed-order perturbative results are most reliable.
We extract the strong coupling constant from the distributions of the five-jet resolution parameter and the five-jet rate, as measured at LEP1 and LEP2 by the ALEPH collaboration. We find α s (M Z ) = 0.1156 +0.0041 −0.0034 , which compares well with other recent determinations of the strong coupling constant, and is somewhat lower than the current world average value. We stress that our treatment of the uncertainties on α s is conservative. A detailed knowledge of the experimental systematics, and a more sophisticated approach to theoretical errors will very likely lead to a higher precision in the determination of α s from five-jet observables at LEP.

Acknowledgments
We are grateful to Hasko Stenzel for the participation in earlier stages of this work. We would like to thank Stefan Höche for providing us with SHERPA results. Useful conversations with Andrea Banfi, Guenther Dissertori, Stefano Forte, Michelangelo Mangano, Gavin Salam, Peter Skands, Roberto Tenchini, Paolo Torrielli, and Bryan Webber are gratefully acknowledged. During the work on this paper, we have benefited from the hospitality extended to us by the Aspen Center for Physics, the Fermilab Theory group, and the Theory Unit at CERN. This research is supported by the NSF under grant PHY-0855365, by the start-up funds provided by Johns Hopkins University, by the British Science and Technology Facilities Council, and by the Swiss National Science Foundation under contract 200020-126691.

Appendix: details of the fit
In this Appendix, the details of the fitting procedure are described. In the fit, we consider two observables -the five-jet resolution parameter distribution and the five-jet rate. These observables are measured at several energies; the results are available in the form of binned distributions. In principle, we can extract the value of the strong coupling constant from any of the bins but, clearly, the availability of many bins helps in decreasing the errors. The problem is that both the experimental and theoretical errors affecting different bins may be correlated, and it is important to take these correlations correctly into account, to obtain a proper estimate of the overall uncertainty in the determination of α s . Our procedure follows closely the approach of the LEP QCD working group described in ref. [68]. A possible way to treat theoretical uncertainties is discussed in ref. [69], but we follow a simpler approach that is described below.
The set of all available bins, for a chosen range of y 45 and y cut , is a set of observables L O from which values of the strong coupling constant can be determined. Each member of L O is represented by three numbers , where X i , and σ stat,syst i are the central value and the statistical and systematic uncertainties respectively, for a given observable in the bin i. If uncertainties are asymmetric, O i is a collection of five numbers, and the discussion below applies to positive and negative errors separately. In what follows, we use α s ± δα s as the shorthand notation instead of the full one α s +δ + αs −δ − αs . As we already mentioned, each of the observables from the list L O can be used to determine the value of the strong coupling constant. This is done by solving the equation where T i is the (parton level) theoretical prediction, H i is the hadronization correction, and E i is the experimental value for the bin i. The theoretical prediction T i depends on α s (M Z ) and the renormalization scale µ. As discussed in the text, for the hadronization correction we can use either the cluster or the Lund string model and we choose the former as our default. We write the value of the strong coupling constant at M Z , obtained by solving eq. (5.1), as

4) δα i,hadr
s is obtained by solving eq. (5.1) for α s with the same settings as in item 1), except that the Lund string model is used for hadronization. We define ±δα i,hadr s = ±|α s − α i s |. We note that this is a conservative choice, since clearly α s − α i s is either positive or negative.
The result of this procedure is a set of values of the strong coupling constants, α i s , with the corresponding errors. They need to be combined to obtain the average value of α s . To this end, it is necessary to construct the covariance matrix, which we define as the sum of the covariance matrices for statistical, systematic, perturbative, and hadronization errors. If we denote generically by δα i s one of these four errors, the corresponding covariance matrix is where C ij is the statistical correlation between α i s and α j s (see later). The covariance matrix is used to calculate the average value of the strong coupling constant and its error in a standard way. We compute the weights (V −1 ) kl , 1 ≤ i ≤ N , N = dim(V ), (5.4) and obtain the estimate of the average of the strong coupling constant and of its error It is easy to see that if two or more errors that enter the definition of the covariance matrix are close numerically, and the absolute values of the corresponding statistical correlations C ij is close to one, eq. (5.3) may lead to pathological results, since the weights w i tend to grow large in absolute value and to have opposite signs. To avoid this, the LEP QCD working group [68] adopts the formula which is essentially equivalent to taking the largest possible C ij in eq. (5.3), that still leads to non-pathological weights w i . Clearly, eq. (5.7) is an overestimate of the correlation if the actual C ij is small. For this reason, in our fit we use a slightly modified formula for the covariance matrix V ij = δ ij δα i s 2 + (1 − δ ij ) min (δα i s ) 2 , (δα j s ) 2 , C ij δα i s δα j s . (5.8) As we explain below, for the observables and errors that we consider, eq. (5.8) coincides with eq. (5.7) in all cases, except for the statistical correlation between σ −1 tot dσ/dy 45 and R 5 .
In summary, we take the covariance matrix to be where each of the terms on the right-hand side of eq. (5.9) is constructed according to eq. (5.8) using δα i,stat s , δα i,syst s , δα i,scale s , and δα i,hadr s respectively. As far as the off-diagonal terms of the various V matrices are concerned, we have assumed what follows: • Statistical errors are uncorrelated between different center-of-mass energies. At a given center-of-mass energy, σ −1 tot dσ/dy 45 data are uncorrelated, R 5 data are fully correlated, and R 5 data are correlated with σ −1 tot dσ/dy 45 ones for y cut ≤ y 45 . We have explicitly computed the coefficient C ij relevant to the σ −1 tot dσ/dy 45 − R 5 correlation, and have used them in eq. (5.8).
• Since we do not have detailed information about correlations of systematic uncertainties in ALEPH data, we assume conservatively that all systematic errors at a given center-of-mass energy are fully correlated. We also assume that systematic errors are completely uncorrelated between LEP1 and LEP2, but that they are fully correlated in the measurements performed at different LEP2 energies.
• Perturbative errors are assumed to be fully correlated, for all observables and energies. See also the main text for a comment on this point.
• Hadronization errors, that we compute only at LEP1, are assumed to be fully correlated.
Finally we point out that in the computation of the central value of α s , according to eq. (5.5), we neglect the off-diagonal entries of V scale and V hadr [68]. On the other hand, we include all the off-diagonal entries in the computation of the standard deviation according to eq. (5.6). In fact, in the presence of errors numerically very close to each other (which is the case for the perturbative and hadronization errors), the result for the average value of α s tends to assume the value of the input with the smallest error, which is again an artifact of the combination procedure. We have checked that the value of the strong coupling constant that we obtain in this way is statistically fully compatible with the result we would have obtained by considering all correlations when determining the central value of α s .