Precise determination of αS(MZ) from a global fit of energy–energy correlation to NNLO+NNLL predictions

We present a comparison of the computation of energy–energy correlation in e+e− collisions in the backto-back region at next-to-next-to-leading logarithmic accuracy matched with the next-to-next-to-leading order perturbative prediction to LEP, PEP, PETRA, SLC and TRISTAN data. With these predictions we perform an extraction of the strong coupling constant taking into account non-perturbative effects modelled with Monte Carlo event generators. The final result at NNLO+NNLL precision is αS(MZ ) = 0.11750 ± 0.00018(exp.) ± 0.00102(hadr.) ± 0.00257(ren.) ± 0.00078(res.).


Introduction
The strong interaction in the Standard Model (SM) is described by quantum chromodynamics (QCD) [1][2][3][4]. The theory successfully models the interactions between quarks and gluons and is a source of numerous predictions. Verifying the predictions of QCD is instrumental for searches for physics beyond the SM at the LHC, since the reliable prediction of SM processes as sources of backgrounds for searches is essential.
Precision measurements of event shape distributions in e + e − annihilation have provided detailed experimental tests of QCD and remain one of the most precise tools used for extracting the strong coupling α S from data [5,6]. Quantities related to three-jet events are particularly well suited for this task.
The state of the art for QCD for event shape observables currently includes exact fixed-order next-to-next-to-leading order (NNLO) corrections for the six standard three-jet event shapes of thrust, heavy jet mass, total and wide jet broadening, C-parameter and the two-to-three jet transition variable a e-mail: andrii.verbytskyi@mpp.mpg.de y 23 [7][8][9] as well as jet cone energy fraction [9], oblateness and energy-energy correlation [10]. The numerical matrix element integration codes described in the references allow the straightforward computation of any suitable, i.e. collinear and infrared safe event shape or jet observable.
However, fixed-order predictions have a limited kinematical range of applicability. For small values of an event shape observable y corresponding to events with two-jetlike topologies the fixed-order predictions do not converge well. This is due to terms where each power of the strong coupling α n S is enhanced by a factor (ln y) n+1 (leading logs), (ln y) n (next-to-leading logs) etc. For three-jet event shapes such logarithmically enhanced terms can be resummed at next-to-next-to-leading logarithmic (NNLL) accuracy [11][12][13][14][15][16][17], i.e. up to terms ∼ (ln y) n−1 . Resummation in nextto-next-to-next-to-leading logarithmic (N 3 LL) accuracy has been achieved for the C-parameter [18] and thrust [19]. A prediction incorporating the complete perturbative knowledge about the observable can be derived by matching the fixed-order and resummed calculations.
For the frequently used event shapes of thrust, heavy jet mass, total and wide jet broadening, C-parameter and y 23 , NNLO predictions matched to NLL resummation were presented in [20]. Predictions at NNLO matched to N 3 LL resummation are also known for thrust [12,19] and the C-parameter [18].
In this paper we consider the energy-energy correlation (EEC) in e + e − annihilation and present NNLO predictions matched to NNLL resummation for the back-to-back region. EEC was the first event shape for which a complete NNLL resummation was performed [11] while the fixedorder NNLO corrections to this observable were computed recently [10]. Moreover, EEC is the first event shape observable for which an analytic fixed-order NLO correction was computed [21].
The agreement between the predictions at NNLO+NNLL accuracy and the measured data is still not perfect. The discrepancy can be attributed mainly to non-perturbative hadronization corrections. We extract these corrections from data by comparison to state-of-the-art Monte Carlo predictions and determine the value of the strong coupling by comparing our results to measurements over a wide range of centre-of-mass energies. Our analysis allows us to target the highest precision of α S determination and we present the first global fit of the strong coupling to EEC at NNLO+NNLL accuracy. Our analysis also represents the first extraction of α S based on Monte Carlo hadronization corrections obtained from NLO Monte Carlo setups at NNLO+NNLL precision.

EEC distribution in perturbation theory
EEC is the normalized energy-weighted cross section defined in terms of the angle between two particles i and j in an event [22]: where E i and E j are the particle energies, Q is the centre-ofmass energy, θ i j = χ is the angle between the two particles and σ t is the total hadronic cross section. The back-to-back region θ i j → 180 • corresponds to χ → π , while the normalization ensures that the integral of the EEC distribution from χ = 0 • to χ = 180 • is unity. 1

Fixed-order and resummed calculations
The differential EEC distribution has been computed numerically at NLO accuracy in perturbation theory some time ago [24][25][26][27][28][29][30][31][32][33][34] and efforts towards obtaining an analytic result at this order [35,36] have culminated in a complete calculation very recently [21]. The NNLO prediction has also been obtained in ref. [10] using the CoLoRFulNNLO method [9,37,38]. At the default renormalization scale 2 of μ = Q the fixed-order prediction reads where A, B and C are the perturbative coefficients at LO, NLO and NNLO, normalized to the LO cross section for 1 Refs. [11] and [23] use the opposite convention of θ i j = 180 • − χ such that the back-to-back region corresponds to χ → 0 • . Here we use θ i j = χ throughout which agrees with the experimental convention. 2 We use the MS renormalization scheme throughout the paper.
e + e − → hadrons, σ 0 . In massless QCD this normalization cancels all electroweak coupling factors, and the dependence on the collision energy enters only through α S (Q). However, experiments measure the distribution normalized to the total hadronic cross section, so physical predictions must be normalized to σ t . The distribution normalized to the total hadronic cross section can be obtained from the expansion in Eq.
(2) through multiplying by σ 0 /σ t . For massless quarks, this ratio is independent of all electroweak couplings and reads The colour factors which appear above are given by while n f denotes the number of light quark flavours. The renormalization scale dependence of the fixed-order prediction can be restored using the renormalization group equation for α S and one finds wherē with x R = μ/Q. Finally, using three-loop running the scale dependence of the strong coupling is given by Here, t = ln(μ 2 /Λ 2 QCD ) and the β i are the MS-scheme coefficients of the QCD beta function, The fixed-order perturbative predictions diverge for both small and large values of χ , due to the presence of large logarithmic contributions of infrared origin. Concentrating on the back-to-back region χ → 180 • , these contributions take the form α n S log 2n−1 y, where y = cos 2 χ 2 .
As y decreases, the logarithms become large and invalidate the use of the fixed-order perturbative expansion. In order to obtain a description of EEC in this limit, the logarithmic contributions must be resummed to all orders. This resummation has been computed at NNLL accuracy in Ref. [11] 3 while in Ref. [40] a factorization theorem for EEC was derived based on soft-collinear effective theory which will allow to preform the resummation at N 3 LL accuracy once the corresponding NNLO jet function is computed. Since the complete jet function is currently not available, we use the NNLL results and formalism of Ref. [11] in the following. The resummed prediction at the default scale of μ = Q can be written as The large logarithmic corrections are exponentiated in the Sudakov form factor, The zeroth order Bessel function J 0 in Eq. (4) and b 0 = 2e −γ E in Eq. (5) have a kinematic origin. The functions A, B (not to be confused with the fixed-order expansion coefficients appearing in Eq. (2)) and H in Eqs. (4) and (5) are free of logarithmic corrections and can be computed as perturbative expansions in α S , Explicit expressions for the expansion coefficients (up to NNLL accuracy) in our normalization conventions can be found in Ref. [23]. It is possible to perform the q 2 integration in Eq. (5) analytically and the Sudakov form factor can be written as where a S = α S (Q)/(4π) and L = ln(Q 2 b 2 /b 2 0 ) corresponds to ln y at large b (the y 1 limit corresponds to Qb 1 through a Fourier transformation). Writing the Sudakov form factor this way clearly shows that S(Q, b) depends on its variables only through the dimensionless combination b Q. The functions g 1 , g 2 and g 3 correspond to the LL, NLL and NNLL contributions. Their explicit expressions can be found in Refs. [11,23].
So far, we have not considered the dependence of the resummed prediction on the renormalization scale. Besides the replacement of α S (Q) by α S (μ) in Eqs. (4) and (9), the resummation functions g i (λ) also acquire renormalization scale dependence, where the prime denotes differentiation with respect to λ. The factorization between the constant and logarithmic terms H (α S ) and S(Q, b) in Eq. (4) also involves some arbitrariness, since the argument of the large logarithm L can always be rescaled as 1. This arbitrariness is parametrized by x L , which plays a role in the resummed computation which is analogous to the role played by the renormalization scale in the fixed-order calculation. This rescaling of the logarithm introduces some modifications of the resummed formulae and the expansion coefficients in Eqs. (6)- (8). We find while the Sudakov form factor in Eq. (9) is also modified as follows

