e+e− angularity distributions at NNLL′ accuracy

We present predictions for the e+e− event shape angularities at NNLL′ resummed and Oαs2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}_s^2\right) $$\end{document} matched accuracy and compare them to LEP data at center-of-mass energies Q = 91.2 GeV and Q = 197 GeV. We perform the resummation within the framework of Soft-Collinear Effective Theory, and make use of recent results for the two-loop angularity soft function. We determine the remaining NNLL′ and Oαs2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}_s^2\right) $$\end{document} ingredients from a fit to the EVENT2 generator, and implement a shape function with a renormalon-free gap parameter to model non-perturbative effects. Using values of the strong coupling αs(mZ) and the universal non-perturbative shift parameter Ω1 that are consistent with those obtained in previous fits to the thrust and C-parameter distributions, we find excellent agreement between our predictions and the LEP data for all angularities with a ∈ [−1, 0.5]. This provides a robust test of the predictions of QCD, factorization, and the universal scaling of the non-perturbative shift across different angularities. Promisingly, our results indicate that current degeneracies in the {αs(mZ), Ω1} parameter space could be alleviated upon fitting these parameters to experimental data for the angularity distributions.


Introduction
Event-shape variables are classic QCD observables that characterize the geometric properties of a hadronic final-state distribution in high-energy collisions [1]. They have been measured to high accuracy at LEP and other e + e − colliders in the past, and are commonly used for determinations of the strong coupling constant α s . Whereas past α s extractions from event shapes were often limited by perturbative uncertainties, various theoretical developments triggered a renewed interest in event-shape data in the past decade.
The theoretical description of event-shape distributions combines elements from QCD in different regimes. For generic values of the event shape, the distributions can be described in fixed-order perturbation theory where, currently, O(α 3 s ) corrections are known for a large class of variables [2][3][4][5]. In phase-space regions in which the QCD radiation is confined to jet-like configurations, the perturbative expansion develops large logarithmic corrections that need to be controlled to all orders. The resummation of these Sudakov logarithms was pioneered in [6] using the coherent branching formalism, and has since been reformulated and pushed to new levels of accuracy with methods from Soft-Collinear Effective Theory (SCET) [7][8][9][10]. To date, resummations for a variety of event-shape variables have been worked out to impressively high accuracies, e.g. to next-to-next-to-leading logarithmic (NNLL) [11][12][13] and even N 3 LL order [14][15][16].
The resummation of the logarithmic corrections tames the singular behaviour of the cross section near the kinematic endpoint and generates a characteristic peak structure whose shape is affected by non-perturbative physics. For the tails of the distributionsthe region commonly used for the α s determinations -the non-perturbative effects typically result in a shift of the perturbative distribution that is driven by a non-perturbative but universal parameter Ω 1 . Current α s extractions from event shapes (see e.g. [17][18][19][20][21]) are therefore organized as two-dimensional {α s (m Z ), Ω 1 } fits to the experimental data. The parameter Ω 1 is universal in the sense that the leading non-perturbative shift to a large number of event shapes depends only on this single parameter, scaled by exactly calculable observable-dependent coefficients, as can be demonstrated using factorization and an operator product expansion [22][23][24]. This universality does, however, require a careful consideration of hadron mass effects [25,26]. It is therefore possible to reduce the correlation between α s (m Z ) and Ω 1 in the two-dimensional fits by including data from different event-shape variables.
The most precise α s extractions from event shapes to date are based on the thrust [18,19] and C-parameter [21] variables. In these works the authors combine an N 3 LL resummation of large logarithmic corrections in the two-jet region together with fixed-order predictions up to O(α 3 s ) accuracy in the multi-jet region, and they use a shape function that reproduces the aforementioned shift in the tail region to model non-perturbative effects. While the values of α s (m Z ) the authors find from these fits are very precise and consistent with each other, they are noticeably below the world average [27] (see also [28]). As can be seen, e.g., in figure 4 of [21], the low value of α s (m Z ) seems to be associated with the implementation of non-perturbative effects. It is therefore important to test the underlying assumptions of the non-perturbative physics with other event-shape variables. . Angularity distributions at NNLL + O(α 2 s ) accuracy, convolved with a renormalon-free non-perturbative shape function, whose calculation is the subject of this paper. We display the predictions for three values of a (for now without uncertainties), illustrating roughly where two-jet and three-or-more-jet events lie in each τ a spectrum. For this illustration, the boundary is drawn at the value of τ a for a four-particle state that is grouped into pairs of jets with opening angle 30 • . As a becomes larger (smaller), the peak region is more (less) dominated by purely two-jet events.
In the present work we analyze a class of event shapes known as angularities, which are defined as [29] where Q is the center-of-mass energy of the collision and the sum runs over all final-state particles i with rapidity η i and transverse momentum p i ⊥ with respect to the thrust axis. The angularities depend on a continuous parameter a, and they include thrust (a = 0) and total jet broadening (a = 1) as special cases. Whereas infrared safety requires that a < 2, we restrict our attention to values of a ≤ 0.5 in this work, since soft recoil effects which complicate the resummation are known to become increasingly more important as a → 1 [30]. It is also possible to define τ a in eq. (1.1) with respect to an axis other than the thrust axis, such as the broadening axis or another soft-recoil-insensitive axis [31]. We stick to the standard thrust-axis-based definition here, to coincide with the available data. See [32] for a recent calculation with an alternative axis.
The phenomenological effect of varying a is to change the proportions of two-jet-like events and three-or-more-jet-like events that populate the peak region of the τ a distributions (see figure 1). The relevant collinear scale that enters the factorization of angularity distributions in the two-jet limit then varies accordingly with a, to properly reflect the transverse size of the jets that are dominating each region of the distributions.
The resummation of Sudakov logarithms for the angularity distributions is based on the factorization theorem [29,[33][34][35] 1 σ 0 dσ dτ a (τ a ) = H(Q 2 , µ) dt a n dt ā n dk s J a n (t a n , µ) J ā n (t ā n , µ) S a (k s , µ) δ τ a − t a n + t ā n Q 2−a − k s Q , (1.2) which arises in the two-jet limit τ a → 0. Here H is a hard function that contains the virtual corrections to e + e − → qq scattering at center-of-mass energy Q (normalised to the Born cross section σ 0 ); J a n,n are quark jet functions that describe the collinear emissions -3 -

JHEP01(2019)147
into the jet directions, and are functions of a variable t a n,n of mass dimension (2 − a); and S a is a soft function that encodes the low-energetic cross talk between the two jets and depends on a variable k s of mass dimension 1. In the language of SCET, the angularities are of SCET I type, since the virtuality of the soft modes µ 2 S ∼ Q 2 τ 2 a is smaller than that of the collinear modes µ 2 J ∼ Q 2 τ 2/(2−a) a for a < 1. In contrast, jet broadening with a = 1 corresponds to a SCET II observable, which requires somewhat different techniques for the resummation of the logarithmic corrections [36,37].
For very small values of τ a , the soft scale becomes non-perturbative and the soft effects are often parametrized by a non-perturbative shape function that is convolved with the perturbative distribution [38,39]. The dominant effect of the shape function in the tail region is a shift of the perturbative distribution: This form of the shift, scaling as 1/Q and with a value of c τa = 2/(1 − a) for the angularitydependent coefficient, was argued for in [40,41] based on the infrared behavior of the resummed perturbative exponent [39] and on dressed gluon exponentation [42], building upon earlier models of an effective infrared coupling leading to the proposed universal form of the shift in eq. (1.3) for numerous event shapes [43][44][45]. This universal form was later proven to all orders in soft emissions by an operator analysis in [23]. In the latter work, a definition of the non-perturbative parameter Ω 1 was obtained as a vacuum matrix element of soft Wilson lines and a transverse energy-flow operator. More importantly, it turns out that Ω 1 is the same non-perturbative parameter that enters the thrust analysis, related through the coefficient c τa . The growth of the shift with larger a reflects the greater dominance of narrow two-jet events in the peaks of the τ a distributions than for smaller a (see figure 1). It is one of the purposes of this work to provide a state-of-the-art description of the angularity distributions in order to verify if the angularity-dependent shift is reflected in the experimental data. To this end we will significantly improve the theoretical description of the angularity distributions for which, up until this year, the NLL +O(α s ) calculation from [34] represented the highest accuracy achieved. However, a recent calculation of the two-loop angularity soft function [46][47][48] allows the distributions to immediately be extended to NNLL accuracy, and indeed NNLL predictions using these two-loop results appeared recently in [32] in the context of multi-differential observables. NNLL predictions were also obtained using the ARES formalism in [49]. In this paper, we will further extend these results to NNLL + O(α 2 s ) accuracy 1 by utilizing the EVENT2 generator to determine the missing singular ingredients (the constant in the two-loop angularity jet function) and the full non-singular part of the fixed-order O(α 2 s ) prediction. Furthermore, in order to treat non-perturbative effects, we convolve the improved perturbative distribution with a shape JHEP01(2019)147 function, 1 σ 0 σ(τ a ) = dk σ PT τ a − k Q ; µ H , µ J , µ S , µ ns e −2δa(µ S ,R) d dk f mod k − 2∆ a (µ S , R) , (1.4) where σ PT represents the perturbative prediction. In our notation, σ and σ PT can represent either the differential or integrated spectra, (1.5) The perturbative cross section σ PT depends on the scales µ H,J,S associated with the factorization into hard, jet, and soft functions in the two-jet region (see eq. (1.2)), and a scale µ ns related to nonsingular terms in τ a which are predicted in fixed-order perturbation theory, but which are not contained in the factorized part of σ PT . That is, one can write the perturbative distribution as σ PT (τ a ; µ H , µ J , µ S , µ ns ) = σ sing (τ a ; µ H , µ J , µ S ) + σ ns (τ a , µ ns ) . (1.6) The non-perturbative shape function f mod contains a gap parameter ∆ a accounting for the minimum value of τ a for a hadronic (as opposed to partonic) spectrum, while R is the scale of subtraction terms δ a (µ S , R) that remove renormalon ambiguities between ∆ a and the perturbative soft function [51]. Formal definitions of all these objects will be given below. Finally, we compare our predictions to experimental data from the L3 Collaboration [52], which provides measurements of the angularity distributions for eight different values of a ∈ {−1.0, −0.75, −0.5, −0.25, 0.0, 0.25, 0.5, 0.75} at center-of-mass energies Q = 91.2 GeV and Q = 197 GeV (though we will stick to seven values in our analysis, leaving out a = 0.75 due to uncontrolled a → 1 corrections). The Q = 91.2 GeV data for two values of a is shown in figure 2 for illustration. The L3 analysis also includes a comparison of this data with different Monte-Carlo event generators and NLL theory predictions from [40,41,53], which include the non-perturbative scaling rule. However, with the higher-order perturbative contributions and more sophisticated treatment of nonperturbative effects now available, we think that the time is ripe for an updated comparison. In particular, our setup in eq. (1.4) allows for a clear separation of perturbative and nonperturbative effects, which is not possible with Monte Carlo hadronization models that were tuned to LEP data and which entered many of the theory comparisons in [52]. We can therefore rigorously assess the impact of the non-perturbative corrections in our framework. This paper is organized as follows: in section 2 we collect the formulae required for the resummation of Sudakov corrections in the two-jet region, which includes the new two-loop ingredients from the soft function calculation in [46][47][48]. In order to achieve NNLL accuracy, one in addition needs to obtain the corresponding two-loop jet function terms, which we determine from a fit to the EVENT2 generator in section 3. In this section we also perform the matching of the resummed distribution to the fixed-order O(α 2 s ) prediction. Then, in section 4, we discuss our implementation of non-perturbative effects and we present the final expressions of our analysis after renormalon subtraction. We further discuss our scale  -47 - Figure 2. Sample data at Q = 91.2 GeV from the L3 collaboration [52] for two angularities: a = −0.25 and a = 0.25. The theory predictions in the top panels are based on [40,41,53], which include the prediction of the non-perturbative scaling rule in eq. (1.3). choices in section 5, and compare our results to the L3 data in section 6. Finally, we conclude and give an outlook about a future α s determination from a fit to the angularity distributions in section 7. Some technical details of our analysis are discussed in the appendix.

NNLL resummation
The formalism for factoring and resumming dijet event shapes within a SCET I factorization framework is well developed and documented in many places (see, e.g., [33,35,54]) and will not be re-derived here. Below we will simply display the final results of these analyses and collect the required ingredients to achieve the NNLL resummation we desire. The precise prescriptions for which parts of eq. (1.4) are needed to which order in α s will be given in table 6 in section 4.3. In particular, to reach NNLL accuracy, we need to know the heretofore unknown two-loop jet and soft anomalous dimensions γ 1 J,S and finite terms of the two-loop jet and soft functions c 2 J,S (in a notation we shall define below). These have recently been determined or can be obtained from results in [46][47][48] and the EVENT2 simulations we report in this paper. The rest of this section details what these ingredients are and how they enter the final cross sections that we use to predict the angularity distributions.

Resummed cross section
The analytic forms for the resummed differential or integrated cross sections in τ a , derived in standard references like [34,35], are given by where the two cases are for σ being the differential or integrated distributions in eq. (1.5), and with the two functions F , G given by .
The Born cross-section contains a sum over massless quark flavours f = {u, d, s, c, b} with Q f being the charge of the associated flavour in units of the electronic charge e, and v f and a f are the vector and axial charges of the flavour: The jet and soft functions J, S appearing in eq. (2.1) are the Laplace transforms of J a n,n , S a from eq. (1.2), with their arguments written in terms of the logarithms on which they naturally depend (we suppress their indices to simplify the notation). The total evolution kernels K, Ω accounting for the running of the hard function H and the jet and soft functions J, S are given by constructed out of the individual evolution kernels which are determined from the anomalous dimensions of the functions F = H, J, S: -7 -

JHEP01(2019)147
The coefficients j F , κ F in eq. (2.6) are given by and RG invariance of the cross section eq. (2.1) imposes two consistency relations on these anomalous dimension coefficients, The appearance of partial derivatives ∂ Ω in the arguments of eq. (2.1) uses the notation of [55,56], and is due to the implementation of the following identity for arbitrary powers of logarithms: such that functions originally dependent on a logarithm L can be rewritten as F (L, µ) → F (∂ n , µ). The arguments are further shifted by logarithms of scale ratios in eq. (2.1) because we have pulled the evolution kernels through the fixed-order H, J, and S functions from the right. Note that in eq. (2.5) Ω is actually independent of the factorization scale µ due to eq. (2.9), whereas K still depends on µ, but this dependence cancels against the ω Fdependent factors on the first line of eq. (2.1). Of course, even the dependence on µ H,J,S cancels in the all-order cross section, but a residual dependence will remain at any finite order of resummed accuracy. The µ-dependence, on the other hand, should cancel exactly because of the consistency relations satisfied by the hard, jet, and soft anomalous dimensions. With standard perturbative expansions of K Γ , η Γ in eq. (2.7), however, this property does not precisely hold in practice. We explain why and show how to restore exact, explicit µ-independence in appendix C. Though the difference is formally subleading and numerically quite small, it is our aesthetic preference to work in this paper with the explicitly µ-invariant form. The result from appendix C for the cross section is (2.12) in terms of the modified cusp evolution kernel,

JHEP01(2019)147
and where K γ is the sum of just the non-cusp evolution kernels in eq. (2.7), (2.14) The reorganization of individual pieces of the evolution kernels that leads to eq. (2.11) restores explicit µ-invariance to every piece of this formula. In what follows we use eq. (2.11) to generate our theory predictions. (The numerical differences with eq. (2.1) will be negligible.)

Evolution kernels
The evolution kernels in eq. (2.7) or eq. (2.13) can be computed explicitly in terms of the coefficients of the perturbative expansions of the anomalous dimensions given by and where the running of α s is given by The values of the first few coefficients Γ n , γ n F , β n are given in appendix B. The formulas for the evolution kernels K Γ , η Γ can be evaluated at any finite N k LL accuracy in closed form by first making the change of integration variables from µ to α s , using eq. (2.16): and similarly for K γ F . Up to N 3 LL accuracy, the expansion of the kernel η Γ becomes, for instance,

JHEP01(2019)147
and K γ F is directly obtained from eq. (2.18) by using the replacement rules, The non-cusp kernels only need to be evaluated to one order lower than the cusp kernels, since they do not multiply an additional logarithm of µ F /Q F in the anomalous dimension or in eq. (2.1). Even though an appropriate choice of µ F will keep these logarithms small, cancellation of the µ-dependence of the K Γ kernels requires K Γ and η to be kept to the same accuracy, see [35] or appendix C. Meanwhile, the expansion of K Γ in eq. (2.7) is quoted in eq. (C.2). As mentioned above, we observe that, when using these standard expansion formulae for K Γ , the sum of evolution factors on the first line of eq. (2.1) does not exhibit exact independence of the factorization scale µ at any truncated order of resummed logarithmic accuracy, but instead has a residual (subleading) numerical dependence on it. This is due to a reliance on the relation used to obtain eq. (C.2), which is not exactly true if β[α] is truncated to finite accuracy. It is true at one-loop order and to infinite order, but not at any other fixed order. The corrections are, of course, subleading compared to the order to which eq. (2.21) is evaluated -see appendix C for details. The expansion for K Γ in eq. (2.11), on the other hand, does exhibit exact numerical independence on µ at any order. Its expansion up to N 3 LL accuracy is given by where r = α s (µ)/α s (µ F ) as in eq. (2.19), and We use these expressions to evaluate eq. (2.11) at a given resummed accuracy. The hard, jet, and soft functions H, J, S in eqs. (2.1) and (2.11) can be expanded to fixed order, as their corresponding logarithms are small near their natural scales µ H,J,S . We can obtain the generic form of their fixed-order expansions by making use of the solutions of the renormalization group equations (RGEs) that they satisfy, where F = H, J, S, and where Q H = Q, and Q F = Q/(e γ E ν a ) 1/j F for J, S. Here ν a is the Laplace-conjugate variable to τ a . From the solution of eq. (2.24), which is in general we can derive the form of each function F order-by-order in their perturbative expansions, 27) and to order α 2 s the coefficients F n are given by The corresponding logarithms are L H ≡ ln(µ H /Q) for F = H and L F = ln(µ j F F e γ E ν a /Q j F ) for F = J, S. In eqs. (2.1) and (2.11), each factor L F gets replaced by the differential operator shown in the argument of J, S for the jet and soft functions. The quantities Γ n F , γ n F are the coefficients of the perturbative expansions of the anomalous dimensions given by eq. (2.15), and the individual cusp parts Γ F of the anomalous dimensions are defined by where j F , κ F were given in eq. (2.8). In evaluating the ∂ Ω derivatives acting on the functions F , G in eqs. (2.1) and (2.11) up to two-loop order, we need to know the first four derivatives.

