Measurement of the strong coupling alpha_S from the three-jet rate in e+e- - annihilation using JADE data

We present a measurement of the strong coupling alpha_S using the three-jet rate measured with the Durham algorithm in e+e- -annihilation using data of the JADE experiment at centre-of-mass energies between 14 and 44 GeV. Recent theoretical improvements provide predictions of the three-jet rate in e+e- -annihilation at next-to-next-to-leading order. In this paper a measurement of the three-jet rate is used to determine the strong coupling alpha_s from a comparison to next-to-next-to-leading order predictions matched with next-to-leading logarithmic approximations and yields a value for the strong coupling alpha_S(MZ) = 0.1199+- 0.0010 (stat.) +- 0.0021 (exp.) +- 0.0054 (had.) +- 0.0007 (theo.) consistent with the world average.


Introduction
The annihilation of an electron-positron pair into a pair of quarks provides an ideal laboratory to test the theory of the strong interaction, Quantum Chromodynamics (QCD) [2][3][4]. The free parameter of QCD, the strong coupling α S , can be determined with events with more than two jets in the final state. To first order perturbation theory the radiation of a gluon from a quark is proportional to α S . To determine this fundamental constant the observed three jet rate is compared to a perturbative expansion which predicts the three jet rate as a function of a single parameter α S [5][6][7].
Recently theoretical progress has been made leading to a significant improvement in the prediction of the three jet rate [8][9][10][11] as a function of α S . Previously e + e − -event shape distributions and the three-jet rate were only known to next-to-leading-order (NLO) accuracy, now QCD calculations to next-to-next-to-leading-order (NNLO) are available. These predictions were used to determine α S based on a single value of the resolution parameter of the three-jet rate using the Durham algorithm [12], with data taken at a centre-of-mass energy at 91 GeV with the ALEPH experiment [13].
In this analysis we use data taken with the JADE experiment located at the PETRA collider at DESY between the years 1979 and 1986. We measure the strong coupling α S at six centre-of-mass-energies in the range between 14 and 44 GeV. Besides using a different energy range we perform a fit in a range of the three-jet resolution parameter. We present a matching scheme to combine NNLO predictions together with next-to-leading logarithm approximations (NLLA). The matched predictions are used to determine the strong coupling α S . The analysis follows closely the determination of α S measuring the four-jet rate using data collected with the JADE experiment [14].

