Interplay between perturbative and non-perturbative effects with the ARES method

We present a new semi-numerical method to compute leading hadronisation corrections to two-jet event shapes in $e^+e^-$ annihilation. The formalism we present utilises the dispersive approach, where the magnitude of power corrections is controlled by suitable moments of an effective strong coupling, but it can be adapted to other methods. We focus on observables where the interplay between perturbative and non-perturbative effects is crucial in determining the power corrections. A naive treatment of power corrections for some of these observables gives rise to an unphysical behaviour in the corresponding distributions for moderate observable values, thus considerably limiting the available range to fit the non-perturbative moments. We present a universal treatment to handle such observables, based on a suitable subtraction procedure, and compare our results to the analytic result in the case of total broadening. Finally, for the first time we present predictions for the thrust major, which cannot be handled with analytic methods.


Contents
Monte Carlo simulation of multiple soft and collinear emissions, accompanied by at most one "special" emission triggering NNLL corrections. In the present case, the special emission is an ultra-soft gluon. This is the first step to having predictions for event-shape distributions, valid in the twoand the three-jet regions, which will be needed for a new and more accurate global simultaneous fit of α s and α 0 . The paper is organised as follows. In section 2 we describe the method and its scope. In section 3 we use the method to compute the NP corrections to known event-shape distributions and to the thrust major. The latter, which gives the main novel result of this paper, does not admit an analytic treatment, and can only be handled with the semi-numerical approach presented here. For the total broadening and the thrust major, the predictions presented in section 3 cannot be naively matched to the three-jet region due to an unphysical divergence arising in the three-jet region. Such divergence can be "cured" in a general way, which is the topic of section 4. After a brief discussion in section 5 on how to treat the mean values of event shapes, in section 6 we perform simultaneous fits of α s (M Z ) and α 0 (µ I ) using NLL resummations matched to NLO and our new determination of 1/Q power corrections. The concluding section 7 contains a summary of our main findings and further research directions.

Leading hadronisation corrections in the two-jet region
In this section we introduce a general method to compute leading power corrections to a large class of event shapes in the two-jet region. Before explaining the method in its full generality, it is instructive to look at the known examples of the thrust and heavy-jet mass, which display some of the theoretical issues involved.

Two examples: the thrust and heavy-jet mass
Starting with kinematics, the thrust axis splits an event in two hemispheres H 1 and H 2 . Let ρ 1 be the invariant mass squared of hemisphere H 1 (normalised to Q 2 , the centre-of-mass energy squared of the e + e − collision) and similarly ρ 2 the invariant mass squared of hemisphere H 2 , in formulae where each q i denotes the momentum of a final-state hadron. From ρ 1 and ρ 2 , we can construct the heavy-jet mass ρ H = max(ρ 1 , ρ 2 ). Also, close to the two-jet limit, where the thrust T is close to one, we have 1 − T ≃ ρ 1 + ρ 2 . The bulk of events contributing to the distributions in the thrust and heavy-jet mass comes from soft and/or collinear partons, still with momenta much larger than Λ QCD , so that their emission probabilities can be safely computed in perturbative QCD. In this regime, one can approximately write the fraction of events such that 1 − T < τ or ρ H < ρ in terms of the distribution J q (k 2 ), the probability that a quark-initiated jet has an invariant mass k 2 [42], as follows (2. 2) The above distributions are better expressed in terms of the Laplace transform of J q (k 2 ), defined asJ This gives where the contour of the ν-integration runs parallel to the imaginary axis, to the right of all the singularities of the integrand. Following the approach of [15], looking at the explicit expression ofJ q (ν) reported there, we observe that it contains a contribution arising from soft radiation, as follows The coupling α s (q) in eq. (2.5) is to be interpreted as a non-perturbative effective coupling. The main result of [15] is that one can single out in the q-integral an ultra-soft contribution with q < µ I , where Λ QCD ≲ µ I < τ Q. This ultra-soft contribution produces a number of corrections suppressed by powers of 1/Q, which arise from expanding e −νuQ 2 in eq. (2.5). Following [15], the first order in the expansion gives the leading 1/Q non-perturbative correction where α NP s (q) is the non-perturbative component of the physical coupling α s (q). The representation of eq. (2.6) assumes that ultra-soft contributions toJ q (ν) have the same form as perturbative contributions, just with a different coupling. One can take a more general approach and model the whole δ lnJ q (ν) as a non-perturbative shape function. In this case, the first term in the expansion in ν will still have the form as in eq. (2.6), but ∆ will not be interpreted as the moment of an effective coupling, but rather as a matrix element of Wilson lines [43]. Let us now write lnJ q (ν) as the sum of a fully perturbative contribution, lnJ PT q (ν), 1 and the leading NP correction −νQ 2 ∆. For the thrust distribution, this gives where (2.9) Keeping the leading power of ν in δ lnJ q (ν) gives rise to a shift in the perturbative thrust distribution. If we analyse more closely the assumptions under which the shift approximation works, we see that it holds mathematically because we have been able to expand the exponential e −νuQ 2 in eq. (2.5). But what is the physical meaning of such expansion? It is known that the values of ν giving the largest contribution in eq. (2.9) are of the order 1/(τ Q 2 ). Also, the dimensionless variable u is just the contribution to the thrust due to the ultra-soft gluon. This is of the order µ I and hence much smaller than τ Q 2 . In this limit, we have ∆ ∼ µ I /Q ≪ τ ∼ 1/(νQ 2 ), hence we can expand Σ(τ ) as follows . (2.10) 1 The actual form of lnJ PT q (ν), which depends on the perturbative order considered, is not relevant for our discussion. The only important aspect is that the coupling appearing in the integrals defining lnJ PT q (ν) is perturbative. This implies that lnJ PT q (ν) produces a renormalon divergence, proportional to 1/Q, which by construction cancels against δ lnJq(ν).
The approximation in the last line of the above equation is not needed for event shapes such as the thrust, heavy-mass, C-parameter [24] and D-parameter [44], whose shift is just a constant value, independent of the Laplace variable ν. The distributions in those event shapes can be expressed in terms of a Laplace transform, where the NP shift appears manifestly in an exponentiated form as in eq. (2.8). There are cases, such as the jet broadenings, see e.g. [39,45], where the shift itself depends on the Laplace variable ν. In that case, under the assumptions that µ I /Q is much less than the observable's value, one only obtains a shift in the approximate way indicated by eq. (2.10).
At this point, an important remark is in order. While in the case of the thrust the nonperturbative shift is completely uncorrelated with perturbative contributions, this is not the case for the heavy-jet mass. In fact, for its integrated distribution, we obtain with (2.12) Why is the shift of the heavy-jet mass half of that of the thrust? This is because the heavy-jet mass collects contributions from all hadrons in the heavier hemisphere. If an ultra-soft gluon is emitted in the lighter hemisphere, as long as ρQ > µ I , it will not be able to contribute to the heavy-jet mass. This is the simplest example of interplay between perturbative and non-perturbative contributions, which is the main topic of this article. Let us comment now on the physical meaning of the variables q and u in eq. (2.5). The variable q is the transverse momentum k t of the soft emission with respect to the emitting parton (quark or antiquark). For the thrust and heavy-jet mass, this is the same as the transverse momentum with respect to the thrust axis. The variable uQ 2 is the contribution of the ultra-soft gluon to the invariant mass of the jet. It can be used to define the rapidity η with respect to the emitting quark or antiquark as η ≡ ln k t uQ . (2.13) Using the variables k t and η, we can express ∆ as follows (2.14) Since the η integral is exponentially damped in rapidity, the rapidity integral can be extended up to infinity instead of the appropriate kinematic limit. The exponential damping implies also that the only ultra-soft emissions that give an appreciable contribution to ∆ are those at η ∼ 0. Note also that, in the assumption that the distribution of ultra-soft gluons is the same as that of perturbative soft gluons except for a different coupling, the rapidity integral (alternatively the u integral in eq. (2.5)) can be performed. In the following, we assume that leading non-perturbative corrections to event-shape distributions have the same origin as for the thrust and the heavy-jet mass, to summarise: • they are due to ultra-soft emissions with k t ≪ v p Q, where v is the value of the considered event shape and p is some positive power, typically one. This is in contrast to perturbative emissions where k t ∼ v p Q.
• they give rise to a shift of the corresponding perturbative distributions; the shift can be computed by considering a single ultra-soft emission (with accompanying virtual corrections).
• the dynamics of these ultra-soft emissions is not known; in the case that non-perturbative corrections are modelled as a shape function, an ultra-soft emission k is produced with an unknown matrix element squared M 2 NP (k); with the additional assumption that the distribution of ultra-soft emissions follows that of soft perturbative gluons, they are emitted uniformly in rapidity and azimuth, with an unknown dependence on transverse momentum, which can be embodied in a soft effective coupling.
As a closing remark, we note that different assumptions on the distribution of ultra-soft emissions result in various degrees of observable dependence for the shift. For instance, if we assume that non-perturbative corrections to the thrust and heavy-jet mass arise from a shape-function modifying the invariant mass distribution of a single jet, the shift in the distributions of these observables will depend on the same non-perturbative quantity ∆. As we will see later, assuming that ultra-soft emissions are produced uniformly in rapidity makes it possible to relate the shifts of many more observables.