JHEP01(2019)147
For G in the integrated distributions, this yields where H(−Ω) ≡ γ E + ψ(1 − Ω) is the harmonic number function, ψ(x) is the digamma function, and ψ (n) (x) is its nth derivative. For F in the differential distributions, the derivatives can be obtained using F (Ω) = −Ω G(Ω), yielding The forms of eqs. (2.30) and (2.31) are convenient as they contain in closed, compact form the results of convolving fixed-order logarithmic plus-distributions in the momentumspace jet and soft functions with the generalized plus-distributions in the momentum-space evolution kernels (see, e.g., [18,34,57]).

NNLL ingredients
In order to achieve NNLL accuracy, one needs to combine the NNLL evolution kernels with the O(α 2 s ) fixed-order expansions of the hard, jet, and soft functions (the required ingredients are listed in table 6 in section 4.3). Whereas the cusp and the non-cusp anomalous dimension for the hard function are currently known to three-loop order, the non-cusp anomalous dimensions for the jet and soft function were previously only known to one-loop order for generic values of a [34]. 2 Recently, some of us developed a generic framework for the computation of two-loop soft anomalous dimensions [47], which we will use to determine γ 1 S . The two-loop jet anomalous dimension can then be determined from the consistency relation γ H + 2γ J + γ S = 0.
According to [47], the two-loop soft anomalous dimension can be written in the form where we have made the a-dependence of the anomalous dimension explicit, and are expressed in terms of the thrust anomalous dimension and a-dependent contributions expressed as integrals: , (2.34) which vanish for a = 0. The integral representations can easily be evaluated numerically to high accuracy for any value of a, and the relevant values for our work are given in table 1.
The constants in eq. (2.28) are known for the hard function at one [58,59], two [55,60], and three loops [18,61]. At NNLL accuracy, we will need these constants to two-loop order, and they are given by These results can be obtained from the two-loop quark form factor [62][63][64][65], also quoted in [55,60]. From the latter results, we replace the logarithm L = ln(Q 2 /µ 2 ) present there for deep-inelastic scattering with L = ln(−Q 2 − i )/µ 2 for e + e − annihilation, keeping the resulting extra π 2 and π 4 terms in the hard function. The result in eq. (2.35) agrees with the expressions in [14]. The expansion of the entire two-loop hard function given by eq. (2. 28) with the values in eq. (2.35) for the c 1,2 H also agrees with [55,60] with the replacement L → ln(−Q 2 − i )/µ 2 , and with the expressions in [14].
The soft function constants are known to one [54] and two loops [66,67] for a = 0, but they were, until recently, only known to one-loop order for generic values of a [34]. The two-loop a = 0 soft function constants have now been determined numerically in [46,48]. In Laplace space, they are where [34]  and c 2 S (a) is given by [46,48] Sample values of c CA 2 and c nf 2 for the values of a we plot in this paper are given in table 2. Finally, the jet function constants are known to one [68,69], two [70,71], and even three loops [72,73] for a = 0, but they were so far only known to one-loop order for generic values of a [34]. In Laplace space, they read where c 1 J is the momentum-space constant computed in [34], and The two loop constant c 2J (a) is thus far unknown and will be derived in the next section.