Observable
To determine a multijet-rate a jet finding algorithm has to be applied to particles observed in the final state. For this analysis we use the Durham algorithm [12], which selects jets according to the jet resolution parameter y cut . The Durham algorithm defines initially each particle as a proto-jet and a resolution variable y ij is calculated for each pair of proto-jets i and j: where E i and E j are the energies of jets i and j, cos θ ij is the cosine of the angle between them and E vis is the sum of the energies of the detected particles in the event (or the partons in a theoretical calculation). If the smallest value of y ij is less than a predefined value y cut , the pair is replaced by a new proto-jet with four-momentum p µ k = p µ i + p µ j , and the clustering starts again. Clustering ends when the smallest value of y ij is larger than y cut , and the remaining proto-jets are counted as final selected jets.
Perturbative QCD calculations predict the fraction of three jet events R 3 (y cut ) as a function of y cut and α S . The NNLO QCD prediction for the three-jet rate can be written as R 3 (y cut ) = σ 3-jet (y cut ) σ tot =α S A 3 (y cut ) +α 2 S B 3 (y cut ) +α 3 S C 3 (y cut ), withα S = α S (µ)/(2π) the only free parameter. For this analysis the coefficients A 3 , B 3 and C 3 are taken from [15]. Equation 2 is shown for renormalisation scale µ = Q, where Q is the physical scale usually identified with the centre-of-mass energy √ s for hadron production in e + e − -annihilation. The terms generated by variation of the renormalisation scale parameter x µ = µ/Q are implemented according to [15]. It is well known that for small values of y cut the fixed order perturbative prediction is not reliable, because the expansion parameter (α S )L 2 , where L = − ln y cut , logarithmically enhances the higher-order corrections. For instance, (α S ) ln 2 (0.01) ≈ O(1). Thus, one has to perform the all-order resummation of the leading and NLLA contributions. This resummation is possible for the Durham algorithm using the coherent branching formalism [12,16]. Matched NLO+NLLA predictions for jet rates have been compared to ALEPH data and good agreement was found, see Fig. 2 of [17]. There an improved resummation formula was used, where part of the subleading logarithms, that can be controlled systematically by including the 'K-term' [18][19][20], are also taken into account (NNLO+NLLA+K). In this paper we use the same resummation formula and extend the matching to the NNLO accuracy, which requires the expansion of the improved resummation formula up to O(α 3 S ). In the resummed prediction we use the one-loop formula for the running coupling. One could also use higher-loop running, but the difference would be in the coefficients of the subleading logarithms (NNLL and higher), which are not controlled systematically in the resummed prediction and therefore are neglected. Expanding the resummation formula in [12], we find with the R (i,j) 3 coefficients are given in Table 1. The dependence on the renormalisation scale µ can be obtained by making the substitution and keeping the terms up to O(α 3 S ) in the expansion. In Eq. (4), with C A and C F are the quadratic Casimir operators of the gauge group in the adjoint and fundamental representation, while n f is the number of light flavours (we use n f =5). The matched NNLO+NLLA+K predictions compared to NNLO, NLO and matched NLO+NLLA+K, NNLO+NLLA predictions are shown in Fig. 1.
As opposed to event-shapes, such as the y 23 -distribution [21,22], jet rates do not obey simple exponentiation (except for the two-jet rate). For an observable that does not exponentiate, the viable matching scheme is the so-called R-matching [23]. To obtain the R-matched predictions, we subtract the expansion of R 3 from the resummation formula and add the corresponding NNLO prediction given by Eq.2, L j and A 3 , B 3 , C 3 as in Eq. 2. Also in the case of jet rates the resummed logarithm is fixed to L = − ln y cut unambiguously, and the theoretical uncertainty of the prediction is estimated by varying the renormalisation scale (see Fig. 1) as in [24]. The theoretical uncertainty estimated by varying the renormalisation scale is found to be larger by using matched NNLO+NLLA+K predictions instead of NNLO+NLLA predictions. x − 10 9 y f , 3 Analysis Procedure

The JADE Detector
A brief description of the JADE detector focusing on the major detector parts used in this analysis is given here. Primarily the momenta of charged and neutral particles are used in this analysis. The trajectories of charged particles are mainly reconstructed with the central tracking detector, consisting mainly of a large scale jet chamber. The tracking detector is located within a 0.48 T solenoidal magnetic field. Charged and neutral particles except muons and neutrinos are reconstructed with an electromagnetic calorimeter which surrounds the magnetic coil. The calorimeter consists of lead glass Cherenkov counters and is separated in a barrel and two endcap regions. A more complete summary of the JADE detector can be found in [25].

Data and Monte Carlo Samples
For this analysis we use the identical data sample as used for previous JADE analyses [14,24,26]. The data were taken between 1979 and 1986, adding up to an integrated luminosity of about 195 pb −1 . It is subdivided into six data taking periods with different average centre-of-mass energies of 14 GeV, 22 GeV, 34.6 GeV, 35 GeV, 38.3 GeV and 43.8 GeV.
The number of selected events ranges between 1403 events at 22 GeV and 20876 events at 35 GeV.
To correct for acceptance and resolution effects as well as for hadronisation effects a large sample of Monte Carlo events is generated. The Monte Carlo generators are tuned to match events taken with the OPAL experiment at a centre-of-mass energy of 91 GeV [27,28]. Using these same parameters, except of setting the appropriate centre-of-mass energy, a good description of the JADE data, down to smallest energies, is achieved [29]. For the correction of acceptance and resolution effects the events are passed through a simulation of the JADE experiment. Events generated with PYTHIA 5.7 [30] are used as default for the correction of acceptance and resolution effects. As a cross check events simulated with the HERWIG 5.9 [31] event generator are utilised. To assess the changes in the three-jet rate distribution originating from the transition from partons to hadrons three different Monte Carlo generators are applied. As default PYTHIA 6.158 is used and events produced with the HERWIG 6.2 and ARIADNE 4.11 [32] generators are used as a consistency check.
Due to the good agreement between data and the predictions from the PYTHIA, HERWIG and ARIADNE Monte Carlo event generators (see section 5.1) we do not utilise event generators which incorporate high-multiplicity matrix elements which were not yet tuned to our purpose. This choice is justified by the fact that generators used in our study use leading-order matrix elements combined with a leading-logarithm parton shower which do provide a satisfactory description of the three parton final states studied in this analysis as shown in [24,33].

