Energy-energy correlation in electron-positron annihilation at NNLL+NNLO accuracy

We present the computation of energy-energy correlation in $e^+e^-$ collisions in the back-to-back region at next-to-next-to-leading logarithmic accuracy matched with the next-to-next-to-leading order perturbative prediction. We study the effect of the fixed higher order corrections in a comparison of our results to LEP and SLC data. The next-to-next-to-leading order correction has a sizable impact on the extracted value of $\alpha_{\mathrm S}(M_Z)$, hence its inclusion is mandatory for a precise measurement of the strong coupling using energy-energy correlation.


Introduction
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. Quantities related to three-jet events are particularly well suited for this task, since the deviation from simple two-jet configurations is directly proportional to α S . Furthermore, since the strong interactions occur only in the final state, non-perturbative QCD corrections are restricted to hadronization and power corrections. These may either be extracted from data by comparison to predictions by Monte Carlo simulations or computed using analytic models. Hence the precision of the theoretical computation is limited mainly by the accuracy of the perturbative expansion in α S .
In this regard, the state of the art 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 y 23 [1][2][3] as well as jet cone energy fraction [3], oblateness and energy-energy correlation [4].
However, fixed-order predictions have a limited kinematical range of applicability. For example when the two-jet limit is approached multiple emissions of soft and collinear gluons give rise to large logarithmic corrections that invalidate the use of fixed-order perturbation theory. In order to obtain a description appropriate to this limit, the logarithms must be resummed to all orders. For three-jet event shapes such logarithmically enhanced terms can be resummed at next-to-next-to-leading logarithmic (NNLL) accuracy [5][6][7][8][9][10][11] and even at next-to-next-to-next-to-leading logarithmic (N 3 LL) accuracy for thrust [12] and the C-parameter [13]. A prediction incorporating the complete perturbative knowledge about the observable can be derived by matching the fixed-order and resummed calculations.
For the standard event shapes of thrust, heavy jet mass, total and wide jet broadening, Cparameter and y 23 , NNLO predictions matched to NLL resummation were presented in ref. [14]. Predictions at NNLO matched to N 3 LL resummation are also known for thrust [6,12] and the C-parameter [13].
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 [5] while the fixed-order NNLO corrections to this observable were computed recently in ref. [4]. We also investigate the numerical impact of our results and perform a comparison of the most accurate theoretical prediction with precise OPAL [15] and SLD [16] data.
The paper is organized as follows. In section 2 we review the ingredients of our calculation, i.e., the fixed-order result as well as the resummation formalism. The matching of the NNLO predictions to the NNLL resummation is not entirely straightforward and we devote section 3 to a careful discussion of our procedure. We compare our results to LEP and SLC measurements in section 4 and in particular perform a fit to OPAL and SLD data. Finally, in section 5 we draw our conclusions.

