Second-order QCD corrections to event shape distributions in deep inelastic scattering

We compute the next-to-next-to-leading order (NNLO) QCD corrections to event shape distributions and their mean values in deep inelastic lepton–nucleon scattering. The magnitude and shape of the corrections varies considerably between different variables. The corrections reduce the renormalization and factorization scale uncertainty of the predictions. Using a dispersive model to describe non-perturbative power corrections, we compare the NNLO QCD predictions with data from the H1 and ZEUS experiments. The newly derived corrections improve the theory description of the distributions and of their mean values.


Introduction
Event shape variables allow various kinematical properties of hadronic final states to be analysed. The resulting event shape distributions were measured extensively in e + e − [1] and ep [2] collisions, enabling a variety of precision QCD studies, including measurements of the strong coupling constant, resummation and parton-shower effects, investigations of non-perturbative power corrections, and tuning of multipurpose event simulation models.
Precision studies of event shapes distributions demand that their theoretical description is of comparable accuracy to the experimental measurements, requiring the calculation of higher order contributions in perturbative QCD. For e + e − event shapes, an appropriate level of theory precision was achieved already some time ago with the calculation of the next-to-next-to-leading order (NNLO) QCD corrections [3][4][5][6][7][8][9] in the form of generic parton-level event generators that allow any infrared-safe event shape distribution to be computed. These fixed-order NNLO results can be combined with resummation of large logarithmic corrections to next-to-nexta e-mail: jmo@physik.uzh.ch to-leading logarithmic level (NNLL) and beyond for specific event shape variables [10][11][12][13][14][15].
For event shapes in deeply inelastic ep scattering (DIS), the currently available level of theoretical accuracy is lower, with fixed-order results only known to next-to-leading order (NLO) [16][17][18] and resummation at next-to-leading logarithmic level (NLL) [19][20][21][22]. The theory uncertainty (as quantified through variation of the renormalization and factorization scales) on these predictions is often comparable to or larger than the experimental errors on the event shape measurements from H1 [23] and ZEUS [24], thereby limiting the extraction of fundamental QCD parameters from these data. To overcome this limitation requires an improvement of the fixed-order predictions to NNLO, which is presented in the following.
This paper is structured as follows. In Sect. 2, we summarise the definitions of the most common DIS event shape variables, and the kinematical ranges covered by the H1 [23] and ZEUS [24] measurements. The calculation of NNLO corrections to event shape distributions is performed in the NNLOJET framework [25] and follows closely the related NNLO calculations of jet production in DIS [26,27] and is documented in Sect. 3. To compare the resulting parton-level NNLO predictions with experimental hadron-level data, we employ a dispersive model [28][29][30], described in Sect. 4, determining the non-perturbative power corrections to the event shape distributions. We perform detailed comparisons of the hadron-level predictions to event shape data from H1 and ZEUS in Sect. 5. Our findings are summarized in Sect. 6.