Event Selection
The measurement of the strong coupling α S is based on the analysis of well measured hadronic events. A detailed description of the hadronic event selection can be found in [14]. Cuts applied to the events are based on the number of charged tracks, total visible energy and momentum imbalance. The cut on the momentum imbalance reduces events emitting a high energetic photon in the initial state, leading to a reduced hadronic centre-of-mass energy. The cuts on the number of charged tracks and total visible energy minimise to an insignificant level the number of events from hadronic tau decays and from hadronic final states originating from two-photon scattering.

Corrections to the Data
The corrections applied to the data follow exactly the same procedure as summarised in [14]. For the calculation of the three-jet rate all charged tracks and electromagnetic clusters are considered. The estimated minimum ionizing energy from tracks associated with electromagnetic calorimeter clusters is subtracted from the cluster energies. For the correction procedure two different categories of jet rate distributions are defined for simulated events, the so-called detector level and the hadron level distribution. The detector level distribution of simulated events is obtained by using all selected tracks of charged particles and electromagnetic clusters. The hadron level distribution is obtained by using the true four-momenta of stable particles, where particles with a lifetime of τ > 300 ps are declared as stable. Simulated events with photon initial state radiation (ISR) leading to the centre-of-mass energy reduced by more than 0.15 GeV are rejected from the hadron level 2 . Thus the correction for experimental effects explained below also takes care of residual effects due to initial state radiation. Only hadronic events originating from the primary production of u,d,s or c quarks are considered.
Before correcting the three-jet rate data distribution the expected contribution from e + e − → bb-events as expected from simulated events at the detector level is subtracted from the three jet-rate. About 1/11 of all qq-events are bb-events and the expected number of bb-events is subtracted from the observed number of data events at each y cut -bin. The hadronic events used in this analysis correspond to events with e + e − annihilating to a pair of u,d,s or c-quarks. The distribution corrected for e + e − → bb-events is then multiplied bin-by-bin with the ratio of the hadron level distribution divided by detector level distribution. A correction method based on a matrix unfolding method returns compatible results within statistical uncertainties [29]. The impact on the measurement due to changes of the b-quark fragmentation in the simulation are covered by the systematic uncertainty (see section 4) assigned to the correction for bbevents [29] . The numerical results of the corrected distributions are summarised in Tables 6 and 7.

Systematic Uncertainties
In order to assess the systematic uncertainty of the measurement of the strong coupling α S we evaluate several possible sources. For each variation the difference to the result with respect to the default analysis procedure is taken as a contribution to the systematic uncertainty. The systematic uncertainty is assumed to be symmetric around the default value. No systematic uncertainty is evaluated related to the fact that massless theoretical predictions are used since the contribution from e + e − → bb-events is subtracted from the data distribution (see section 3.4). The uncertainty originating from correcting for ISR effects is small and no systematic uncertainty is assigned.

Experimental Uncertainties
The assessment of the experimental uncertainty follows exactly the procedure described in detail in [14]. The analysis is repeated with a slightly modified event and track selection, using a different reconstruction software, using different MC models for the correction of detector effects and modified fit ranges. In addition the amount of subtracted e + e − → bbevents is modified by ±5%. The differences between the results obtained from the modified fits and the default fits are added in quadrature and taken as the combined experimental systematic uncertainty. The main contribution originates from using a different version of the reconstruction software and using HERWIG instead of PYTHIA for the correction of detector effects.