Numerical extraction of fixed-order contributions
In this section we determine the remaining fixed-order ingredients needed to predict the angularity distributions to NNLL +O(α 2 s ) accuracy: the constant c 2J in the two-loop jet function (as defined in eq. (2.28)) and the O(α 2 s ) part of the nonsingular distribution in eq. (1.6) predicted by full QCD. We will determine these in section 3.2 and 3.3, respectively, using the EVENT2 generator [74,75]. In section 3.1 we first derive the analytic relation between c 2J and the singular part of the cross section that will allow us to perform this extraction.

Fixed-order expansion
To extract the unknown constant c 2J in the jet function, we need to know the fixed-order expansion of the singular part of the cross section in eq. (2.1) to O(α 2 s ). We can do this -14 -

JHEP01(2019)147
by setting all the scales equal, µ = µ H = µ J = µ S , and multiplying out the fixed-order expansions of the individual H, J, S functions. In Laplace space, The expansions of each function are given in eq. (2.28), with Γ n F given by eqs. (2.29) and (2.8). The explicit expressions for the anomalous dimension coefficients and the constants c 1,2 F for F = H, J, S can be found in section 2.4 and appendix B. The two-loop jet function constant c 2J is the only missing ingredient at NNLL accuracy, and its numerical extraction is the main goal of this and the next subsection.
Writing out the individual expansions for H, J, and S explicitly and multiplying them together gives the Laplace-space cross section: where we have used the consistency relation γ H + 2γ J + γ S = 0 to eliminate γ H , along with γ 0 S (a) = 0. We have further expressed eq. (3.2) in terms of α s (Q), using eq. (B.3) to one-loop order, and we introduced the notationν a ≡ e γ E ν a .
We can now use eq. (A.5) to inverse Laplace transform eq. (3.2) back to momentum space. Specifically, we obtain the integrated distribution via where γ is such that the integration path lies to the right of all the singularities of the integrand. After performing the inverse Laplace transform and plugging in explicit values for the anomalous dimensions and fixed-order constants, we obtain c,sing (τ a ) = where we recall that the function f (a) was defined in eq. (2.41). Explicit values for γ 1 H n.A. and γ 1 S can then be inserted from eqs. (B.7) and (2.32). The total two-loop constant c (2) is given in terms of the Laplace-space constants by This formula immediately gives us c 2J as soon as we determine c (2) (which, we recall from eq. (3.6), is in momentum space), whose extraction from EVENT2 will be described in the next subsection.

Two-loop jet function constant
The program EVENT2 [74,75] gives numerical results for partonic QCD observables in e + e − collisions to O(α 2 s ). Using the method described by Hoang and Kluth [76], we can extract the singular constant c (2) in eq. (3.6), and thus the unknown jet function constant c 2J via eq. (3.7). For pedagogical purposes, we will give our own description of this method in the language of continuous distributions, which we find more intuitive to understand, rather than the language of discrete bins, which we encourage the reader to study in [76], as in practice one implements the discrete method.
The integrated (cumulative) angularity distribution in full QCD has a fixed-order expansion of the form: to O(α 2 s ). The c nm coefficients should agree with the SCET prediction in eq. (3.4) for the singular terms. The r n c functions are the nonsingular remainders that vanish as τ a → 0 and which are not predicted by the leading power expansion in SCET. Since SCET predicts the singular coefficients correctly, the difference of the QCD and SCET results is simply given by these remainders: Figure 3. Top: comparison of the EVENT2 prediction for two angularity distributions at O(α s ) with the analytically-known singular terms that follow from eq. (3.5). Bottom: comparison of their difference with the remainder function r 1 (τ a ) from [34].
which we will use in the next subsection to obtain the nonsingular remainder functions r n c from the difference of the EVENT2 output and the SCET prediction. To do this, however, we must know all the c nm coefficients in eq. (3.8), including the constants in c 20 ≡ c (2) in eq. (3.7). But we do not yet know c 2J .
In the limit of zero bin size, EVENT2 is generating an approximation to the differential distribution, which takes the form: where A is the constant coefficient of the delta function, B is a singular function, turned into an integrable plus-distribution, and r = dr c /dτ a is nonsingular, that is, directly integrable as τ a → 0. However, in practice we only obtain the distribution away from τ a = 0 from EVENT2, thus obtaining a histogram of the function: The delta function coefficient A, which gives the constants c n0 in the integrated distribution eq. (3.8), does not appear. The SCET prediction, on the other hand, gives just the singular parts: Figure 4. Integral eq. (3.14) of the O(α s ) remainder function from EVENT2 as a function of the lower integration limit τ a . As τ a → 0, a suitable fitting region is indicated that gives the constant r 1 c (1), for the case a = 0.
so the difference between the EVENT2 and SCET predictions (away from τ a = 0) gives the nonsingular part: Integrating this difference over 0 < τ a ≤ 1 then gives the expansion of the nonsingular part in terms of the r n c functions in eq. (3.8): Now, the total hadronic cross section, normalized to σ 0 , is simply the integrated distribution eq. (3.8) evaluated at τ a = 1. As the plus-distributions integrate to zero over this region, we use eq. (3.10) to obtain The total cross section is known to the considered order from [77][78][79]: From the perturbative expansion of eq. (3.15), we then obtain the relations Figure 5. Comparison of the EVENT2 prediction for two angularity distributions at O(α 2 s ) with the analytically-known singular terms that follow from eq. (3.6). The plots show the coefficients of the color structures C 2 The coefficient c 20 is precisely the constant c (2) in eq. (3.6), so we can use eq. (3.18) to determine c (2) -and thus c 2J through eq. (3.7) -from the known σ (2) tot in eq. (3.17) and the EVENT2 results for r 2 c (1). Essentially, the method relies on the fact that the total cross section is the sum of the singular constant A and the integral over the nonsingular distribution r c (1). Summing the EVENT2 bins, with the singular terms subtracted off, between a small τ a ∼ 0 and τ a = 1 gives the latter. Its difference from the known total cross section then gives the singular constants c n0 . This is illustrated for the O(α s ) parts of the differential distributions in figures 3 and 4. The EVENT2 data for these plots is based on 1 · 10 11 events with a cutoff parameter of 10 −12 .
In the upper plots of figure 3 we show the output of EVENT2 for two values of the angularity parameter, a = −0.5 and a = 0.5, compared to the known predictions for the singular terms. The difference is shown in the lower plots of figure 3 together with the remainder function r 1 (τ a ) from [34]. One sees that the difference between the EVENT2 output and the singular terms vanishes as τ a → 0 as expected. In figure 4 we show the result of computing the integral eq. (3.14) as a sum over EVENT2 bins, for a = 0. The value of this integral as the lower limit τ a → 0 gives the numerical result for the constant r 1 c (1), which through eq. (3.18) gives an extraction of the singular constant c 10 in the fixed-order distribution. The O(α s ) terms are, of course, already known; we show this just to illustrate the logic of our procedure.
Our EVENT2 results for the O(α 2 s ) parts of the differential distributions are displayed for a = −0.5 and a = 0.5 in figure 5. The functions plotted are the coefficients of the color structures C 2 F , C F C A , C F T F n f , and the EVENT2 output is again shown in comparison to the singular terms predicted by eq. (3.4). The agreement looks very good down to small values of τ a , until one gets so low that cutoff effects in EVENT2 play a role. For these plots we have again set the cutoff parameter to 10 −12 , and the results are based this time Figure 6. Integral eq. (3.14) of the O(α 2 s ) remainder function from EVENT2 as a function of the lower integration limit τ a . Shown are the coefficients of the C 2 F , C F C A , C F T F n f colour structures. Our best fit values for the constants r CF 2 , r CA 2 , r nf 2 defined in eq. (3.19) are shown in red.