Event shape variables
Event shapes in DIS are measured in the Breit frame, defined by the momentum directions of the virtual photon (current axis) and the proton (remnant axis), and boosted such that the energy component of the virtual photon momentum vanishes. The Breit frame provides a separation in pseudorapidity η between the proton remnant (remnant hemisphere, η > 0) and the hard scattering process (current hemisphere, η < 0). The event shape variables are dimensionless quantities that are determined from the four-momenta p h = (E h , p h ) of all particles in the current hemisphere. The different variables [22], which are generically denoted as F, are defined as follows.
The thrust τ γ measures the longitudinal momentum components projected onto the current axis: Thrust τ T is the thrust with respect to the thrust axis in the direction n T which maximizes the longitudinal momentum components projected onto this axis: This is analogous to the definition of thrust in e + e − collisions.
The jet mass parameter ρ is the squared invariant mass in the current hemisphere, normalized to four times the total energy squared: The jet broadening B γ measures the sum of the transverse momenta with respect to the current axis: As with thrust, the jet broadening can also be defined with respect to the thrust axis: Finally, the C-parameter is derived from the linear momentum tensor i j : with the eigenvalues λ 1 , λ 2 , λ 3 of i j yielding Equivalently, it can be expressed as where θ hh is the angle between particles h and h . In the experimental analysis, the event shapes are computed from the hadron momenta in the current hemisphere, while the theoretical calculation uses the parton momenta. For the Born-level contribution to inclusive DIS, leptonquark scattering, only the final state quark is produced in the current hemisphere, with thrust axis and current axis coinciding. Consequently, all event shape variables defined above become trivially zero. The first non-trivial contribution to the event shape distributions arises from two-parton final states: eq → eqg or eg → eqq, such that the leading-order (LO) perturbative contribution is O(α s ). The event shape distributions are thus closely related to DIS two-jet production in the Breit frame.
In higher-multiplicity final states, it is possible that all partons scatter into the remnant hemisphere, leaving the current hemisphere empty. To ensure infrared safety of the observables, these events are not accepted by demanding that the total energy in the current hemisphere of an event exceeds some minimum value lim Event shapes in deep inelastic scattering have been measured at HERA by the H1 [23] and ZEUS [24] experiments, based on the analysis of electron-proton scattering data taken at a centre-of-mass energy of √ s = 319 GeV (the H1 data set also contains a small fraction of data taken at √ s = 301 GeV). The DIS kinematics in the process e(k) + p( p) → e(k ) + X ( p X ), with momentum transfer q = k − k is described by the variables Q 2 = −q 2 , x = Q 2 /(2q · p) and y = Q 2 /(xs).
The H1 analysis [23] selects events with which are then classified into bins in Q = Q 2 , as listed in Table 1. For the event shape determination, lim = Q/10 is used. The ZEUS analysis [24] covers the kinematic range 0.0024 < x < 0.6 , 0.04 < y < 0.9, with events binned into in (Q 2 , x), described in Table 2.
The energy cut in the current hemisphere used by ZEUS is lim = Q/4. Both experiments normalize the event shape distributions to the DIS cross section integrated over the kinematical bin under consideration, which is determined without applying the lim cut.
Both experiments performed measurements [23,24] of the event shape distributions for F = τ γ , τ T , ρ, B γ , C. In addition, they also measured the mean values F for these variables, supplemented in the ZEUS study by a measurement of the mean value B T of the jet broadening with respect to the thrust axis. The measurements of the mean values are done for the same kinematical bins, Tables 1 and 2, as used for the distributions.

QCD corrections to event shapes
The event shape variables defined above assume non-trivial values only for final states containing two or more partons. Consequently, the event shape distributions in DIS receive the same parton-level contributions as two-jet production in DIS. Higher-order QCD corrections to event shape distributions can thus be obtained from the corresponding calculation for di-jet production by replacing the jet reconstruction algorithm by computations of the event shape variables.
We calculate the differential distributions and mean values for the DIS event shapes with the parton-level Monte Carlo event generator NNLOJET, by extending the existing calculation of NNLO corrections to di-jet production in DIS [26,27].
It combines the contributions from four-parton production at tree-level [31][32][33], three-parton production at one loop [34][35][36][37] and two-parton production at two loops [38][39][40][41], using the antenna subtraction method [42][43][44] to isolate infrared singular terms from the different contributions, which are then combined to yield numerically finite predictions for arbitrary infrared-safe observables constructed from the parton momenta. Besides for di-jet production at NNLO, the same ingredients and setup have been used previously in the computation of N 3 LO corrections to single jet production in DIS [45], in extractions of the strong coupling constant from DIS jet data [46,47], and in studies of diffractive di-jet production [48]. The calculations have also been extended to jet production in charged current DIS [49,50] at the same perturbative orders.
We compute the event shapes for electron-proton collisions with √ s = 319 GeV, using the NNPDF3.1 parton distributions with α s (M Z ) = 0.118 and for N F = 5 massless quark flavours. Central renormalization and factorization scales are fixed to μ F = μ R = Q, and theory uncertainties are estimated by the envelope of varying these scales independently by a factor two up and down, avoiding the pairings of variations in opposite directions (seven-point scale variation). Event selection cuts on the lepton variables and on h E h are applied according to the H1 [23] and ZEUS [24] analyses, and events are then classified into the different kinematical bins of Tables 1 and 2. The total hadronic DIS cross section for each kinematical bin (required for the normalization of the event shape distributions and mean values) is obtained to NNLO from NNLOJET, based on the one-jet calculation to this order [45]. Central renormalization and factorization scales are used for the normalization.