Matching the fixed-order and resummed predictions
In order to obtain a prediction which is valid over a wide kinematical range 4 the fixed-order and resummed calculations must be matched. Here we employ the log-R matching scheme as worked out for EEC in Ref. [23], and limit ourselves to recalling the final results.
In the log-R matching scheme for EEC we consider the cumulative distribution 4 We note that another resummation in the forward limit would be required to describe EEC over the full angular range.
The differential EEC distribution is easily recovered from Σ(χ, μ), The particular linear combination of moments introduced in Eq. (11) has the property that the divergence of the differential EEC distribution in the forward region (χ → 0) is suppressed by the factor of (1 − cos χ). Hence, in contrast to EEC itself, the fixed-order cumulative coefficients of Σ(χ, μ) can be computed reliably. Furthermore, one can show that in massless QCD this cumulative distribution is unity when χ = 180 • . Hence, we can integrate the fixedorder differential distribution in Eq. (3) and use the unitarity constraint Σ(π, μ)/σ t = 1 to all orders in α S to fix the constants of integration, Moreover, starting from Eq. (4) and using the definition of Σ, Eq.(11), we obtain the following expression for the resummed prediction 5 : where the Sudakov form factor S(Q, b) is the one given in Eq. (9). The final expression for the matched prediction was derived in Ref. [23] and reads 5 Note a misprint in Eq. (3.12) of Ref. [23] where an overall factor of 1/2 appears erroneously.
is the fixed-order expansion of the resummed result, The expansion coefficientsĀ res. ,B res. andC res. can be found in Ref. [23]. Notice that the function H (α S ) does not appear in Eq. (14) at all. In the log-R matching scheme such nonlogarithmically enhanced contributions should not be exponentiated, instead these terms, as well as subdominant logarithmic contributions, are all implicit in the unsubtracted parts of the fixed-order coefficientsĀ,B andC [41]. Thus the log-R matched prediction can be computed without the explicit knowledge of H (n) .
Finally, we comment on our implementation of the unitarity constraint Σ(π, μ)/σ t = 1. It can be shown that this constraint can be satisfied by modifying the resummation formula in Eq. (4) such that in the kinematical limit y = 1 the Sudakov form factor is unity. This may be achieved in several ways and here we choose a very simple solution and modify the resummation coefficientsÃ (n) andB (n) according tõ where p is a positive number. 6 This modification is fully legitimate since it does not modify the logarithmic structure of the result and introduces only power-suppressed terms. In practice, we set p = 1 and quantify the impact of this modification by comparing the results to those obtained with p = 2.

Finite b-quark mass corrections
The theoretical prediction presented above was computed in massless QCD. However, the assumption of vanishing quark masses is not fully justified, especially at lower energies, where b-quark mass effects are relevant at the percent level [43]. 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 con-tribution. Hence, we include mass effects directly at the level of matched distributions, Here 1 is the NNLO+NNLL matched distribution, computed in the log-R matching scheme in massless QCD as outlined above, while 1 is the fixed-order massive distribution. As the complete NNLO correction to this distribution is currently unknown, we model it by supplementing the massive NLO prediction of the parton level Monte Carlo generator Zbb4 [44], with the NNLO coefficient of the massless fixed-order result.
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 at NNLO accuracy (O(α 2 S ) in the strong coupling) including the exact bquark mass corrections at O(α S ) and the leading mass terms up [45]. We note that the electroweak coupling factors do not cancel in this ratio and the summation over quark flavours has to be carried out explicitly when computing σ massive (e + e − → hadrons).
were generated for each of the considered energies using a pole b-quark mass of m b = 4.75 GeV, which is consistent with world average estimations of pole mass 4.78 ± 0.06 GeV [46].
In order to assess the uncertainty associated to the modelling of b-quark mass corrections, we have investigated two alternative approaches for including them in our predictions. In approach A Eq. (17) is modified to i.e., we simply subtract the massless fixed-order NLO prediction multiplied by the fraction r b (Q) of b-quark events and add back the corresponding massive NLO distribution. Approach B is defined in a way very similar to our baseline, Eq. (17), but we do not include any NNLO corrections to the massive distribution, is simply the prediction obtained with Zbb4.
The criteria to include the data were high precision of differential distributions obtained with charged and neutral final state particles in the full χ range, presence of corrections for detector effects, correction for initial state photon radiation and sufficient amount of supplementary information. Therefore, data sets without supplementary information [59], with large uncertainties [60], superseded datasets [61,62] and measurements unfolded only to charged particles in the final state [63] are not included in the analysis.
The data sets selected for the extraction procedure have high precision and the measurements from different experiments performed at close energy points are consistent. 7 This justifies their use in the extraction procedure in a wide centreof-mass energy interval, similarly to studies of thrust [43] and C parameter [64] and allows us to target the highest precision of α S determination with available theoretical predictions.

Monte Carlo generation setup
In a previous study [23] the non-perturbative effects for the EEC distribution were modelled with an analytic approach. In this paper the non-perturbative effects in the e + e − → hadrons process are modelled using state-of-the-art particlelevel Monte Carlo (MC) generators. The non-perturbative corrections of the energy-energy correlation distributions were extracted as ratios of energy-energy correlation distributions at hadron and parton level in the simulated samples.
In GeV. In all cases, the simulation of initial state radiation was disabled and generator settings were defaults if the opposite is not stated explicitly. The value of the strong coupling used for the hard process was set to α S (M Z ) = 0.1181 [46].
To test the fragmentation and hadronization model dependence, the events generated with SHERPA2.2.4 were hadronized using the Lund string fragmentation model [72] or the cluster fragmentation model [73]. The first setup is labelled below as S L and the second as S C .
To assure proper fragmentation of heavy quarks and heavy hadron decays the cluster fragmentation model was adjusted. The value of SPLIT_LEADEXPONENT parameter was set to 1.0, the parameter M_DIQUARK_OFFSET was set to 0.55, 7 Some observed differences between the measurements performed at √ s = 91.2 GeV are not statistically significant once the systematical uncertainties and correlations are taken into account. 8 Partially updated to version 2.2.5. the production of charm and beauty baryons was enhanced by factors 0.8 and 1.7.
For the cross-check of SHERPA2.2.4 samples, the Herwig7.1.1 generator was used. The Herwig7.1.1 samples were generated with improved unitarized merging [74] using the MadGraph5 [75] matrix element generator and the GoSam [71] one-loop library to produce matrix elements of the e + e − → Z /γ → 2, 3, 4, 5 partons processes. The 2−parton final state processes again had NLO accuracy in perturbative QCD and the matrix elements were calculated assuming massive b-quarks. The merging parameter was set to √ s × 10 −1.25 √ s × 5.623 × 10 −2 . For the modelling of the hadronization process the default implementation of the cluster fragmentation model [76] was used. To improve the modelling of beauty production at the lowest energies, the b-quark nominal mass was changed from the default value of 5.3 GeV to 5.1 GeV. This setup is labelled below as H M .

Estimation of hadronization effects from MC models
Estimation of hadronization corrections is an integral part of comparing the parton-level QCD predictions to the data measured on hadron (particle) level. Despite the fact that under certain conditions the local parton-hadron duality leads to close values of quantities on parton and hadron level, the difference between them is not negligible and should be taken into account in precise analyses. One way to do so is to apply correction factors estimated from MC simulations to the perturbative QCD prediction. The factors, called hadronization corrections H/P, are defined as ratios of the corresponding quantities at parton level to the same quantities at hadron level at every point of the considered distribution.
To obtain the EEC distributions, the generated MC samples were processed in the same way as data (see e.g. Ref. [50]), using partons before hadronization for parton-level calculations and undecayed/stable particles for hadron-level calculations. For the parton-level calculations the parton energies were used as provided by the MC generators.
The predictions obtained with all setups describe the data well for all ranges of χ with the exception of the regions near χ = 0 • and χ = 90 • , for all values of √ s. For √ s < 29 GeV the H M setup is sensitive to the b-quark mass and the corresponding predictions are not reliable.
However, to assure an even better description of data, a reweighting procedure was applied to the simulated samples. The samples were reweighted at hadron level on an event-byevent basis to describe the data and the corresponding event weights were propagated to the parton level. The resulting distributions are shown in Fig. 1 for the SHERPA2.2.4 setups and in Fig. 2 for the Herwig7.1.1 setup.
As the SHERPA2.2.4 setups give the most stable and physically reliable predictions these are used in the analysis for reference hadronization corrections (S L ) and for system- atic studies (S C ). The corresponding hadronization corrections together with parametrizations are shown in Fig. 3 for reweighted samples.

Estimation of statistical correlations between measurements from MC models
To perform an accurate extraction procedure, the available data and uncertainties were examined and for every measured set of data a covariance matrix was built. The procedure consisted of multiple steps.
In the first step the systematic uncertainties were recalculated and separated from statistical uncertainties when this was possible. For the measurements with the uncertainties rounded to one significant digit [50,53,54] the uncertainties were expanded assuming maximal uncertainty before rounding. The measurements of TASSO [53] were converted from the dΣ/d cos χ form to dΣ/dχ using values of cos χ on the bin edges.
For the data from TOPAZ [52] the systematic uncertainty was calculated from an estimated relative systematic uncertainty of ±4% [52].
Taking into account the uncertainties of α S extraction analysis [58] the same was done for PLUTO [58] data. The systematic uncertainty estimation of ±5% for TASSO [53] is based on the upper limit of 10% for the total uncertainty mentioned in the paper [53]. The systematic uncertainties from DELPHI [49] and SLD [47] were used as provided in the original papers.
For all remaining data sets the published combined uncertainty was treated as statistical.
The measurements of Σ are provided in the original publications without correlations between the individual points. The correlation matrix was estimated from the Monte Carlo samples in terms of Fisher correlation coefficients [77,78]. Some of the obtained correlation matrices are shown in Fig. 4.
The obtained correlation coefficients are sizeable, up to 0.5 for the closest points, which highlights the importance of properly taking into account the correlations between measured points in the fits. The obtained correlation matrix together with statistical uncertainties was used to build a statistical covariance matrix for every data set.
To construct the systematic covariance matrix, the systematic uncertainties from the original publications were used with an assumption that these are positively correlated with correlation coefficient ρ = 0.5 between closest points. The correlations between the uncertainties of data from different experiments or different beam energies were neglected. The final covariance matrix used in the fit for every data set was a sum of statistical and systematic covariance matrices.

Fit procedure and estimation of uncertainties
The strong coupling extraction procedure is based on the comparison of data to the perturbative QCD prediction combined with non-perturbative (hadronization) corrections. The perturbative part of the predictions was calculated in every bin as described in previous sections. To tame the statistical fluctuations present in the obtained binned hadronization correction distributions, these were parametrized with analytic functions, expressed as a sum of polynomials of χ − 90 • . The value of the fitted function at the bin centre was used as the correction factor.
To find the optimal value of α S , the MINUIT2 [79,80] program was used to minimize where the χ 2 (α S ) value was calculated for each data set as with D standing for the vector of data points, P(α S ) for the vector of calculated predictions and V for the covariance matrix for D. The default scale used in the fit procedure was μ = Q = √ s. The fit ranges were chosen to avoid regions where resummed predictions or hadronization correction calculations are not reliable. The selected fit ranges were 117 • −165 • , 60 • −165 • and 60 • −160 • . The uncertainty on the fit result was estimated with the χ 2 + 1 criterion as implemented in the MINUIT2 program. The results of the fits are given in Tab. 2 for each fit range. In order to assess the impact of the NNLO corrections, in Tab. 2 we also present the results obtained using NLO+NNLL predictions. The NNLO corrections affect the fit in a moderate but non-negligible way. The obtained values of χ 2 divided by the number of degrees of freedom in the fit are of order unity for all cases. The corresponding distributions obtained from the fits for different √ s points are shown in Figs. 5, 6, 7 and 8. The systematic uncertainties of the obtained results were estimated with procedures used in previous studies [81]. To estimate the bias of the obtained result caused by the absence of higher-order terms in the perturbative predictions, the renormalization scale variation procedure was performed. In this procedure the fits were repeated, with variation of the renormalization scale in the range between x R = 1/2 and The bias of hadronization model selection is studied using the S L and S C setups of hadronization corrections, see results in Fig. 9. The bias related to the ambiguity of resummation scale choice was estimated by varying x L in the range between x L = 1/2 and x L = 2. To estimate the bias related to the ambiguity of our prescription implementing the unitarity constraint in the resummed calculation (see Eq. (16)), two values of p were used: p = 1 and p = 2. The difference between results obtained with two options is negligible. In all cases above the sizes of the biases were estimated numerically as half of the difference between the maximal and minimal α S value obtained in the corresponding set of  Fig. 9 it is seen that the estimated biases are relatively independent and, therefore combined in the final result as such.
Besides the estimations, several cross-checks of the obtained results were performed. First, the datasets were Here a 1 and a 2 are non-perturbative parameters that can be related in the dispersive approach to certain momentsᾱ q, p of the strong coupling α S [82]. These moments are the fit parameters of the analytic model. The results obtained from the fits with this setup are listed in Tab. 2. They show a high degree of dependence on the selected fit range, but are close to results obtained with the Monte Carlo based hadronization corrections in the range 117 • −165 • , see Tab. 2. Hence we conclude that away from the back-to-back region, the analytic model cannot fully account for hadronization effects.

Results and discussions
In this paper we presented the first combined analysis and extraction of α S at NNLO+NNLL precision from energyenergy correlation in electron-positron annihilation. Moreover, our analysis is the first extraction of the strong coupling based on Monte Carlo hadronization corrections obtained from NLO Monte Carlo setups at NNLO+NNLL precision. For the central value of the final result we quote the results obtained from the fits with the S L hadronization model in the range 60 • −160 • with uncertainties and estimations of biases obtained as described above. At NNLO+NNLL accuracy we obtain the best fit value of α S (M Z ) = 0.11750 ± 0.00018(ex p.) ± 0.00102(hadr.) ±0.00257(ren.) ± 0.00078(res.).
We see that the inclusion of the NNLO corrections has a moderate but non-negligible effect on the extracted value of α S . It has been explicitly checked that there are no correlations between estimated biases, therefore, the combined values with combined estimations of bias at NNLO+NNLL accuracy are: The value obtained from the analysis in NNLO+NNLL approximation is in agreement with the world average as of 2017 [83], 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 [83]. The estimated uncertainties are dominated by the uncertainty on the theoretical predictions. The results obtained in this study can be compared to those described in the original publications with NLO+NLL precision as well as the results obtained with analytic hadronization model in the sister paper [23].