Fixed-order and resummed predictions
EEC is the normalized energy-weighted cross section defined in terms of the angle between two particles i and j in an event [17]: where E i and E j are particle energies, Q is the center-of-mass energy, θ ij = π − χ is the angle between the two particles and σ t is the total hadronic cross section. Notice that the back-to-back region, θ ij → π corresponds to χ → 0, while the normalization ensures that the integral of the EEC distribution from χ = 0 • to χ = 180 • is unity.
The fixed-order prediction for EEC has been known in QCD perturbation theory up to NLO accuracy for some time [18][19][20][21][22][23][24][25][26][27][28] and has been computed at NNLO accuracy recently in ref. [4] using the CoLoRFulNNLO method [3,29,30]. At the renormalization scale µ 1 the result can be written as: whereĀ,B andC are the perturbative coefficients at LO, NLO and NNLO, normalized to the total hadronic cross section. In practice, our numerical program computes this distribution normalized to σ 0 , the LO cross section for e + e − → hadrons and at the fixed scale of µ = Q, At the default renormalization scale, the distribution normalized to the total hadronic cross section can be obtained from the expansion in eq. (2.3) by multiplying with with the color factors (2.6) 1 We use the MS renormalization scheme throughout the paper with the number of light quark flavors set to n f = 5. Furthermore, we use the two-loop running of αS for all predictions that incorporate a fixed-order NLO result, while predictions involving a fixed-order NNLO result are obtained using three-loop running.
Scale dependence can be restored using the renormalization group equation, (2.8) Finally, using three-loop running the scale dependence of the strong coupling is given by S and α S represent the one-, two-and three-loop contributions: . (2.11) The physical predictions for EEC up to NNLO accuracy are presented in figure 1 where the data measured by the OPAL collaboration is also shown. The bands represent the effect of varying the renormalization scale by a factor of two around its central value of µ = Q in both directions. Including the higher order corrections reduces the discrepancy between the predictions and data, although sizable differences remain. However, examining the region of intermediate χ (χ 30 • ) i.e., the region of validity of the fixed-order expansion, we observe that the LO scale variation band does not overlap with the NLO one, while the overlap of the NLO band with the NNLO one is marginal up to around χ ∼ 60 • , beyond which the two bands no longer touch. This behavior indicates that up to NLO the customary prescription for scale variation is not a reliable estimate of the size of the higher order corrections and casts some doubts also on the reliability of the NNLO band to estimate the perturbative uncertainty of the calculation. This phenomenon, however, is not unique to EEC and in fact very similar comments apply also to other three-jet event shapes in e + e − annihilation [1][2][3]. Nevertheless, one could argue that a more realistic estimate of the perturbative uncertainty could be obtained by considering a wider range for scale variation, see Ref. [31] for a careful discussion. This observation could explain, at least partially, the difference between the NNLO predictions and experimental data.
The fixed-order predictions clearly diverge for both small and large values of χ. As discussed above, this is the result of large logarithmic contributions of infrared origin. Concentrating on the back-to-back region (θ ij → π, i.e., χ → 0), these contributions take the form α n S ln 2n−1 y, where y = sin 2 χ 2 . (2.12) 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 the small angle limit, these logarithmic contributions must be resummed to all orders. The appropriate resummation formalism has been developed in refs. [32][33][34][35][36] and the coefficients which control this resummation are known completely at NNLL accuracy [5] 2 . At a center-of-mass energy of Q and renormalization scale µ the resummed prediction can be written as where the large logarithmic corrections are exponentiated in the Sudakov form factor, . (2.14) The Bessel function in eq. (2.13) and b 0 = 2e −γ E in eq. (2.14) have a kinematic origin. The functions A, B and H in eqs. (2.13) and (2.14) are free of logarithmic corrections and can be computed as perturbative expansions in α S 3 , It is possible to perform the q 2 integration in eq. (2.14) analytically and the Sudakov form factor can be written as where a S = α S (µ)/(4π) and L = ln(Q 2 b 2 /b 2 0 ) corresponds to ln y at large b (the y 0 limit corresponds to Qb 1 through a Fourier transform). The g i functions read 4 (See also ref. [11].) 3 Notice that our normalization conventions for A (n) , B (n) and H (n) , as well as βn in eq. (2.7) differ from ref. [5]. We follow the conventions of ref. [37]. 4 We note that the form of the gi functions is not affected by our different choice of normalization as compared to ref. [5], hence our expressions agree with those in ref. [5].
− ln The functions g 1 , g 2 and g 3 correspond to the LL, NLL and NNLL contributions. The expansion coefficients A (n) and B (n) appearing in eqs. (2.19)-(2.21) above were obtained in ref. [5] (see also refs. [11,37] for the complete NNLL A (3) coefficient). In our normalization conventions they read For B (1) and B (2) we have Finally, H (1) reads while the values of the higher loop coefficients H (n) (n > 1) are currently unknown.
Notice that the g i functions are singular as λ → 1. This singularity is related to the presence of the Landau pole in the QCD running coupling. Thus a prescription must be introduced to properly define the integral over b in eq. (2.13). Here we follow the same procedure as in ref. [5] and deform the contour of integration to the complex b-space as explained in refs. [38][39][40].
In figure 2 we present the resummed predictions for EEC up to NNLL accuracy together with OPAL data. Clearly, the resummed predictions are finite for χ → 0 and capture the general trends in the data in this limit. However, the purely resummed result badly underestimates the measured data already for moderate angles.