Event shape distributions
The event shape distributions are computed as histograms in the event shape variables. We use a considerably finer bin resolution than in the experimental analyses [23,24], which will subsequently allow us to apply hadronization corrections that result in a dynamical shift of the event shape variables. The histograms are defined in terms of variable ranges and number of equal-sized bins: The fixed-order calculation for an event shape F diverges in the limit F → 0, where all-orders resummation of large log(F)-terms is required. In this limit, the fixed-order expressions become meaningless, and we accordingly apply cuts on Fig. 1 Fixed-order predictions for the event shape distribution for H1 kinematics [23] in Q = 30-50 GeV bin: LO (green), NLO (blue) and NNLO (red), for H1 kinematics [23]. The lower frames display the ratio to the NLO predictions for the central scale μ 2 = Q 2 the minimum values of each shape variable, which set the first few bins of the distributions to zero: These cuts are typically within the first bin of the experimental analysis, which should anyhow be discarded in the comparison of fixed-order theory and experimental data. Figure 1 displays the fixed-order predictions (re-binned from the initial histograms by combining four adjacent bins each) for the H1 kinematics in the Q = 30 − 50 GeV bin. Since the qualitative behaviour of the higher-order corrections to the distributions is similar for all kinematical bins, we show the fixed-order distributions without power corrections only for one representative bin. The quantitative size of the corrections and of their uncertainties decreases with increasing Q, mainly due to the decrease in the running coupling constant α s .
In general, we observe that the NNLO corrections in the bulk of all distributions are typically positive (up to + 20%), often displaying only a marginal or no overlap of the uncertainty bands at NLO and NNLO. The scale uncertainty decreases from NLO (∼10%) to NNLO (∼5%). Even in the bulk, the higher-order corrections are not uniform between the distributions, each displaying a non-trivial shape in the NNLO/NLO ratio.
Towards the kinematical edges F → 0 and F → F max , the higher-order corrections behave differently for each distribution, often displaying large effects well beyond the scale uncertainty estimates. For F → F max , these features are caused by two different but related issues. For some of the shape variables, F max can not yet be realised in the Born process, owing to its low multiplicity. This is the case for the C-parameter which has a Born-level upper limit of 3/4 and for τ T with an upper limit of 0.293. Higher order real radiation corrections allow to attain larger values of F, thereby resulting in a kinematical mismatch between real and virtual contributions (Sudakov shoulder, [51]), which (although finite) produces large perturbative corrections in the vicinity of the Born-level kinematical limit.
In the case of DIS event shape variables, the kinematical constraints of the Born process produce further structures that narrow down the dimensionality of the final state phase space for specific values of different variables. These ridges in the multi-dimensional phase space were investigated in detail in [21] and produce kinks and spikes in the one-dimensional event shape distributions. These are sometimes already present at leading order, and go along with large and unstable higher order corrections in the immediate vicinity of the exceptional points, which are visible in particular in the distributions in C and τ T in Fig. 1. These features are typically localised in small patches of the phase space. For sufficiently large bin sizes, their impact is diluted to an invisible level. High-resolution measurements of event shape distributions, for example at a future electron-ion collider [52] or at the LHeC [53] will be able to resolve these features, thereby potentially necessitating resummation of large corrections associated with them.
At low values of F, the fixed-order predictions contain logarithmic terms log F at each order in perturbation theory, which spoil the convergence of the fixed-order perturbative expansion. In Fig. 1, the onset of these effects is visible in particular in the B γ distribution, while its onset takes place only at lower values of F in all other distributions. A description of the event shape distributions over the full kinematical range, and extending towards lower values of F than probed by currently available measurements [23,24] will need to include the resummation of these log F terms, which is currently known to next-to-leading logarithmic level [19][20][21][22] for all distributions.