Hadronisation
The standard analysis uses PYTHIA to evaluate the change of the three-jet rate distribution originating from the transition from partons to stable hadrons (see section 3.2). Only the process e + e − to a pair of u,d,s or c-quarks is simulated. To assess the uncertainty the fit is repeated with alternative Monte Carlo models. For this we use HERWIG and ARIADNE instead of PYTHIA and the difference to the standard PYTHIA correction is taken as systematic uncertainty. In all cases the larger difference is seen by using the HERWIG Monte Carlo simulation instead of PYTHIA.
A variation of the parameters describing the hadronization model leads to significant smaller changes on the measurement of the strong coupling α S than applying a different hadronization scheme, like HERWIG [29]. For this reason we quote only the largest uncertainty coming from using a different hadronization model. As described in 3.2 we use different versions of the Monte Carlo generator for detector and hadronisation corrections. However, the hadronisation corrections for each generator are consistent within their statistical uncertainty.

Theoretical Uncertainties
The theoretical prediction of the three-jet rate distribution is a truncated asymptotic power series. The uncertainties originating from missing higher order terms are assessed by varying the renormalisation scale parameter x µ = µ/ √ s. For this x µ is set to two and to 0.5. The larger deviation from the fit using the default setting x µ = 1 is taken as systematic uncertainty. We assess the theoretical uncertainties by applying NNLO+NLLA+K QCD predictions in the fit, since the uncertainties using NNLO+NLLA predictions without K-term are found to be smaller. Effects from electroweak corrections are neglected.

Three-Jet Rate Distributions
The three-jet rates as a function of y cut at

Determination of α S
The value of the strong coupling is determined by a minimum-χ 2 fit of the matched NNLO+NLLA+K predictions to the corrected data distributions, separately for each centre-of-mass energy. The reduced renormalisation scale parameter x µ = µ/ √ s in the theoretical predictions is set to the natural choice x µ =1. The QCD predictions describe the three-jet rate at parton level only. To correct for hadronisation effects the matched QCD predictions are multiplied at each y cut point by the ratio of the hadron level distribution divided by the parton level distribution obtained from simulated events. The parton level distribution in simulation is obtained from the final state partons after the parton shower has terminated, i.e. just before the hadronisation step. The hadron level is defined in an identical way as described in section 3.4, containing only hadronic events with e + e − annihilating to a pair of u,d,s or c-quarks. The ratio between the hadron level divided by the parton level estimated with simulated events is shown in the appendix in Fig. 8. Using the QCD-prediction corrected for hadronisation effects and the corrected data distribution a χ 2 -value is calculated with i and j being the y cut points in the chosen fit range and the R(α S ) theo 3,i values are the predicted values of the three-jet rate. Each event can contribute to several points in the three-jet rate distribution leading to correlation between different y cut points. For this reason the covariance matrix V ij is not diagonal and the off-diagonal elements have to be computed. We follow the approach described in [14] using 1000 subsamples with 1000 simulated events each.
When choosing the fit range, we took the following considerations into account. The corrections applied to the data (see section 3.4) reverting the imperfectness of detector and the correction of the QCD-predictions due to hadronisation effects are required to be small. In addition we require the leading log contribution in the low y cut region of the three-jet rate distribution to be well below unity to ensure that the NLLA is valid. The leading log term is proportional toα S · ln 2 (1/y cut ) and requiring this term to be well below unity leads us to a lower limit of y cut = 0.01 assuming a value of the strong coupling α S (M Z 0 )=0.118. The upper limit is determined by requiring the leading order contribution A(y cut ) to be larger than zero. The corrections are small in this range with the detector corrections being less than 30% and the hadronisation corrections being less than 30%, apart from the corrections at the centre-of-mass energy of 14 GeV, which are up to 70%. The considerations described above lead to a fit range from 0.01 to 0.2, which is identical to the fit range used for the determination of α S using the differential y 23 event shape distribution [24]. The results of the fits to the three-jet rate distribution at the various centre-of-mass energies are shown in Fig. 3. The numerical results of the fits are summarised in Table 2. The statistical uncertainty is obtained from the fit and the systematic uncertainty is evaluated as described in section 4.
The fitted theory generally describes the data at hadron level well within the fit ranges, as seen in Fig. 3 and confirmed by the χ 2 /d.o.f. values in Table 2. The χ 2 /d.o.f. values are based only on the statistical uncertainties of the data and thus it is reasonable that they are larger than unity. The extrapolations outside of the fit regions also provide a reasonable description of the data.
For comparison the numerical results for fits using NLO+NLLA+K predictions are compiled in Table 3. The values of α S are consistently smaller and the theoretical uncertainty is about a factor of three larger compared to the fit using NNLO+NLLA+K calculations. The fit quality reflected by the χ 2 /d.o.f.-value is similar for both QCD predictions, but slightly better for the fit based on NNLO+NLLA+K predictions. The difference between the α S value returned with matched NLO and matched NNLO predictions is of similar size as the theoretical uncertainty evaluated for the fit using NLO+NLLA+K calculations. The reduction of the theoretical uncertainty using NNLO+NLLA+K QCD predictions is associated to the higher-order terms available in the new calculations.
The result obtained here cannot be compared directly to the measurement of α S using the identical data set and the y 23 event shape observable [24]. While for the α S measurement in [24] matched NNLO+NLLA QCD predictions are used, in this analysis NNLO+NLLA+K predictions are applied, which take subleading logarithms into account.

