Determination of the Strong Coupling \boldmath{\as} from hadronic Event Shapes and NNLO QCD predictions using JADE Data

Event Shape Data from $e^+e^-$ annihilation into hadrons collected by the JADE experiment at centre-of-mass energies between 14 GeV and 44 GeV are used to determine the strong coupling $\alpha_S$. QCD predictions complete to next-to-next-to-leading order (NNLO), alternatively combined with resummed next-to-leading-log-approximation (NNLO+NLLA) calculations, are used. The combined value from six different event shape observables at the six JADE centre-of-mass energies using the NNLO calculations is $\alpha_S(M_Z)$= 0.1210 +/- 0.0007(stat.) +/- 0.0021(expt.) +/- 0.0044(had.) +/- 0.0036(theo.) and with the NNLO+NLLA calculations the combined value is $\alpha_S$= 0.1172 +/- 0.0006(stat.) +/- 0.0020(expt.) +/- 0.0035(had.) +/- 0.0030(theo.) . The stability of the NNLO and NNLO+NLLA results with respect to missing higher order contributions, studied by variations of the renormalisation scale, is improved compared to previous results obtained with NLO+NLLA or with NLO predictions only. The observed energy dependence of $\alpha_S$ agrees with the QCD prediction of asymptotic freedom and excludes absence of running with 99% confidence level.


Introduction
Analyses of events originating from e + e − annihilation into hadrons allow studies [2][3][4][5] of Quantum Chromodynamics (QCD), the theory of the strong interaction [6][7][8][9]. Comparison of observables like jet production rates or event shapes with theoretical predictions provides access to the determination of the strong coupling α S . Recently significant progress in the theoretical calculations of event shape observables has been made and next-to-next-to-leading order (NNLO) predictions are now available [10] as well as matching with resummed calculations in the next-toleading-log-approximation (NLLA) [11]. As a first application, measurements of α S at centre-of-mass-system (cms) energies between √ s = 91 GeV and √ s = 206 GeV were presented [12]. The same theoretical predictions are used in this paper to determine the strong coupling α S from JADE 1 data recorded at lower cms energies. As in [12] and the previous standard LEP and JADE analyses [5] we use Monte Carlo simulations to treat hadronisation effects. In [13] data for thrust at cms energies 14 ≤ √ s ≤ 207 GeV are analysed with combined NNLO+NLLA calculations and an analytic model for non-perturbative physics.
tion of α S [15]. The data correspond to a total integrated luminosity of 195/pb taken at average cms energies between 14.0 GeV and 43.8 GeV. The breakdown of the data samples, including the average cms energies, the energy ranges, data taking periods, integrated luminosities and the overall numbers of selected hadronic events are summarised in table 1. Monte Carlo events are generated in large numbers to correct the data for experimental acceptance, resolution effects and background. Events are simulated using PYTHIA 5.7 [16], and for systematic studies with HER-WIG 5.9 [17]. Subsequently the generated events are processed through a full simulation of the JADE detector and are reconstructed in the same way and with the same program chain as the data. For comparison with the corrected data and for the correction of hadronisation effects large samples of Monte Carlo events have been produced using PYTHIA 6.158, HERWIG 6.2 and ARIADNE 4.11 [18]. We use the model parameters as determined at √ s=91 GeV by the OPAL experiment at the LEP e + e − collider [19,20].
4 Experimental Procedure

Event Selection
The selection of identified and well measured hadronic event candidates follows the procedure outlined in [15]. Events with a large momentum imbalance due to photons emitted in the initial state are rejected. The event selection is based on minimal requirements for charged particle multiplicity, visible energy and longitudinal momentum imbalance. The dominating backgrounds from hadronic τ decays and two-photon interactions with hadronic final states are supressed to negligble levels.

Event Shape Distributions
The properties of hadronic events can be described by event shape observables. Event shape observables used for this analysis are thrust (1 − T ) [21,22], heavy jet mass (M H ) [23], wide and total jet broadening (B W and B T ) [24], C-Parameter (C) [25][26][27] and the transition value between 2 and 3 jet configurations [28,29] defined by the Durham jet algorithm (y D 23 ) [30]. Whenever we refer to a generic event shape observable 1 − T , M H , B T , B W , C or y D 23 we use the symbol y.
The event shape observables are calculated from selected charged particle tracks and calorimeter clusters after correcting for double counting of energy as described in [15]. We compared the data with the predictions from Monte Carlo simulations as described above and found good agreement at all energy points. Similar observations were made in [15] for related observables.

Corrections to the Data
The corrections to the data for limited experimental resolution, acceptance and bb background follows exactly the treatment described in [15]. Selected charged particle tracks as well as electromagnetic clusters are used to calculate the event shape observables.
From simulated events two different distributions are built: the detector-level distribution and the hadron-level distribution. The detector-level distributions are calculated exactly in the same way as for data using measured charged particle tracks and calorimeter clusters. The hadron-level distributions use the true four-momenta of the stable particles 4 in events where the centre-of-mass energy is reduced due to initial state radiation (ISR) by less than 0.15 GeV. The bin-by-bin corrections for the data distributions are derived from the ratio of hadron-level to detector-level distributions for simulated events with u, d, s or c primary quarks. Contributions from B hadron decays bias the measurement of event shape observables and therefore the expected contribution from e + e − → bb events is subtracted from the detector-level distributions before the corrections are applied. The simulations were optimised by OPAL to describe production and decays of B hadrons [19,20]. The good description of our uncorrected data by the simulations confirms that using the simulations to subtract the e + e − → bb background is justified.
In order to study systematic uncertainties the selection and correction procedures are modified and the whole analysis is repeated. The evaluation of the systematic uncertainties follows identically the procedure described in [15].

QCD Calculations
The distributions of the event shape observables are predicted by O(α 3 S ) (NNLO) perturbative QCD calculations [10]: withα S = α S (µ)/(2π). Equation (1) is shown for renormalisation scale µ = Q, where Q is the physical scale usually identified with the cms energy √ s for hadron production in e + e − annihilation. The coefficient distributions for leading order (LO) dA/dy, next-to-leading order (NLO) dB/dy and NNLO dC/dy were kindly provided by the authors of [10]. In [31] a problem at small values of y with the NNLO terms calculated in [10] was shown, but it does not affect the kinematic regions selected in our fits. The normalisation to the total hadronic cross section and the terms generated by variation of the renormalisation scale parameter x µ = µ/Q are implemented according to [10]. The prediction in equation (1) may be combined with resummed NLLA calculations [11] using the ln R-matching scheme; we refer to these predictions as NNLO+NLLA. The ln R-matching procedure ensures that in the combination of the fixed order NNLO calculations and the resummed NLLA calculations no double counting of common terms occurs. The NNLO+NLLA predictions are here compared with experimental data for the first time.
The theoretical predictions provide distributions at the level of quarks and gluons, the so-called parton-level. The distributions calculated using the final state partons after termination of the parton showering in the models are also said to be at the parton-level. To compare the QCD predictions with measured hadron-level event shape distributions the predictions are corrected for hadronisation effects. These corrections are obtained by calculating in the Monte Carlo models the ratio of the cumulative distributions at hadron-level and parton-level. The corrections are applied to the cumulative prediction R(y) = y 0 1/σdσ/dy ′ dy ′ as in [32]. It was shown in [12] that the event shape observable distributions derived from the parton-level of the Monte Carlo generators are described reasonably well by the NNLO calculation in their fit ranges. We compared the partonlevel predictions of the Monte Carlo generators with the QCD predictions in NNLO or NNLO+NLLA with α S (m Z 0 ) = 0.118 and x µ = 0.5, 1.0, 2.0 at √ s = 14, 22, 35 and 44 GeV. We study the quantity r(y) theory,MC = dσ/dy theory /dσ/dy MC − 1. In addition we compute the corresponding quantities r(y) MCi,MCj for any pair i, j of Monte Carlo predictions and r(y) xµ=0.5;1;2,xµ=0.5;1;2 for theory predictions with different renormalisation scale values. The largest values of the abs(r(y) MCi ,MCj ) and the abs(r(y) xµ=0.5;1;2,xµ=0.5;1;2 ) at each y are added in quadrature to define the uncertainty ∆r(y) of r(y) theory,MC . The average valuesr of the abs(r(y)) over the fit ranges (see below) are taken as a measure of the consistency between between theory and Monte Carlo predictions. The ratios ofr theory,MC with the average error ∆r(y) are generally smaller than or about equal to unity and reach values of about two for C at √ s = 44 GeV. The model dependence of the hadronisation correction and the renormalisation scale dependence of the theory will be studied as systematic uncertainties below. Our studies show that systematic uncertainties introduced by discrepancies between the theory predictions and the Monte Carlo parton-level predictions will be cov-ered by the combined hadronisation and theory systematic variations.

Measurements of α S
The strong coupling α S is determined by a χ 2 fit to each of the measured event shape distributions at the hadronlevel, i.e. corrected for experimental effects. A χ 2 value is calculated at each cms energy: where i, j count the bins within the fit range of the event shape distribution, d i is the measured value in the ith bin, t i (α S ) is the QCD prediction for the ith bin corrected for hadronisation effects, and V ij is the covariance matrix of the d i . The final prediction is obtained by integrating the QCD predictions in equation (1) over the bin width after application of the hadronisation correction as explained above in section 5.1. The χ 2 value is minimised with respect to α S while the renormalisation scale factor is set to The evolution of the strong coupling α S as a function of the renormalisation scale is implemented in three loops as shown in [5]. Since the evolution of α S in the cms energy range considered here does not involve the crossing of flavour thresholds it does not introduce significant uncertainties. In order to quantify the uncertainty from the evolution procedure we evolve α S (m Z 0 ) = 0.118 from m Z 0 to 14 GeV in three loops and two loops and find a relative difference of 0.1%.
In order to take the correlations between different bins into account the covariance matrix V ij is computed following the approach described in [33]: number of events after subtraction of background b i from bb events and multiplying by a correction α i for detector effects, P i is the normalised hadron-level distribution at bin i and N = k N k . The fit ranges are determined by several considerations. We require the leading log terms to be less than unity, because we also use a fixed order expansion without resummation of log-enhanced terms, see e.g. [34], and for consistency we use the same fit ranges in the NNLO and NNLO+NLLA analyses. The leading log term of dA/dy is ln y/y and we requireα S ln y/y < 1 for α S (m Z 0 ) = 0.118, 14 ≤ √ s ≤ 200 GeV and y > 0.1. The upper limit is given by the requirement that all three orders of the NNLO calculations contribute, i.e. the fit range extends to the kinematic limit of the LO coefficients dA/dy. The fit range for M H should be compared with the related fit range for 1−T after squaring, because M 2 H ∼ 1 − T in LO. The fit ranges for C and 1 − T are related by a factor of (ln 6)/6 ≃ 0.3. The resulting fit ranges are shown in table 2. The detector corrections (see section 4.3) are generally ±20% or less within the fit ranges. The hadronisation correction factors are maximally 3.5 for B T at √ s = 14 GeV but are generally smaller than 2 for the larger cms energy points. The hadronisation corrections are smallest and have the least variations over the fit ranges for y D 23 . The evaluation of the systematic errors of the α S measurements takes into account experimental effects, the hadronisation correction procedure and uncertainties of the theory. The three sources of systematic uncertainty are added in quadrature to the statistical error taken from the fits to obtain the total errors. Below we describe how we find the systematic uncertainties: Experimental Uncertainties The analysis is repeated with slightly varied event and track selection cuts and a systematic uncertainty from variation of the fit ranges is studied [15]. The cross section used in the subtraction of bb events is varied by ±5% which takes account of possible differences in the efficiency determination using the simulations of e + e − → qq (q=u,d,s,c) and e + e − → bb events. For each experimental variation the value of α S is determined and compared to the central (default) value. The quadratic sum of the differences and the fitrange uncertainty is taken as the experimental systematic uncertainty. Hadronisation For the default analysis, PYTHIA is used to estimate the corrections originating from hadronisation effects (section 5.1). As a systematic variation HERWIG and ARIADNE are used to evaluate the effects of hadronisation. The larger of the deviations is taken as systematic hadronisation uncertainty. It was observed in [35,36] that systematic uncertainties between the PYTHIA, HERWIG and ARIADNE models are generally much larger than systematic uncertainties from varying the parameters of a given model. Theoretical Uncertainties The theoretical prediction of event shape observables is a finite power series in α S . The uncertainties originating from missing higher order terms are assessed by changing the renormalisation scale factor to x µ = 0.5 and x µ = 2.0. The larger deviation from the default value of α S is taken as the systematic uncertainty.

Results from NNLO Fits
The results of the NNLO fits are summarised in table 3.
In figure 1  only while the combined experimental and hadronisation uncertainties are at least a factor of two larger than the statistical errors leading to a reduction of χ 2 /d.o.f. by a factor of at least 4 if these uncertainties were taken into account in the fits. We conclude that there is no significant disagreement between the event shape data and the QCD fits.
The results at each cms energy are remarkably consistent with each other, we find root-mean-square (rms) values for α S (

Results from NNLO+NLLA Fits
The results of the NNLO+NLLA fits are given in GeV, i.e. the scatter of individual results is essentially the same as for the NNLO analysis. The pattern of statistical errors and experimental and hadronisation uncertainties is the same as for the NNLO fits discussed above. Compared with the NNLO analysis the theoretical uncertainties are reduced by 10 − 20% and the values of α S are lower by 4% on average. The hadronisation uncertainties of the NNLO+NLLA fits are also smaller in most cases. The difference in α S between NNLO and NNLO+NLLA calculations is smaller than the difference in α S between NLO and NLO+NLLA calculations as expected in [11]. As discussed above in section 5.3 for NNLO fits the sometimes large χ 2 /d.o.f. values can be explained by the small statistical errors in some data sets.

Combination of Results
The results obtained at each energy point for the six event shape observables are combined using error weighted averaging as in [5,33,39,40]. The statistical correlations between the six event shape observables are estimated at each energy point from fits to hadron-level distributions derived from 50 statistically independent Monte Carlo samples. The experimental uncertainties are determined assuming that the smaller of a pair of correlated experimental errors gives the size of the fully correlated error (partial correlation). The hadronisation and theory systematic uncertainties are found by repeating the combination with changed input values, i.e. using a different hadronisation model or a different value of x µ . The results are given in table 5 and shown for the NNLO analysis in figure 2, be-cause this allows a direct comparison with the results of the NNLO analysis of ALEPH event shape data [12]. The statistical uncertainties of the combined results are reduced as expected. The systematic uncertainties of the combined results tend to be close to the best values from individual observables, because the systematic uncertainties are not completely correlated and because observables with smaller uncertainties have larger weights in the combination procedure.
A combination of the combined results at the six JADE energy points shown in table 5 after running to a common reference scale m Z 0 using the combination procedure described above results in α S (m Z 0 ) = 0.1210 ±0.0007(stat.) ± 0.0021(exp.) (4) ±0.0044(had.) ± 0.0036(theo.)  (α S (m Z 0 ) = 0.1172 ± 0.0051) for the NNLO+NLLA analysis. The NNLO+NLLA result has smaller hadronisation and theory uncertainties compared with the values in the NNLA analysis. We choose the latter result from NNLO+NLLA fits as our final result, because it is based on the most complete theory predictions and it has smaller theory uncertainties. It is consistent with the world average of α S (m Z 0 ) = 0.119±0.001 [41], the recent NNLO analysis of event shape data from the ALEPH experiment α S (m Z 0 ) = 0.1240 ± 0.0033 [12] as well as with the related average of α S (m Z 0 ) = 0.120 ± 0.005 from the analyses of the LEP experiments using NLO+NLLA QCD predictions [5]. The total error for α S (m Z 0 ) of 4% is among the most precise determinations of α S (m Z 0 ) currently available.
After running the fit results for α S ( √ s) for each observable to the common reference scale m Z 0 we combine the results for a given observable to a single value. We use the same method as above and obtain the results for α S (m Z 0 ) shown in table 6. The rms values of the results for α S (m Z 0 ) are 0.0029 for the NNLO analysis and 0.0026 for the NNLO+NLLA analysis; both values are consistent with the errors of the corresponding combined results shown in equations (4) and (5). Figure 3 shows the combined results of α S (m Z 0 ) for each observable together with results from alternative analyses discussed below. Combining the combined results for each observable or combining all individual results after evolution to the common scale m Z 0 yields results consistent with equation (5)   The results from the NNLO analysis of ALEPH data [12] are shown as well. In order to study the compatibility of our data with the QCD prediction for the evolution of the strong coupling with cms energy we repeat the combinations with or without evolution of the combined results to the common scale. We set the theory uncertainties to zero since these uncertainties are highly correlated between energy points. We conservatively assume the hadronisation uncertainties to be partially correlated, because these uncertainties depend strongly on the cms energy. The χ 2 probabilities of the averages for running (not running) with NNLO+NLLA fits then become 0.39 (9.9 · 10 −3 ). With the NNLO fits the χ 2 probabilities for running (not running) are 0.48 (1.2·10 −3 ). We interpret this as strong evidence for the dependence of the strong coupling on cms energy as predicted by QCD from JADE data alone.

Comparison with NLO and NLO+NLLA
For a comparison of our results with previous α S measurements the fits to the event shape distributions are repeated with NLO predictions and with NLO predictions combined with resummed NLLA with the modified ln Rmatching scheme (NLO+NLLA), both with x µ = 1. The NLO+NLLA predictions with the modified ln R-matching scheme were the standard of the final analysis of the LEP experiments [33,40,42,43]. The fit ranges as well as the procedures for evaluation of the systematic uncertainties are identical to the ones in our NNLO and NNLO+NLLA analyses.
The combination of the fits using NLO predictions returns α S (m Z 0 ) = 0.1301 ± 0.0009(stat.) ± 0.0029(exp.) ± 0.0054(had.)±0.0086(theo.), the fits using combined NLO+ NLLA predictions yield α S (m Z 0 ) = 0.1172±0.0007(stat.)± 0.0022(exp.) ± 0.0039(had.) ± 0.0054(theo.) and these results are shown in figure 3. The result obtained with the NLO+NLLA prediction is consistent with our NNLO and NNLO+NLLA analyses, but the theory uncertainties are larger by about 60%. The analysis using NLO predictions gives theoretical uncertainties larger by a factor of 2.6 and the value for α S (m Z 0 ) is larger compared to the NNLO or NNLO+NLLA results. It has been observed previously that values for α S from NLO analysis with x µ = 1 are large in comparison with most other analyses [44]. Both the NLO+NLLA or NNLO+NLLA analyses yield a smaller value of α S (m Z 0 ) compared to the respective NLO or NNLO results. The difference between NNLO+NLLA and NNLO is smaller than the difference between NLO+NLLA and NLO, since a larger part of the NLLA terms is included in the NNLO predictions.

Renormalisation Scale Dependence
The theoretical uncertainty due to missing higher order terms is evaluated by setting the renormalisation scale parameter x µ to 0.5 or 2.0. In order to assess the dependence of α S on the renormalisation scale the fits are repeated using NNLO, NNLO+NLLA, NLO and NLO+NLLA predictions with 0.1 < x µ < 10. The strong coupling α S as well as the χ 2 /d.o.f. as a function of x µ for 1 − T at √ s = 35 GeV are shown in figure 4. The χ 2 /d.o.f. curves for the NLO+NLLA and NNLO+NLLA fits show no local minimum in the x µ range studied. The α S (m Z 0 ) values from NLO predictions are the largest for x µ > 0.2. The α S (m Z 0 ) values using NLO+NLLA and NNLO calculation almost cross at the natural choice of the renormalisation scale x µ = 1 while the α S (m Z 0 ) value from the NNLO+NLLA fit is slightly lower. The NLLA terms at x µ = 1 are almost identical to the O(α 3 S )-terms in the NNLO calculation. A similar behaviour can be observed for B T and y D 23 . The slopes of the α S (m Z 0 ) curves of the NNLO and NNLO+NLLA fits around the default choice x µ = 1 are smaller than the slopes for the NLO and NLO+NLLA fits leading to the decreased theoretical uncertainties in our analyses.

Conclusion
In this paper we present measurements of the strong coupling α S using event shape observable distributions at cms energies 14 < √ s < 44 GeV. To determine α S fits using NNLO and combined NNLO+NLLA predictions were used. Combining the results from NNLO+NLLA fits to the six event shape observables 1 − T , M H , B W , B T , C and y D 23 at the six JADE energy points returns α S (m Z 0 ) = 0.1172±0.0006(stat.)±0.0020(exp.)±0.0035(had.)±0.0030(theo.), with a total error on α S (m Z 0 ) of 4%. The investigation of the renormalisation scale dependence of α S (m Z 0 ) shows a reduced dependence on x µ when NNLO or NNLO+NLLA predictions are used, compared to analyses with NLO or NLO+NLLA predictions. The more complete NNLO or NNLO+NLLA QCD predictions thus lead to smaller theoretical uncertainties in our analysis. The combined results for α S at each cms energy are consistent with the running of α S as predicted by QCD and exclude absence of running with a confidence level of 99%.