Kinematics of ultra-soft emissions
In many cases, event shapes are defined using the thrust axis as a reference axis to project momenta upon. However, the thrust axis is generally not a direction that corresponds to a singularity of the QCD matrix elements, which are better parameterised in terms of the directions of the actual finalstate particles. It is therefore useful to relate different parameterisations of soft-gluon momenta. To be precise, we consider 2-jet events in e + e − annihilation. The thrust axis defines two light-like vectors (p 1 , p 2 ), in terms of which the final-state qq pair, (p 1 ,p 2 ), perturbative emissions, k i , and the ultra-soft gluon, k, have the following Sudakov decomposition where k ⊥ etc. are space-like vectors with purely transverse components, i.e. k ⊥ = (0, ⃗ k t , 0), p ⊥1 = (0, ⃗ p t,1 , 0) and so on. The rapidity of the emission with respect to the thrust axis, η, is given by Let us consider the case in which emission k is collinear top 1 . The collinear singularity is encoded in the following transverse momentum Here κ (1) ≡ |⃗ κ (1) | defines the transverse momentum with respect top 1 . From the properties of the thrust axis, we have that the total transverse momentum vanishes in each hemisphere and thus we have where ⃗ p t,1 denotes the recoil in the hemisphere due to all perturbative emissions. More precisely, for each leg ℓ we can define Using the fact that these are soft, we can approximate z (1) p ≃ 1 in eq. (2.17) up to corrections suppressed by powers of v. Neglecting contributions quadratic in κ (1) , we obtain Given κ (1) , we can define the rapidity η (1) with respect to thep 1 direction as follows Similarly, we define suitable kinematic variables for the case in which k is collinear top 2 The boundaries for η (ℓ) (ℓ = 1, 2) are a bit involved. Nevertheless, identifying the region collinear top ℓ with hemisphere H ℓ , we have that the limits on the Sudakov variable k t /Q < z (ℓ) < 1 can be recast in terms of η (ℓ) and κ (ℓ) as follows where ϕ is the azimuthal angle with respect to the direction of ⃗ p t,ℓ . In the presence of soft emissions only, p t,ℓ ≪ Q and thus the lower limit of η (ℓ) vanishes.

General treatment of leading non-perturbative corrections
We now consider a generic recursive infrared and collinear safe (rIRC) safe [46] observable, denoted by V ({p}, k 1 , . . . , k n ), in e + e − annihilation. Here {p} = {p 1 ,p 2 } are the momenta of a hard quarkantiquark pair and k 1 , . . . , k n subsequent emissions. Without any additional emissions {p 1 ,p 2 } are the momenta of a back-to-back quark-antiquark pair. In this case V ({p}) = 0, whereas in general V ({p}, k 1 , . . . , k n ) ≥ 0. We consider the region in which V ({p}, k 1 , . . . , k n ) = v ≪ 1, i.e. we are close to the Born limit, but v is not too small, so that the value of the observable is determined by perturbative QCD emissions. To enforce this condition, we restrict ourselves to the region v ≫ Λ QCD /Q, with Q ∼ √ s the typical hard scale of the process. In this region, the observable cumulant Σ(v), the fraction of events such that V ({p}, k 1 , . . . , k n ) < v, can be in largest part computed in perturbative QCD, with small non-perturbative corrections, as follows: Here dZ[{k i }] is an integration measure, associated with multiple emissions and the corresponding virtual corrections, with the normalisation In particular, if we consider only soft and collinear emissions widely separated in angle, the relevant configurations giving NLL accuracy [46], we have where [dk] is the Lorentz-invariant measure for a gluon of four-momentum k = (ω, ⃗ k) . (2.28) In terms of these variables, the matrix element squared for a soft emission off the qq dipole, M 2 (k), is given by where α CMW s is the strong coupling in the physical CMW scheme [42]. Here κ is the invariant transverse momentum of emission k with respect to the hard qq pair, defined as For emissions collinear to leg ℓ, κ ≃ κ (ℓ) introduced in section 2.2. Last, for such configurations, for any function G({p}, k 1 , . . . , k n ), we define As for the cases of the thrust and heavy-jet mass, we model leading NP corrections as the contribution of an ultra-soft emission k, whose transverse momentum is much smaller than the typical transverse momentum of soft emissions contributing to Σ PT (v). The latter is of the order v p Q, with p some positive power. The ultra-soft emission is produced with an unknown matrix element squared M 2 NP (k) and gives the following correction to Σ(v) where virtual corrections are implemented via unitarity, i.e. by imposing that this contribution is zero if there are no constraints on any emissions. 2 The difference between the observable with an additional ultra-soft emission and without it is in general much smaller than v in the region we are interested in. Therefore, we can approximate the difference between the two step functions in eq. (2.32) as follows where we have introduced the change in the observable due to a NP emission k, as follows (2.34) Substituting the approximation in eq. (2.33) into eq. (2.32), we obtain where we have introduced the average of δV NP over all perturbative configurations, defined as Observing also that we can write This means that, whenever v ≫ ⟨δV NP ⟩, NP corrections amount to a shift of the perturbative distribution by an amount given by the average ⟨δV NP ⟩.
In principle, ⟨δV NP ⟩ is different for each observable and might depend non-trivially on unknown NP dynamics. However, here we make the assumption that the production rate of ultra-soft emissions is driven by PT dynamics, i.e. it is uniform in rapidity, defined with respect to final-state partons, and azimuth, as follows where η (1) (η (2) ) have been defined in Sect. 2.2. In general, the ultra-soft emission k is accompanied by soft and/or collinear emissions. In this paper, we assume these accompanying emissions are soft, collinear and widely separated in angle. This is appropriate to achieve NLL accuracy for Σ PT (v).
We also restrict ourselves to event-shape variables (i.e. we exclude jet-resolution parameters) and furthermore only to those event shapes for which the NP contribution to the observable is linear in κ, as follows This implies and (2.43) Furthermore, we restrict ourselves to observables where the rapidity integral in eq. (2.43) is convergent when the upper bound is pushed to infinity. For such observables, setting the rapidity boundary to the actual kinematic limit ln(Q/κ) would give a contribution that has a further suppression in κ, i.e. a sub-leading power correction. In eq. (2.43), V sc ({p}, k 1 , . . . , k n ) is the value that V ({p}, k 1 , . . . , k n ) assumes when all emissions are soft and collinear. Formally, in the presence of only soft and collinear emissions k 1 , · · · , k n , An important remark is in order here. Since V ({p}, k 1 , . . . , k n ) is a rIRC safe observable, all soft and collinear emissions k 1 , . . . , k n have momenta in the region V ({p}, k i ) ∼ v. This is in contrast with the contribution of the NP emission, δV NP ∼ κ ≪ vQ.
In the dispersive approach of [14] the relation between ⟨κ⟩ NP and the phenomenological parameter α 0 is given by Before moving forward, we comment on the relevance of the result in eq. (2.41). It means that, for all event shapes for which δV NP has the property in eq. (2.40), the NP shift to the corresponding distributions is given by the product of a genuine NP quantity ⟨κ⟩ NP , the average transverse momentum of ultra-soft emissions, and a calculable coefficient ⟨h V ⟩. The aim of this paper is precisely to devise a semi-numerical procedure to compute ⟨h V ⟩ for a generic rIRC safe event shape whose leading NP corrections are given as the shift in eq. (2.41). The general procedure is the same as used to compute NLL and NNLL corrections to event-shape distributions. First, we rewrite the measure dZ sc [{k i }] as where the "radiator" R(v) is given by For perturbative emissions which are soft and collinear, it is convenient to express any soft momentum k in terms of the leg ℓ it is collinear to, the variable ζ = V sc ({p}, k)/v and its azimuth ϕ with respect to a suitably chosen plane. For event shapes only, which is what we consider here, we can freely integrate over the rapidity fraction of emission k with respect to the total available rapidity for fixed ζ and ϕ [46]. This gives where R ′ ℓ (ζv) is the Jacobian arising from the integration over the rapidity fraction. Therefore, the soft-collinear phase space measure reads where we have introduced a modified soft-collinear measure normalised in such a way that The quantity ϵ in eq. (2.51) acts as a cutoff for the integration over the variables ζ i . Considering again the calculation of ⟨h V ⟩, the factor exp[−R(v)] cancels out between numerator and denominator, and we obtain (2.52) We know that, at NLL accuracy and keeping only soft contributions, we have 3 One can also show (see appendix A) that We can then recast the denominator of eq. (2.52) in terms of F(R ′ ) and obtain (2.56) This expression is suitable for both analytic calculations and numerical determinations according to the method explained in appendix A. The idea, developed originally in [47], is to label k 1 the emission such that ζ 1 is the largest of the ζ i , neglect all emissions with ζ i < ϵζ 1 and use the constraint on the observable to perform the integration over ζ 1 . This gives