Combination of α S Measurements
The results of the measurements of the strong coupling α S at the various centre-of-mass energies are combined to a single value of α S (M Z 0 ). For this the values of α S ( √ s) obtained are evolved to a common energy scale M Z 0 and combined using a weighted mean. The theoretical uncertainty as well as the uncertainty originating from modelling the hadronisation process are likewise determined by calculating the weighted mean of the uncertainty  evaluated at each single energy point. The difficulties arise in estimating the correlations between the systematic uncertainties obtained at the different energy points. The identical problem is present for the combination of α S results from the LEP-collaborations and for this reason we use the same method as outlined in [14,24,34].

Simultaneous Variation of α S and the Renormalisation Scale
Besides fixing the renormalisation scale parameter to the natural choice x µ =1 we repeat the fit as a cross-check with both, α S and x µ , being varied within the minimisation procedure, thus choosing the so called optimal renormalisation scheme [36]. The results of the fits are summarised in Table 4. The is not used as input for the combined value of α S . The result shown at a centre-of-mass energy of 91.2 GeV (triangle) is obtained from a fit to the three-jet rate using NNLO predictions only with data taken by the ALEPH detector [13].
to a single value using the method described in 5.3. The combined value is:

Measurements of α S using NNLO Predictions only
To compare our result with results obtained with NNLO predictions only, we repeat the fit without NLLA matching. The fit range is chosen to be identical to the fit range used in the default fit. We perform two different fits using different choices for the renormalisation scale parameter: a fixed value of x µ =1 and the optimised renormalisation scheme, where x µ is an additional free parameter in the fit. The numbers in Table 5 show the result obtained at a centre-of-mass energy of 35.0 GeV, the energy point with the largest number of selected hadronic events. The results obtained at all energy points are summarised in the appendix in Table 8 and Table 9. Fig. 6 shows the fit result using NNLO predictions with x µ =1 and x µ =x opt µ at a centre-of-mass energy of 35.0 GeV. The shape of the three-jet distribution is not well matched by using NNLO predictions only, even at large y cut values. The χ 2 /d.o.f.-value obtained with x µ =1 increases significantly compared to the corresponding NNLO+NLLA+K fit. The fit using the optimised renormalisation returns reasonable χ 2 /d.o.f.-values and smaller x µ -values. A similar behavior was already observed using NLO predictions only with the renormalisation factor being set to x µ =1 and x µ =x opt µ [36,37] where variation of the renormalisation scale factor led to an improved description of the data as well. For both renormalisation scale schemes the scale uncertainty is considerably increased compared to the measurement using matched NNLO+NLLA+K predictions. In addition the difference between the fit using the natural and the optimised renormalisation scale is increased compared to a comparison using matched NNLO+NLLA+K predictions. A large sensitivity to the choice of the fit range is observed for a fit using NNLO predictions only. Again, the measurement of α S with data taken at a centre-of-mass energy of 14 GeV returns by far the largest We conclude, contrary to [13] which performs a fit to a single y cut -bin only and therefore being insensitive to the shape of the distribution, that resummation affects the fit significantly by decreasing the scale dependence and making the result of the fit more reliable. For this reason we consider the result based only on NNLO predictions as a cross-check only.