JHEP01(2019)147
on 4 · 10 10 events. 3 The difference between the EVENT2 output and the singular terms is the nonsingular remainder function r 2 (τ a ), which we will show in section 3.3. Integrating this difference according to eq. (3.14) gives the numerical result for the integral r 2 c (1). The relevant plots for computing these integrals as a sum over EVENT2 bins are shown in figure 6, and the values we extract for r 2 c (1) are listed in table 3. The tabulated results are the coefficients of each color structure: The errors we report on the extracted central values are determined by the values of the fitted parameters at which the χ 2 per degree of freedom increases by 100% relative to that at the best fit value over the fit range of log 10 τ a ∈ (−5 + 4a/10, −4 + 4a/10) for C 2 F , log 10 τ a ∈ (−5, −4) for C F C A , and log 10 τ a ∈ (−4.5, −3.5) for C F T F n f , which we determined empirically by looking for stable plateau regions in figure 6 from which to perform the fits. These plateaus are very stable for the C F C A and C F T F n f color structures for most values of a, but less so for many of the C 2 F results -stability issues with C 2 F singular terms in EVENT2 have been encountered in prior analyses as well (e.g. [15,80]). This is possibly due to an undersampling of C 2 F contributions in the infrared region. Also from figure 6 it can be seen that we did not achieve stable plateaus for a = 0.5 (except arguably for C F T F n f ) before hitting the cutoff-affected region. However we went ahead and applied the same fit procedure as described above, and the result in table 4 for a = 0.5 correspondingly has a much larger uncertainty. Despite these issues with C 2 F and a = 0.5 contributions, for our present analysis we were nevertheless able to obtain sufficiently reliable results for the final cross sections we present in section 6. Namely, even doubling the uncertainties on c 2J shown in table 4, the total NNLL uncertainty bands on the cross sections in section 6 are unchanged for all a ≤ 0.25, and only change a few percent for a = 0.5 in the resummation region, with the methods for uncertainty estimation described in section 5.2. Of course, when other theoretical uncertainties are pushed sufficiently low to make the uncertainties on r CF 2 and for a = 0.5 more prominent (and before making a robust uncertainty estimate on any definitive extraction of α s ), one should revisit the methods used to compute these numbers.
The constant part of the singular O(α 2 s ) cross section is finally given by eq. (3.18), which we can plug into eq. (3.7) to obtain the so-far unknown two-loop jet function constant c 2J . Our results are shown in table 4, where we have set n f = 5 and we have added the uncertainties of the individual r 2 c (1) coefficients in table 3 linearly. In our phenomenological analysis below, we will vary c 2J over the uncertainties shown in table 4, and we will account for its contribution to the overall uncertainties presented in section 6. This particular contribution to the total uncertainty turns out, however, to be almost negligible. Our fit result for a = 0 agrees well with the analytically known result from [70], c 2J = −36.3. 3 We also have 1 · 10 11 events with a cutoff parameter set to 10 −15 , which are included in our results for the nonsingular remainder function, but which we do not use to extract the singular constants in the limit τa → 0 here. This is because the extracted values for the C 2 F constant departs substantially, and the CF CA constant marginally, from the known values at a = 0 when such a small cutoff parameter is used. We thus stick to the 10 −12 data, where the a = 0 constants come out correct, for extracting the singular constants.   Table 3. Fit values for the coefficients of the integral r 2 c (1) of the nonsingular QCD distribution as defined in eq. (3.19). The central values and their uncertainties have been extracted from the plots in figure 6 as described in the text. Table 4. Extracted values of the two-loop jet function constants c 2 J , determined by eq. (3.7) and the value for c (2) implied by eq. (3.20) and the numerical results for r 2 c (1) in table 3.

Remainder functions
We next move on to the integrable remainder function r 2 (τ a ) in eq. (3.13), and consequently its nonsingular integrated version r 2 c (τ a ) in eq. (3.8). A prior issue to sort out in this context is the kinematic endpoint of the τ a distribution in full QCD, which the singular predictions from SCET do not "see" since they assume τ a τ max a ∼ 1. The full QCD distribution dσ/dτ a in eq. (3.10) vanishes at O(α k s ) above the maximum kinematic value of τ a for (k + 2) particles in the final state, which we will call τ k,max a . The SCET distribution dσ sing /dτ a , however, continues above τ a > τ k,max a at this order. The remainder distribution r(τ a ) in eq. (3.13) is therefore just the negative of the SCET distribution above τ k,max a : where B k is the O(α k s ) coefficient of the singular distribution B(τ a ) in eq. (3.10) (in units of (α s (Q)/(2π)) k ), and we have slightly abused notation in using the same symbol r k for the remainder function for τ a < τ k,max a .
At O(α s ) there are up to three particles in the final state. The maximum value of τ a occurs for the symmetric "Mercedes-star" configuration with three particles in a plane, with equal energy E 1,2,3 = Q/3 and an angle θ ij = 2π/3 between any two of them (for the values of a we consider, but see [34] for exceptions). The thrust axis may then be taken to be along any of the three particles, and one easily derives τ 1,max s ) there are up to four particles in the final state, and it is similarly straightforward to compute the angularity τ a for the symmetric four-particle configuration, which we derive in appendix D. Assuming that the symmetric configuration again determines the maximum for the values of a relevant for our analysis, we obtain Figure 7. Left: O(α s ) remainder function for a = 0.25, which is known from [34]. Right: the corresponding O(α 2 s ) remainder function from EVENT2 (data points) with a suitable interpolation as described in the text.
We next turn to the functional form of the remainder functions r 1 (τ a ) and r 2 (τ a ), which are displayed in figure 7 for a = 0.25. For these results we collect EVENT2 data from all events at cutoffs of 10 −12 and 10 −15 as mentioned above. In the left panel the predicted curve from [34] for the O(α s ) remainder function is shown, whereas our EVENT2 results at O(α 2 s ) are given by the data points in the right panel. In this panel we also illustrate how, following [18], we divide the τ a domain into two regions τ a < 0.1 and τ a > 0.1, where we bin linearly in log 10 τ a or in τ a itself, respectively, to accurately capture the shape of the remainder function across the whole domain. In the linearly binned region, the EVENT2 uncertainties are negligible and we take a direct interpolation of the data points to be our working remainder function. In the low τ a region, on the other hand, we perform a fit to a predetermined basis of functions and use the result as our functional prediction. We then join these predictions smoothly to determine the final remainder function entering our cross section predictions: with f (z, z 0 , ) = 1/(1 + e −(z−z 0 )/ ). We actually transition to the direct interpolation of the EVENT2 data r 2 lin at τ a = 0.05, which is a bit below the boundary τ a = 0.1 between the logarithmically and linearly binned regions, since the EVENT2 data is so precise well below this point anyway.
The basis of functions we use to fit the data in the low τ a region is given by which is meant to be a representative but not necessarily complete (or even fully accurate) set of functions that can appear in the fixed-order O(α 2 s ) distribution. This basis is simply chosen since it gives a good fit to the growth of r 2 (τ a ) at small values of τ a to the accuracy of the EVENT2 data that we have. The fit is performed subject to the constraint that the total area r 2 c (1) under the remainder function r 2 (τ a ) is consistent with the values extracted in  Table 5. Results for the fit coefficients in eq. (3.24) obtained from the EVENT2 data displayed in figure 8. The uncertainties on these fits are illustrated in the plots on the right panel of figure 8.
shown gives a negligible shift in the extracted fit parameters in eq. (3.24). The errors on our remainder function r 2 (τ a ) actually come from the 1-σ error bands determined by the NonlinearModelFit routine in Mathematica, with the weight of each data point inversely proportional to the square of its EVENT2 error bar, and with a 0 determined by the area condition, leaving five free parameters in the fit. We give the results of these fits in table 5.
Our final results for the O(α s ) and O(α 2 s ) remainder functions for all seven values of the angularity a we consider in this work are shown in figure 8. The uncertainties of the aforementioned fit procedure are hardly visible in the plots in the middle panel of figure 8, but they become noticeable in the low τ a region, which is magnified in the right panel. Although the statistical errors from EVENT2 at low τ a grow substantially, the fits are constrained by the highly precise larger τ a data and the constraint on the total area under the remainder function. More precise results for r 2 (τ a ) at low τ a , and a more rigorous derivation of the fit functions that should appear in eq. (3.24), are certainly desirable for precision phenomenology in this region, but such exercises lie outside the scope of this paper. Here we seek only to obtain sufficiently smooth and stable visual results for very low τ a values, and focus on precision comparisons to data for intermediate-to-large τ a .
The remainder functions r 1 (τ a ) and r 2 (τ a ) have been extracted at a scale µ = Q, cf. eq. (3.8). The singular parts predicted by the factorization theorem eq. (2.1) and these non-singular parts must separately add up to be scale independent, if summed to all orders in α s (µ). Thus we can probe the perturbative uncertainty of the non-singular parts by varying them as a function of a separate scale µ ns : where and similarly for the integrated version r c . These functions r(τ a , µ ns ), r c (τ a , µ ns ) enter our final NLL +O(α s ), NNLL+O(α s ), and NNLL +O(α 2 s ) resummed and matched predictions for the cross sections we will present in section 6. We have now assembled all the pertur- bative ingredients that we need to make these predictions. But before we do so, we turn our attention to the implementation of non-perturbative effects, whose proper treatment can even improve the behaviour of the perturbative convergence.