The Milan factor
Due to the non-inclusive nature of event shapes, it was initially suspected [19] that the leading hadronisation correction is not universal, which cast doubts about the utility of the dispersive approach in predicting the NP parameter α 0 . Using the thrust variable, [23] considered the hadronisation corrections due to the decay of the NP gluon to a gluon (or qq) pair. The non-inclusiveness of the thrust variable turned out to be encoded in a calculable multiplicative constant, called the Milan factor, multiplying the leading moment of the effective coupling. In this section we describe, within our approach, how the Milan factor arises. The Milan factor is historically given by the sum of two contributions. One is the 'inclusive' contribution, where one upgrades the NP corrections due to one emission by dressing the ultra-soft gluon with its inclusive splittings, as well as virtual corrections. This contribution, computed in [23,24], is constructed in such a way as to give a multiplicative factor times the contribution of a single ultra-soft gluon. The remaining contribution is the so-called 'non-inclusive' correction, arising from the difference of the equivalent of eq. (2.33) in the presence of two NP emissions (k a , k b ), and of their inclusive contribution.
The introduction of an inclusive contribution is just a convenience and, as shown in [23,24], it is possible to obtain the Milan factor for the thrust and C-parameter by combining the contribution of single and double soft emission and the corresponding virtual corrections. We are not interested in computing the Milan factor here, but want to show that for certain observables the contribution of two ultra-soft emissions gives rise to the same ⟨h V ⟩ as one ultra-soft emission. As the contribution of two ultra-soft emissions has a collinear divergence, we employ an inclusive subtraction as a technical convenience, in a similar way as subtraction terms are employed in fixed-order calculations. The so-regularised contribution of two ultra-soft emissions to the observable constraint reads We can now express the non-inclusive correction in our approach as follows .
(2.60) The crucial assumption now is that the ultra-soft matrix-element squared M 2 NP (k a , k b ) has the same form as the perturbative one, where the physical, i.e. infrared finite, coupling is assumed to be a function of the invariant mass m 2 = (k a + k b ) 2 . In particular, M 2 NP (k a , k b ) is independent of the rapidity of the parent ultra-soft gluon, k = k a + k b , while it depends on the azimuthal angle difference ϕ a − ϕ b solely through the dependence on the invariant mass. The two-body phase space can be cast in terms of the following standard variables [23] [ The Dirac delta function is used above to keep the invariant mass fixed. One can easily use the delta function to integrate over the azimuthal angle difference ϕ a − ϕ b , and we find where ϕ is a residual azimuthal angle, i.e. ϕ a , ϕ b or the azimuth of k a + k b , which the squared matrix-element does not depend upon. We do not need the explicit form of M 2 NP (k a , k b ), but only its generic dependence on the phase space variables (2.64) In terms of these variables we define the 'inclusive' observable one can integrate freely over η (ℓ) and ϕ.
a and η (ℓ) b , then we must express the observable function h V in terms of η (ℓ) which requires the following relations (2.67) We notice that the relation between η (ℓ) a , η (ℓ) b and η (ℓ) corresponds to a boost that depends on u a , u b and z. Last, we need to discuss the integration limits on η (ℓ) . As for a single massless gluon, the limits on the individual rapidities are as follows 0 < η (ℓ) a,b < ∞, therefore, we integrate the functions and thus obtain The 'residual' two-body phase space is given by It has been shown for the first time in [23] that, in the dispersive approach of [25], we have with M (n.i.) the non-inclusive part of the Milan factor whose correct numerical form is given in [24,48], and was later determined analytically in [49]. The analysis we presented for the Milan factor confirms that, in our approach, any observable satisfying eq. (2.66) exhibits a non-inclusive contribution to the leading hadronisation correction proportional to the contribution of a single ultra-soft gluon multiplied by the non-inclusive contribution to the Milan factor, viz.
Therefore, considering also the inclusive contribution to the Milan factor, the total shift in our approach is given by where M is the full Milan factor.
3 NP shifts for "popular" e + e − event shapes The function h V (η (ℓ) , ϕ, {p}, {k i }) needs to be computed separately for each observable. In this section we perform the calculation for rIRC safe event shapes that have been routinely measured at LEP, namely thrust, C-parameter, heavy-jet mass, jet broadenings and thrust major. Among those, only the thrust major has not been computed so far. This is due to the fact that this observable does not admit an analytical treatment in the two-jet limit. Therefore, only a semi-numerical approach such as the one presented here is viable to compute the NP shift to its distribution. From the discussion below it will be evident that event shapes can be grouped into classes showing a similar structure for h V (η (ℓ) , ϕ, {p}, {k i }). However, it is not clear to us whether a general classification can be formulated and hence we leave this question for future work. In the following we derive the expressions for the NP shifts for the above observables, as provided by our method.

Recovery of known results
We first show that with our method we can compute the shift to all observables that have been studied analytically so far and obtain perfect agreement.
Thrust and C-parameter: The function h V (η (ℓ) , ϕ, {p}, {k i }) for one minus the thrust, 1−T , and C-parameter does not depend explicitly on the momenta of PT soft-collinear emissions In this case, one has The lower limit on η (ℓ) can be set to zero because p tℓ ∼ √ vQ for such observables. Also, which implies that the upper limit on η (ℓ) can be set to infinity due to exponential damping. This gives Heavy-jet mass: In the case of the heavy-jet mass there is an interplay between NP and PT radiation since, in the presence of perturbative emissions, a non-zero NP correction to the heavyjet mass arises only when the ultra-soft emission is in the heavier hemisphere. This is due to the fact that the contribution to the invariant mass of either hemisphere due to a NP emission is of order κ, which is much less than the contribution of all PT emissions, which is of order ρ ℓ . As a consequence, in the presence of multiple soft and collinear PT emissions, a NP emission can never determine which hemisphere is heavier. Therefore where ρ 1 and ρ 2 are the invariant masses of the two hemispheres. This gives Total and wide-jet broadening: The interplay between PT and NP effects is more complicated for the total and wide-jet broadening than it is for the heavy-jet mass. To better understand the issues involved, we consider an ultra-soft emission k in the hemisphere containing leg ℓ. The broadening of the corresponding hemisphere (which we denote by B ℓ ) is then given by From the definition of the thrust axis, we havẽ with ⃗ p t,ℓ defined in eq. (2.19). The hadronisation contribution to B ℓ is then given by We remark that k t is the transverse momentum of the NP emission with respect to the thrust axis and not with respect to its emitter. In fact, k t is different from κ and depends also on η (ℓ) , ϕ and the recoiled momentum p t,ℓ . More precisely, if we express all quantities in eq. (3.10) in terms of κ, η (ℓ) and ϕ, we find where in the last line we have omitted a term that vanishes after integration over the azimuthal angle ϕ. This gives, for the total broadening For the wide-jet broadening, where the ultra-soft emission gives a contribution only if it is emitted in the broader hemisphere, we get (3.14) Notice that when integrating h B T and h B W we can take the limits on η (ℓ) to be between 0 and infinity since the integral is convergent. Thus we obtain the following integral, with p > 0, Since p is the transverse recoil induced by perturbative emissions, all terms that vanish as a power of p/Q give contributions that are at most suppressed by powers of B T , B W and hence can be neglected. With this approximation, our final result for ⟨h B T ⟩ and ⟨h B W ⟩ reads Remarkably, the fact that the hard legs are displaced from the thrust axis has profound implications for the NP shift to the jet broadenings. In fact, from eq. (3.15), one can see that the rapidity of an ultra-soft emission collinear to leg ℓ is effectively cut at ln(Q/p t,ℓ ), which is far away from the collinear limit ln(Q/κ). This is consistent with our treatment of NP corrections, which requires that the ultra-soft emission is at large angles.
To explicitly compute ⟨h B T ⟩ and ⟨h B W ⟩, we further manipulate eq. (3.17) by rescaling p t,ℓ with respect to B T , B W as follows and by introducing the two-dimensional vectors ⃗ ζ i ≡ ζ i (cos ϕ i , sin ϕ i ). This gives where (3.20) As has been noted already in [39], for small B T , B W the shift depends logarithmically on the event shape's value, with corrections given by η (B) 0 and the functions χ T (R ′ ) and χ W (R ′ ). The latter have been computed analytically in appendices B.2 and B.1 respectively and agree with the results of [39].
We also compute the NP shift fully numerically and compare with the analytic calculation. This is shown in Figs. 1 and 2. We observe a perfect agreement between the analytic calculation and the output of our numerical program. We also plot the limiting behaviour of the shift for This behaviour is also reflected in numerical issues when calculating ⟨h B T ⟩ with the procedure described in appendix A. In fact, it may transpire, especially when R ′ is small, that only one hemisphere is populated by a small number of emissions with k t ∼ B T Q, whereas all emissions in the other hemisphere (which we can refer to as the "empty" hemisphere) fall below the cutoff of the Monte Carlo integration. As a consequence, one would find a zero value of p t,ℓ in the empty hemisphere and eq. (3.20), which requires the calculation of ln(Q/p t,ℓ ), would give a floating point exception. This issue is ultimately due to the fact that the Monte Carlo procedure of appendix A assumes that only perturbative emissions with comparable transverse momenta contribute to the shift, which is clearly not the case for p t,ℓ of the empty hemisphere. In order to obtain finite numerical predictions we are forced to decrease the cutoff ϵ more and more as R ′ approaches zero.
We reiterate that this problem is of physical and not technical nature. Its ultimate solution requires upgrading the probability of soft-gluon emissions in the empty hemisphere as was done in [39]. This leads to a finite value for ⟨h B T ⟩ even for R ′ → 0. The generalisation of this procedure to an arbitrary event shape will be discussed in section 4.

NP shift for the thrust major
Another observable that has a similar behaviour to the jet broadenings is the thrust major. In this case, only PT soft and collinear emissions determine the thrust-major axis ⃗ n M , which we conventionally set to give the y-direction for all momenta. In fact, when varying the trial direction for the thrust-major axis, one finds a set of local maxima which are the candidates for the magnitude of T M . In the presence of soft and collinear PT emissions that contribute to the observable in the two-jet region, all these local maxima are of order T M as all emissions have comparable transverse momenta. The addition of an ultra-soft emission, with κ ≪ T M Q, will not change the hierarchy of the local maxima and thus will not change the value of T M up to sub-leading power corrections. This is similar to the argument that is used to compute NP corrections to the heavy-jet mass. Ultimately this implies that the thrust-major axis is not altered by a NP emission. This gives Using the transverse momentum component p y,ℓ defined in eq. (2.19), the change in the thrust major due to a NP gluon k collinear to leg ℓ is given by We now need to recast the quantities in square bracket in terms of the transverse momentum with respect to the emitter κ y = κ sin ϕ. Using the symmetry properties of the ϕ-integral, we obtain that δT M depends only on |p y,ℓ |. This gives

(3.23)
In the last equality we have dropped a term that vanishes upon integration over ϕ. We then obtain Therefore, the η and ϕ integrals now give, with p > 0 The constant ln(2e −2 ) is the same as found for the thrust minor [45]. This is expected because, in both cases, it is harder emissions that fix the thrust-major axis and here we are considering the magnitude of one of the two components of ⃗ p t,ℓ . Inserting the above expression into the general formula for ⟨h V ⟩ gives In a similar way as for the broadenings, we can introduce the rescaled variables x ℓ = 2p y,ℓ /(T M Q) and the two-dimensional vectors ⃗ ζ i ≡ ζ i (cos ϕ i , sin ϕ i ) and obtain (3.28) We first observe that, similar to the broadenings, the NP shift for the thrust major has an explicit logarithmic dependence on T M . The function χ M (R ′ ) can only be computed numerically with the procedure outlined in appendix A. This calculation constitutes the first new result of this paper.
The resulting value of ⟨h T M ⟩ is plotted as a function of T M in Fig. 3, and compared to its limiting behaviour for R ′ → 0 (which has been computed analytically in eq. (C.37)). Figure 3: The numerical calculation of the NP shift for T M , compared to its limiting behaviour for R ′ → 0 (with α s = 0.118).
As was the case for B T , the shift for T M has a minimum and increases for large T M , or equivalently small values of R ′ , where a 1/R ′ singular behaviour dominates the shift. The calculation of χ M (R ′ ) has the same problem as that of χ T (R ′ ). In fact, when R ′ → 0, there will be one emission with k t ∼ T M Q in one hemisphere and all emissions in the other hemisphere will have much smaller transverse momenta. From eq. (C.37) we find that in this limit χ M (R) = 4 . This implies that a naive implementation of eq. (3.26) will result in p y,ℓ of the empty hemisphere to be equal to zero, and hence in a floating point exception, unless one decreases the cutoff for the integration over the ζ i . In the next section, we discuss how to deal with this problem in full generality and devise a procedure to obtain finite NP shifts down to R ′ = 0.
We observe that, for soft and collinear PT emissions, the scaling of V sc ({p}, k i ) with respect to η (ℓ) is very different for the different observables. In particular, there is a rapidity suppression for thrust, C-parameter and heavy-jet mass, but no rapidity dependence for the jet broadenings and thrust major. This is not the case for a NP emission in the presence of multiple soft and collinear PT emissions. In this case the contribution to δV NP of the collinear regions η (ℓ) → ∞ is always suppressed for the event shapes considered here, see Fig. 4. Also, the suppression is such that the rapidity integral in eq. (2.52) is convergent.
The analytic expressions for ⟨h B W ⟩ and ⟨h B T ⟩ as a function of R ′ can be found in eqs. (B.13) and (B.23) respectively. There one can see that, while ⟨h B W ⟩ is finite for all values of R ′ , ⟨h B T ⟩ has a divergence for R ′ → 0. As already anticipated, this divergence originates from the fact that, when R ′ → 0, one hemisphere contains a few emissions with k t ∼ B T Q, while emissions in the other hemisphere have much lower transverse momenta. Mathematically, this corresponds to the situation where R ′ is much smaller than all higher derivatives of R(v), which can no longer be neglected. Numerically, this creates a problem for a Monte Carlo implementation of ⟨h B T ⟩ which assumes all transverse momenta of perturbative emissions to be of the same order. For B T , the problem is solved in [39] by retaining higher derivatives of the radiator, e.g. R ′′ . While this improved evaluation can be carried out fully analytically for ⟨h B T ⟩, the same is not true for ⟨h T M ⟩ because the direction of the thrust-major axis depends on all perturbative emissions. Therefore, we need to devise a procedure to compute ⟨h V ⟩ that is suitable for all observables and gives a finite result for R ′ → 0. What we propose is to add and subtract to ⟨h V ⟩ a counterterm that displays the appropriate 1/R ′ behaviour. The counterterm is designed in such a way as to cancel the divergence of ⟨h V ⟩ for R ′ → 0 at the integrand level. It must also be simple enough to be computed fully analytically for all values of V . This procedure ensures that ⟨h V ⟩ is finite for all values of R ′ .

Counterterm for ⟨h B T ⟩
To see how a subtraction procedure might work, we consider the case of the total broadening for which we have full analytic control. In the limit R ′ → 0, one of the hemispheres will contain a single emission, which we denote k 1 , with a transverse momentum of order B T Q. In the limit R ′ → 0, there exists δ with ϵ ≪ δ ≪ 1 such that all other emissions have transverse momenta less than δB T Q, and thus Using eq. (2.51), invoking the symmetry between the two hemispheres defined by the thrust axis (H 1 ↔ H 2 symmetry) and the fact that We rescale ζ i and introduce a two-dimensional Fourier transform, noting that we must apply a cut on large values of the transverse momenta of the recoil in H 2 due to the small transverse momenta of emissions in this hemisphere. We achieve this by imposing an upper-bound |⃗ x 2 | < x max (where x max ∼ δ) and a corresponding lower-bound on the conjugate parameter This may be written in a simplified exponential form to give As emissions in H 2 do not contribute to the observable, as per eq. (4.1), we notice that there is no exponential damping factor in the ζ-integral. Using the fact that, for large b, we obtain Performing the x 2 -and b 2 -integrations and extracting the singular R ′ → 0 behaviour (setting We note that we do not control terms of O(1) in this calculation as they may arise from the interplay of unknown terms of order R ′ multiplied by the leading 1/R ′ . From the full analytical evaluation of ⟨h B T ⟩ performed in appendix B.2, the leading R ′ → 0 behaviour is presented in eq. (B.24) and we find that the uncontrolled O(1) term is in fact equal to zero in this limit. Therefore in the limit At this point we face essentially two possibilities. One is to relax this assumption and try computing ⟨h V ⟩ from scratch, a framework that we do not adopt here because of its limited scope. The other is to devise a local counterterm, which we call ⟨h c.t. V ⟩, that exhibits the same singularity for R ′ → 0 as ⟨h V ⟩. We then subtract (add) this counterterm from (to) the full shift. The subtraction of the singularity for R ′ → 0, within ⟨h V ⟩ − ⟨h c.t.
V ⟩ which will be carried out via a general Monte Carlo procedure, takes place at the integrand level and thus we obtain a finite result as R ′ → 0. The counterterm is constructed so as to have full analytic control; in particular we must be able to obtain all constant terms appearing in the shift of any observable as R ′ → 0.
Let us return to the case of the total broadening. Fully generally, for R ′ → 0 one emission, emitted in hemisphere H ℓ , has a value of transverse momentum that is much larger than all of the other emissions. In this situation B T is determined by this emission and the hemisphere that does not contain this emission, Hl, is automatically the hemisphere with the smaller broadening. Motivated by the above discussion, a good counterterm is then constructed as follows where ⃗ k t,max = max i {k ti }. In a similar way as for ⟨h B T ⟩, we can introduce the rescaled vari- with f T (0) = 1, the function σ(R ′ /2) defined in eq. (B.8) and where xl ≡ |⃗ xl|. Such counterterm has the same 1/R ′ singularity as ⟨h B T ⟩ and is constructed so that it may be computed analytically, as is performed in appendix C.1. Therefore, the combination The function χ c.t. T (R ′ ) is computed analytically in appendix C.1 and its expression can be found in eq. (C.5). We may also rewrite χ c.t.
T (R ′ ) in a way that is suitable for Monte Carlo integration, following the strategy described in appendix A. We find where, as explained in appendix A, now B T,sc ({p}, is finite for all values of R ′ . When the calculation involves the recoil due to emissions in the same hemisphere as emission k 1 , i.e. hemisphere H ℓ1 , the counterterm χ c.t. T (R ′ ) has no effect whatsoever. When the calculation involves the recoil due to emissions in the other hemisphere, the counterterm makes sure that the Monte Carlo gives a finite result even for R ′ → 0. From a numerical point of view, this regime is problematic when we have no emissions other than k 1 above the cutoff ϵ and are required to compute ln i / ∈H ℓ 1 ⃗ ζ i . However, in this limit B T,sc {p}, k 1 , . . . , k n+1 /B T → 1 and the theta-constraint of the counterterm will be trivially satisfied, thus the contributions of χ T (R ′ ) and χ c.t.
T (R ′ ) will cancel perfectly as shown in Fig. 5.
Analytic Results T (R ′ ) cancel perfectly in the limit R ′ → 0.

Numerical Results
We now describe how to compute the counterterm in eq. (4.9). If we do this naively we of course obtain a 1/R ′ divergence from eq. (4.12). This is because in devising the measure dZ[{R ′ ℓi , k i }] we have neglected all higher derivatives of the radiator. In particular, the second derivative regularises the 1/R ′ divergence. As a result the shift now behaves as 1/ √ α s and the product of such a contribution with a finite correction of order α s , which is beyond our nominal accuracy, gives a √ α s contribution. Therefore, in the improved version of χ c.t.
T we can account for all contributions up to order √ α s excluded. In order to do so, we also need to upgrade the expression for χ c.t.
T to take into account hard-collinear real and virtual corrections. The exact procedure to perform this upgrade is technically involved and is explained in appendix C. The outcome of this procedure is the following improved counterterm with where The final expression for the shift for the total broadening is then where The behaviour of the shift for the total broadening can be appreciated by considering two separate regimes, R ′ ≫ √ R ′′ and R ′ ≪ √ R ′′ . As such, we observe which gives: • R ′ ≫ √ R ′′ : this regime describes the small-B T region where many emissions populate both hemispheres and we have We see that ⟨h c.t. B T ⟩ imp. cancels almost completely against ⟨h c.t. B T ⟩, leaving a finite contribution of order α s which is beyond our control. Indeed, the precise form of the latter crucially depends on the choice of counterterm, arising from ⟨h c.t.
• R ′ ≪ √ R ′′ : this regime describes the large-B T region where the hemisphere that does not contain the emission setting the broadening is almost empty and we obtain We re-emphasise that, as anticipated, we do not control terms of order √ α s as they could arise from the product of terms of order α s beyond our nominal accuracy, multiplied by the prefactor 1/ √ α s .
The procedure outlined above to obtain the shift for B T is similar to the one originally proposed in [39]. There, a counterterm was introduced at the level of the integral transforms that were employed to compute the shift. Our procedure is designed to handle a generic observable, in a way that is amenable to an efficient numerical implementation. We plot our results for ⟨h B T ⟩ in Fig. 6 and compare them both to ⟨h B T ⟩ and the analytic expressions in [39] (which we refer to as ⟨h B T ⟩ DMS ). We first remark that in the region of large B T , corresponding to small R ′ , we have specifically checked that ⟨h B T ⟩ and ⟨h B T ⟩ DMS agree up to corrections of order √ α s , which are beyond our nominal accuracy. In the region of small B T , corresponding to large R ′ , we note that our result differs both from [39] and the analytic behaviour of ⟨h B T ⟩. This deviation is due to the residual O(α s ) term in eq (4.21). This occurs because our method relies on cancelling the 1/R ′ singularity locally at the integrand level, which unavoidably leaves a residual O(α s ) contribution in our final shift. It would be desirable to devise a counterterm which automatically switches off these spurious contributions for B T → 0, but we leave this for future work.
At the phenomenological level, the results that we present are valid up to values of B T that are not too small, which contain the range in which simultaneous fits of α s and α 0 take place and is thus appropriate for the study performed in this paper.

Counterterm for ⟨h T M ⟩
Once we have validated our procedure we are in a position to compute the shift for the thrust major. Fully generally, for R ′ → 0, one emission, emitted in hemisphere H ℓ , has a value of transverse momentum that is much larger than all of the other emissions. In this situation T M is determined by this emission only, and it is this emission that sets the thrust-major axis. This reduces enormously the complexity of the calculation of the shift for the thrust major and makes it amenable to an analytic calculation.
A suitable counterterm that may be computed analytically is constructed as follows where again ⃗ k t,max = max i {k ti } and the y-direction is along ⃗ k t,max . In a similar way as for ⟨h c.t. B T ⟩, we can introduce the rescaled variables x ℓ = 2p y,ℓ /(T M Q) and the two-dimensional vectors ⃗ ζ i ≡ ζ i (cos ϕ i , sin ϕ i ) and obtain with f M (0) = 1, the function ρ 1 (R ′ /2) defined in eq. (C.31) and where ϕ i,max denotes the angle between ⃗ k t,i and ⃗ k t,max . Such counterterm has the same 1/R ′ singularity as ⟨h T M ⟩ and is constructed so that it may be computed analytically, as is performed in appendix C.3. Therefore, the combination ⟨h T M ⟩ − ⟨h c.t. T M ⟩ is free of singularities for all values of R ′ with the cancellation of singularities occurring locally, i.e. for each Monte Carlo configuration used to compute ⟨h T M ⟩ − ⟨h c.t.
T M ⟩. What we obtain is The function χ c.t. M (R ′ ) is computed analytically in appendix C.3 and its expression can be found in eq. (C.34). We may also rewrite χ c.t.
M (R ′ ) in a way that is suitable for Monte Carlo integration, following the strategy described in appendix A. We find where, as explained in appendix A, now T M,sc ({p}, k 1 ) = T M . We observe that χ M (R ′ ) − χ c.t. M (R ′ ) is finite for all values of R ′ . When the calculation involves the recoil due to emissions in the same hemisphere as emission k 1 , i.e. hemisphere H ℓ1 , the counterterm χ c.t.
M (R ′ ) has no effect whatsoever. When the calculation involves the recoil due to emissions in the other hemisphere, the counterterm makes sure that the Monte Carlo gives a finite result even for R ′ → 0. From a numerical point of view, without a counterterm this regime would be problematic when we have no emissions other than k 1 above the cutoff ϵ and are required to compute ln i / ∈H ℓ 1 ζ i sin ϕ i . However, in this limit T M,sc {p}, k 1 , . . . , k n+1 /T M → 1, ϕ i → ϕ i,max and the theta-constraint of the counterterm will be trivially satisfied, thus the contributions of χ M (R ′ ) and χ c.t.
M (R ′ ) will cancel perfectly as shown in Fig. 7. Performing a similar upgrade as for the total broadening we obtain the needed, improved expression for ⟨h c.t.

Numerical Results
T M ⟩ imp. , computed in appendix C.3 and repeated here with where The final expression of the shift for the thrust major is given by where The behaviour of the shift for the thrust major can be appreciated by considering two separate regimes using eq. (4.20): • R ′ ≫ √ R ′′ : this regime describes the small-T M region, where many emissions populate both hemispheres and we have We see that ⟨h c.t. T M ⟩ imp. cancels almost completely against ⟨h c.t. T M ⟩, leaving a finite contribution of order α s which is beyond our control. Indeed, the precise form of the latter crucially depends on the choice of counterterm, arising from ⟨h c.t.
T M ⟩ imp. . • R ′ ≪ √ R ′′ : this regime describes the large-T M region where the hemisphere that does not contain the emission that sets the thrust-major axis is almost empty and we obtain Note that we do not control terms of order √ α s as they could arise from the product of unknown terms of order α s multiplied by the prefactor 1/ √ α s .
The function ⟨h T M ⟩ is plotted in Fig. 8. As for the total broadening we notice that the 1/R ′ divergence has been replaced with a constant value for large values of T M , while at very small values of T M discrepancies between ⟨h T M ⟩ and ⟨h T M ⟩ will be of order α s and beyond our control. We conclude by commenting on the fact that the exact behaviour of the shift at small values of T M (the black dots in Fig. 8) is perfectly under control. The fact that ⟨h T M ⟩ does not tend to ⟨h T M ⟩ for T M → 0 is due to a residual O(α s ) term, as was explained in the case for the total broadening.

Mean values
With the formalism we have developed so far we can also compute the non-perturbative corrections to the mean values of event shapes. The mean value of an event shape is defined as , from eq. (2.24), into the above equation, we obtain with δΣ NP given in eq. (2.32). The perturbative element of the mean value, ⟨v⟩ PT , can be computed as an expansion in powers of α s . The non-perturbative element, ⟨v⟩ NP , can be simplified by performing an integration by parts We can now use the fact that v max is the same for both perturbative and non-perturbative distributions to set the first term on the right-hand side of the above equality to zero. In the second term we can substitute the expression for δΣ NP (v) in eq. (2.32) and approximate the resulting expression in terms of δV NP defined in eq. (2.34) Using the assumption in eq. (2.39), inserting the expression for δV NP ({p}, k, {k i }) in eq. (2.40) and assuming that the ultra-soft emission k is accompanied by perturbative emissions that are soft, collinear and widely separated in angle, we obtain (5.6) We can now compute the NP corrections to the mean values of the event-shape variables considered in section 3.

(5.7)
Jet broadenings and thrust major. Starting from eq. (5.6), the NP corrections to the mean values of the jet broadenings and thrust major (unlike the shift for the thrust major distribution) may be computed fully analytically. Since no novel numerical procedure is required, we leave the details of the calculation to appendix D. Here we report the final results: We note that we do not control terms of order √ α s and to this accuracy eqs. (5.8) and (5.9) agree with the corresponding results in [39]. Comparing results explicitly, the sole difference is the replacement of the α s (Q) term in the denominator of eqs. (5.8) and (5.9) with α s,CMW (Q), withQ = Qe −3/4 , in [39]. This is highlighted in [39] and results in a difference of order √ α s which is not formally under control.

Mean Values
We take the non-perturbative corrections to the mean, ⟨v⟩ NP , computed in section 5 for the various event-shape observables of interest and add these to ⟨v⟩ PT , the perturbative element of the mean value. This provides the full mean value, ⟨v⟩ = ⟨v⟩ PT + ⟨v⟩ NP , for which we perform simultaneous fits for α s (m Z ) and α 0 (2 GeV) to data. The PT element of the mean value is given by with A V and B V the O (α s ) and O α 2 s fixed-order perturbative coefficients. A V and B V are known analytically for C [50] and have otherwise been determined using the program EVENT2 [51]. The values used have been set out in Table 1.  Table 1: Fixed-order perturbative coefficients of the event-shape observables taken from [50] for C and otherwise determined using EVENT2 [51]. These are used in performing the simultaneous fits for α s and α 0 .
The latest ALEPH QCD publication [35] provides references to the data [33,36,[52][53][54][55][56][57][58][59][60][61][62] used to perform the simultaneous fits for the means of 1 − T , C, ρ H , B W and B T . Where possible we have identified the exact same data, however it has not been possible to obtain 13 of the data points used for the fits for 1 − T and 22 of the data points used for the fits for ρ H . For T M , as this observable was not considered in [35], we have identified all available data for T M that we were able to find in the literature [29,30,36,59,60,[63][64][65]. We note that there is considerably less data available for T M than was used in the fits for the other event-shape observables in [35]. With this data we have performed the simultaneous fits for α s and α 0 with the results set out in Table 2. Using the results from Table 2 we have plotted the 95% CL contours for the event-shape means in Fig. 9 and the energy dependence of the mean values of the event-shape observables, comparing LO perturbative, NLO perturbative and NLO+NP predictions with data in Fig. 10.
For 1−T (despite the 13 data points that we have been unable to obtain), C and B T (despite the O( √ α s ) difference between our result for ⟨B T ⟩ and that of [39]) our simultaneous fitted values for α s and α 0 lie within the 95% confidence level contours deduced from the fitted values and associated errors calculated by the ALEPH collaboration [35]. For ρ H and B W however, our simultaneous fitted values for α s and α 0 do not lie within the 95% confidence level contours deduced from [35].  Table 2: Results for the simultaneous fits for α s and α 0 to experimental data for the mean values of event-shape observables.
For B W , this is due to the O( √ α s ) difference between our result for ⟨B W ⟩ NP and that of [39].
For ρ H , this is due to the 22 data points (at low centre-of-mass energies) that we have been unable to find in public repositories. This has resulted in a noticeably large 95% confidence level contour in Fig. 9 which is particularly driven by the uncertainty in α 0 . We note that the α s and α 0 values for ⟨ρ H ⟩ computed in [35] lie within our large 95% confidence level contour. It is our belief that including the 22 data points would both reduce the size of our contour and also produce simultaneous fitted values for α s and α 0 that lie within the deduced ALEPH 95% confidence level contours. Values of α s between 0.1117 and 0.1230 are compatible with the Particle Data Group's world average value of 0.118 [1]. The values of α 0 lie close to 0.5 apart from for ρ H which is noticeably higher. The large α 0 value for ρ H is in part due to the low centre-of-mass energy data points that have not been obtained (cf. [35] which included these data points and found α 0 = 0.627). This may also be due to the effect of hadron masses, as discussed in [66], which we neglect in this preliminary study.
In Fig. 9, the 95% confidence level contours for all 6 observables show a strong negative correlation between the fitted values for α s and α 0 . We note that the contours for 1 − T , C, B W and B T  Table 4 and compared with experimental data.
lie close together with the contour for T M lying a little below that for B W . We stress that far fewer data points were identified for T M than for the other observables and that, for all event shapes, very few measurements were available at low centre-of-mass energies.

Distributions
For the event-shape observables considered in section 3, we performed an NLL resummation using the CAESAR program [46]. These NLL resummed distributions were then matched with fixed order distributions at NLO, obtained from EVENT2, using the log-R scheme [67]. To the resummed, matched distributions (NLL+NLO) we then apply the NP shifts that we have computed in sections 3 and 4.
The latest ALEPH QCD publication [35] provides experimental data covering centre-of-mass energies between 91.2 GeV and 209 GeV, as well as the fit ranges that were employed for the simultaneous fits for α s (m Z ) and α 0 (2 GeV) for each of the event-shape observables. As a test of our method we have attempted to repeat these simultaneous fits. Where possible we have used the same fit ranges, however for the C-parameter and jet broadenings it was necessary to slightly reduce the upper bound of the fit ranges. It was noted in [35] that the fit ranges for all observables have been selected in the central region of three-jet production.
Using our method, we obtained simultaneous fit values for α s and α 0 for 1 − T , C and B T that are consistent with the one-sigma confidence level contours calculated by the ALEPH collaboration [35]. For ρ H and B W it was not possible to obtain reliable and consistent fit values for α s and α 0 . ρ H and B W exhibit a known property whereby the resummed, matched distributions (NLL+NLO) must be squeezed to smaller observable values to ensure that the non-perturbative shift is positive over the whole fit region and thus enabling an accurate fit to experimental data. This requires a significant reduction in the value of α s from the Particle Data Group's world average value of 0.118 [1]. In addition, fit ranges in the central region of three-jet production cover a region where this squeeze requirement is particularly pronounced for ρ H and B W .
While the fit ranges in [35] have been selected in the central region of three-jet production, the NP shifts are computed in the two-jet region. As a result we have performed a new simultaneous fit for α s and α 0 where we have extended the lower bound of the fit ranges in [35] to the two-jet peak. This criterion has been used to determine appropriate fit ranges for the thrust major. The fit ranges that we have used in this study are set out in Table 3. Using the experimental data in [35] and the fit ranges detailed in Table 3 we have performed simultaneous fits for α s and α 0 with the result set out in Table 4.  Table 3: Fit ranges for the simultaneous fits for α s and α 0 to ALEPH data from 91.2 GeV to 206 GeV for the distributions of event-shape observables.
Using the results from Table 4 we have plotted the 95% CL contours for the event-shape distributions in Fig. 11 and the comparison between the resummed, matched distributions (NLO+NLL) with α s = 0.118, our fitted distributions (NLO+NLL and NLO+NLL+1/Q) using the results in Table 4 and experimental data at a centre-of-mass energy of 91.2 GeV in Fig. 12.   Figure 11: 95% confidence level contours for the fitted values of α s and α 0 to experimental data for the distributions of event-shape observables.

Event Shape
As for the 95% confidence level contours for the mean values, the 95% confidence level contours for the distributions of each of the observables in Fig. 11 show a strong negative correlation between the fitted values of α s and α 0 . However, we notice that the 95% confidence level contours for the distributions lie noticeably further apart from one another than was the case for the mean values in Fig. 9.
In Fig. 12 we observe that T M exhibits the same squeeze requirement as for ρ H and B W . It is clear therefore that an extension of our calculation to the three-jet region is particularly needed for these observables. Nevertheless, we have still performed the simultaneous fits with these observables. We find that this squeeze requirement leads to noticeably smaller fitted values of α s for ρ H , B W and T M than for T , C and B T . In particular, for T M we observe a fitted value of α s that is comparable with that for B W and a fitted value of α 0 located between that of B T and the grouping of C and  Table 4 and compared with experimental data. The grey shading indicates the regions excluded from the fit range.

Conclusions
In this paper we have presented a general method to compute the leading non-perturbative corrections to event-shape distributions in the two-jet region. The power of our method relies on a numerical algorithm suitable to treat any observable, including those leading to logarithmic enhancement in the linear power correction. Indeed, prior to our work, only few examples of such observables were studied, specifically those who could be handled analytically. The crucial point is that leading power-suppressed corrections can be modelled in terms of the emission of an ultra-soft gluon, accompanied by an ensemble of perturbative soft and collinear emissions. The only incalculable quantity is a moment, known as α 0 , of the soft effective coupling which determines ultra-soft emission probability. The integration over the remaining kinematic variables, i.e. rapidity and azimuth, is assumed to follow the perturbative soft emission matrix element squared. Therefore, the ultra-soft gluon can be considered as a "special" emission, which together with an arbitrary number of soft and collinear emissions, can be simulated numerically with the ARES method. In fact, similar contributions appear when computing NNLL corrections to two-jet event-shape distributions. First, we have validated the method to reproduce leading hadronisation corrections to all known event-shape distributions. Then, we have been able for the first time to compute the hadronisation correction to the distribution of the thrust major, which does not allow an analytic treatment, but for which data has existed for a long time. In order to do so, we had to tackle the problem of unphysical divergences of hadronisation corrections occurring for large values of recoil-sensitive event shapes, e.g. the total broadening and the thrust major itself. This is done by performing a subtraction procedure, so that all numerical integrations are finite, and unphysical divergences are treated fully analytically in a general way.
Finally, we have performed new simultaneous fits of the strong coupling α s and α 0 using a selection of experimental data obtained by the ALEPH collaboration. We have been able to obtain consistent results for all event shapes for which hadronisation corrections were known. Moreover, we have performed a similar fit for the thrust major distributions and mean values. The resulting values of α s and α 0 are in the same ballpark as for other event shapes.
The method we have devised is just our own first attempt to improve the phenomenology of event-shape distributions and means at LEP. A further step is to extend our calculation to the three-jet region, along the lines of what was done in [40,41]. Such extension is particularly needed for the heavy-jet mass, wide-jet broadening and thrust major. Another direction of improvement is more sophisticated procedures to subtract unphysical divergences for recoil-sensitive event shapes. Last, we remark that although here for simplicity we have restricted our method to selected event shapes, our procedure could be extended to the two-jet rate, which is particularly important for precise α s determinations. This requires including information on the rapidity of subsequent soft and collinear emissions, in a similar way as is done for NLL [47] and NNLL [6] resummations and of course including the interplay between PT and NP emissions. We remark that for jet rates, the treatment of multiple ultra-soft emissions does not lead to the Milan factor and needs to be computed along the lines of [68]. We leave both issues for future work.

A Monte Carlo determination of non-perturbative shifts
This appendix describes the Monte Carlo procedure that we adopt to compute the various nonperturbative shifts. As an illustrative example, let us first consider the denominator of eq. (2.52) Here, the contribution with zero emissions does not satisfy the observable constraint δ(v−V sc ({p}, {k i })). We therefore select emission k 1 , such that V sc ({p}, k 1 ) is the largest of all V sc ({p}, k i ), and neglect all emissions such that V sc ({p}, k i ) < ϵV sc ({p}, k 1 ). This gives We introduce the rescaling ζ i → ζ 1 ζ i and rescale the corresponding momenta k i such that V sc ({p}, k i ) → ζ 1 V sc ({p}, k i ). As V is rIRC safe we have the rescaling We can therefore write eq. (A.2) as The ζ 1 -integration can now be performed using the delta-function constraint to obtain where F(R ′ ) is the NLL multiple-emission function, whose expression in terms of {k i } reads The above expression can be evaluated numerically using a Monte Carlo procedure. One can apply the same approach to the numerator of eq. (2.52) and we find where now V sc ({p}, k 1 ) = v.

B Analytical determination of non-perturbative shifts
This appendix presents analytical computations of the non-perturbative shifts for the wide-jet and total broadenings within our framework. For the thrust major, the analytical computation of the non-perturbative shift in the limit R ′ → 0 is also presented.

B.1 Wide-Jet Broadening
For the wide-jet broadening we wish to analytically compute ⟨h B W ⟩, thus enabling us to determine ⟨δB W, NP ⟩ via eq. (2.41). We start with eqs. (3.19) and (3.20) which we repeat here for clarity Using eq. (2.51), invoking the symmetry between the two hemispheres defined by the thrust axis (H 1 ↔ H 2 symmetry) and performing the rescaling where x ℓ ≡ |⃗ x ℓ |. Introducing the Fourier and Mellin transforms of the kinematic and observable constraints respectively, this may be written as We make use of the following integral 4 We highlight that R ′ ℓ i is the radiator corresponding to hemisphere H ℓ i and that R ′ = R ′ 1 + R ′ 2 (where we will eventually set R ′ 1 = R ′ 2 = R ′ /2 due to the H 1 ↔ H 2 symmetry). and find It is convenient to redefine µ → µ + ν, rescale b 1 → µb 1 , b 2 → νb 2 , x 1 → x 1 /µ and x 2 → x 2 /ν, after which we perform the x 1 -and x 2 -integrals directly and find where y i ≡ 1 + b 2 i . To evaluate the y i -integrals we make use of the following auxiliary functions taken from [39] For the wide-jet broadening, we know that with the polygamma function ψ(a) = d ln Γ(a)/da. In agreement with [39] we obtain It is important to inspect the leading R ′ → 0 behaviour of ⟨h B W ⟩. To do so we require the following limits such that

B.2 Total Broadening
The analytical determination of the non-perturbative shift for the total broadening follows a similar process to that for the wide-jet broadening in appendix B.1. As before, we start with eqs. (3.19) and (3.20) which we repeat here for clarity We follow similar steps as for the wide-jet broadening in appendix B.1. Using eq. (2.51), performing the rescaling ζ i → 2ζ i , introducing the Fourier and Mellin transforms of the kinematic and observable constraints respectively and utilising the H 1 ↔ H 2 symmetry, we write χ T (R ′ ) as where x ℓ ≡ |⃗ x ℓ |. Using eq. (B.5) we find We rescale b i → νb i and x i → x i /ν, after which we perform the x 1 -and x 2 -integrals directly and find where y i ≡ 1 + b 2 i . We note that compared to the equivalent expression for B W , eq. (B.7), the key difference is that for B T there is only one Laplace variable whereas for B W there are two. We will see below that this will have the effect of leaving a residual 1/R ′ singularity which is not cancelled. For the total broadening we know that Therefore, using eqs. (B.8), (B.9) and (B.10) and setting R ′ 1 = R ′ 2 = R ′ /2 we find 22) and in agreement with [39] we obtain Using eq. (B.14) we deduce the leading R ′ → 0 behaviour of ⟨h B T ⟩

B.3 Thrust Major
As explained in section 4, the shift for the thrust major cannot be obtained in closed form. Nevertheless we can analytically determine its behaviour in the limit R ′ → 0. We start with eqs. (3.27) and (3.28) which we repeat here for clarity In the limit R ′ → 0, one of the hemispheres will contain a single emission, which we denote k 1 , with a transverse momentum of order T M Q/2 and which sets the thrust-major axis. All other emissions have transverse momenta less than δ T M Q 2 , where ϵ ≪ δ ≪ 1, and thus Using eq. (2.51), the H 1 ↔ H 2 symmetry and the fact that R ′ = R ′ 1 + R ′ 2 , we write χ M (R ′ ) in the limit R ′ → 0 as We rescale ζ i and introduce a one-dimensional Fourier transform, noting that we must apply a cut on large values of the transverse momenta of the recoil in H 2 due to the small transverse momenta of emissions in this hemisphere. We achieve this by imposing an upper-bound |x 2 | < x max (where x max ∼ δ) and a corresponding lower-bound on the conjugate parameter |b 2 | > b min (where b min ∼ 1/δ). We obtain This may be written in a simplified exponential form to give As emissions in H 2 do not contribute to the observable we notice that there is no damping factor in the ζ-integral. From eq. (4.5) we obtain Performing the x 2 -and b 2 -integrations and extracting the singular R ′ → 0 behaviour, setting We note that we do not control terms of O(1) in this calculation as they may arise from the interplay of unknown terms of order R ′ multiplied by the leading 1/R ′ .

C Analytical determination of the counterterms
This appendix presents analytical computations of the counterterms for the total broadening and thrust major.

C.1 Counterterm for the Total Broadening
In section 4 we introduced the counterterm for the total broadening which we wish to determine analytically. We start with eq. (4.12) which we repeat here where Hl denotes the hemisphere that does not contain the emission with the largest transverse momentum and xl ≡ |⃗ xl|. We use eq. (2.51), the H 1 ↔ H 2 symmetry and the fact that R ′ = R ′ 1 +R ′ 2 . Introducing the Fourier and Mellin transforms of the kinematic and observable constraints respectively, we may write χ c.t.
T (R ′ ) in a simplified exponential form to give By construction, as discussed in section 4, χ c.t. T (R ′ ) must reproduce the same 1/R ′ behaviour of ⟨h B T ⟩ such that it can act as a suitable local counterterm in the numerical procedure. To see this, we take the NLL approximation of the Sudakov radiator, i.e. R ′ 2 (ζv) ≃ R ′ 2 (v), in eq. (C.2) and using the H 1 ↔ H 2 symmetry we find (C.3) We rescale b 2 → νb 2 and x 2 → x 2 /ν, after which we perform the x 2 -integrals directly and find where f T (R ′ ) is defined in eq. (4.11). Using eq. (B.14) we deduce the leading R ′ → 0 behaviour of This is the same 1/R ′ divergence as for the leading behaviour of χ T (R ′ ) in the limit R ′ → 0 (in eq. (B.24)) and therefore χ c.t.
is regular (and in fact tends to 0) in the limit R ′ → 0. As discussed in section 4, in the limit R ′ → 0 one can no longer neglect the higher derivatives of the radiator in eq. (C.2) which, using the H 1 ↔ H 2 symmetry, now reads This can be evaluated exactly to give where R (n) ≡ R (n) (v) denotes the n-th logarithmic derivative of the radiator. Rescaling b 2 → νb 2 and x 2 → x 2 /ν and introducing a change of variables y = 1 + b 2 2 , the x 2 -integrals can be performed exactly to give In eq. (C.9), the accuracy sought is to have full control over all constant terms in the limit R ′ → 0. Therefore, in all the terms on the first line of eq. (C.9) that converge as y → ∞ we can neglect contributions that involve two or higher derivatives of the radiator. This leaves the dy/y term, which originally triggered the R ′ → 0 divergence, where higher derivatives of the radiator are now required to be taken into account. As will become apparent below, this requires the exponent in (C.9) to be expanded to O(R (3) ) as follows where we have introduced the change of variable, z = ln[(1 + y)/2]. We notice that it is enough to retain R ′′ 2 in the exponent to guarantee the convergence of the z-integration as z → ∞. Also, in the last line we have neglected any terms beyond our accuracy. Setting R where the functional form ofH T is given bỹ (C.12) The above discussion highlights the fact that terms which would be subleading from the point of view of logarithmic resummations, such as those containing R ′′ and R (3) , do contribute to the shift at an accuracy that is within the control we claim. This suggests the need to include in χ c.t. imp T (R (n) ) further sub-leading effects such as those induced by hard-collinear radiation. This is what we do in the next section.

C.2 Improved Sudakov and Hard Collinear Element
In the limit of R ′ → 0, our goal is to have control over terms up to and including constants. To do so we need to consider a perturbative configuration with one hard-collinear emission on top of an ensemble of soft-collinear emissions. On the perturbative level, a single hard-collinear emission gives rise to an NNLL contribution which is of the same order as R ′′ .
To include the effect of a hard-collinear emission, say to leg ℓ, we start with the emission probability where the vector ⃗ κ ≡ ⃗ κ (ℓ) is defined in eq. (2.17). In the collinear limit, we also have that z = z (ℓ) The change in the observable due to the hard-collinear emission leads to the following contribution The above expression is manifestly finite in the singular limits ζ → 0 and z → 0. To understand the contribution of the above term in the limit R ′ → 0 we merely need to Taylor expand the first step function around small recoil, xl → 0. Thanks to the fact that the difference between the step functions vanishes at xl = 0, the Taylor expansion at leading order is linear in xl. This leads us to conclude that eq. (C.21), in the limit xl → 0, yields an O(α s ) contribution to the shift which is beyond our accuracy.
Notice that in eq. (C.21) the double-counting is properly subtracted, albeit that the subtraction is only needed for the 2/z part of the splitting function. Therefore, for the regular portion of the splitting function we need to consider the following contribution where the term entering with negative weight is nothing but the hard-collinear virtual corrections. It should be understood that the soft-collinear phase space measure is upgraded by including higher derivatives of the radiator, as was done in the previous subsection. Once again, the above expression is manifestly finite in the singular limit ζ → 0. As we mentioned before, our goal is to have full control over any constants in the limit R ′ → 0. Therefore, we end up with the following expression (setting R which easily evaluates to π s e s 2 Erfc(s) .
where H T is given by

C.3 Counterterm for the Thrust Major
In section 4 we introduced the counterterm for the thrust major which we wish to determine analytically. We start with eq. (4.26) where ϕ i,max denotes the angle between ⃗ k t,i and ⃗ k t,max . We use eq. (2.51), the H 1 ↔ H 2 symmetry and the fact that Introducing the Fourier and Mellin transforms of the kinematic and observable constraints, respectively, we may write χ c.t.
M (R ′ ) in a simplified exponential form to give As was the case for χ c.t. T (R ′ ), by construction χ c.t. M (R ′ ) must exactly reproduce the leading 1/R ′ behaviour of ⟨h T M ⟩ such that it can act as a local counterterm in the numerical procedure. To demonstrate this, we take the NLL approximation of the Sudakov radiator, i.e. R ′ 2 (ζv) ≃ R ′ 2 (v), in eq. (C.28) and using the H 1 ↔ H 2 symmetry we find (C.29) We rescale b 2 → νb 2 and x 2 → x 2 /ν, after which we perform the x 2 -integral directly and find To evaluate the y-integrals we make use of the following (C.34) To determine the leading R ′ → 0 behaviour of χ c.t.
M (R ′ ) we observe the following limits Therefore we find We observe that this is the same 1/R ′ divergence as for the leading behaviour of χ M (R ′ ) in the limit R ′ → 0 in eq. (B.32). We note from section 4.2 that in the limit ) → 0 and we may therefore deduce the following As discussed in section 4, in the limit R ′ → 0 one can no longer neglect the higher derivatives of the radiator in eq. (C.28) which, using the H 1 ↔ H 2 symmetry, now reads This may be evaluated exactly to give where R (n) denotes the n-th logarithmic derivative of the radiator. Rescaling b 2 → νb 2 and x 2 → x 2 /ν and introducing a change of variables y = 1 + b 2 2 , the x 2 -integral can be performed exactly to give (C.40) -49 -As was the case for χ c.t. imp T (R (n) ), the accuracy sought in eq. (C.40) is to have full control over all constant terms in the limit R ′ → 0. Therefore, in all the terms on the first line of eq. (C.40) that converge as y → ∞ we can neglect contributions containing two or higher derivatives of the radiator. This leaves the dy/y term, which originally triggered the R ′ → 0 divergence, where higher derivatives of the radiator are now required to be taken into account. As will become apparent below, this requires the exponent in (C.40) to be expanded to O(R (3) ), thus which we write as where we have introduced the change of variable, z = ln[(1 + y)/2]. We notice that it is enough to retain R ′′ 2 in the exponent to guarantee convergence of the z-integration as z → ∞. In addition, in the final line we have neglected any terms beyond our accuracy. Setting R where the functional form ofH M is given bỹ As forH T (R ′ , R ′′ , R (3) ), we must also consider the contribution of one hard-collinear emission on top of an ensemble of soft-collinear emissions. The calculation of this contribution follows the same steps as is set-out in appendix C.2 and we find analogously that π s e s 2 Erfc(s) .
where H M is given by

D Analytical determination of NP corrections to the means
This appendix presents analytical computations of the non-perturbative corrections to the mean values for the jet broadenings and the thrust major.

D.1 Single-Jet Broadening
Without loss of generality, for the single-jet broadening we shall consider B 1 . We start with eq. (5.6) (D.1) Following the approach of section 3, from eq. (3.15) Therefore we obtain We note the presence of B 1,max in the denominator of the observable constraint. As B 1,max corresponds to R ′ = 0 we may write our above equation as It is possible to solve eq. (D.4) by evaluating the square brackets first, for general R ′ , before taking the limit R ′ → 0. However, this approach will only be amenable for observables that allow an analytic computation. Instead, we may consider the limit R ′ → 0 at first. In the limit R ′ → 0, our hemisphere of interest will contain a single emission, which we denote k 1 , with a transverse momentum of order B 1 Q. We consider the observable when all other emissions, k i , have transverse momenta less than δB 1 Q where ϵ ≪ δ ≪ 1 thus We use eqs. (2.49) and (2.51) and, in a similar way as for ⟨h B T ⟩, introduce the rescaled variables ⃗ x 1 = ⃗ p t,1 /(B 1 Q) and the two-dimensional vectors ⃗ ζ i ≡ ζ i (cos ϕ i , sin ϕ i ) and find We note that the factor of e −R(v) in eq. (2.49) has been dropped as we are in the limit R ′ → 0, where B 1 → B 1,max and thus e −R(vmax) → 1. We also drop the ln(1/B 1 ) term as this will vanish in the limit R ′ → 0 as B 1,max ≃ 1 (with any resultant discrepancies of order √ α s and thus outside of our desired accuracy).
To identify the precise form of any 1/R ′ behaviour we take the NLL approximation of the Sudakov radiator, i.e. R ′ 1 (ζ 1 v) ≃ R ′ 1 (v), in eq. (D.6). Performing the ζ 1 -integration gives As discussed in section 4, in the limit R ′ → 0 one can no longer neglect the higher derivatives of the radiator in eq. (D.6). The accuracy that we require is to have full control over all constant terms in the limit R ′ → 0. We therefore return to eq. (D.6) and, where appropriate, perform the logarithmic expansion of R ′ 1 (ζ 1 v), expanding to O(R (3) ) to achieve the required accuracy. We find Performing the ζ 1 -integration by parts and setting R We must also include the contribution of one hard-collinear emission on top of an ensemble of soft-collinear emissions, the calculation of which is set-out in appendix C.2. We therefore obtain (D. 10) In the limit R ′ → 0, R ′ ≪ √ R ′′ and from eq. (4.20) we find, in agreement with [39] up to terms of order √ α s , that

D.2 Wide-Jet Broadening
The analytic determination of the non-perturbative correction to the mean value for the wide-jet broadening follows a similar approach to that for the single-jet broadening in appendix D.1, but with one important difference that we shall see below. As before we start with eq. (5.6) (D.12) Following the approach of section 3, from eq. (3.15) (D.13) We find (D.14) In the limit R ′ → 0, our wide-hemisphere will contain a single emission, which we denote k 1 , with a transverse momentum of order B W Q. We consider the observable when all other emissions, k i , have transverse momenta less than δB W Q where ϵ ≪ δ ≪ 1 thus We use eqs. (2.49) and (2.51) and introduce the rescaled variables ⃗ x ℓ = ⃗ p t,ℓ /(B W Q) and the twodimensional vectors ⃗ ζ i ≡ ζ i (cos ϕ i , sin ϕ i ). We perform the rescaling ζ i → ζ 1 ζ i and, as R ′ = R ′ 1 +R ′ 2 , we find 16) We note that the factor of e −R(v) in eq. (2.49) has been dropped as we are in the limit R ′ → 0 where B W → B W,max and thus e −R(vmax) → 1. We also drop the ln(1/B W ) term as this will vanish in the limit R ′ → 0 as B W,max ≃ 1 as before.
We note a factorisation of the integrals over emissions in H 1 and H 2 . The second set of square brackets, containing the integrals over H 2 , evaluates to 1 + O (R ′ 2 ) 2 in the limit R ′ → 0 and may therefore be dropped. For the first set of square brackets, containing the integrals over H 1 , we note that the rescaling of the emissions in the non-wide hemisphere produces an exponent of R ′ (ζ 1 v) rather than R ′ 1 (ζ 1 v) as was the case for ⟨B 1 ⟩ NP . To identify the precise form of any 1/R ′ behaviour we take the NLL approximation of the Sudakov radiator, i.e. R ′ ℓ (ζ 1 v) ≃ R ′ ℓ (v), in eq. (D.16). Performing the ζ 1 -integration gives As discussed in section 4, in the limit R ′ → 0 one can no longer neglect the higher derivatives of the radiator in eq. (D. 16). The accuracy that we require is to have full control over all constant terms in the limit R ′ → 0. We therefore return to eq. (D. 16) and, where appropriate, perform the logarithmic expansion of R ′ ℓ (ζ 1 v), expanding to O(R (3) ) to achieve the required accuracy. Using the H 1 ↔ H 2 symmetry and setting R = R (n) /2, we find where we highlight that crucially here s = R ′ / √ 2R ′′ . We must also include the contribution of one hard-collinear emission on top of an ensemble of soft-collinear emissions, the calculation of which is set-out in appendix C.2. We therefore obtain √ π s R ′ e s 2 Erfc(s) + 1 3 √ π s e s 2 Erfc(s) − 1 .
(D. 19) In the limit R ′ → 0, R ′ ≪ √ R ′′ and from eq. (4.20) we find, in agreement with [39] up to terms of order √ α s , that We note that eq. (D.20) takes a form similar to that in eq. (D.11) but with C F → 2C F due to the rescaling of the emissions in the non-wide hemisphere producing an exponent of R ′ (ζ 1 v) rather than R ′ 1 (ζ 1 v).

D.3 Total Broadening
The analytic determination of the non-perturbative correction to the mean value for the total broadening follows a similar approach to that for the single-jet broadening in appendix D.1. As before we start with eq. In the limit R ′ → 0, one hemisphere, say H 1 , will contain a single emission, which we denote k 1 , with a transverse momentum of order B T Q. We consider the observable when all other emissions, k i , have transverse momenta less than δB T Q where ϵ ≪ δ ≪ 1 thus We use eqs. (2.49) and (2.51) and introduce the rescaled variables ⃗ x ℓ = ⃗ p t,ℓ /(B T Q) and the twodimensional vectors ⃗ ζ i ≡ ζ i (cos ϕ i , sin ϕ i ). We perform the rescaling ζ i → ζ 1 ζ i and, as R ′ = R ′ 1 +R ′ 2 , we find (D. 24) We note that the factor of e −R(v) in eq. (2.49) has been dropped as we are in the limit R ′ → 0 where B T → B T,max and thus e −R(vmax) → 1. We also drop the ln(1/B T ) term as this will vanish in the limit R ′ → 0 as B T,max ≃ 1 as before.
We consider the three terms in the final pair of brackets: • For the ln(1/ζ 1 ) and η (B) 0 terms there is no dependence on emissions in H 2 and therefore the integrals over H 2 will trivially evaluate to 1.
• For the ln(1/|⃗ x 2 |) term, we of course have a dependence on emissions in H 2 but we note that the integrals over H 1 and H 2 will still factorise. We set the upper-bound of the ζ iintegration to δ (to reflect the imposed restriction on transverse momenta of emissions in H 2 ) and recognise the integrals over H 2 as precisely that in eq. (4.2). We recall from eq. (4.7) that this will give 1/(2R ′ 2 ) + O(R ′ ).
For the integrals over H 1 we note that, as was the case for ⟨B W ⟩ NP , rescaling the emissions in H 2 produces an exponent of R ′ (ζ 1 v) rather than R ′ 1 (ζ 1 v) (as was present for ⟨B 1 ⟩ NP ). To identify the precise form of any 1/R ′ behaviour we take the NLL approximation of the Sudakov radiator, i.e. R ′ ℓ (ζ 1 v) ≃ R ′ ℓ (v), in eq. (D.24). Performing the ζ 1 -integration we find which simplifies to give As discussed in section 4, in the limit R ′ → 0 one can no longer neglect the higher derivatives of the radiator in eq. (D.24). The accuracy that we require is to have full control over all constant terms in the limit R ′ → 0. We therefore return to eq. (D.24) which, using the H 1 ↔ H 2 symmetry, we write in a simplified manner as To achieve the required accuracy we will perform the logarithmic expansion of R ′ ℓ (ζ 1 v), expanding to O(R (3) ). To enable convenient analytic computation we may introduce the same theta-constraint as we had for χ c.t.
T (R ′ ) and χ c.t. imp T (R ′ ). This constraint is trivially satisfied in the limit R ′ → 0 and thus will not affect the result to the given accuracy. We therefore write noting that the theta-constraint allows us to send the upper-bound of the ζ i -integration to infinity. Setting R (D. 29) In the limit R ′ → 0, R ′ ≪ √ R ′′ and from eq. (4.20) we find, in agreement with [39] up to terms of order √ α s , that We note, as expected, that ⟨B T ⟩ NP = 2⟨B 1 ⟩ NP .

D.4 Thrust Major
The analytic determination of the non-perturbative correction to the mean value for the thrust major follows a similar approach to that for the total broadening in appendix D.3. As before we start with eq. (5.6) (D.32) In the limit R ′ → 0, one hemisphere, say H 1 , will contain a single emission, which we denote k 1 , with a transverse momentum of order T M Q/2 that will set the thrust-major axis. We consider the observable when all other emissions, k i , have transverse momenta less than δT M Q/2 where ϵ ≪ δ ≪ 1 thus We use eqs. (2.49) and (2.51) and introduce the rescaled variables x ℓ = 2p y,ℓ /(T M Q) and the twodimensional vectors ⃗ ζ i ≡ ζ i (cos ϕ i , sin ϕ i ). We perform the rescaling ζ i → ζ 1 ζ i and, as R ′ = R ′ 1 +R ′ 2 , we find (D. 34) We note that the factor of e −R(v) in eq. (2.49) has been dropped as we are in the limit R ′ → 0 where T M → T M,max and thus e −R(vmax) → 1. We also drop the ln(1/T M ) term as this will vanish in the limit R ′ → 0 as T M,max ≃ 1 as before.
We consider the three terms in the final pair of brackets: • For the ln(1/ζ 1 ) and (2 ln 2−2) terms there is no dependence on emissions in H 2 and therefore the integrals over H 2 will trivially evaluate to 1.
• For the ln(1/|x 2 |) term, we of course have a dependence on emissions in H 2 but we note that the integrals over H 1 and H 2 will still factorise. We set the upper-bound of the ζ iintegration to δ (to reflect the imposed restriction on transverse momenta of emissions in H 2 ) and recognise the integrals over H 2 as precisely that in eq. (B.28). We recall from eq. (C.37) that this will give (2/π)(1/R ′ 2 + ln 2) + O(R ′ ).
For the integrals over H 1 we note that, as was the case for ⟨B W ⟩ NP and ⟨B T ⟩ NP , rescaling the emissions in H 2 produces an exponent of R ′ (ζ 1 v) rather than R ′ 1 (ζ 1 v) (as was present for ⟨B 1 ⟩ NP ). To identify the precise form of any 1/R ′ divergence we take the NLL approximation of the Sudakov radiator, i.e. R ′ ℓ (ζ 1 v) ≃ R ′ ℓ (v), in eq. (D.34). Performing the ζ 1 -integration we find