Mean values
The mean values of the different event shapes variables are computed using NNLOJET by weighting each event with the reconstructed value of the event shape variable under consideration. The phase space integrations are performed by imposing only a very low technical cut-off of F min = 10 −4 on the event shape variables, since the weighting with the shape variable regulates the divergent behaviour of the integrals for F → 0, rendering the mean value integrals finite. The numerical stability of the mean value integrals has been checked with an even lower technical cut-off of F min = 10 −5 , from which we conclude that F min = 10 −4 is already sufficient for a stable result. The mean values are also normalized to the inclusive hadronic cross sections.
The fixed-order predictions for the mean values are displayed in Figs. 2 and 3, for the H1 [23]and ZEUS [24] kinematics. With the exception of the broadenings B γ and B T , the NNLO corrections to the mean values are positive for all event shapes, and decrease in magnitude with increasing Q 2 . The NNLO predictions are often at the upper boundary of the NLO theory uncertainty band, for the lowest Q 2 bins they are even outside the NLO band. For the broadenings, the NNLO corrections to B γ are positive at low Q 2 , and Fixed-order predictions for the mean value of the event shapes at LO (green), NLO (blue) and NNLO (red) compared to ZEUS data [24] become negative at large Q 2 , to B T they display the opposite behaviour, and are smaller in absolute magnitude. For all mean values, inclusion of the NNLO corrections leads to a reduction of the scale uncertainty compared to NLO, which is most pronounced for the broadenings, whereas being more modest for the other shape variables. For large values of Q 2 > 2500 GeV 2 , the NNLO theory uncertainty is limited to below 5%.
Comparing the fixed-order predictions to the measurements of the mean values from H1 and ZEUS, we observe that the data are considerably above the theory predictions throughout all shape variables and for all values of Q 2 , although the discrepancy is most pronounced at low Q 2 . This behaviour indicates the relevance of power corrections from hadronization effects, which can have large effects on the mean values [28][29][30].

Hadronization effects
In the previous section, we computed higher-order corrections to the DIS event shape distributions and mean values at parton level. To compare these predictions with hadron-level data requires accounting for the impact of the partonhadron transition, which is a non-perturbative process. Consequently, these hadronization effects cannot be computed in perturbation theory, but require a non-perturbative model description. The hadronization corrections are expected to be suppressed by positive powers of /Q, such that their relative numerical impact is decreasing with increasing Q. In the following, we employ the dispersive model [28] to estimate the leading power corrections at order ( /Q) to event shape distributions. This model has been worked out in detail for the DIS event shapes in Ref. [29], and its implications are briefly summarized in the following.
In the dispersive model, an effective coupling α eff is introduced at low scales, which is matched to the running QCD coupling α s (μ) at a scale μ I = 2 GeV. This gives a constant α 0 which is defined as the first moment of the effective coupling below the scale μ I , The power corrections are suppressed by powers of 1/Q, and result in a shift of the perturbative differential distribution where the power correction P is universal for all the event shape variables. The perturbative ingredients to the dispersive model are the running of the coupling constant and the relation between the MS-coupling and the effective coupling, whose definition [54] absorbs universal correction terms from the cusp anomalous dimension. It can be expanded in α s (Q), and its expression up to NNLO is given by [55,56]: with M = 1.49 a constant normalization factor (Milan factor [30,57]) accounting for higher-order contributions. In our numerical results, we use α 0 (μ I ) = 0.5 at μ I = 2 GeV, which has been estimated from fits to event shape moments in DIS [23,24,58] and e + e − annihilation [56]. The betafunction coefficient and cusp anomalous dimension [59] in the above expression are: with C A = 3, C F = 4/3, T F = 1/2. The coefficients a F depend on the event shape variable, they were computed in [29] and are tabulated in [58]. Their values are repeated here: where the shift of the jet broadening has an additional enhancement [60] given by with α s evaluated at the scale μ = e − 3 4 μ R and η 0 = −0.614. For this analysis a B varies from 1.6 to 2.3.
The dispersive model is based on an analytic treatment of hadronization effects on a two-parton correlator [28], which corresponds to the mean value integral of each event shape. The effect of the dispersive power correction on the mean values is additive: where F pert. is the mean value obtained in fixed-order perturbation theory, described in Sect. 3.2 above. When applied to differential event shape distributions, the power correction P in the shift (15) can in principle depend on the numerical value of F. Using a constant shift P for the full distribution amounts to an approximation, which may be overcome by an improved treatment of the non-perturbative corrections.
In combining the fixed-order predictions derived in the previous section with the power corrections, we truncate the factor P in (16) to α 2 s (Q) for NLO, and to α 3 s (Q) for NNLO. Inclusion of the α 3 s (Q) terms leads to a substantial reduction of P, with P NNLO ≈ 0.60 P NLO at Q = 15 GeV and P NNLO ≈ 0.75 P NLO at Q = 100 GeV. Fig. 4 Event shape distribution for thrust with respect to boson axis: τ γ fixed-order predictions at NLO (dashed cyan), NNLO (dashed brown), and corrected for hadronization effects at NLO (blue) and NNLO (red), compared to H1 data [23]. The lower frames display the ratio to the NLO prediction for the central scale μ 2 = Q 2