Matching
Fixed-order results are valid for moderate to large y (α S ln 2 y 1), while resummed results apply to small y (y 1). In order to obtain predictions over a wide kinematical range 5 the two computations must be matched. A number of different matching procedures have been proposed in the literature (see for example [41] for a review), but conceptually they all involve adding the two computations and subtracting the doubly counted terms: Here, the first term on the right hand side is the resummed result of eq. (2.13), the second term is the fixed-order prediction of eq. (2.2), while the last term is obtained by expanding the resummed component to the same order in α S as was used to compute the fixed-order result. Nevertheless, the subtraction of doubly counted terms alone is in general not sufficient to produce a physically sensible matched prediction. Indeed, in eq. (3.1) the resummed contribution on the right hand side is assumed to contain all logarithmically enhanced terms, hence the difference of the second and third terms, should be free of such contributions. However, unless the order of the logarithmic approximation is high enough to correctly reproduce the complete singular behavior of the fixed-order result as χ → 0, the difference in eq. (3.2) will contain non-exponentiated subleading logarithmic terms which make the matched distribution divergent at small χ. In contrast, the physical requirement is that the distribution should vanish at least as fast as a positive power of χ. The matching procedure is thus in general more involved than a simple subtraction of the terms that have been doubly counted.
For EEC the NNLL approximation is sufficient to reproduce all singular terms in the LO and NLO fixed-order differential distributions. This is no longer true at NNLO accuracy. (Similarly the NLL approximation will only reproduce the complete singularity structure of the LO fixed-order result.) Hence eq. (3.1) may be used to define a matched result at NNLL+NLO accuracy but not at NNLL+NNLO accuracy.
In order to obtain a matched prediction at NNLL+NNLO accuracy which behaves physically for small χ, the prescription of eq. (3.1) must be refined. This refinement has been worked out for EEC at NLL+NLO accuracy explicitly in ref. [42] and corresponds to what is commonly referred to as 'R matching' for event shapes [43]. However, this matching scheme has the drawback that some matching coefficients must be extracted from the behavior of the fixed-order result around χ → 0. Since the fixed-order calculation is particularly challenging in this region, the matching coefficients can only be extracted with large numerical uncertainties. This issue becomes more severe as we go to higher orders in the fixed-order computation. Hence we will not develop R matching for EEC at NNLL+NNLO. Nevertheless, the complete R matching formula does simplify to just eq. (3.1) when applied at NNLL+NLO accuracy, thus we will refer to NNLL+NLO predictions obtained with eq. (3.1) as R matched predictions below.
An alternative procedure for combining fixed-order and resummed results is log-R matching [43]. An attractive feature of this scheme is that all matching coefficients can be extracted analytically from the resummed calculation. In the log-R scheme one considers the cumulative event shape distribution, which we denote by the generic variable R(y, µ) for some given event shape y: This quantity has the following fixed-order expansion: where the fixed-order coefficientsĀ,B andC are obtained by integrating the corresponding differential distribution (e.g., eq. (2.2) for EEC) and using the constraint R(y max , µ) = 1 to all orders in α S in order to fix the constants of integration. (y max is the kinematically allowed maximum value of the variable y, for EEC χ max = 180 • .) The specific formulas for log-R matching in the literature [43] pertain to event shapes where the resummed prediction for the cumulative distribution can be written in a fully exponentiated form, (with L = ln y) where the C n and g n are known constants and functions. The g n functions can be expanded in powers of α S and L such that The log-R matching scheme simply amounts to taking the logarithm of eq. (3.4) and expanding it as a power series in α S yielding: and similarly rewriting eq. (3.5) as: (3.8) Removing the terms up to O(α 3 S ) from eq. (3.8) and replacing them by the O(α 3 S ) terms from eq. (3.7) yields the final expression for the log-R matched prediction at NNLL+NNLO: Notice that the constants C n do not enter eq. (3.9), since they are removed from eq. (3.8) and replaced by the fixed-order coefficients. Indeed, constant terms of the form C n α n S must be factorized with respect to the form factor [43] and should not be exponentiated.
For the case of EEC, the straightforward application of eq. (3.9) faces two difficulties. First, the resummed expression is not directly in the form of eq. (3.5). Second, EEC exhibits a particular problem because the fixed-order differential distribution diverges at both small and large χ, so that the cumulative coefficientsĀ,B andC cannot be reliably determined [15].
The latter complication can be conveniently solved by focusing on the following cumulative distribution: where we used the definition of y in eq. (2.12) to obtain the second equality. Hence, Σ(χ, µ)/σ t is just a linear combination of the zeroth and first moments of the differential EEC distribution. It is straightforward to reconstruct the original differential EEC distribution from the quantity defined in eq. (3.10): In addition, Σ(χ, µ) has the following properties. First, the singularity in the forward region (χ → π or y → 1) of the differential EEC distribution present in the fixed-order perturbative predictions is regularized by the factor of (1 + cos χ) which goes to zero in this limit. Second, it is not difficult to show that in massless QCD the value of this cumulative distribution is unity when χ = 180 • , i.e., Σ(χ max , µ)/σ t = 1. These properties together ensure that the fixed-order cumulative coefficients for Σ (defined in eq. (3.4) for a generic observable R) can be computed accurately by integrating the corresponding differential distribution and using Σ(χ max , µ)/σ t = 1 to all orders in α S to fix the constants of integration.
Furthermore, using eq. (2.13) in the definition of Σ, eq. (3.10), we find that the integration over y is straightforward to perform analytically and we obtain the following expression for the resummed prediction: where the Sudakov form factor S(Q, b) is unchanged and given in eq. (2.14).
Now let us turn to the issue of extending the definition of log-R matching, eq. (3.9), to the case of EEC. Although the resummed prediction for Σ is not directly in the form of eq. (3.5), one can repeat the constructions in eqs. (3.7)-(3.9) to arrive at the analog of eq. (3.9) for Σ. However, one must pay attention to the treatment of the non-logarithmically enhanced constant terms of the form H (n) α n S . These terms must not be exponentiated [43] and thus should not appear in the formula for the log-R matched expression, just as the C n constants are absent in eq. (3.8). Hence, we find the following expression for EEC in the log-R matching scheme 6 whereĀ res. ,B res. andC res. are the coefficients obtained by expanding the resummed component in the curly brackets above in a power series in α S : (3.14) The explicit expressions for these coefficients are somewhat long and we present them in appendix A. Eq. (3.13) is our final result for the log-R matched prediction for Σ at NNLL+NNLO accuracy. The log-R matched prediction at NNLL+NLO accuracy is obtained simply by dropping the O(α 3 S ) term in eq. (3.13). We emphasize that the quantities H (n) do not appear in eq. (3.13) at all. In the log-R matching scheme these terms, as well as subdominant logarithmic contributions are all implicit in the unsubtracted parts of the fixed-order coefficientsĀ,B andC [43]. Hence the log-R matched prediction can be computed without the explicit knowledge of H (n) .

Phenomenological results
In the following we investigate the numerical impact of the NNLO corrections and show quantitative predictions for the differential EEC distribution at NNLL+NLO and NNLL+NNLO accuracy.
We start by presenting the NNLL+NLO predictions. As discussed above, at this accuracy both R matching (eq. (3.1)) as well as log-R matching (eq. (3.13)) may be used to define physically sensible matched predictions. The results obtained with both the R matching and log-R matching schemes are shown in figure 3, together with the fixed-order NLO result. Throughout we set the center-of-mass energy to the Z-boson mass, Q = M Z = 91.2 GeV while the strong coupling is fixed to α S (M Z ) = 0.118. The fixed-order prediction is seen to diverge to −∞ as χ → 0. On We see that the two matching schemes give very similar results with the relative difference of the R matched prediction to the log-R matched prediction changing from about −2% at χ ∼ 0 • to 0% for χ ∼ 120 • . Around χ ∼ 180 • , the relative difference is about +0.5%.
Next, we include the NNLO corrections and present our NNLL+NNLO results. At this accuracy only the log-R matching scheme gives a prediction which behaves physically at small χ, and this prediction is shown in figure 4. Here too, the center-of-mass energy was set to Q = 91.2 GeV and we used α S (M Z ) = 0.118. The fixed-order prediction at NNLO accuracy diverges to +∞ as χ → 0. (The divergence to +∞ is not visible on the plot where the NNLO result seems to diverge to −∞.) As previously, the matched prediction is well-behaved down to very small values of χ. The bottom panel shows the ratio of the fixed-order prediction to the matched result. The band again represents the effect of varying the renormalization scale by a factor of two in either direction around its central value of µ = Q. In this case, three-loop running of the strong coupling is used.
To better appreciate the impact of the NNLO corrections on the matched prediction, in figure 5 we compare the NNLL+NLO and NNLL+NNLO results obtained in the log-R matching scheme.  The lower panel shows the ratio of the NNLL+NLO prediction to the NNLL+NNLO one. We see that the inclusion of the NNLO corrections slightly lowers the prediction in and below the peak region by from −5% to −2%, while the prediction is enhanced for medium and high values of χ. This enhancement goes from around +7% at χ = 60 • to around +14% at χ = 120 • and up to +20% to +25% for values of χ near 180 • . Hence the inclusion of NNLO corrections has a sizable impact on the shape of the distribution.
Next, we compare our predictions to precise OPAL [15] and SLD [16] data. In particular, we perform a fit of our most accurate NNLL+NNLO prediction to the experimental data with the strong coupling α S as a free parameter. We use a χ 2 analysis for the fitting procedure. In general, both statistical and systematic errors are correlated between bins, but unfortunately the experimental publications provide practically no information on the correlations. Therefore, we simply add statistical and systematic uncertainties in quadrature and treat them as uncorrelated between all data points. In order to quantify the impact of the NNLO corrections, we perform the same fit also for the NNLL+NLO prediction computed in both the R and the log-R matching schemes. Since we cannot take into account correlations properly, we do not aim to produce the most accurate extraction of α S (M Z ) here, but rather to assess the impact of the NNLO corrections on the extraction.
In a first attempt, we neglect hadronization corrections, however, we come back to this point Fit range below. Obviously, the results obtained in this way must be interpreted with care.
In order to ease the comparison of our results to previous work, we choose the fit ranges of ref. [5] for our analyses. In the first case we include data in the range 0 • < χ < 63 • where the effects of resummation are rather pronounced. However, given that the low χ region is particularly sensitive to non-perturbative corrections, we also investigate the range 15 • < χ < 63 • where the lower cut is expected to mitigate the effects of these contributions. Finally, we perform fits including data in the 15 • < χ < 120 • interval. Since for large χ the matched prediction is controlled by the fixed-order result, we expect the effects of the NNLO correction to be most prominent here. The upper limit is chosen to cut the forward region where another resummation would be required.
We collect the results of these fits in table 1 where we show the best fit value of α S and the χ 2 /d.o.f. for each fit. The quoted uncertainty on the extracted value of α S is obtained by adding the fit uncertainty and the theoretical uncertainty from missing higher-order contributions in quadrature. We assess the latter by repeating the fit with several values of µ in the range Q/2 ≤ µ ≤ 2Q and taking the envelope of the obtained results. This theoretical contribution is dominant in the total uncertainty.
In the first case, when we include data in the 0 • < χ < 63 • interval, we observe that the quality of the fit is actually better for the NNLL+NLO predictions than the theoretically most accurate NNLL+NNLO prediction, with the latter fit being rather poor as evidenced by the high value of χ 2 /d.o.f. = 206.4/50 = 4.13. The extracted values of α S (M Z ) are generally quite high compared to the world average α S (M Z ) = 0.1181 ± 0.0011 [44], in the range α S (M Z ) = 0.129 -0.133, with the fit using the NNLL+NNLO prediction giving the smallest value. However, as mentioned above, non-perturbative corrections are the largest in the low-χ region, hence the results of fits using purely perturbative predictions should be interpreted with great care, especially in this fit range.
Indeed, repeating the fit in the 15 • < χ < 63 • interval, we observe that the χ 2 /d.o.f. decreases in each case and the change is particularly significant for the fit using the NNLL+NNLO prediction, where we find a much more reasonable value of Our extracted values of α S (M Z ) based on the NNLL+NLO predictions using R matching are quite close to the values obtained in ref. [5] for all three fit ranges, although our results are marginally higher. We have checked that these differences are due to the fact that the determinations in ref. [5] used the incomplete A (3) NNLL resummation coefficient.
Overall, we observe that the inclusion of the fixed-order NNLO corrections reduces the extracted value of α S (M Z ). This reduction is about −2% to −3% when data in the range 0 • < χ < 63 • are taken into account, about −2% to −4% for the range 15 • < χ < 63 • and between −5% to −7% when 15 • < χ < 120 • , depending on the matching prescription used for the NNLL+NLO prediction. Hence, these corrections must be included in a precise determination of α S using EEC.
In our analysis so far, we have neglected hadronization corrections. However, non-perturbative contributions are expected to be relevant, especially at small angles [17,[32][33][34]45], and indeed the OPAL analysis of ref. [15] found hadron-parton correction factors from around 1.5 for very small χ to around 0.9 for large χ 7 . Hence it is important to account for these non-perturbative contributions. As already mentioned, these can be determined either by extracting them from data by comparison to Monte Carlo predictions, or by performing analytic model calculations. Here, we follow the latter option and use the non-perturbative model of ref. [42] to describe the hadronization contributions. Thus we multiply the Sudakov form factor of eq. (2.14) with a correction of the form and treat a 1 and a 2 as free parameters of the non-perturbative model to be fitted from data.
We have performed a three-parameter fit including data in the 0 • < χ < 63 • range using our NNLL+NNLO prediction, as well as the predictions obtained at NNLL+NLO with both R matching and log-R matching. In the R matching scheme at NNLL+NLO accuracy, we extract the following parameters: 2) with χ 2 /d.o.f. = 38.7/48 = 0.81. All uncertainties are again obtained by adding the fit uncertainties and the theoretical uncertainties in quadrature. The theoretical uncertainties are assessed by varying the renormalization scale µ between Q/2 and 2Q and repeating the fit. The total uncertainties are mostly dominated by the theoretical uncertainty with the exception of the upper limit of strong coupling. In this case, we find that the maximal best fit value of α S is obtained for µ Q, hence the upper limit is controlled by the fit uncertainty. We also report the correlation matrix of Evidently the strong coupling α S is highly anti-correlated with the non-perturbative parameter a 2 .
The analysis of ref. [5] performed on the same data gave |a 2 | 0.002 GeV, a very small value compatible with a 2 = 0. After fixing the parameter a 2 to zero, a two-parameter fit to the strong coupling and the remaining non-perturbative parameter a 1 produced the best fit values of α S (M Z ) = 0.130 +0.002 −0.004 and a 1 = 1.5 +3.2 −0.5 GeV 2 with χ 2 /d.o.f. = 0.99. Our results in eq. (4.2) are compatible with these values within uncertainties. We have nevertheless verified that the source of the discrepancy between the two extractions is, again, due to the fact that ref. [5]  The strong coupling α S and the a 2 non-perturbative parameter is even more strongly anti-correlated than in the R matching scheme. As before, the uncertainties in eq. (4.4) include the fit and theoretical uncertainties added in quadrature. We observe that the quality of the fits as measured by χ 2 /d.o.f. is very similar in the two matching schemes and the fit results are compatible between the two schemes within uncertainties. The extracted value of the strong coupling is reduced by about −5% in the log-R scheme compared to the R scheme, however, it remains high compared to the world average in both schemes.
We present the comparison of the best fit NNLL+NLO predictions in the R and log-R matching schemes to the data in figure 6. The figure shows a nice overall agreement between the predictions and experiment and it is clear that the calculations can reproduce the measurements up to the smallest measured values of χ. Nevertheless, we observe a small but systematic deviation of the prediction from data in the region of medium χ (from about χ 30 • ) and it is clear that the shape of the measured distribution is not fully reproduced. The bottom panel shows the ratio of the data and the R matched prediction to the log-R matched result, with the bands representing scale uncertainty.
Finally, we investigate the impact of NNLO corrections and repeat the three-parameter fit in the same range of 0 • < χ < 63 • , but using our most accurate NNLL+NNLO theoretical prediction. We see that the quality of the fit improves drastically compared to the purely perturbative fit reported in table 1. Moreover, the extracted value of α S (M Z ) is sizably reduced compared to the fits based on NNLL+NLO predictions and is indeed compatible with the world average within uncertainties. Figure 7 shows the comparison of the best fit NNLL+NNLO result to the measured data. We again observe that the measurement is very well described by the theoretical prediction and, in particular, the impact of the NNLO correction is clearly visible in the medium χ range, where the agreement between the data and the prediction is now excellent. The systematic deviation which is present in the NNLL+NLO predictions in this range is completely erased when the NNLO correction is taken into account. At the same time the best fit value of α S (M Z ) is shifted by about −6%. We conclude that the inclusion of the fixed-order NNLO correction is essential for a precise determination of α S from EEC. Finally, the three-parameter fits show that in this approach to hadronization corrections, the non-perturbative parameter a 1 is more important than a 2 . As stressed already in ref. [5], this indicates that the parametrization in eq. (4.1) is not able to fully describe the non-perturbative corrections, especially at medium and large χ. Hence, part of the hadronization effects are absorbed into the strong coupling. This is also apparent from the very strong anti-correlation in the fits between α S and the non-perturbative parameter a 2 . Thus, it would be very interesting to repeat our analysis with hadronization corrections extracted from data by comparison to Monte Carlo simulations. The results of such an analysis will appear elsewhere [46].

Conclusions
In this paper we presented precise QCD predictions for the energy-energy correlation in e + e − collisions. Our computation includes fixed-order perturbative corrections up to NNLO accuracy, as well as a resummation of the logarithmically enhanced terms in the back-to-back region at NNLL accuracy. In order to obtain a description which incorporates the complete perturbative knowledge about the observable and is valid over a wide kinematical range, the fixed-order and resummed predictions must be matched. We have implemented this matching in the R scheme at NNLL+NLO and also, for the first time, in the log-R scheme at both NNLL+NLO and NNLL+NNLO accuracy.
All of our matched results satisfy the physical requirement that the EEC distribution should vanish as χ → 0.
We also presented perturbative predictions at NNLL+NLO and NNLL+NNLO accuracy and compared these to precise OPAL and SLD data. In particular, we have performed a fit of our results to the data with the strong coupling α S as a free parameter. Using an analytic model to account for hadronization corrections, we obtain a very good description of the data down to the smallest measured angles. We observe that the inclusion of the NNLO corrections has a significant impact on the extracted value of α S (M Z ), shifting the best fit value by around −6% compared to the NNLL+NLO computation. Hence, the inclusion of these corrections in a precise measurement of α S from EEC is mandatory.
Using our most accurate NNLL+NNLO theoretical prediction and eq. (4.1) to model the nonperturbative contributions, we obtain our best fit value of α S (M Z ) = 0.121 +0.001 −0.003 which is compatible with the world average within uncertainties. It would be very interesting to perform a more comprehensive phenomenological analysis and a precise measurement of α S from EEC, using modern Monte Carlo tools to extract the hadronization corrections from data. This work is in progress.