Non-perturbative corrections
Having resummed the angularity distributions in section 2 and matched them to the fixedorder predictions in section 3, we will now include a treatment of non-perturbative corrections that will influence the overall shape and position of the distributions.

Non-perturbative shape function
As with any hadronic observable, event shapes are sensitive to low energetic QCD radiation and effects of confinement. The importance of these non-perturbative effects depends on the domain of τ a considered. For angularities with a < 1, power corrections from the collinear sector are suppressed with respect to those from the soft sector [23,24]. The non-perturbative effects can then be parameterized into a soft shape function f mod (k) that is convolved with the perturbative distribution [38,39,51]: which ultimately leads to eq. (1.4) for the cross section. Here S PT is the soft function computed in perturbation theory, and ∆ a is a gap parameter, which we will address in the next subsection. The shape function f mod (k) is positive definite and normalized. We follow previous approaches and expand the shape function in a complete set of orthonormal basis functions [57]: and P n are Legendre polynomials. The normalization of the shape function implies that the coefficients b n satisfy ∞ n=0 b 2 n = 1. In practice, we only keep one term in the sum (4.2), setting b n = 0 for n > 0 (cf. [16,18,81]). The parameter λ is then constrained by the first moment of the shape function as explained in the next subsection. More terms can in principle be included in eq. (4.2) if one wishes to study higher non-perturbative moments beyond the first one.
This function, when convolved with the perturbative distribution from the previous sections, reproduces the known shift in the tail region [43,44,82], which can be shown rigorously via an operator product expansion (OPE) [23] to be the dominant non-perturbative Here Ω 1 is a universal non-perturbative parameter that is defined as a vacuum matrix element of soft Wilson lines and a transverse energy-flow operator (for details, see [23]).

Gap parameter and renormalon subtraction
As advocated in [51], we use a soft function with a non-perturbative gap parameter ∆ a , as already displayed in eq. (4.1). The gap parameter accounts for the minimum value of τ a for a hadronic spectrum (the distribution can go down to zero for massless partons, but not for massive hadrons). In the tail region of the distributions, eq. (4.1) then leads to a shift of the perturbative cross section, eq. (4.4), with Since the first moment is shifted linearly by ∆ a , this parameter was rescaled in [34] from its default definition in [51] via where ∆ ∼ Λ QCD is an a-independent parameter. Note that this determines that, with eq. (4.2) truncated at n = 0, the model function parameter λ = 2(Ω 1 − ∆)/(1 − a). Up to this point, the barred quantities Ω 1 , ∆ (a) are taken to be defined in a perturbative scheme like MS in which S PT has been calculated. In [51] it was pointed out that such a definition of the gap parameter ∆ a has a renormalon ambiguity, shared by the perturbative soft function S PT in eq. (4.1). This is similar but not identical to the renormalon in the pole mass for heavy quarks (see, e.g. [83]). To obtain stable predictions, it is necessary to cancel the ambiguity from both S PT and ∆ a in eq. (4.1). This can be done by redefining the gap parameter as ∆ a = ∆ a (µ) + δ a (µ) , where δ a has a perturbative expansion with the same renormalon ambiguity as S PT (but opposite sign). The remainder ∆ a is then renormalon free, but its definition depends on the scheme and the scale of the subtraction term δ a . We adopt here the prescription chosen in [76] and also later implemented by [16,18,34], which is based on the position-space subtraction for the heavy-quark pole-mass renormalon introduced in [84,85]. We will 4 In the peak region, the OPE does not apply and the full shape function f mod (k) is required to capture the non-perturbative effects. Furthermore, the result in (4.4) is not only leading order in the OPE, it is also subject to other corrections like finite hadron masses and perturbative renormalization effects on the quantity Ω1, as described in [26]. 5 The expression for cτ a diverges in the broadening limit a → 1, where the SCETI factorization theorem we use breaks down. A careful analysis revealed that the non-perturbative effects to the broadening distributions are enhanced by a rapidity logarithm, cB T = ln Q/BT [24].

JHEP01(2019)147
translate this notation to the Laplace-space soft function we have been using in this paper -see eq. (2.28). The two formulations are completely equivalent with ν ↔ ix. The subtraction and its evolution equations are easier to formulate in Laplace space than in momentum space. For the Laplace-space soft function 6 S(ν, µ) = ∞ 0 dk e −νk S(k, µ) , (4.8) the convolution in eq. (4.1) becomes where we have grouped the renormalon-free gap parameter ∆ a with the shape function f mod , and the perturbative subtraction term δ a with S PT , rendering each group of terms in brackets renormalon free. There are a number of schemes one can choose to define δ a that achieve cancellation of the leading renormalon. A particularly convenient one found in [76,84,85] is a condition that fixes the derivative of the soft function at some value of ν to have an unambiguous value: where S PT (ν, µ) = e −2νδa(µ) S PT (ν, µ), which is sufficient to render S PT , and thus ∆ a , to be renormalon-free. This is known as the "Rgap" scheme, and its condition determines δ a as a function of a new, arbitrary subtraction scale R, which should be taken to be perturbative, but small enough to describe the characteristic fluctuations in the soft function. Explicitly, eq. (4.10) defines the subtraction term as and we see that δ a (and thus ∆ a ) depends on two perturbative scales, µ and R. Expanding the subtraction terms as we obtain, for the MS Laplace-space soft function S PT , using its expansion in eq. (2.28), where from eq. (2.29), Γ n S = −2Γ n /(1 − a), and γ 1 S and c 1 S are given by eq. (2.32) and eq. (2.36).

JHEP01(2019)147
Eq. (4.13) exhibits logarithms of µ/R that appear in the subtraction term δ a and thus the renormalon-free gap parameter ∆ a . Since µ = µ S will be chosen in the next section to be a function of τ a , which varies over a large range between µ 0 ∼ 1 GeV and Q, a fixed value of R can only minimize these logarithms in one region of τ a , but not everywhere. We therefore need to allow R to vary as well to track µ S , and so we need to know the evolution of ∆ a in both µ and R. The µ-RGE is simple to derive. Since ∆ a in eq. (4.7) is µ-independent, we obtain where from the perturbative expansion of δ a in eq. (4.13), we can determine explicitly to O(α 2 s ), and which can be shown to hold to all orders [76]. The solution of this RGE is given by with an initial condition at some scale µ ∆ , and where κ S = 4/(1 − a) was given in eq. (2.8) and the kernel η Γ was defined in eq. (2.7). The evolution of the gap parameter ∆ a (µ, R) in R is a bit more involved, and was solved in [85] for quark masses and applied to the soft gap parameter in [76]. We follow this derivation here (in our own notation). Since from eq. (4.16) we know how to evolve ∆ a (µ, R) in µ, we just need to derive the evolution of ∆ a (R, R) in R. Since ∆ a in eq. (4.7) is also R-independent, we can derive from the perturbative expansion of δ a in eq. (4.13) the "R-evolution" equation: where γ R has a perturbative expansion, whose first two coefficients we read off from eqs. (4.12) and (4.13), Even though γ 0 R = 0 for the soft gap parameter (since γ 0 S (a) = 0), we will keep it symbolically in the solution below for generality (and for direct comparison with the quark mass case in [85]).

JHEP01(2019)147
multiplying and dividing by R in the integrand, anticipating using eq. (2.21) to change integration variables to α s . But first we need to invert α s (R) to express R. To this end, we write eq. (2.21) in the form where G[α] is the antiderivative of 1/β[α], This determines G up to a constant of integration (we effectively choose it such that G[α] → 0 as α → ∞). If R, R ∆ are scales for which α s is perturbative, we can determine G explicitly order by order, and we can use this relation in eq. (4.20) to achieve the change of variables from R to α s , Keeping the explicit expression for G in eq. (4.23) up to the B 2 term, and expanding out the integrand in powers of α s where we can, we obtain This series of integrals is conveniently evaluated by changing the integration variable to t = −2π/(β 0 α), upon which These integrals can then be expressed in terms of the incomplete gamma function Γ(c,

JHEP01(2019)147
which is consistent with the solution of the R-evolution equation given in [76,85]. The solution of both µ-and R-evolution equations in eqs. (4.16) and (4.28) finally determines the renormalon-subtracted gap parameter ∆ a (R, µ) at any perturbative scales R, µ in terms of an input value ∆ a (R ∆ , R ∆ ) at a reference scale R ∆ : This expression is real, and we recall that γ 0 R = 0 in our case. In the solution of the µ-and R-evolution of ∆ a in eq. (4.29), we truncate at each order of logarithmic accuracy as follows: at NLL( ) we keep γ µ ∆ , γ R to O(α s ) (i.e. up to the second line of eq. (4.29)), and at NNLL( ) we keep γ µ ∆ , γ R to O(α 2 s ) (i.e up to the last line). These rules are summarized with all other truncation rules in table 6 below. Note that this means that η Γ is actually kept to one order of accuracy lower than indicated by eq. (2.18). This is because the µ-evolution of ∆ a in eq. (4.16) is not multiplied by an extra logarithm as for the hard, jet, and soft functions in the full factorized cross section. 7 Transforming the renormalon-free soft function in eq. (4.9) back to momentum space, we obtain the shifted version of eq. (4.1), (4.30) Then the parameter Ω 1 in eq. (4.5), describing the non-perturbative shift of the perturbative cross section induced by the shape function, turns into a renormalon-free shift: We will take the input gap parameter at a reference scale R ∆ = 1.5 GeV to be ∆(R ∆ , R ∆ ) = 0.1 GeV in our phenomenological analysis below. The exact value of this parameter is not particularly relevant to the tail region in which we focus our comparisons to data [18]. The shift in the perturbative part of eq. (4.30) can also be expressed in terms of a differential translation operator that acts on the perturbative soft function S PT : In comparing to the R-evolution for quark masses in [85], it may also appear that we keep one fewer order at N k LL accuracy than in that paper. But the counting for the logarithms is different in the two cases, since the logarithms appear as single logarithms for quark masses, but as double logarithms for event shapes; so the terms we call N k LL correspond to terms that are called N k−1 LL in [85]. Also, our truncation scheme seems to differ from the one applied for the gap parameter in [16,18], as described in eq. (A31) of [18] or eq. (56) of [16]. However, it is consistent with the corresponding tables in these papers and with the actual numerical implementations used by these authors in their results [86].

JHEP01(2019)147
which, after integrating by parts, gives In the final cross section, the renormalon-subtracted shape function then enters as a convolution against the perturbative distribution, which is the expression we anticipated in eq. (1.4). This implies we convolve both the singular and nonsingular parts of the cross section eq. (1.6) with the same, renormalonsubtracted shape function. In doing this we follow [18] and ensure a smooth transition from the resummation to fixed-order regime even after non-perturbative effects are included.
In practice we expand out the shape function terms to the order we work in α s , with δ 1,2 a (µ S , R) from eq. (4.13). The order to which these terms are kept at each accuracy are included in table 6.

Final resummed, matched, and renormalon-subtracted cross section
We now collect all pieces described above, giving our final expressions for the resummed cross section, matched to fixed-order and convolved with a renormalon-free shape function.
In evaluating the convolution in eq. (4.34), we must truncate the product of the fixed-order perturbative pieces contained in eqs. (2.11) and (3.13) along with the nonperturbative pieces in eq. (4.35) to the appropriate order in α s for N k LL ( ) accuracy. Namely, starting with eq. (2.11) for the integrated distribution, we expand the fixed-order coefficients in powers of α s : Table 6. Ingredients we include at various orders of unprimed N k LL (Left), primed N k LL (Middle), and matched (Right) accuracies, up to a given fixed order O(α n s ). The tables apply to the integrated distribution in eq. (4.38) and the Laplace-transformed distribution, but not, for unprimed accuracies, directly to the differential form in eq. (2.1) -see [35] for details. We have included a counting for the renormalon subtractions terms δ a in eq. (4.12) and the µ-and R-evolution anomalous dimensions γ µ ∆ , γ R in eqs. (4.15) and (4.18) as described in the text.
with the fixed-order operators F k at each order given by where α F ≡ α s (µ F ), H n , J n , S n are given by eq. (2.28), and r In the convolved cross section eq. (4.34), the proper fixed-order expansion up to NNLL( ) accuracy is then given by: where Eq. (4.40) is to be truncated at the fixed order demanded by table 6 (first term up to NLL, second term for NLL and NNLL, and up to the last term for NNLL ), while the evolution kernels contained in each expression are to be evaluated to the appropriate resummed logarithmic accuracy as described in section 2.
Eqs. (4.40) and (4.41) represent our final expression for the renormalon-free resummed and matched cross section that is convolved with a non-perturbative shape function. We should point out that we will perform the convolution in k prior to choosing particular values for the scales µ H,J,S,ns and R. We have clearly exhibited the dependence on these scales versus the explicit dependence on τ a appearing in eqs. (4.38) and (4.39). In the next section we will describe how we choose these scales, but for now it suffices to say that they are functions of the measured τ a in the cross section, and not functions of the convolution variable τ a − k/Q. Thus the τ a dependence inside the scales µ H,J,S,ns and R are not convolved over in eq. (4.41). In our numerical code, the convolution between the explicit τ a − k/Q dependence from eqs. (4.38) and (4.39) and the k dependence in f mod is then computed for a given set of scales µ i (τ a ), R(τ a ).
-34 -JHEP01(2019)147 Figure 9 illustrates the practical effect of implementing the renormalon subtraction and the convolution with the non-perturbative shape function as described above. We present the central theory curves at NNLL accuracy for the cumulant cross sections at seven values of the parameter a, all in the low-τ a domain. The curves in the left panel reflect the purely perturbative calculation, which exhibit unphysical (negative) values for the cross section. The curves in the right panel, on the other hand, represent the calculation performed with eq. (4.41). It is clear that the renormalon cancellation is successful and no unphysical behavior is observed.

Profile functions
From the arguments of the logarithms in the fixed-order hard, jet, and soft functions appearing in eq. (2.11), one can identify the natural scales at which these logarithms are minimized, finding In the tail region, where the resummation is critical, we want to evaluate the distributions near these scales. However, we ultimately predict the distributions over a domain of τ a that can be roughly broken into three regions where the comparative scalings differ (see, e.g. [18]): • Far-tail Region: µ H = µ J = µ S Λ QCD .
In the peak region the soft scale is non-perturbative, and it is here that the full model shape function described in section 4 becomes necessary for making reliable predictions. In this region we will adjust the scales to plateau at a constant value just above Λ QCD . On the other hand, the scales are well separated in the tail region where the resummation is most important. We want to minimize the logarithms in the resummed distributions, and hence the scalings are close to the natural values in eq. (5.1). Finally, our predictions should match onto fixed-order perturbation theory in the far-tail region. The resummations should therefore be switched off, and the scales should merge at µ H,J,S ∼ Q.
Getting the scales to merge near µ H,J,S ∼ Q in the far-tail region will require µ J,S to rise faster with τ a than the natural scales in eq. (5.1), since the physical maximum value of τ a is less than 1. We will achieve this below by defining a smooth function to transition between the resummation and fixed-order regions. But the transition can be made less sudden by increasing the rate of change of µ J,S even in the resummation region. Such an increased slope was used for the C-parameter and thrust distributions in [16]. For thrust, i.e. a = 0, the authors used the central values with r s = 2 in the resummation region. We will follow this strategy here and give a physical interpretation to the slope parameter r s . The maximal value for thrust is τ 0 = 1/2, which is achieved for a perfectly spherically symmetric distribution of particles in the final state. The slope r s = 2 thus ensures that µ J,S merge with µ H at this maximum value τ sph 0 , instead of at τ a = 1 as the natural scales eq. (5.1) do. For arbitrary a, the angularity of the spherically symmetric configuration is which ranges from τ sph −1 ≈ 0.356 to τ sph 1/2 ≈ 0.616 for the values of a we consider in this work. These may be compared to the maximum values of a three-and four-particle configuration in figure 20 in appendix D. We will then choose a default slope r s = 1/τ sph a in eq. (5.2) that ensures that the scales µ H,J,S meet at τ a = τ sph a . We actually want the scales to merge a bit before τ sph a , so that there is a non-vanishing region where the predicted distributions are purely fixed order.
We have designed a set of profile functions (see, e.g. [18,57,87,88]) that fulfill all of the criteria discussed above while smoothly interpolating between the various regions. The precise form of our profiles depends on a running scale defined by The function ζ ensures that µ run and its first derivative are smooth. Specifically, we adopt the functional form from [16], which connects a straight line in the region τ a < t 0 with slope r 0 and intercept y 0 with another straight line in the region τ a > t 1 with slope r 1 and intercept y 1 via where the coefficients of the polynomials are determined by continuity of the function and its first derivative: The parameters t i control the transitions between the non-perturbative, resummation, and fixed-order regions of the distributions, and can be varied as well as part of the estimation  Table 7. Central values and/or parameter ranges we use in the band and scan methods described in section 5.2. For the band method, the lists of parameters are discrete variations whose results are added in quadrature. For the scan method, all parameters are chosen randomly within the ranges shown. It is important to remember that the range of values chosen for any given parameter need not be identical in the two independent approaches we use to estimate the theoretical uncertainties.
of the theoretical uncertainties. We will set these parameters to with coefficients n i that we can adjust. The design of profile functions is somewhat of an art. The chosen a-dependence of t 0,1,2 in eq. (5.7) is based on some empirical observations about the theory distributions ultimately predicted. The first two, t 0,1 , control the transition between the non-perturbative and resummation regions, and we have chosen them to track the location of the peak of the differential τ a distributions. Very roughly, this location scales like 3 a . The parameter t 2 was determined as a numerical approximation to the point where singular and nonsingular contributions become equal in magnitude, since the resummation should be turned off once the latter become as sizable as the former. The formula for t 2 in eq. (5.7) is our rough empirical fit to this crossing point in τ a . Finally, the parameter t 3 is chosen so that our predictions for the distributions revert to those of fixed-order perturbation theory a bit below the maximum value τ sph a , as described above, which we will achieve by choosing n 3 1.
The parameters n i in eq. (5.7), and r and µ 0 in eq. (5.4), will be treated differently in the discussion of section 5.2, where we probe the theoretical uncertainty of our predictions in two different ways: the band and scan methods. In the former, certain parameters including n i , r and µ 0 will be taken as constants, whereas in the latter these will be varied over. The central values and/or scan ranges for all parameters we consider are given for both methods in table 7, and the details of the two procedures are described in section 5.2.

JHEP01(2019)147
Finally then, we use the following forms for the hard, jet, and soft profiles: where varying the parameter e H will adjust the overall magnitude for all the scales together (since µ H enters µ run in eq. (5.4)), and varying e S,J controls the width of the respective soft and jet bands, thereby allowing a variation about the default shape and the canonical relation µ 2−a J = µ 1−a H µ S . In addition to eq. (5.8), our predictions depend on the scale R associated with the renormalon subtraction in the soft function, see eqs. (4.33) and (4.34). For this scale, following [16,18], we use a profile that mimics the soft scale µ S in eq. (5.8b) but that starts at an initial value of R 0 < µ 0 for τ a < t 0 : where our choice for R 0 is given in table 7. As pointed out in [18], a choice like this ensures that the hierarchy between R and µ S is never large enough to generate large logarithms in the subtraction terms δ n a in eq. (4.12), while also keeping a nonzero O(α s ) subtraction δ 1 a (with the right sign) below t 1 where the effect is most pronounced, by strictly taking R < µ S . Finally, our predictions depend on the scale µ ns that enters the nonsingular fixed-order contribution in eq. (4.38). A variation of this scale probes missing higher-order terms in the fixed-order prediction. For small values of τ a , this includes subleading logarithms of the form τ a ln k τ a (in the integrated distribution) (see, e.g., [89]) which are not resummed by our leading-power factorization theorem in eq. (2.11) (see [90] for a resummation of subleading logarithms of thrust in H → gg). To this end, we again adopt the strategy of [16,18] (cf. also [81]) and consider three choices of µ ns as a function of τ a : where n s is a discrete parameter, which we will vary to select between the three different choices of µ ns . In order to obtain a comprehensive theory error estimate, we must consider uncertainties associated with all of the scale and profile parameters that contribute to our calculation. Far from adding arbitrariness to our predictions, these scales and parameters have meanings that reflect the physics contributing to the different regions of the event shape distributions, while the observation of reduced dependence on these parameters from one -38 -

JHEP01(2019)147
order of accuracy to the next provides a stringent test that we have organized the perturbative expansion correctly. We have carried out these variations with two different methods, which we discuss next.

Scale variations
There exist numerous approaches to varying the different scale and profile parameters that enter the theoretical predictions in order to gauge their respective uncertainties. Two methods were contrasted in [18], the so-called band and scan methods. In the first, one parameter is varied at a time and the envelope of the resulting predictions is taken to be the total uncertainty band. In the second, many random choices for all the parameters are taken at once, within predefined ranges for each, and the resulting envelope of the predictions is again taken to be the total uncertainty band. The analysis in [18] concluded that the band method implemented this way tends to underestimate uncertainties, and the authors therefore advocated the use of the scan method. We will apply both methods to all seven angularity distributions treated in this study, but we will purposefully make the band method more conservative by adding uncertainties from each individual variation in quadrature (cf. also [81]). We furthermore adjust the ranges of the parameter variation in each method to achieve convergence plots for the final cross section that exhibit both good convergence and display sufficiently conservative error estimates. In the end, the two methods can be made to give similar results. We undertake this comparison not to advocate one or the other method, but more as a conservative test of the reliability of our displayed uncertainties.
Band method. A simple and intuitive procedure for obtaining uncertainty estimates is to execute independent variations about the relevant scales µ i (i = H, J, S, R with µ R ≡ R) and µ ns , with each variation designed to probe missing orders of perturbative accuracy not calculated in our study. As mentioned above, we then add the uncertainties from the individual variations in quadrature to obtain our final theory error estimate. We have performed six independent variations, with the particular scale settings in each one given by:

JHEP01(2019)147
with the soft scale), we have maintained the relative splitting and hierarchy between µ S and R (and therefore the size of the logarithms of their ratios) by defining R high and R low as 11) and similarly for the functions R R,low and R R,high .
In the left panel of figure 10, we illustrate the width of the jet-and soft-scale variations for a = {−1.0, −0.5, 0.0, 0.5} and Q = m Z . One clearly observes the leveling off in the low τ a region, the natural behavior in the mid τ a region, and finally the convergence of all three scales in the far-tail region. Again, the coloured bands represent independent variations of the jet and soft scales through e J,S ∈ {−0.5, 0.0, 0.5}. The third hard-scale variation then corresponds to shifting the overall scale of the plots in figure 10 up and down, and we do so through Q/2 ≤ µ H ≤ 2Q.
The last three variations are not illustrated in figure 10. The nonsingular scale variation is designed to probe missing fixed-order terms not accounted for in our calculation, and the variation of c 2J and r 2 estimate the systematic uncertainty of our Scan method. For the scan method we take 64 random selections of profile and scale parameters within the ranges shown in table 7. Note that in this method we have chosen not to vary e S away from 0, thus allowing the variations of the other relevant parameters (e.g. n 0 , n 1 , µ 0 , r) to probe the relevant range of the soft scale µ S . Note also that we fix R 0 to always be (µ 0 − 0.4 GeV) so that while µ 0 varies, the difference between µ 0 and R 0 remains fixed, to preserve the scheme choice for the gap parameter ∆ a (cf. [16,18]). The parameters δc 2J and δr 2 are defined to probe the uncertainties coming from scanning over the 1-σ ranges of c 2J and r 2 (τ a ) determined in section 3. That is, we take   methods are of a similar form and they exhibit all of the qualities we demand of the profile functions across the whole τ a domain. Given the current scan ranges and parameter choices, the errors in the scan method indeed appear to be larger than those of the band method, but this is somewhat compensated for by taking the simple envelope of resulting variations rather than adding in quadrature. In the next section we will apply both methods to our final predictions for the resummed and matched cross sections, and we will compare the quality of their convergence both for integrated and differential distributions.
6 Results and data comparison  [21] for α s (m Z ) (to two signficant digits) and Ω 1 (R ∆ , R ∆ ) (to one signficant digit) at NNLL accuracy. Some discussion on the phenomenological impact of choosing different values for these input parameters will be given in section 7. We first show in figure 11 our predictions for the integrated distributions σ c (τ a ) given by eqs. (4.40) and (4.41) for four values of a at Q = m Z up to NNLL + O(α 2 s ) accuracy, including renormalon subtractions. At the same time we compare the two methods discussed in section 5.2 to estimate the overall uncertainties of our analysis. The band method has been applied in the left panel and the scan method in the right panel. One clearly observes that moving to primed accuracies not only dramatically reduces the scale uncertainties, but that also the variations converge across the entire spectra as we move to higher accuracies. One also sees that the two methods that have been applied to estimate the theory uncertainties are consistent with one another, given the parameter ranges and variations chosen. However, when computing differential distributions by taking the derivative of eq. (4.40), we notice a slight improvement in numerical stability when using the scan method. This is partially due to the envelope of the many (64) variations smoothing out wiggles coming from derivatives of profile functions µ i (τ a ), especially in transition regions. Figure 12 illustrates this convergence of the differential distribution (multiplied by τ a ) in the central τ a domain for four values of the angularity parameter a. For this reason, and to allow for a more direct comparison of our results to those of [16,18,21], we implement theory uncertainties as obtained with the scan method in the following.
In order to demonstrate the improvement in precision that we achieve for the differential distributions in moving from NLL + O(α s ) [34] to NNLL + O(α 2 s ) accuracy, we present our (renormalon-subtracted) predictions for a = {−0.5, 0.25} and Q = m Z across the entire τ a domain in figure 13. This figure also illustrates the differences between taking the derivative of the integrated distributions (which we call dσ c /dτ a here) in the left panel, and directly evaluating the resummed ("naïve") differential distributions (which we call dσ n /dτ a here) in the right panel. We see that the former give better convergence and that they better  Figure 11. preserve the total integral under the distributions. These issues with the naïve distributions were extensively discussed in [35]. In fact, as shown there, at unprimed orders the naïve formulas do not even preserve the correct order of accuracy, and even the primed orders suffer from the illustrated mismatch with the total integral under the curve. These issues can be remedied by supplementing the naïve formula with additional terms that both preserve accuracy (at unprimed orders) and maintain agreement with the total integral under the curve (at any order). See [32,35] for such strategies, and [91] for a beautiful mathematical solution to this problem. Here, for simplicity, we have not implemented such strategies, and we simply stick to the integrated distributions from eq. (4.40) as the basis for all of our predictions.
In figure 14 we illustrate the effect of choosing a slope parameter r s = 1/τ sph a for the running scale in eq. (5.4) in the resummation region. In this figure we compare predictions for integrated distribution for a = {−0.5, 0.25} and Q = m Z at different logarithmic accuracies. One observes that the convergence is slightly better for the choice r s = 1/τ sph a than for r s = 1, although the differences between both choices become smaller at higher logarithmic accuracies. This is consistent with the findings of [16], who made a similar comparison for thrust (a = 0), and this also validates our interpretation of the slope parameter r s as being related to the maximum value of the angularity τ sph a .  Figure 13. Differential angularity distributions for a = −0.5 and a = 0.25 at Q = m Z over the entire τ a domain. The distributions correspond to NLL + O(α s ) (orange) and NNLL + O(α 2 s ) (purple) accuracy with renormalon subtractions, and they are either obtained as the derivative of the integrated distributions (left) or directly resummed as differential distributions (right). As expected from the analysis in [35], the former show a bit better convergence than the latter. To obtain the binning, we have integrated the upper and lower bounds of the differentiated cumulative distribution over each bin, which we then divided by the bin width. In both cases we compare to the data from the L3 collaboration [52], and we explicitly illustrate the effect of including the non-perturbative effects from eq. (4.41) (red bins). In the figures we focus on the central τ a domain where the effect of the resummation is most relevant. It is clear that -given the parameters we have applied in our analysis -the data is consistent with the non-perturbative shifted curves for the entirety of the angularity spectra. On the other hand, the bins obtained using only the perturbative predictions (blue bins) systematically undershoot the available data throughout the resummation-sensitive region. 9 These observations are less obvious with the less precise 197 GeV data in figure 16, which nevertheless provides a decent test of the Q dependence of α s and the 1/Q scaling of the leading nonperturbative correction.

JHEP01(2019)147
Slope  These results serve as an impressive visual confirmation of the leading non-perturbative correction discussed in section 4, which clearly indicates a scaling with the angularity parameter a. To further emphasize that this scaling is actually reflected in the data, in figure 17 we have drawn analogous plots to those in figure 15, but we have purposefully altered the non-perturbative scaling coefficient such that c τa → 2.0, instead of the predicted c τa = 2/(1 − a). 10 The associated bins for the predictions convolved with a renormalonfree shape function, now shown in orange, clearly overshoot the data at a = {−1.0, −0.5}, obviously remain unchanged for thrust (a = 0), and are beginning to slightly undershoot some bins at a = 0.5. This indicates that the scaling c τa = 2/(1 − a) is indeed a good theoretical prediction.
Note that the leading bin depicted in all panels of figures 15-17 is actually the first bin after the peak of the distribution, as determined by the experimental data. It is unsurprising that in most cases our theory predictions do not successfully describe this data point, as it is particularly sensitive to the form of the non-perturbative shape function. As described in section 4.1, we have only included the first term in its expansion in a series of basis functions, and one would need to include additional terms in this expansion in order to accurately describe the data in the peak region. The blue bins represent the purely perturbative prediction and the red bins include a convolution with a gapped and renormalon-subtracted shape function, with a first moment set to Ω 1 (R ∆ , R ∆ ) = 0.4 GeV. Overlaid is the experimental data from [52].  The above results serve as a powerful visual confirmation of the predictions of factorization and resummation in QCD using the SCET framework, the evolution of α s , and the scaling and universality of leading non-perturbative corrections for e + e − event shapes. We have not yet performed a full, robust statistical fit to the available L3 data for α s (m Z ) and Ω 1 (R ∆ , R ∆ ), but visually the fit values from [21] appear consistent within our calculations when compared to the data in figure 15. This motivates carrying out such a fit to the available LEP angularity data in the future.

Summary and outlook
We have presented NNLL resummed and O(α 2 s ) matched distributions for the e + e − event shape angularities at Q = 91.2 GeV and Q = 197 GeV. Our results are the most precise achieved for this observable to date, and they are made possible by a recent calculation of the two-loop soft anomalous dimension γ 1 S [46,47] and finite constants c 2 S [46,48]. We determined the remaining unknown NNLL ingredient c 2J from a fit to the EVENT2 generator, from which we also obtained the fixed-order matching to QCD at O(α 2 s ). We further modeled non-perturbative effects with a gapped shape function whose renormalon ambiguity was cancelled using the "Rgap" scheme to define the appropriate subtractions Figure 18. Top: a binned data comparison at a = {−1.0, 0.5} for distributions generated with α s (m Z ) = 0.118 and Ω 1 (R ∆ , R ∆ ) = 0.2 GeV. All other settings are the same as in figure 15. As is clear, the final renormalon-corrected (red) curves overshoot the data at a = −1.0 but, thanks to the small shift parameter, successfully describe the a = 0.5 data. Bottom: difference between central curves and curves evaluated with single variations of either Ω 1 (dashed, purple) or α s (m Z ) (solid, green) at the same values of a and Q = m Z . Note that the lower plots are generated with canonical scales, are unmatched, and the variation about Ω 1 is illustrated using the shift (1.3), rather than the fully shape and renormalon-corrected formulae. and obtain the "R-evolution" of the non-perturbative moment Ω 1 (µ, R). Finally, reliable theory uncertainty estimates have been obtained via variations about a suitable profile function describing the transition of non-perturbative and perturbative scales throughout the relevant domain.
We have compared the predictions resulting from this analysis to data from the L3 collaboration, and find excellent agreement for all angularities considered in this study. In particular, we have visually confirmed the predicted sensitivity of the distributions to leading non-perturbative scaling effects that are simultaneously encoded in the angularity parameter a and the universal shift parameter Ω 1 -the distributions indeed seem to require shifts proportional to 2/(1 − a) in order to accurately describe the data.
This characteristic non-perturbative scaling strongly motivates an extraction of both α s and Ω 1 from a fit to experimental data on angularities. As previous extractions of α s from e + e − event shapes like thrust and C-parameter have shown [18,19,21], the dis--50 -

JHEP01(2019)147
entangling of perturbative and non-perturbative effects is crucial and has a large impact on the determined value of α s itself. Thus angularities, with their tunable relative contributions of these effects, can play an important role in confirming (or not) the values coming from these other extractions, which as mentioned before, are considerably lower than the PDG world average [27]. While we have not performed a dedicated extraction in this study, in the top panels of figure 18 we have shown the analogous plots to those in figure 15 at a = {−1.0, 0.5}, but with α s (m Z ) = 0.118 (consistent with the PDG value) and Ω 1 (R ∆ , R ∆ ) = 0.2 GeV, just large enough to consistently implement our gap parameter and provide a good fit to the data for the angularity a = 0.5. As can be seen, however, at a = −1.0 the final theory bins are then well above the data in the tail region. As the practical effect of lowering (raising) α s (m Z ) (Ω 1 ) away from the best values used in figure 15 is to the lower (raise) the theory prediction, this naïve exercise leads one to provisionally expect that such parameter values will be difficult to reconcile with the available data.
In the bottom panels of figure 18 we finally show the difference (dσ/dτ a ) central −dσ/dτ a over the range 0.085 ≤ τ a ≤ 0.35 for a = {−1.0, 0.5}, where (dσ/dτ a ) central is an (unmatched, non-renormalon corrected) NNLL resummed curve. For (dσ/dτ a ) we have varied 2Ω 1 by ± 0.1 GeV and α s (m Z ) by ± 0.001, corresponding to the purple and green curves, respectively. Our plots are to be compared to figure 10 in [18], where the same variations were made at different center-of-mass energies Q. We find that varying a (Q) down (up) from high (low) values leads to an enhanced sensitivity to the relative effects of α s and Ω 1 variation. This suggests that extractions from angularity data at a single centre-of-mass energy Q, but different values of a, will be able to discriminate between α s and Ω 1 in a similar way to varying Q. In fact, angularities for −2 ≤ a ≤ 0.5 exhibit a factor of six variance in the overall non-perturbative shift, which represents an essentially equivalent sensitivity to measurements made between Q = 35 GeV and Q = 207 GeV, as analyzed for thrust in [18]. Unfortunately, angularity data is currently only available down to a = −1.0 and only at two values of Q. Thus a reanalysis of additional LEP data would be highly desirable. Of course, a high-precision determination of α s and Ω 1 should ideally also include further theory improvements beyond those we have already implemented here, such as hadron mass effects or higher-precision numerical (or even analytical) determinations of c 2J and r 2 (τ a ). New methods and calculations for O(α 3 s ) jet and soft functions at a = 0 would be needed to go to N 3 LL( ) accuracy, while O(α 3 s ) fixed-order matching to QCD could be obtained from a program like EERAD3 [4]. In conclusion, we are optimistic that the a-dependence of the angularities may ultimately be helpful in lifting degeneracies between α s and Ω 1 in the two-parameter event-shape fits. We leave it to future work to bring such an analysis to full culmination.

JHEP01(2019)147
the EVENT2 data. We also thank Massimilano Procura for discussions. JT thanks the T-2 Group at LANL, and CL in particular, for hospitality and support while much of this work was completed. JT also acknowledges research and travel support from DESY and the Senior Scholarship Trust of Hertford College, Oxford. The work of AH and CL was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under Contract DE-AC52-06NA25396 and an Early Career Research Award, and by the LANL/LDRD Program.
Finally, we would like to dedicate this work to the memory of Andrew Hornig (1982-2018) who passed away suddenly before its final publication.

A Laplace transforms
In this appendix we collect results for the Laplace transforms and inverse Laplace transforms between the logarithms L ≡ ln 1 τ ,L ≡ ln(νe γ E ) . where γ lies to the right of all the poles of F in the complex plane, are given by

B Anomalous dimensions
The coefficients of the beta function up to four-loop order in the MS scheme are given by [92][93][94]  The running coupling α s (µ) up to four-loop accuracy is given by the formula, The four-loop term in eq. (B.3) agrees with the numerical form given in [18], and is rederived below in appendix C.
The cusp anomalous dimension coefficients are given up to 3-loop order [95,96] Γ q 0 = 4C F , (B.5) T F n f , The MS non-cusp anomalous dimension for the hard function can be obtained [55,60] from the infrared divergences of the on-shell massless quark form factor, which are known evaluating eq. (C.4a), so the canceling logs are treated symmetrically. But it is straightforward to show that this actually leads to a final expression for the cross section still plagued by mismatches like eq. (C.13) and consequently a (subleading) violation of µ-invariance: the problem is that ln(1/τ  similar to eq. (2.17). It is now easy to show that the combination is by itself exactly µ-invariant, at any order of resummed accuracy, not just LL. This leads us to reorganize the pieces of the cross section eq. (2.1), where K γ is the sum of just the non-cusp evolution kernels in eq. (2.7): K γ (µ H , µ J , µ S ) ≡ K γ H (µ, µ H ) + 2K γ J (µ, µ J ) + K γ S (µ, µ S ) , (C. 19) which is exactly µ-invariant by itself at any order of accuracy, and where we will evaluate K, Ω using closed-form expressions like eqs. (C.2) and (2.18), replacing eq. (C.2) for K Γ with the expansion for K Γ we have given in eq. (2.22). Note the expressions in eq. (2.22) for K Γ reduce to those for K Γ in eq. (C.2) for Q = µ F or r Q = 1, as they must. The departure of r Q from 1 is an O(α s ) effect, so the difference between K Γ and K Γ at each order is always subleading, consistent with our observations above. However, the r Q terms precisely compensate for the subleading terms that violate the µ invariance arising from using eq. (C.5), restoring exact µ invariance of K in eq. (C.17).
To see the small numerical effect of using the µ-invariant eq. (C.18) in place of the original form eq. (2.1), we show in table 8 the predictions of these two formulas for the purely perturbative resummed integrated distributions σ c,sing /σ 0 given by eqs. The maximum values of τ a for the three-particle (blue), four-particle (orange) and spherical distribution of particles (green) as a function of a.
For some extreme values of a, it is possible that the maximum τ a configuration is not the symmetric one in figure 19, as occurs for the three-particle configuration for a −2.6 [34]. We assume that this does not occur in the four-particle configuration for the values of a we consider, −1 ≤ a ≤ 0.5. We in fact did not notice any anomalous deviation in the endpoints of our O(α 2 s ) distributions away from this formula. The maximum three-and four-particle values of τ a are illustrated in figure 20.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.