Results
With the inclusion of power corrections, the fixed-order parton-level predictions can now be compared to hadronlevel data from the HERA experiments on event shape distributions and mean values.

Event shape distributions
Figures 4, 5, 6, 7, 8, 9, 10, 11, 12 and 13 display the theory predictions obtained by combining the fixed-order predictions up to NNLO with power corrections computed using the dispersive model as described in the previous section to the experimental data from H1 [23] and ZEUS [24]. To illustrate the magnitude of the power corrections, the uncorrected fixed-order predictions for central scales μ = Q are indicated by blue lines at NLO and brown lines at NNLO. The shift (15) is applied on the high-resolution histograms (12) which were computed with a lower cut-off affecting their first bin (where all-order resummation of large logarithmic corrections [19][20][21][22] is required to obtain a finite prediction). The shifted high-resolution histograms are then Owing to the interplay of the lower cut-off on the distributions and the power correction shift, the prediction for the left-most non-vanishing bin of each distribution is unreliable, and should not be taken into account when comparing the experimental data with the theory predictions. A prediction for F → 0 will have to include resummation in order to become meaningful. This limitation should be kept in mind in the following comparisons.
For the thrust distribution τ γ , Figs. 4 and 5, we observe that the NNLO corrections at low and moderate values of Q 2 are leading to an increase of the distribution in the bulk, and a slight decrease at high and low τ γ . At the highest values of Q 2 , the NNLO corrections become very small and nega- Fig. 6 Event shape distribution for jet mass: ρ fixed-order predictions at NLO (dashed cyan), NNLO (dashed brown), and corrected for hadronization effects at NLO (blue) and NNLO (red), compared to H1 data [23]. The lower frames display the ratio to the NLO prediction for the central scale μ 2 = Q 2 tive even in the bulk. Overall, the NNLO corrections improve the description of the data. Compared to NLO, inclusion of the NNLO correction leads to a reduction of the scale uncertainty. This reduction is only moderate at the lowest values of Q 2 , where the NNLO scale uncertainty remains at the 6% level. At higher Q 2 , the reduction of scale uncertainty at NNLO is more pronounced, leading to predictions with residual uncertainty below 4%. These uncertainties should be compared to the experimental errors. The ZEUS data [24] are slightly more precise than the H1 data [23], and also reach to lower values of Q 2 . In the low-Q 2 bins, the NNLO scale uncertainty remains larger than the experimental errors, as was also observed [61] for jet production in DIS at low Q 2 . For moderate and high values of Q 2 , the scale uncertainty is now well below the experimental errors, thereby allowing for the use of the event shape distributions in precision QCD studies.
A similar pattern is also observed in the jet mass distribution, Figs. 6 and 7: positive NNLO corrections in the bulk at moderate Q 2 , which turn negative when going to large values of ρ or to large Q 2 . At the lowest values of Q 2 , the NNLO corrections are negative throughout the distribution (despite Fig. 7 Event shape distribution for jet mass: ρ fixed-order predictions at NLO (dashed cyan) and NNLO (dashed brown), and corrected for hadronization effects at NLO (blue) and NNLO (red), compared to ZEUS data [24]. The lower frames display the ratio to the NLO prediction for the central scale μ 2 = Q 2 positive corrections at parton-level) due to the reduced size of the power correction at NNLO. The NNLO corrections lead to an improved description of the shape of the experimental data. This improvement is particularly visible at large ρ for all values of Q 2 , where negative NNLO contributions lead to a considerably better description of the kinematical shape of the data. Overall, the agreement between data and theory is however somewhat worse for ρ than it was for τ γ . The NNLO scale uncertainties also follow a similar pattern as for τ γ : compared to NLO only a modest reduction at low Q 2 and a substantial reduction to the level of a few per cent at high Q 2 . Again, the experimental errors are larger than the scale uncertainty for moderate and high Q 2 , thus enabling precision QCD studies.
In the C-parameter, Figs. 8 and 9, we must distinguish the region below and above the Sudakov shoulder, which is located at C = 0.75 in the perturbative parton-level expression. The dispersive power corrections shift the location of this shoulder to higher values of C. This shift is largest at low Q 2 , and decreases in magnitude towards higher Q 2 . The region above the Sudakov shoulder is kinematically forbidden at LO, and receives contributions only from NLO Fig. 8 Event shape distribution for C-parameter: C fixed-order predictions at NLO (dashed cyan), NNLO (dashed brown), and corrected for hadronization effects at NLO (blue) and NNLO (red), compared to H1 data [23]. The lower frames display the ratio to the NLO prediction for the central scale μ 2 = Q 2 onwards. Already for values of C below the Sudakov shoulder, the pattern of NNLO corrections is more intricate than what was observed in τ γ and ρ. The observed structure is due to presence of a kinematical ridge [21] in the perturbative expressions at C ≈ 0.515, which destabilizes the perturbative convergence of the distribution in its vicinity, clearly visible in the high-resolution C-parameter distribution, Fig. 1. The perturbative predictions for the C-parameter distributions become quite precise at NNLO, with scale uncertainties of typically less than 8% below the Sudakov shoulder and away from the kinematical ridge. They are however affected by large hadronization corrections, which shift the whole distribution by more than two bins in C at low Q 2 . Compared to all other event shape distributions, these power corrections are particularly large in the C-parameter distributions, see (18). For the lower values of Q 2 , we also observe that the shape of the distribution is poorly described. For medium and large values of Q 2 , the power corrections are much smaller, NNLO corrections are relatively small and uniform, and a satisfactory description of the experimental data is observed.
The thrust distribution with respect to the thrust axis τ T , Figs. 10 and 11, displays a similar pattern, with non- Fig. 9 Event shape distribution for C-parameter: C fixed-order predictions at NLO (dashed cyan) and NNLO (dashed brown), and corrected for hadronization effects at NLO (blue) and NNLO (red), compared to ZEUS data [24]. The lower frames display the ratio to the NLO prediction for the central scale μ 2 = Q 2 trivial structures in its perturbative expressions around the Sudakov shoulder at τ T = 0.293 and the kinematical ridge around τ T ≈ 0.13, which are both nicely visible in the τ T -distribution at high resolution, Fig. 1. These exceptional points are shifted to larger values of τ T by the power corrections. The NNLO corrections are negative and small throughout almost the whole distribution for all values of Q 2 , except in the bin with the highest τ T , which is already well above the Sudakov shoulder, and where the cross section is already very small. As for τ γ , we note that the smallness of the NNLO effect is mainly due to a cancellation between positive parton-level corrections and a decrease in the power corrections. The corrections are typically within the NLO scale uncertainty band. The overall agreement between experimental data and theory predictions is satisfactory only at higher Q 2 . Substantial discrepancies at low Q 2 are observed especially in the vicinity of the kinematical ridge, where the theory predictions are systematically and considerably below the data. This behaviour may be indicative of the need of an re-consideration of the perturbative expansion and of the hadronization corrections in regions around the kinematical ridges. Fig. 10 Event shape distribution for thrust with respect to thrust axis: τ T fixed-order predictions at NLO (dashed cyan), NNLO (dashed brown), and corrected for hadronization effects at NLO (blue) and NNLO (red), compared to H1 data [23]. The lower frames display the ratio to the NLO prediction for the central scale μ 2 = Q 2 Finally, for the jet broadening with respect to the boson axis B γ , Figs. 12 and 13, the NNLO corrections assume a non-trivial shape, changing from negative at small B γ to positive at large B γ , thereby leading to a considerably better description of the data. In these distributions, the onset of large logarithmic terms at low B γ is well visible, indicating the need for their resummation. The NNLO corrections lead to a considerable reduction of the scale uncertainty in the bulk of the distributions, which is more pronounced at low Q 2 than in most other event shape distributions. The NNLO scale uncertainty ranges from 8% at low Q 2 to 3% at high Q 2 , which is comparable to or below the experimental uncertainties throughout.
Across the different event shape distributions, several common features are observed. The NNLO corrections are typically moderate, and fall usually within the NLO scale uncertainty bands. This is particularly remarkable since the NLO corrections were typically large (often comparable in size to the LO predictions), and well outside the LO scale uncertainty bands. The numerical smallness of the NNLO effect is often due to a partial cancellation between the parton- Fig. 11 Event shape distribution for thrust with respect to thrust axis: τ T fixed-order predictions at NLO (dashed cyan) and NNLO (dashed brown), and corrected for hadronization effects at NLO (blue) and NNLO (red), compared to ZEUS data [24]. The lower frames display the ratio to the NLO prediction for the central scale μ 2 = Q 2 level correction and modification of the power corrections at this order. Except in the low-F region, where large logarithmic corrections require an all-order resummation, and in the vicinity of Sudakov shoulders and kinematical ridges, we observe the onset of a good convergence of the perturbative fixed-order expansion. The corrections at low values of Q 2 are inevitably larger (due to the larger expansion parameter), which also translates in a sizeable scale uncertainty remaining at NNLO of about 10%. At larger values of Q 2 , this scale uncertainty improves considerably to the typically 4% or below, clearly highlighting the potential of precision QCD studies with event shapes based on existing HERA data [23,24], or with much larger data sets for hadronic final states that could be obtained at a future electron-ion collider [52] or at the LHeC [53]. The non-perturbative power corrections that we obtained in the dispersive model induce large shifts in some of the distributions, especially at low Q 2 , where the statistical quality of the data is largest. Moreover, their application to the distributions in the form of a constant shift is only an approximation, which should be revisited carefully as soon as more precise data are becoming available. Fig. 12 Event shape distribution for jet broadening with respect to thrust axis: B γ fixed-order predictions at NLO (dashed cyan), NNLO (dashed brown), and corrected for hadronization effects at NLO (blue) and NNLO (red), compared to H1 data [23]. The lower frames display the ratio to the NLO prediction for the central scale μ 2 = Q 2