Summary
In this paper we present a measurement of the strong coupling α S using the threejet rate taken with the JADE experiment at centre-of-mass energies between 14 and 43.8 GeV. The three-jet rate is compared to simulated events obtained with PYTHIA, HERWIG and ARIADNE. All Monte Carlo models reproduce the measured distributions well. The strong coupling α S is measured by a fit to matched NNLO+NLLA+K predictions. The combined value using the measurements between 22 and 43.8 GeV results in α S (M Z 0 ) = 0.1199 ± 0.0060 (total error). The value is consistent with the world average [35]. The theory uncertainty has shrunk considerably by using matched NNLO+NLLA+K calculations and among the uncertainties it is the smallest. The dominant uncertainty originates from applying different Monte Carlo models to estimate the transition from partons to hadrons. A fit using NNLO predictions only with the renormalisation scale parameter x µ set to one cannot describe the shape of the three-jet rate distribution well. The shape can only be described well if α S and x µ are fitted simultaneously. Fig. 7 shows a comparison of the α S measurement presented in this paper with previous α S measurements using higher QCD predictions. The results obtained with a fit to the y 23 -distribution using data taken with the JADE [24] and the OPAL experiments [33] Fig. 7: Comparison of the α S measurement obtained in this analysis with previous α S measurements using the three-jet rate, the four-jet rate and the y 23 -differential distribution in e + e − -annihilation between 14 and 91 GeV using matched and unmatched QCD predictions. The 1 st measurement (top down) is the result of the present analysis, the 2 nd is the measurement using the differential two-jet rate distribution y 23 with NNLO+NLLA calculations and data taken between 14 and 44 GeV [24], the 3 rd is the measurement using the four-jet rate using matched NLO+NLLA+K predictions and data between 14 and 44 GeV [14], the 4 th is the measurement using the differential two-jet rate distribution y 23 with NNLO+NLLA calculations and data taken at 91 GeV [33] and the 5 th is a measurement using the three-jet rate at a centre-of-mass energy of 91 GeV using NNLO predictions only [13].
expected to be strongly correlated with the three-jet rate and are considered as a good cross check. The measured mean values, similar to the result obtained with this analysis, are slightly above the world average value. The measurements using the four-jet rate and data taken with the JADE experiment [14] and the three-jet rate using data taken with the ALEPH experiment [13] are consistent within the uncertainties. All measurements are in agreement with the world average value of α S [35].  0.140 ± 0.003 ± 0.008 0.152 ± 0.010 ± 0.015 0.126 ± 0.005 ± 0.007 −1.18 0.106 ± 0.002 ± 0.006 0.117 ± 0.009 ± 0.012 0.097 ± 0.005 ± 0.004 −1.05 0.077 ± 0.002 ± 0.006 0.088 ± 0.008 ± 0.009 0.072 ± 0.004 ± 0.005 −0.93 0.049 ± 0.002 ± 0.001 0.060 ± 0.006 ± 0.010 0.043 ± 0.003 ± 0.004 −0.81 0.028 ± 0.001 ± 0.003 0.034 ± 0.005 ± 0.007 0.023 ± 0.003 ± 0.001 −0.68 0.011 ± 0.001 ± 0.002 0.007 ± 0.002 ± 0.002 0.008 ± 0.002 ± 0.002 −0.56 0.002 ± 0.000 ± 0.000 0.006 ± 0.002 ± 0.009 0.003 ± 0.001 ± 0.002 Table 7: Hadron-level value of the three-jet fraction using the Durham algorithm at 35, 38.3 and 43.8 GeV. The value is corrected for contributions from e + e − → bb-events. In all cases the first quoted error indicates the statistical error while the second quoted error corresponds to the total experimental uncertainty. Uncertainties consistent with zero indicate that the corresponding value is smaller than the precision shown in the Table. B Hadronisation correction estimated with simulated events     Table 9: The value of α S and the statistical, experimental, hadronisation, renormalisation scale parameter x opt µ , the correlation between α S and the renormalisation scale parameter and the χ 2 /d.o.f. value of the fit using NNLO predictions. The fit is performed using the optimised renormalisation scheme for all energy points.