High precision determination of $\alpha_s$ from a global fit of jet rates

We present state-of-the-art extractions of the strong coupling based on N$^3$LO+NNLL accurate predictions for the two-jet rate in the Durham clustering algorithm at $e^+e^-$ collisions, as well as a simultaneous fit of the two- and three-jet rates taking into account correlations between the two observables. The fits are performed on a large range of data sets collected at LEP and PETRA colliders, with energies spanning from $35$ GeV to $207$ GeV. Owing to the high accuracy of the predictions used, the perturbative uncertainty is considerably smaller than that due to hadronization. Our best determination at the $Z$ mass is $\alpha_s(M_Z) = 0.11881 \pm 0.00063(exp.) \pm 0.00101(hadr.) \pm 0.00045(ren.) \pm 0.00034(res.)$, which is in agreement with the latest world average and has a comparable total uncertainty.


Introduction
Measurements using hadronic final states in e + e − annihilation provide a unique opportunity to study the strong interaction in a well controlled environment without strongly interacting particles in the initial state. Accordingly, Quantum Chromodynamics (QCD) was tested extensively at LEP (see e.g. [1][2][3][4]).
QCD is a well-established theory by now, and its coupling constant α s has been measured in a variety of different processes at different energies. Still, the last PDG average of α s has an uncertainty of about 1% [5], which is considerably larger than errors in other gauge couplings. This uncertainty has an important impact on the current LHC precision physics program. Furthermore, a number of outlier fits of α s exist, hence any further precise determination of the value of the strong coupling constant is very valuable.
Many measurements of the strong coupling α s at e + e − colliders are based on comparisons of differential distributions of event shapes or jet rates to perturbative predictions. Presently progress in such measurements depends solely on the improvement in the accuracy of theoretical predictions as new data are not foreseen in the near future.
Compared to event shapes, jet-rates are known to be less sensitive to hadronization corrections 1 , hence more suited for precise determinations of the strong coupling constant. Fully differential predictions for the process e + e − → 3 partonic jets are available to NNLO (α 3 s ) accuracy [7][8][9][10][11]. Using the predictions for the three-jet rate at NNLO and the total cross-section at N 3 LO [12], the two-jet rate can be deduced at N 3 LO accuracy.
For a jet rate with a resolution parameter y, the fixed order predictions have a limited range of validity and are not reliable near the boundaries of the phase space, dominated by soft and collinear QCD radiation. In particular, for y → 0, the perturbative prediction at order O(α n s ) features logarithmic terms of the type α n s L m with L = ln(y), and m ≤ 2n. When L is large, such terms have to be summed up to all orders in perturbation theory.
The accuracy of the resummation for observables for which double logarithmic terms O(α n s L 2n ) exponentiate, is usually defined in terms of the logarithm of the cumulative distribution Σ. For such exponentiating observables, we define leading logarithms (LL) as terms of the form α n s L n+1 , next-to-leading logarithms (NLL) as α n s L n , next-to-next-toleading logarithms (NNLL) as α n s L n−1 , in ln Σ. The state of the art for most event shapes and jet rates measured at LEP is either NNLL or even N 3 LL [13][14][15][16][17][18][19][20][21][22][23][24]. These resummed predictions matched to fixed-order ones were used for precise extractions of α s [13,14,17,25,26].
Since the structure of logarithmic terms of predictions for jet rates is commonly more involved than that of most event shapes, until very recently, only next-to-double logarithmic corrections α n s L 2n−1 to these observables were known [27]. In this study, we focus on jet rates obtained with the Durham clustering algorithm [27] and will use for the first time NNLL predictions for the two-jet rate R 2 (y), which became available in Ref. [20].
The main result of this paper is an extraction of the strong coupling from the Durham two-jet rate to data measured at the LEP and PETRA colliders which relies on N 3 LO+NNLL accurate theoretical predictions. As an additional result, we also present a fit based simultaneously on the two-and three-jet rates, where the latter is computed at NNLO accuracy. For the first time, a fit of the coupling based on the two-jet rate features perturbative uncertainties that are considerably smaller than those of hadronization. This is due to the accurate N 3 LO+NNLL prediction adopted in the extraction.