Mean values
The power corrections to the mean values result in an additive shift of the perturbative predictions, see (20). As for the event shape distributions, we truncate P in (16) to order α 2 s (Q) for the power corrections to the NLO fixed order predictions and to α 3 s (Q) for power corrections applied to the NNLO predictions. Applying this shift to the perturbative results of Sect. 3.2, we obtain hadron-level predictions for the mean values, which are compared to the data from H1 [23] and ZEUS [24] in Figs. 14 and 15. The fixed-order predictions for central scales μ = Q are indicated by blue lines at NLO and brown lines at NNLO, showing that the power corrections are sizeable for all mean values. Their inclusion eliminates the tension between data and purely perturbative results seen in Figs. 2 and 3 above.
Comparing the mean values with and without power corrections, we observe that the large positive NNLO corrections at low Q 2 that are seen in Figs. 2 and 3 are more than compensated by the decrease in the numerical magnitude of the power correction in going from NLO to NNLO in P. The combined effect of the NNLO contributions to the fixed Fig. 13 Event shape distribution for jet broadening with respect to thrust axis: B γ fixed-order predictions at NLO (dashed cyan) and NNLO (dashed brown), and corrected for hadronization effects at NLO (blue) and NNLO (red), compared to ZEUS data [24]. The lower frames display the ratio to the NLO prediction for the central scale μ 2 = Q 2 order predictions and the power corrections is typically a small reduction of the mean values at low Q 2 , thereby leading to an improved description of the H1 and ZEUS data. At larger Q 2 , this combined effect results in a very small change of the predictions from NLO to NNLO, which comes with a substantial reduction of the perturbative scale uncertainty, which is almost halved.
Both the H1 [23] and ZEUS [24] studies used their measurements of event shape distributions for a simultaneous fit of the QCD coupling constant α s (M Z ) and the effective coupling α 0 that appears in the power correction, performed using NLO fixed-order results. While ZEUS lists only the results obtained for the individual event shape variables (displaying a substantial scatter), H1 also performed a combined fit, resulting in α s (M Z ) = 0.1198 ± 0.0012(exp) +0.0056 −0.0043 (th) and α 0 = 0.476 ± 0.008(exp) +0.018 −0.059 (th). Especially the theory error on α s (M Z ) is largely dominated by the NLO scale uncertainty.
A similar NNLO study has been performed previously on event shape moments in e + e − annihilation [56], resulting in decreased scatter between different shape variables, and in Fig. 14 Mean value of the event shapes at NLO (blue) and NNLO (red) including power corrections, compared to H1 data [23]. The fixed-order NLO and NNLO predictions (dashed cyan and brown lines) are included to illustrate the magnitude of the power corrections lower theory uncertainties. Also, a slight increase in the fitted central value of α 0 from NLO to NNLO was observed in e + e − , resulting in an NNLO value of α 0 = 0.5132 ± 0.0115(exp) ± 0.0381(th). Our predictions in Figs. 14 and 15 use α 0 = 0.5 throughout. They result in a very good description of the data at NNLO, and we notice that the NLO curves could be brought into better agreement with the data by a slight lowering of α 0 , towards its H1 fit value.
With the newly derived NNLO corrections, the combined fit of α s (M Z ) and α 0 to event shape distributions and their mean values can now be repeated fully consistently to NNLO accuracy. It can be anticipated that the main effect of the NNLO corrections will be in a reduction of the theoryinduced uncertainty on the extracted value of α s (M Z ), which was found to be about 5% in the NLO-based study by H1 [23]. This reanalysis of the experimental data will require our fixed-order results to be re-cast into convolution grids [47] that enable an efficient re-evaluation for multiple parameter values and parton distributions, which is beyond the scope of the present paper. Fig. 15 Mean value of the event shapes at NLO (blue) and NNLO (red) including power corrections, compared to ZEUS data [24]. The fixed-order NLO and NNLO predictions (dashed cyan and brown lines) are included to illustrate the magnitude of the power corrections

