1 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) [24]. 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 [57].

Recently theoretical progress has been made leading to a significant improvement in the prediction of the three jet rate [811] as a function of α S. Previously \(\mathrm{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].

2 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:

$$ y_{ij}=\frac{2\mathrm{min}(E_i^2,E_j^2)}{E_{\mathrm{vis}}^2}(1-\cos \theta_{ij}), $$
(1)

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}^{\mu}= p_{i}^{\mu}+ p_{j}^{\mu}\), 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

(2)

with \(\hat{\alpha}_{\mathrm{S}} = \alpha _{\mathrm{S}}(\mu) / (2\pi)\) 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 \(\sqrt {s} \) for hadron production in \(\mathrm{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 \(( \hat{\alpha}_{\mathrm{S}} ) L^{2}\), where L=−lny cut, logarithmically enhances the higher-order corrections. For instance, \(( \hat{\alpha }_{\mathrm{S}} ) \ln^{2}(0.01)\approx \mathrm{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’ [1820], 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(\(\hat{\alpha}_{\mathrm{S}} ^{3}\)). 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

$$ R_3 = \sum_{i=1}^3 \sum _{j=i}^{2i} ( C_{\mathrm {F}} \hat{\alpha}_\mathrm{S} )^i R_3^{(i,j)} L^j + {\rm O}\bigl( \hat{\alpha}_\mathrm{S} ^4 \bigr)\,, $$
(3)

with the \(R_{3}^{(i,j)}\) coefficients are given in Table 1. The dependence on the renormalisation scale μ can be obtained by making the substitution

$$ \hat{\alpha}_\mathrm{S} \to \hat{\alpha }_\mathrm{S} \sum_{n=0}^2 \biggl( \beta_{0} \hat{\alpha }_\mathrm{S} \log \frac{\mu}{Q} \biggr)^n\,, $$
(4)

and keeping the terms up to O(\(\hat{\alpha}_{\mathrm{S}} ^{3}\)) in the expansion. In Eq. (4),

$$ \beta_{0} = \frac{11}{3} C_{\mathrm{A}} - \frac{2}{3} n_{\mathrm{f}} $$
(5)

with C A and C F being 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.

Fig. 1
figure 1

The QCD predictions are shown for the Durham three-jet rate calculated with a value of \(\alpha_{\mathrm{S}}(\mathrm{M}_{\mathrm {Z^{0}}})=0.1180\) and at a centre-of-mass energy of 35 GeV. In (a) the NNLO prediction is compared to the matched NNLO+NLLA and NNLO+NLLA+K predictions. The insert shows the ratio between the different predictions. In (b) the matched and unmatched NLO predictions are compared to the matched NNLO prediction with the insert presenting the ratio to the fixed order and the matched NLO predictions. For all calculations the uncertainty band reflects the uncertainty originating from setting the renormalisation scale factor x μ =μ/Q to 0.5 and 2

Table 1 Coefficients \(R_{3}^{(i,j)}\) of the expansion in Eq. (3) with \(K= (\frac{67}{18}-\frac{\pi^{2}}{6} )x-\frac {10}{9} y_{\mathrm{f}} \), \(x = \frac{ C_{\mathrm{A}} }{ C_{\mathrm{F}} }\), \(y_{\mathrm{f}} = \frac{ n_{\mathrm{f}} }{2C_{\mathrm{F}}}\) and b 0=β 0/C F

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),

(6)

where \(A^{{\rm NLL}}_{3} = C_{\mathrm{F}} \sum_{j=1}^{2} R_{3}^{(1,j)} L^{j}\), \(B^{{\rm NLL}}_{3} = C_{\mathrm{F}} ^{2} \sum_{j=2}^{4} R_{3}^{(2,j)}\) L j, \(C^{{\rm NLL}}_{3} = C_{\mathrm{F}} ^{3} \sum_{j=3}^{6} R_{3}^{(3,j)} 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=−lny 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.

3 Analysis procedure

3.1 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].

3.2 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 Sect. 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].

3.3 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.

3.4 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.Footnote 1 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 \(\mathrm{e^{+}e^{-}} \to \mathrm{b}\bar{\mathrm{b}} \)-events as expected from simulated events at the detector level is subtracted from the three jet-rate. About 1/11 of all \(\mathrm{q}\bar{\mathrm{q}}\)-events are \(\mathrm{b}\bar{\mathrm{b}}\)-events and the expected number of \(\mathrm{b}\bar{\mathrm{b}}\)-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 \(\mathrm{e^{+}e^{-}}\) annihilating to a pair of u, d, s or c-quarks. The distribution corrected for \(\mathrm{e^{+}e^{-}} \to \mathrm{b}\bar{\mathrm{b}} \)-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 Sect. 4) assigned to the correction for \(\mathrm{b}\bar{\mathrm{b}}\)events [29] . The numerical results of the corrected distributions are summarised in Tables 6 and 7.

4 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 \(\mathrm{e^{+}e^{-}} \to \mathrm{b}\bar{\mathrm{b}} \)-events is subtracted from the data distribution (see Sect. 3.4). The uncertainty originating from correcting for ISR effects is small and no systematic uncertainty is assigned.

4.1 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 \(\mathrm{e^{+}e^{-}} \to \mathrm{b}\bar{\mathrm{b}} \)-events 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.

4.2 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 Sect. 3.2). Only the process \(\mathrm{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 hadronisation model leads to significant smaller changes on the measurement of the strong coupling α S than applying a different hadronisation scheme, like HERWIG [29]. For this reason we quote only the largest uncertainty coming from using a different hadronisation model. As described in Sect. 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.

4.3 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_{\mu} = \mu/ \sqrt{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.

5 Results

5.1 Three-jet rate distributions

The three-jet rates as a function of y cut at \(\sqrt{s}= 14\), 22, 34.6, 35, 38.3 and 43.8 GeV are shown in Fig. 2. The rates are corrected for resolution and acceptance effects. The estimated contributions from \(\mathrm{e^{+}e^{-}} \to \mathrm{b}\bar{\mathrm{b}} \)-events are subtracted. The distributions are compared with the predictions obtained from Monte Carlo models used in this analysis. All models reproduce the data distribution well. Most of the data points are within the one sigma uncertainty band and a couple of points show a deviation of up to two sigma uncertainty. There is an apparent deviation of the simulation from the data in Fig. 2 (bottom left) for \(\sqrt{s}=38.3~\mathrm{GeV}\). This is due to the positive correlations between the data points which implies that a fluctuation in one data point is also visible in the neighbouring data points. For additional information inserts show the deviation of the three-jet rate obtained from Monte Carlo with respect to the data, normalised to the statistical and experimental uncertainty added in quadrature.Footnote 2

Fig. 2
figure 2

The three-jet rate as a function of the resolution parameter y cut is shown for centre-of-mass energies \(\sqrt{s}\) of 14.0 to 43.8 GeV. The data distributions are corrected for resolution and acceptance effects from the detector and the contributions from \(\mathrm{e^{+}e^{-}} \to \mathrm{b}\bar{\mathrm{b}} \)-events are subtracted. Data points are shown with statistical uncertainty (inner part) and combined statistical and experimental uncertainty. The histograms show the comparison to the predictions obtained with Monte Carlos simulation using PYTHIA, HERWIG and ARIADNE. The inserts show the deviation from the simulated distribution normalised to the combined statistical and experimental uncertainty

5.2 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_{\mu} =\mu / \sqrt{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 Sect. 3.4, containing only hadronic events with \(\mathrm{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

(7)

with i and j being the y cut points in the chosen fit range and the \(R( \alpha_{\mathrm{S}} )_{3,i}^{\mathrm {theo}}\) 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 Sect. 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 \(\hat{\alpha}_{\mathrm{S}} \cdot\ln^{2}(1 / y_{\mathrm{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 \(\alpha_{\mathrm{S}}(\mathrm{M}_{\mathrm {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 Sect. 4.

Fig. 3
figure 3

The fit results at centre-of-mass energies from 14.0 to 43.8 GeV are shown. The inserts show the deviation of the data points from the QCD-prediction with the α S-value obtained from the fit, normalised to the combined statistical and experimental error

Table 2 The value of α S using matched NNLO+NLLA+K predictions together with the statistical, experimental, hadronisation and theoretical uncertainties as described in Sect. 4. The last column shows the \(\chi^{2}/\mathrm{d.o.f.}\) value of the fit obtained with statistical uncertainties only at the respective energy points

The fitted theory generally describes the data at hadron level well within the fit ranges, as seen in Fig. 3 and confirmed by the \(\chi^{2}/\mathrm{d.o.f.}\) values in Table 2. The \(\chi^{2}/\mathrm{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 \(\chi^{2}/\mathrm {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.

Table 3 The value of α S using NLO+NLLA+K predictions together with the statistical, experimental, hadronisation and theoretical uncertainty as described in Sect. 4. The last column shows the \(\chi^{2}/\mathrm{d.o.f.}\) value of the fit obtained with statistical uncertainties only at the respective energy points

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.

5.3 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 \(\alpha_{\mathrm{S}}(\mathrm{M}_{\mathrm{Z^{0}}})\). For this the values of \(\alpha_{\mathrm{S}}(\sqrt{s})\) obtained are evolved to a common energy scale \({\mathrm{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]. For this analysis we combine the values obtained at centre-of-mass energies between 22 and 43.8 GeV, excluding the result obtained at 14 GeV because of the large hadronisation corrections and the large \(\chi^{2}/\mathrm {d.o.f.}\) returned by the fit. The combination of the results obtained between 22 and 43.8 GeV results in

consistent with the world average of α S=0.1184±0.0007 [35]. To the combined result the value of α S measured at 22 GeV contributes with a weight of 0.13, at 34.6 GeV with 0.23, at 35.0 GeV with 0.16, at 38 GeV with 0.07 and at 43.8 GeV with a weight of 0.41. The results at the various centre-of-mass energies are visualised in Fig. 4. Here, the values obtained at 34.6 and 35.0 GeV are combined to a single value using the same method as described above

The combination of the α S measurements at centre-of-mass energies between 22 and 43.8 GeV determined with NLO+NLLA+K predictions leads to

Fig. 4
figure 4

The open and solid points show the measurements of α S from the three-jet rate at the various centre-of-mass energies. The error bars indicate the statistical (inner part) and the total error. The lines show the world average of \(\alpha_{\mathrm{S}}(\mathrm{M}_{\mathrm{Z^{0}}})\) [35]. The results at 34.6 and 35.0 GeV are combined to a single value. The result obtained at 14 GeV (open point) 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]

5.4 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 χ 2-values obtained are almost identical to that from the default fit. The results for the scale \(x_{\mu}=x_{\mu}^{\mathrm{opt}}\) at the minimal \(\chi^{2}/\mathrm{d.o.f.}\)-value vary between 0.31 and 2.40, all being within the errors consistent with the x μ -range used in the systematic variation of the default fit. To estimate the theoretical uncertainty the fit is repeated with x μ being set to twice and half of the value obtained for the optimal scale \(x_{\mu}^{\mathrm{opt}}\). The α S-results at the various energy points between 22 and 43.8 GeV are combined to a single value using the method described in Sect. 5.3. The combined value is:

leading to almost the identical result as obtained with the default fit. The variation of α S and \(\chi^{2}/\mathrm{d.o.f.}\) with respect to the renormalisation scale is summarised in Fig. 5. The changes for both, the strong coupling α S as well as the \(\chi^{2}/\mathrm{d.o.f.}\)-value, are small within the x μ -range considered, indicating that our results are only moderately sensitive to missing higher order terms.

Fig. 5
figure 5

The variation of the strong coupling α S as a function of the renormalisation scale parameter x μ and the corresponding \(\chi^{2}/\mathrm{d.o.f.}\)-value for the data taken at a centre-of-mass energies from 14 to 43.8 GeV using the matched NNLO+NLLA+K predictions. The arrows indicate the variation of the renormalisation scale parameter used to evaluate the theoretical systematic uncertainty

Table 4 The value of α S using matched NNLO+NLLA+K predictions applying the optimised renormalisation scheme. The result is shown together with the statistical, experimental, hadronisation, theoretical uncertainties, the renormalisation scale parameter \(x_{\mu}^{\mathrm{opt}}\), the correlation between α S and \(x_{\mu}^{\mathrm{opt}}\) as well as the \(\chi^{2}/\mathrm{d.o.f.}\) value of the fit

5.5 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 Appendix C in Table 8 and Table 9. Figure 6 shows the fit result using NNLO predictions with x μ =1 and \(x_{\mu}=x_{\mu }^{\mathrm{opt}}\) at a centre-of-mass energy of 35.0 GeV.

Fig. 6
figure 6

The three-jet rate distribution at a centre-of-mass energy of 35 GeV is shown together with the NNLO QCD prediction evaluated at the α S-value obtained from the least-χ 2 fit (solid line). In addition the QCD-predictions obtained with the optimised renormalisation scale parameter (dash-dotted line) is shown. The QCD predictions are corrected for hadronisation effects. The insert shows the deviation of the data points from the QCD-prediction with the α S-value obtained from the fit, normalised to the combined statistical and experimental error

Table 5 The value of α S determined with NNLO predictions only together with the statistical, experimental, hadronisation and theoretical uncertainty as described in Sect. 4. The last columns show the value for x μ and the \(\chi^{2}/\mathrm {d.o.f.}\) The first row shows the fit result for a fixed renormalisation scale parameter x μ =1, the second row for x μ being varied within the fit. The χ 2-value is obtained with the statistical uncertainty only, taking bin-to-bin correlations into account

The shape of the three-jet distribution is not well matched by using NNLO predictions only, even at large y cut values. The \(\chi^{2}/\mathrm{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 \(\chi^{2}/\mathrm{d.o.f.}\)-values and smaller x μ -values. A similar behaviour was already observed using NLO predictions only with the renormalisation factor being set to x μ =1 and \(x_{\mu }=x_{\mu}^{\mathrm{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 \(\chi^{2}/\mathrm{d.o.f.}\)

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.

6 Summary

In this paper we present a measurement of the strong coupling α S using the three-jet 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 \(\alpha_{\mathrm{S}}(\mathrm{M}_{\mathrm {Z^{0}}}) =0.1199\pm0.0060~(\mathrm{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.

Figure 7 shows a comparison of the α S measurement presented in this paper with previous α S measurements using higher order 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] are 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].

Fig. 7
figure 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 \(\mathrm {e^{+}e^{-}}\)-annihilation between 14 and 91 GeV using matched and unmatched QCD predictions. The 1st measurement (top down) is the result of the present analysis, the 2nd 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 3rd is the measurement using the four-jet rate using matched NLO+NLLA+K predictions and data between 14 and 44 GeV [14], the 4th 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 5th is a measurement using the three-jet rate at a centre-of-mass energy of 91 GeV using NNLO predictions only [13]