Jet rates
A jet clustering algorithm is a procedure to classify final-state events into different jet multiplicities. This categorisation depends on the underlying algorithm used. In this paper we adopt the Durham clustering [27]. This is a sequential recombination algorithm, which requires a distance measure in phase space between momenta 2 p µ i and p µ j , where θ ij is the angle between the spatial components of the pair and E vis the visible energy in the event. If the smallest of these, y min = min({y ij }) is below a pre-defined number, y cut , then the corresponding pair of momenta is recombined into a single one. The procedure continues until all distance measures become larger than y cut . The momenta are recombined using some recombination scheme. Here we adopt the E-scheme [27], according to which the four-momenta of the two clustering particles are simply added together. The n-jet rate is then defined as where σ n-jet is the cross section for n-jet production in hadronic final states obtained with the above algorithm and σ tot is the total hadronic cross section. In the following sections we briefly review the fixed-order and resummed predictions used in this study. The two results can be combined by means of a matching procedure to obtain the predictions that we ultimately use in the fit.

Fixed-order predictions
In this work we use predictions obtained with the CoLoRFulNNLO method [11,28,29]. The perturbative expansion of the n-jet rate R n as function of y cut at the default renormalization scale µ ren = Q reads

3)
2 These momenta belong to either particles, or pseudo-jets obtained during the recombination process.
where A, B and C are the perturbative coefficients computed by the MCCSM code [30]. For massless quarks, these coefficients are independent of Q. The renormalization scale dependence of the fixed-order prediction can be restored using the renormalization group equation for α s , The numerical values for fixed-order coefficients for R 3 , R 4 and R 5 are reported in Appendix A, Tab. 3. These can be used to build up the necessary fixed-order predictions for R 2 that we use in the following.

Resummed predictions
The resummation technique adopted here was formulated in Refs. [19,20,23], hence we present only the main features of these results. The essence of the procedure described in Ref. [19] is that the NLL cross section is given by all-order configurations made of partons independently emitted off the Born legs and widely separated in angle [31]. The NNLL corrections are obtained by correcting a single parton of the above ensemble to account for all kinematic configurations that give rise to NNLL effects [20]. The two-jet rate at NNLL can be written as where the Sudakov radiator R NNLL and the coefficients H (1) , C hc are defined in Ref. [23] 3 , while the functions F NLL and δF NNLL are given in Refs. [19,20].
For the two-jet rate R 2 the resummation is performed with the ARES program [20] and the matching to fixed-order is done according to the ln R scheme [32].
The resummation of the three-jet rate is much more involved due to the extra number of emitting particles. Accordingly, the state-of-the-art resummed predictions have a much lower logarithmic accuracy and include only terms O(α n s L 2n ) and O(α n s L 2n−1 ) in R 3 (y) [27], in contrast to the logarithmic counting that we introduced for R 2 which refers to the logarithm of the jet rate. While R 3 is more sensitive to α s than R 2 , the low theoretical accuracy of the resummation does not guarantee a good theoretical control in the region where the logarithms are large. Therefore, for the present analysis, we do not perform any resummation of R 3 and limit the fit to a range where the fixed order is reliable.

Effects of quark masses
The effect of the non-vanishing b-quark mass on the predictions has also been considered in the literature. In particular, predictions for e + e − → partons are known including α 2 s corrections with massive b-quarks [33]. The resummed predictions for the Durham R 2 and R 3 observables with non-zero b-quark masses are only known at next-to-double logarithmic accuracy [34] (i.e. α n s L 2n−1 in the jet-rates) and are not used in this analysis as this does not match our target accuracy needed to guarantee a robust theoretical control even in the region where the logarithms become large.
In order to take b-quark mass corrections into account, we subtract the fraction of b-quark events r b (Q) from the massless result and add back the corresponding massive contribution computed at fixed order. Hence, we include mass effects directly at the level of the final distributions according to the formulae (2.7) The latter quantities are obtained by combining the total cross section at NNLO including mass corrections as obtained from Ref. [35], and the three-and four-jet rate O(α 2 s ) predictions as computed with the Zbb4 program [33,36]. The strong coupling used in the prediction of Zbb4 is then converted into a 5-flavour coupling by means of the matching relation for α s (m b ) [37,38].
We define the fraction of b-quark events as the ratio of the total b-quark production cross section divided by the total hadronic cross section, We evaluate the ratio of these cross sections to approximate O(α 3 s ) according to Ref. [35]. For the bottom quark we used a pole mass of m b = 4.78 GeV, which is consistent with the corresponding world average [39].

Determination of α s
To extract the strong coupling we compare the theory predictions described above to the available data, taking into account the non-perturbative corrections from Monte Carlo (MC) models, as described in detail in the following.

Data sets
We select experimental data that satisfy the following basic requirements: (i) measurements are obtained with both charged and neutral final state particles, (ii) corrections for detector effects have been taken into account, and (iii) corrections for initial-state QED radiation have been taken into account. We found that the data from Refs. [40][41][42][43][44][45][46][47][48] satisfy these basic criteria.
However, the measurements of Ref. [41] and Ref. [42] are superseded by the measurements in Ref. [43] and hence not included in our analysis. The analysis in Ref. [44] was excluded as the provided combined uncertainties are much lower than those from later refined analyses and are close to expected statistical uncertainties estimated from earlier analysis in Ref. [49]. We also excluded the measurements from Ref. [45] as these contain explicit corrections for the contributions of the process e + e − → bb from MC simulations and cannot be treated in the same way as the data without such corrections.
From the selected analyses we also excluded the measurements at √ s = 172 GeV from Ref. [46] as the background subtraction procedure performed with data of limited statistics has biased the measurement, e.g. introduced a non-monotonous behaviour of R 2 (y), see Ref. [50] for details. For similar reasons we exclude data from Ref. [47]. The sum of the rates for the data sets in Ref. [48] deviates from unity, and the largest deviation reaches almost 0.03 for the data set at √ s = 200 GeV. For this reason we excluded this data from the fits as well.
We summarise the information on the selected data sets in Tab. 1. For some data sets the uncertainties were updated before the fit, as follows. We added in quadrature the two available systematic uncertainties for the measurements at √ s = 91 GeV from Ref. [48]. The data from JADE [43] does not include systematic uncertainties related to the choice of MC samples used for the calculations of detector corrections. Later studies with the same data [45] indicated that such uncertainties can be at least as large as the statistical uncertainties. Therefore, in this case we added a relative uncertainty of 1.5% as an estimation of missing systematic uncertainties.
The measurements of the jet rates selected for the analysis are provided in the original publications without correlations between the individual points and without decomposition of the total systematic uncertainties. To perform an accurate extraction procedure, we examined the available data and uncertainties and built a covariance matrix for all measured sets of data. This procedure consists of multiple steps.
In the first step we estimated the statistical correlation matrix of the individual points in the fit range as described in Ref. [51] from the MC generated samples. We find that these are in good agreement with statistical correlation matrices obtained from data in unpublished Ref. [52]. In the second step we built the full covariance matrix for each measurement from the correlation matrix, statistical and systematic uncertainties.
we make the following assumptions on the correlation coefficients of the systematic uncertainties: We selected ρ = 0.006 in order to mimic patterns of systematic uncertainties observed in Ref. [52]. This corresponds to corr syst [R n (y 1 ), R m (y 2 )] ≈ 0.5K corr stat [R n (y 1 ), R m (y 2 )] for measurements with | log(y 1 ) − log(y 2 )| = 0.125, i.e. in the neighbouring bins for typical binning used in the measurements. K is equal to 1 if n = m and 0.5 otherwise. This approach approximates correlations between observables R 2 and R 3 .

Monte Carlo event generation setup
In our analysis we model non-perturbative effects in the e + e − → hadrons process using state-of-the-art particle level MC event generators. As usual, we estimate the nonperturbative corrections of the jet rate distributions by comparing distributions at hadron and parton level in the simulated samples. In particular, we used the Herwig7.1.4 [53] MC event generator to deliver our final results, and the Sherpa2.2.6 [54] MC event generator for cross-checks. We generated event samples for the process e + e − → hadrons at the centre-of-mass energies listed in Tab. 1. In all cases, we switched off the simulation of initial state radiation and used default generator settings, unless stated otherwise. We used α s (M Z ) = 0.1181 [39] as input for the generation. Moreover, we adopted the G µ scheme with input parameters M Z = 91.1876 GeV, M W = 80.379 GeV and G F = 1.1663787 × 10 −5 GeV −2 . The pole masses of b-and t-quarks were set to 4.78 GeV and 173 GeV respectively. The Herwig7.1.4 samples were generated with the unitarised MENLOPS method [55] using the MadGraph5 [56] matrix element generator and the OpenLoops [57] one-loop library to produce matrix elements for the 2, 3, 4 and 5 parton final states in the hard process. The two-and three-parton final state matrix elements have again NLO accuracy in perturbative QCD and the QCD matrix elements were also calculated using massive b-quarks. The merging parameter was set to √ s × 10 −1.25 . To test the fragmentation and hadronization model dependence, the events generated with Herwig7.1.4 were hadronized using either the cluster fragmentation model [58] or the Lund string fragmentation model [59]. The cluster fragmentation model is natively implemented in Herwig7.1.4. The Lund string fragmentation model is implemented in Pythia8.2.35 [60] and was used in Herwig7.1.4 via the ThePEG1.6.1++ [61] toolkit. For both setups the hadron decays were performed by the EvtGen1.7.0 [62] package. We label the first setup H C and the second one H L .
The Sherpa2.2.6 samples were generated with the MENLOPS method using the matrix element generators AMEGIC [63], COMIX [64] and the OpenLoops [57] one-loop library to produce matrix elements of the processes e + e − → Z/γ * → 2, 3, 4, 5 partons. The two parton final state matrix elements again have NLO QCD accuracy. The QCD matrix elements were calculated assuming massive b-quarks. The merging parameter Y cut was set to 10 −2.5 . The events were hadronized using the cluster fragmentation model [65] as implemented in Sherpa2.2.6. We label this setup S C .

Estimation of hadronization effects from MC models
The estimation of hadronization corrections is an integral part of comparing the partonlevel QCD predictions to the data measured at hadron (particle) level. While the principle of local parton-hadron duality leads to close values of observable quantities at parton and hadron level, the difference between them is not negligible and must be taken into account if one aims at a precise determination of α s . One possible way to do this is to apply correction factors estimated from MC simulations to the perturbative predictions.
To obtain the jet rates at parton and hadron level, we processed the MC generated samples in the same way as for the data, using undecayed/stable particles for hadron level calculations. To define jets with the Durham clustering, we used the implementation of the FastJet/fjcore3.3.0 package [66]. The selected resulting distributions are shown in Fig. 8 in Appendix B.
The predictions obtained with all setups describe the data well for all ranges of y, with the exception of the regions of small y at all values of √ s (see Figs. 9, 10 and 11). Among all considered MC setups, H L was selected to be the reference. The selected setup uses the very well tested Lund hadronization model that provides stable and physically reliable predictions throughout a wide range of centre-of-mass energies in e + e − collisions.
Moreover, precisely this hadronization model was used for the modelling of e + e − collisions in the original publications [43,46,48].
To estimate the uncertainty related to the hadronization modelling, we consider half the difference between the H L and H C setups, as explained in more detail in the following.
MC generators were already used to model hadronization effects in the previous QCD analyses of e + e − data [26,67]. Typically, hadron level predictions for every observable were obtained by multiplying the perturbative predictions by some factor derived from the analysis of MC generated events. In the present analysis the approach was amended to take into account physical constraints on R n , namely that jet rates are positive and that their sum should be one. These constraints are implemented by introducing the variables ξ 1 and ξ 2 such that at parton level The corresponding relations at hadron level are The functions δξ 1 (y) and δξ 2 (y) take into account the non-perturbative corrections. We obtained the numerical values for δξ 1 (y) and δξ 2 (y) from the MC simulated samples. In order to extract the hadronization corrections we proceed as follows. For a given y bin, we extract ξ 1 from the parton level prediction for the two-jet rate, and ξ 2 from the same prediction for the three-jet rate. The shifts δξ 1 and δξ 2 are then extracted in the same fashion from the hadron level predictions. The extracted values and corresponding interpolated functions (splines) are shown in Fig. 12 in App. C for selected centre-of-mass energies. The size of the derived hadronization corrections can be read off Fig. 13. As expected, we see that hadronization corrections increase at small values of the two-and three-jet rates and that the corrections become less important at higher energies.

Fit procedure
To find the optimal value of α s , we used the MINUIT2 [68] package to minimize χ 2 = data sets χ 2 (α s ) data set , where χ 2 (α s ) was calculated for each data set as where D stands for the vector of data points, P (α s ) is the vector of theoretical predictions, and V is the covariance matrix for the experimental data D.
We choose the fit range as follows. In order to assure that the implementation of the hadronization corrections is unitary, i.e. the constraint of eq. (3.2) is satisfied, we set the upper bound of the fit range below the kinematical limit for four jet production, log 10 (y) = log 10 (1/6) −0.8. We therefore choose log 10 (y) = −1 as an upper bound. Moreover, we adapt the lower bound to the centre-of-mass energy in order to take into account that hadronization corrections become more important at lower energies. Accordingly, we fix the lower bound log 10 y min (Q) of the fit range as log 10 y min (Q) = log 10 y min (M Z ) + L with L = log 10 (M 2 Z /Q 2 ). Different values used for log 10 (y min (M Z )) for R 2 and R 3 are indicated in the first columns in Tabs. 2, 4 and 5.
4.1 Fit of the coupling with the two-jet rate R 2 To obtain the most precise results, we first concentrate on fits that include the two-jet rate R 2 solely. The results of the fits, together with the used fit ranges, are given in Tab. 2. We show the results obtained via a fit of R 2 both at N 3 LO and N 3 LO+NNLL, with different hadronization models. The corresponding values of χ 2 divided by the number of degrees of freedom (ndof) in the global fit is also reported. The values χ 2 /ndof in the global fit as well as the values for every data set (not shown) are all of order unity, which can be viewed as a support for our correlation model. These values can be compared to the ones obtained in similar analyses. For instance, in Ref. [69], the χ 2 /ndof for fits with statistical uncertainties only varies for different observables between 0.5 and 60.
From the results given in the table one notices that the effect of the resummation is to move the fitted α s (M Z ) to slightly lower values. The quality of the fit, with or without resummation, is very similar. The benefit of including the resummation will become evident when perturbative uncertainties of the results are discussed.
As our best fit we quote the result obtained from the fits of the  Table 2. Fit of α s (M Z ) from experimental data for R 2 obtained using N 3 LO and N 3 LO+NNLL predictions, three different hadronization models and four different choices of the fit range, as given in the brackets, with L = log 10 (M 2 Z /Q 2 ). The reported uncertainty is the fit uncertainty as given by MINUIT2.

Estimation of uncertainties
The systematic uncertainties in α s are determined following the procedure of [70]. To estimate the size of higher-order terms in the perturbative prediction, we vary the renormalization scale in the range µ ren = Q/2 and µ ren = 2Q. Moreover, while keeping µ ren = Q fixed, we vary the resummation scale in the range µ res = Q/2 and µ res = 2Q. The effects of the individual variations are displayed in Fig. 4, where different hadronization models are compared. We notice that, when resummation is included, a much reduced dependence on the renormalization scale is observed.
The bias due to the selection of the hadronization model is studied by means of the difference between the H L and H C setups, see Fig. 4. In particular, considering the results from the Lund string and cluster hadronization models, the desired systematic uncertainty is obtained as half of the difference between the α S (M Z ) results obtained in nominal fits with H L and H C setups. The obtained numerical value is close to 0.001 (i.e. slightly below 1%). This can be briefly compared to previous estimations obtained with Monte Carlo event generator models in similar analyses. Namely, the values 0.001 [67,69] and 0.0005 [71] obtained previously allow us to validate our estimation. We stress that we have not performed any tuning of the adopted hadronization models to the data in order to artificially reduce the related uncertainties. This leads to a more conservative, and thus more robust estimate of hadronization uncertainties.
The final uncertainty is obtained by combining each of the above uncertainties in quadrature.

Validation of the procedure and further fits
In this section we perform some consistency checks to validate the fitting procedure used above. Moreover, we present an extraction of α s (M Z ) from a simultaneous fit of R 2 at N 3 LO+NNLL and R 3 at NNLO.

Fit consistency tests
We have performed a number of consistency tests, as outlined in the following: 1. We repeat the nominal fits in different ranges of √ s separately, instead of simultaneously. The results of these fits are shown in Fig. 5. We do not observe any significant dependence of the results on the centre-of-mass energy. This value is close to the reference result, but the hadronization uncertainty estimated from this setup is more sensitive to changes in the fit range.
3. To test the reliability of the correlation model used for the systematic uncertainties in the reference fit, we use the OPAL data and systematic shifts (uncertainties) from Ref. [52]. With this data we perform the fits using the χ 2 definition from eq.
The result is 6. We perform variations of the renormalization and resummation scales simultaneously, as done in Ref. [45]. The results are reported in Fig. 7 and show smaller uncertainty than in the case of independent variations. A similar behaviour was observed in Ref. [26].

5.2
Simultaneous fit of the coupling with the two-and three-jet rates R 2 and R 3 An alternative fit was performed with R 2 and R 3 observables simultaneously, with Unlike the results of the reference fits, the obtained result is sensitive to the selected fit range, see Tab. 4 in App. D. Taking this effect into account would result in another uncertainty of order 0.001 that is not included in the uncertainties given above.
As a final cross-check we perform a fit for a single point of ALEPH R 3 data at y = 0.02 without resummation. The obtained result (with statistical uncertainties only), α s (M Z ) = 0.11905±0.00251(exp.), can be compared to the results from Ref. [73], α s (M Z ) = 0.1175 ± 0.0020(exp.) ± 0.0015(theo.). The results agree well despite the differences in the implementation of the hadronization corrections.

Discussion
The value obtained from the analysis relying on N 3 LO+NNLL predictions for R 2 is in agreement with the world average as of 2017 [74], however it is visibly lower than the results from measurements performed for other e + e − observables using NNLO perturbative QCD predictions and MC hadronization models [74]. The estimated uncertainties are on the other hand approximately of the same sizes.

Summary
The main result of this paper is a first fit of the strong coupling for the two-jet rate that relies on very accurate (N 3 LO+NNLL) predictions. Our main result reads α s (M Z ) = 0.11881 ± 0.00063(exp.) ± 0.00101(hadr.) ± 0.00045(ren.) ± 0.00034(res.).
The uncertainty on α s induced by scale variations is now considerably smaller than that related to hadronization modelling. This is not the case if the fit is performed using only N 3 LO predictions, see Fig. 4. Furthermore the experimental uncertainty is now comparable to the perturbative one.
After combining the uncertainties in quadrature we obtain where the largest estimation bias comes from the used hadronization model. Our results agrees with the last PDG average with an uncertainty that is of the same size.
We have also performed a combined fit of the two-and three-jet rate, taking for the first time the correlation between these observables into account. The results of the two fits are fully compatible. However, the fit including R 3 shows a stronger dependence on the fit range than the results of our reference fit based on R 2 only. An accurate resummation for the R 3 observable could potentially reduce the sensitivity to the fit range selection and lead to an even more precise determination of α s (M Z ).

A Perturbative ingredients
We report in Tab. 3 the numerical results for the fixed-order coefficients introduced in eq. (2.5).  Table 3. The perturbative coefficients entering eq. (2.5) as calculated by the MCCSM program.

B MC simulations at hadron and parton levels
The figures in this section contain data measurements and predictions for jet rates calculated in very fine binning at parton (Fig. 8) and hadron (Figs. 9, 10 and 11) levels using different MC setups. The plots with ratios of hadron level predictions to data measurements are calculated with the binning used in the measurements. The error bars in the plots are calculated from the total uncertainty of the data points. The vertical lines in the plots show the fit ranges used in this work for fits of the R 2 and R 3 observables.     Figure 11. Predictions obtained with S C , H C and H L MC setups at hadron level. Fig. 12 shows the δξ 1 and δξ 2 distributions used to model hadronization corrections of jet rates at different centre-of-mass energies. As before, the vertical lines in the plots show the fit range for reference fits of the R 2 and R 3 observables. Fig. 13 shows the size of hadronization corrections obtained from the δξ 1 and δξ 2 values using eqs.

D Additional fits
In this section we present results of additional fits that we have performed. Table 4 shows a simultaneous fit of α s (M Z ) from experimental data for R 2 and R 3 obtained using N 3 LO and N 3 LO+NNLL predictions for R 2 and NNLO predictions for R 3 . Three different hadronization models are used and four different choices of the fit range are shown. The reported uncertainty is only the statistical uncertainty of the fit, as given by MINUIT2. We note that these results show a significant dependence on the fit range used.  Table 4. Simultaneous fit of α s (M Z ) from experimental data for R 2 and R 3 obtained using N 3 LO and N 3 LO+NNLL predictions for R 2 and NNLO predictions for R 3 . Three different hadronization models are used and four different choices of the fit range are shown, as given in the brackets, with L = log 10 (M 2 Z /Q 2 ). The fit range for R 2 is given in the first pair of brackets and the the fit range for R 3 in the second. The reported uncertainty is the fit uncertainty only, as given by MINUIT2.
We also show in Tab. 5 a fit of α s (M Z ) from experimental data for R 3 only obtained using NNLO predictions, three different hadronization models and four different choices of the fit range. As before, the reported uncertainty is the fit uncertainty only, as given by MINUIT2. These results can be compared directly with results of similar analyses, such as the study of Ref. [73]. Altogether we find good agreement between the fitted values reported in Tab. 5 and the results of Ref. [73].  Table 5. Fit of α s (M Z ) from experimental data for R 3 obtained using NNLO predictions, three different hadronization models and four different choices of the fit range, as given in the brackets, where L = log 10 (M 2 Z /Q 2 ). The reported uncertainty is the fit uncertainty only, as given by MINUIT2.