Conclusions
In this paper, we computed the NNLO QCD corrections to event shape distributions and their mean values in deep inelastic lepton-proton scattering. Our calculation was performed in the NNLOJET framework, and is largely based on the NNLO corrections [26,27] to di-jet production in DIS. The NNLO corrections to the distributions are not uniform, although some general trends are observed: positive corrections in the bulk of the distributions at low and medium Q 2 , negative corrections in the bulk at high Q 2 and at the upper kinematical boundaries of the shape variables for all Q 2 . Several perturbative instabilities due to Sudakov shoulders [51] or kinematical ridges [21] were observed in C and τ T . Predictions in the kinematic vicinity of these exceptional points will require novel resummation approaches to overcome the associated instabilities of the fixed-order predictions. Moreover, at low values of the event shape variables, we observed the onset of large logarithmic corrections at each order in perturbation theory. These were particularly pronounced in B γ . Resummation of these corrections is currently understood to next-to-leading logarithmic accuracy [22]. Aiming for a matching between fixed order and resummation in a form where the fixed-order expansion of the resummation formula reproduces all logarithmically enhanced terms up to NNLO (as was done for e + e − event shapes [10,12,15]) will require two more logarithmic orders in the resummation. To compare our parton-level predictions with hadron-level data, we used the dispersive model [28] to estimate the leading power correction effects from hadronization. The model is based on the study of two-point correlators which relate to the mean values of the event shape distributions. On the event shape distributions, additional assumptions must be made concerning the kinematical dependence of the power corrections. The power correction factors receive higher order contributions in the strong coupling constant, which we truncate to the same level as used in the fixed-order parton-level predictions.
Our resulting hadron-level predictions were compared to data from the H1 [23] and ZEUS [24] experiments. On the event shape distributions, we observe that inclusion of the NNLO corrections leads in general to an improved description of their kinematical shapes. Especially at medium and high Q 2 , the NNLO corrections result in a substantial reduction of the scale uncertainties of the predictions, to the level of a few per cent. A similar reduction of scale uncertainty is also observed on the mean values. On these mean values, we observe a compensation between the positive NNLO corrections to the fixed-order parton-level predictions and the negative NNLO contributions to the power corrections, resulting in a relatively small net effect at NNLO. Our newly derived NNLO results yield predictions with scale uncertainties that are typically below the experimental errors of the available HERA data on event shape distributions. They motivate a full NNLO-based reanalysis of event shape distributions and mean values. This should be leading to an improved determination of α s (M Z ) and α 0 , which was previously limited by the uncertainty on the NLO theory.
With high-resolution measurements of event shape distributions in deep inelastic scattering at a future electronion collider [52] or at the LHeC [53], our results will enable a broad spectrum of precision QCD studies, aiming for an improved understanding of its perturbative and nonperturbative aspects.
Data Availability Statement This manuscript has associated data in a data repository. [Authors' comment: The manuscript has associated data (annotated data tables for all predictions, with and without power corrections) which are stored as supplementary material by EPJC.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .