Towards massless and massive event shapes in the large-$\beta_0$ limit

We present results for SCET and bHQET matching coefficients and jet functions in the large-$\beta_0$ limit. Our computations exactly predict all terms of the form $\alpha_s^{n+1} n_f^n$ for any $n\geq 0$, and we find full agreement with the coefficients computed in the full theory up to $\mathcal{O}(\alpha_s^4)$. We obtain all-order closed expressions for the cusp and non-cusp anomalous dimensions (which turn out to be unambiguous) as well as matrix elements (with ambiguities) in this limit, which can be easily expanded to arbitrarily high powers of $\alpha_s$ using recursive algorithms to obtain the corresponding fixed-order coefficients. Examining the poles laying on the positive real axis of the Borel-transform variable $u$ we quantify the perturbative convergence of a series and estimate the size of non-perturbative corrections. We find a so far unknown $u=1/2$ renormalon in the bHQET hard factor $H_m$ that affects the normalization of the peak differential cross section for boosted top quark pair production. For ambiguous series the so-called Borel sum is defined with the principal value prescription. Furthermore, one can assign an ambiguity based on the arbitrariness of avoiding the poles by contour deformation into the positive or negative imaginary half-plane. Finally, we compute the relation between the pole mass and four low-scale short distance masses in the large-$\beta_0$ approximation (MSR, RS and two versions of the jet mass), work out their $\mu$- and $R$-evolution in this limit, and study how their implementation improves the convergence of the position-space bHQET jet function, whose three-loop coefficient in full QCD is numerically estimated.

for any n ≥ 0, and we find full agreement with the coefficients computed in the full theory up to O(α 4 s ). We obtain all-order closed expressions for the cusp and non-cusp anomalous dimensions (which turn out to be unambiguous) as well as matrix elements (with ambiguities) in this limit, which can be easily expanded to arbitrarily high powers of α s using recursive algorithms to obtain the corresponding fixed-order coefficients. Examining the poles laying on the positive real axis of the Borel-transform variable u we quantify the perturbative convergence of a series and estimate the size of non-perturbative corrections. We find a so far unknown u = 1/2 renormalon in the bHQET hard factor H m that affects the normalization of the peak differential cross section for boosted top quark pair production. For ambiguous series the so-called Borel sum is defined with the principal value prescription. Furthermore, one can assign an ambiguity based on the arbitrariness of avoiding the poles by contour deformation into the positive or negative imaginary half-plane. Finally, we compute the relation between the pole mass and four low-scale short distance masses in the large-β 0 approximation (MSR, RS and two versions of the jet mass), work out their µ-and R-evolution in this limit, and study how their implementation improves the convergence of the position-space bHQET jet function, whose three-loop coefficient in full QCD is numerically estimated.

Introduction
In the LHC era, the scientific interest in jet physics has experienced a significant boost. It has been noticed by Cacciari and Salam (see e.g. [1]) that roughly 60% of experimental articles by ATLAS and CMS make use of jets (sprays of energetic particles traveling in nearly the same direction). Some of the theoretical efforts to better understand their dynamics seek to tame so-called hadronization effects, that is, corrections that arise from the fact that, even at high energies, what we observe are not quarks and gluons but hadrons. These corrections are enhanced since jet cross sections are not completely inclusive observables. The two most prominent approaches to pursue such goal are either parametrizing those corrections from first principles, or designing observables that are less afflicted by them. Since the most relevant non-perturbative effects come from low-energy particles, so far only soft hadronization has been dealt with. In the first approach, one can show within factorization theorems derived the context of Soft-Collinear Effective Theory (SCET) [2][3][4][5][6] (or the equivalent CSS formalism [7][8][9][10][11]) that the soft function (which describes largeangle soft radiation) can be re-factorized in a partonic part and a non-perturbative shape function, see Refs. [10,12,13] (other theoretical approaches can be found in [14][15][16]), whose first moment encompasses in a single parameter the leading hadronization effects in the tail of the distribution (the effect of hadron masses on power corrections has been investigated in Refs. [17,18]). In the second approach one defines grooming algorithms that consistently remove soft radiation in a way in which theoretical computations are not overly complicated, such as Soft-Drop [19] of its recursive version [20].
Taking a different perspective, one can investigate power corrections through so-called renormalons, that is, the asymptotic behavior of the associated perturbative series. This approach was applied to the soft function appearing in the factorized cross section for the doubly differential hemisphere mass distribution (which can be projected onto thrust and heavy jet mass in a straightforward manner) in Ref. [21], while a generalization for any other SCET-I type observable was worked out in [18]. The same soft function appears in cross sections for massless or massive quarks [22]. Earlier studies of renormalons in the context of e + e − event shapes can be found in Refs. [23][24][25]. In this article we investigate non-perturbative corrections coming from the remaining pieces of the SCET and a boosted version of Heavy Quark Effective Theory (bHQET) [26][27][28][29][30] factorization theorems, which are therefore associated to different physical scales. Power corrections take the generic form (Λ QCD /Q) n with n ≥ 1 and Q a magnitude with dimensions of energy appearing in the corresponding matrix element. As we shall discuss, a double pole implies a logarithmic enhancement in the power correction, which signals the presence of an anomalous dimension in the associated operator, with n f dependence at leading order within an operator product expansion (OPE).
While a thorough assessment of the large-order behavior of a series in full QCD is only possible in a few cases (such as the pole mass renormalon), one can carry out the relevant computations in the large-β 0 limit of QCD. In order not to violate important properties of the strong interactions such as confinement, this limit formally corresponds to making the number of active flavors n f tend to −∞. Since for phenomenological applications one takes a finite (and positive) number of flavors, in practice one considers the (equivalent) β 0 → ∞ limit. 1 In this setup one reorders the perturbative series by considering α s β 0 ∼ O (1) such that the leading contribution -beyond tree level -contains all terms of the form α s × (α s β 0 ) n ∼ O(α s ) for any n ≥ 0. The series might have an O(α 0 s ) term, such that treelevel and one-loop terms are exactly reproduced. The infinitely many terms that contribute at leading order are equally important and should be summed up. While this sum can be carried out unambiguously for convergent series, for asymptotic ones (that is, with null convergence radius), the sum is ill-defined. Using the Borel resummation method one can, however, define the sum of the series through the principal value prescription and assign an ambiguity to this procedure. It is believed this Borel sum is close to the "true answer" within the estimated ambiguity. For an excellent introduction to renormalons the reader is refereed to Ref. [31], while a very useful technical review can be found in [32].
In this article we compute in the large-β 0 approximation the following renormalized matrix elements: a) the hard function H Q for quark-initiated processes appearing in many SCET factorization theorems, b) the hemisphere jet function J n for massless quarks describing collinear radiation in the thrust, heavy-jet-mass and C-parameter event-shape cross sections, c) the SCET to bHQET matching condition C m showing up in factorization theorems for any massive-quark event-shape distribution in the peak region, in the form of the hard factor H m , and d) the bHQET hemisphere jet function B n , describing radiation which is soft in the rest frame of the boosted heavy quark, for the set of event shapes just mentioned. As a byproduct of these computation we will also obtain closed expressions for the universal cusp anomalous dimensions and their non-cusp counterparts for the SCET hard, jet, and soft functions, as well as the bHQET matrix element and jet function. 2 We find that all matrix elements are ambiguous, while the cusp and non-cusp anomalous dimensions have the same non-zero convergence radius. Additionally, we provide recursive algorithms to expand all matrix elements and anomalous dimensions to any order in perturbation theory.
The jet function(s) calculation's complexity is similar to that of the matching conditions, since the collinear thrust measurement is inclusive and the computation can be cast as the discontinuity of a 1-loop forward-scattering matrix element with a modified gluon propagator. This property does not hold for non-boosted momenta, and therefore to determine the large-β 0 soft-function it is required to compute Feynman diagrams of two-loop complexity. 3 We therefore relegate this calculation to the future (note that in Refs. [18,21] only the asymptotic behavior of the soft function was computed, which does not exactly reproduce the terms leading in 1/β 0 at low orders).
To carry out these computations we generalize the procedure presented in Ref. [32] to series with cusp anomalous dimension [ that is, for cases with 1/ε 2 poles at O(α s ) in dimensional regularization ]. We solve the corresponding RGE equation exactly at 1/β 0 accuracy. The closed form for the renormalized series can be written as the sum of a µ-independent inverse-Borel integral, where all ambiguities are confined, plus other µ-dependent terms with finite convergence radius. The former depends explicitly on Λ QCD , the non-perturbative QCD scale.
We also discuss the relation between the pole and MSR [34,35] masses in the large-β 0 approximation, and compute a closed form for its R-anomalous dimension γ R . We show that, even though the Borel transform of the series defining the MSR mass has a singularity at u = 1/2, the poles for γ R appear only for u ≥ 3/2. We analytically solve the fixed-order R-evolution [36] equation in this limit. From our result for the position-space bHQET jet function we compute the large-β 0 expression for two versions of the jet mass: as defined in Ref. [37] and an alternative scheme which we define in analogy to the soft-gap subtraction introduced in Ref. [38]. We compute their µ-(only in the original definition) and Ranomalous dimensions in a closed form.
The results of our computations are used in three phenomenological applications: i) the associated leading ambiguity is used to quantify the dominant hadronization correction, ii) the convergence radius of all anomalous dimensions is found, and iii) we study under which circumstances the various matrix elements converge best. In particular, we asses the effect of π-resummation and the various low-scale short-distance masses in the convergence of H Q and B n , respectively. We present a tentative procedure to remove the O(Λ QCD ) renormalon present in H m by expressing the pole mass in terms of the MS scheme. Phenomenological studies of this cancellation beyond 1/β 0 are relegated to future work.
This article is organized as follows: In Sec. 2 we provide an introduction to the largeβ 0 expansion and the general treatment of the leading 1/β 0 correction; in Sec. 3 we work out closed expressions for series which do not have a cusp anomalous dimension, providing a solution for their renormalization group equation, while in Sec. 4 we adapt the formalism for series with cusp anomalous dimension. Our prescription to compute the inverse Borel transform when there are poles in the integration path is presented in Sec. 5. The massive-quark self-energy is dealt with in Sec. 6, including the computation of closed forms for the relation between the pole and MS masses, the wave-function renormalization in the on-shell scheme, the MS mass anomalous dimension and the MSR R-anomalous dimension. The computation of the SCET hard and jet functions is carried out in Sec. 7. Analogous computations for bHQET are contained in Sec. 8, where it is shown how the leading renormalon found in each one cancels when expressing the pole mass in terms of a short-distance mass. In Sec. 9 we discuss two alternative versions for the so-called jet mass in the large-β 0 approximation, and analyze how effectively they remove the renormalon of the bHQET jet function. In Sec. 10 we provide a numerical estimate of the unknown 3-loop position-space jet function for a boosted heavy quark. Our conclusions can be found in Sec. 11. 2 General remarks about the large-β 0 expansion Let us consider a renormalized perturbative series starting with a non-vanishing tree-level coefficient, which for simplicity is taken to be 1 (we also assume that the lowest-order Feynman diagrams contain no gluon propagators), written in terms of the renormalized coupling α s (for simplicity we do not show the dependence of α s on µ, but we mark with a superscript that of the fixer-order coefficients c µ i ): c µ i,n n n f , (2.1) where the sub-index i counts the loops while n matches the powers to which the number of quark flavors n f is raised to. We can now rewrite the series in terms of powers of β 0 = 11 − 2n f /3 as follows: For a n,n−1 (that is, the leading terms in the large-β 0 expansion) only one term contributes to the sum in the second equality. In the limit we are interested in one considers α s β 0 ∼ O (1), and the series is reorganized accordingly: where the sum defining f µ j starts contributing at O(α j s ). The parameter β is order one and should not be mistaken with the QCD β-function, which we denote β QCD . The small parameter upon which we reorganize our perturbative series is obviously 1/β 0 . Since all terms that enter the definition of f j are equally important, the sum has to be carried out all the way to infinity. In this article we will be concerned with f µ 1 (β) ≡ β 0 δA(µ) only. In this counting one has β n α n+1 s ∼ O(α s ) for n ≥ 1. 4 Therefore, in d = 4 − 2ε dimensions the QCD β-function takes a simple form, β QCD = −2α s (ε + β), such that one can write dlog(β)/dlog(µ) = −2(ε+β). This result can be used to write down a closed form for the strong coupling renormalization factor in this approximation. To derive this result one has to recall the relation between the renormalized strong coupling α s (µ) and the (scaleindependent) bare coupling g 0 : α s (µ) = Z −1 α (µ 2 e γ E ) −ε g 2 0 /(4π) 1−ε , where we already adopt the MS convention for the renormalization scale µ 2 → µ 2 e γ E /(4π). Taking logarithms on both sides and applying a derivative with respect to log(µ) one gets dlog(Z α )/dlog(µ) = 2β. Since the µ dependence of Z α comes through β only, using the tree-level boundary condition Z α (β = 0) = 1, we obtain a unique solution for the renormalization factor The RGE solution for the strong coupling in the large-β 0 limit simply corresponds to the full QCD leading-logarithm (LL) result, which for the β parameter can be written as where the µ-independent quantity with dimensions of energy Λ QCD = µ e −1/[2β(µ)] determines the non-perturbative regime of QCD and specifies by itself a boundary condition for the strong coupling. Let us consider an unrenormalized perturbative series at O(1/β 0 ) expressed in powers of the bare coupling g 0 , whose (µ-independent) coefficients a (0) n,n−1 contain 1/ε m divergences with m ≥ n if the series does not have a cusp anomalous dimension, or m ≥ n + 1 if it has. 5 If the series carries an anomalous dimension of either kind, even after expanding it in powers of the renormalized coupling α s , additional (multiplicative) renormalization will be necessary: where in the second equality of the first line g 0 has been expressed in terms of β, and in the second line the Pochhammer symbol (a) n = Γ(a + n)/Γ(a) has been used to re-expand the series in powers of β.
It can be shown that when inserting the vacuum polarization function n times in the gluon propagator, and after performing a naïve non-abelianization [31,39,40] with the replacement n f → −3β 0 /2, one gets the following effective gluon propagator: The transverse part of the effective propagator does not contribute to the computations carried out in this article, therefore in practice we use the Feynman gauge with a modified gluon propagator. 6 The gauge dependence of course cancels after summing up all diagrams. 5 One has m ≥ n − 1 if the series has no anomalous dimension at all, making the one-loop term finite. 6 If the quantity of interest has no anomalous dimension one can renormalize the one-loop fermion bubble in the MS scheme, which implies the renormalized coupling is being used. When inserted n times one gets When summing up all possible insertions of the fermion bubble one gets an inverse Borel transform Therefore, also in this case the gluon propagator has a shifted power (−p 2 ) 1+u . The modified 1-loop computation is carried out directly in d = 4 and the following result for the renormalized is found The effective propagator for n ≥ 0 is O(1) in the large-β 0 counting, and therefore any term contributing at leading order can be obtained from gluon bubbles. In practice it is enough to compute the one-loop contribution to the series under study with a shifted power (−p 2 ) 1+h in the gluon's propagator denominator. This computation defines an ε-and h-dependent bare one-loop result which can be written as where Q represents a parameter with dimensions of energy that renders the computation non-vanishing in dimensional regularization. For the computations carried out in this article Q will stand for Q (center of mass energy), m (heavy quark mass), R (scale of the MSR, RS and jet masses) Q 2 = s, the (squared) invariant mass of a jet, andŝ = (s 2 − m 2 )/m. The bare coefficient appearing in Eq. (2.7) can be expressed in terms of the modified 1-loop computation of Eq. (2.11), When inserted in Eq. (2.6) the bare series can be expressed as follows: The F function is standard in renormalon calculations, see e.g. Refs. [24,31,32,40]. In the second line, a factor of u is included such that F (ε, u) is regular (and in general non-zero) when u → 0 for series without cusp anomalous dimension. Moreover, F µ (ε, 0) = F (ε, 0) does not depend on µ. Computations in the context of effective field theories for jets have 1/ε 2 poles at one loop and, when using a shifted gluon denominator, this translates into F (ε, u) ∼ 1/u when u → 0. In Sec. 4 we shall adapt the procedure to account for this fact.

Series without cusp anomalous dimension
In this section we review the standard procedure to find closed expressions for series without cusp anomalous dimension. In Sec. 6 we apply the formalism to the case of the massive quark self-energy. The presentation follows closely Ref. [32], but we provide explicit proofs for all intermediate steps. The regularity condition allows us to expand the function F (ε, u) in powers of β and ε Since the coefficients F µ j,0 do not depend on µ we denote them simply by F j,0 . Inserting this expansion in Eq. (2.13), shifting some indices, manipulating the order of sums and defining A 0 = 1 + δA 0 we obtain: where to get the second line we shift i → i − n and push the summation over n to the right-most position; to obtain the third line we shift j → j − m and swap the order between the sums over j and m; the last line takes its form after shifting j → j + i and interchanging the resulting sums over j and i. Since after renormalization one takes ε → 0, in the last line the j upper summation limit is set to zero. In Appendix A it is shown that the rightmost sum over n in the last line vanishes unless m = 0 or m i, and in the former case it equals −1/i. Therefore we isolate the m = 0 contribution and restrict the remaining sum to values of m larger or equal than i. Shifting in this last sum m → m + i one gets For negative values of j the first sum vanishes and the lower limit of the sum over i becomes −j. Reversing the sing of j and shifting i → i+j one gets for the coefficients of 1/ε j , defined in the first line of Eq. (3.2), the result For j = 0 the lower limit in the sum over i in Eq. (3.3) becomes 1, and the sum over m has only one term, m = 0. As shown in Appendix A, this sum can be carried out explicitly such that the second piece of the finite term takes also a simple form. Adding the divergent terms the following form for the bare series in terms of the renormalized coupling is obtained: From this equation one can read off the expression for the renormalized fixed-order coefficients in the large-β 0 limit. Setting µ = Q the non-logarithmic coefficients, which are factorially divergent, can be obtained (the logarithmic pieces can be easily obtained, as discussed next). Following the notation introduced in Eq. (2.3) and defining a i,i−1 ≡ a Q i,i−1 one obtains: whereγ i A will be defined in Eq. (3.15) below. At this point one can figure out closed forms for the various sums appearing in Eq. (3.6) in terms of the function F defined in the second line of Eq. (2.13). The second term of the finite contribution [ at leading order in 1/β 0 the finite term coincides with the renormalized series β 0 δA(µ) ] involves a UV subtraction related to renormalization, that makes sure the integrand is convergent when ε → 0. For the terms proportional to 1/ε j with j > 0 subtractions are not necessary (as we shall see, the 1/ε term fully determines the anomalous dimension of the series in the large-β 0 limit), and all of them can be worked out at once: The closed expression for the first contribution to the renormalized series takes, as expected, the form of an inverse Borel transform: where again the UV subtraction makes the integrand converge as u → 0. With this result we can write a closed form for the renormalized series. Defining A(µ) = 1 + δA(µ) we get The ambiguities appear in the first integral [ the second is not ambiguous because F (ε, 0) has poles for ε > 0 only ], and written in this way it seems they might depend on µ. Later in this section it will shown that they do not depend on any renormalization scale. The series is renormalized multiplicatively A 0 = ZA = (1 + δZ)(1 + δA) 1 + δZ + δA where A(µ) is finite in the ε → 0 limit and δZ ≡ Z −1 in the MS scheme has only 1/ε j poles. We have discarded the crossed term δZ × δA since it is O(1/β 2 0 ), making renormalization effectively additive at O(1/β 0 ). It is straightforward to find a closed expression for δZ For series with no anomalous dimension F (ε, 0) = 0 and therefore Z = 1. Since F i,0 and, consequently F (−β, 0), are µ independent, the Z factor has no dependence on µ beyond that of β. Since in this limit one also has log(Z) δZ and log(A) δA, the anomalous dimension γ A ≡ dlogA/dlog µ for the renormalized quantity A(µ) can be computed as the negative of the logarithmic µ-derivative acting on Z. We use the chain rule in the large-β 0 approximation for d = 4 − 2ε effectively replacing Taking into account that γ A is finite as ε → 0 it is clear that the anomalous dimension must come entirely from the ε term in the previous equation acting on the 1/ε contribution to Z. The remaining contributions are individually O(1/ε k ) with k > 0 and therefore must cancel when added up to render a finite γ A . Let us show this fact with an explicit computation: 7 where one can shift j → j − 1 and i → i + 1 in the first term of the second line to show that all divergent terms cancel, obtaining the first equality in the last line. Therefore, the corresponding fixed-order coefficients have the following form: By the same reasons exposed above,γ i A are µ independent, and the dependence of γ A on the renormalization scale happens through β only. The RGE for A(µ) can be easily solved, and in the large-β 0 takes the form [ β i ≡ β(µ i ) ]: where in the second line we have pulled out the term which is not regular when either β i tend to zero. In the last equation we use the expanded version of the anomalous dimension, 7 The same result can be obtained applying the derivative directly to the integral form of δZ: where one can shift j → j − 1 in the first term of the second equality to see that all divergent terms cancel. and given the fact that γ A is not ambiguous, provided the values of β are not too large, the sum can be carried out to arbitrarily high orders. This approach, employed also in Ref. [41], provides an effective way of computing the numerical integral in the second line to arbitrary precision. Furthermore, the first term in the last line corresponds to LL accuracy, while adding n additional terms from the sum yields N n LL accuracy. This result shows that the difference of two series with different values of µ is renormalon free. Taking the values µ 2 = µ, µ 1 = µ 0 and using Eq. (3.11) to obtain the boundary condition A(µ 0 ), from Eq. (3.17) we obtain the following closed resummed form for δA(µ): . Written in this way, we have explicitly singled out the dependence on µ -as expected, renormalon-free -from the ambiguous piece which now depends on the matching scale µ 0 . However, since the matching scale is arbitrary, the ambiguities must be also independent of µ 0 . To proof this explicitly we need the following relations where in the first line we have simply added and subtracted F (0, 0)(µ 0 /Q) 2u , and noted that F 0,0 = γ 0 /2. The solution to the RGE equation for β in Eq. (2.5) has been used to re-write the argument of the logarithm in the second line and to show the independence of Λ QCD on µ in the last. The integral in the second line is solved by Taylor expanding the term in square brackets around u = 0, integrating term by term and summing up the resulting series: Introducing the results in Eq. (3.19) into Eq. (3.18) and defining β Q ≡ β(Q) one finds (3.21) The same answer can be obtained using Eq. (3.19) with µ 0 = µ directly in Eq. (3.11). The term with potential ambiguities depends only on the value of Λ QCD , and once this parameter has been fixed it can be computed once and used for any value of µ. 8 We finish this section computing the dependence of the fixed-order coefficients on ≡ log(µ/Q). Making this logarithmic dependence explicit the series can be written as with a i,i−1 ≡ a i,i−1,0 . Applying a derivative with respect to log(µ) (which acts both on β and ) and imposing that the series has the right anomalous dimension we obtain where in the first and second terms of the top line we have shifted j → j + 1 and i → i − 1, respectively, as compared to the summation in Eq.

Series with cusp anomalous dimension
In this section we present a novel computation, adapting the derivation presented in Sec. 3 to series with cusp anomalous dimension. We will use our result to quantities relevant for effective field theories applied to jets initiated by massless and massive quarks.

Derivation
As already anticipated, the function F (ε, u) defined in the second line of Eq. (2.13) is now singular when u → 0. This is related to 1/ε 2 poles in one-loop amplitudes, which signals as u → 0. Taking into account that F (ε, 0) = 0 one has that the anomalous dimension is zero, and from Eq. (3.21) the closed from for the series can be readily inferred: H(0, u) . the presence of Sudakov (double) logarithms at this order (the function is however regular at ε = 0). Therefore we define The bare series A 0 = 1 + δA 0 , after performing manipulations completely analogous to those used in the derivation of Eq. (3.2), takes the following form: where again the sum in j has been restricted to non-positive values, since those are the only ones necessary for renormalizing the series at O(1/β 0 ). As already discussed, the sum over n is zero unless m = 0, 1 or m ≥ i + 1. The result for m = 1 has been already presented in Eq. (3.5), while for m = 0 the sum is worked out as follows where we have expressed the binomial in terms of factorials and have used (i − 1)! = i!/i, n 2 (n − 1)! = n × n! in the second equality, such that a binomial can be written in the third. The last form is obtained with Eq. (A.9), derived in Appendix A. All in all we find If j < 0 the sum over m vanishes and only the terms in the second line contribute. Given the "max" condition in the lower limit of the sum over i we treat j = −1 separately. On the other hand, j = 0 implies m = 0, leaving a single term in the sum: From this result we obtain the following expression for the non-logarithmic term of the renormalized fixed-order coefficients: where we have anticipated the results for the anomalous dimensions to be found in Eqs. (4.14) and (4.16). We derive next closed expressions in terms of G(ε, u) for the various sums. We start with the divergent terms, and work out the following sums which are valid for any value of j: (4.7) This time, the coefficients A 1/ε j exhibit explicit µ dependence, which however happens to be linear in log µ. Using G µ (ε, u) = (µ/Q) 2u G(ε, u) we find G µ i,1 = G i,1 + 2 log(µ/Q)G i,0 and the related result For series which involve harmonic numbers we use identity for which a proof is provided in Eq. (A.10) of Appendix A. With this result we obtain a closed form necessary to find a resumed expression for A 1/ε , while for higher-power divergences we obtain integrals which do not involve UV subtractions at ε = 0: (4.10) Finally we need the following closed forms for sums related to the renormalized series: The inverse Borel transform integral involves now two subtractions around u = 0, which are necessary to render the integrand finite at this value given the overall 1/u 2 factor. The contribution in the first line involves a u-derivative of the once subtracted G(ε, u) function around ε = 0 that ensures a smooth behavior towards the upper integration limit. On the other hand, the u derivative makes this term dependent on µ. We can also explicitly display the µ dependence of the derivative in the first line of Eq. (4.11) and the subtraction in the inverse Borel transform (last line) We can combine the various results to write down a closed form for the bare series, which we split as δA 0 = δZ +δA(µ), with renormalized series A(µ) = 1+β 0 δA(µ) and renormalization factor δZ = Z − 1 = δZ nc + 2 log(µ/Q)δZ cusp : 9 .
Ambiguities are again confined to the inverse Borel transform integral, and it will be shown in Sec. 4.2 that they do not depend on the renormalization scale. We finish this section deriving a closed expression for the cusp (Γ A ) and non-cusp (γ A ) anomalous dimensions, which we define as dlogA/dlogµ = γ A + log(µ/Q)Γ A . Following the same steps employed to find the anomalous dimension in Sec. 3 we find for Γ A the following closed form and fixed-order expansion: 14) making clear that the cusp anomalous dimension is not afflicted by renormalons. For a series without cusp G(ε, 0) = 0 and therefore Γ A = 0. Likewise, the non-cusp anomalous dimension is derived as follows where we have pulled out of each sum the 1/ε term and have shifted i and j by one unit wherever necessary to write all sums in terms of ε j and (−β) i (except for some terms involving harmonic numbers which have an additional explicit 1/β factor). It is not hard to realize that all terms cancel out except for the very first sum, from which the fixed-order coefficients can be read off. The sum can be expressed in a closed form: Since, as we shall see, the function G(β, 0) has singularities on the positive real axis only, no pole is crossed when carrying out the integration in the second line, making thus clear that the non-cusp anomalous dimension has no ambiguities. Furthermore, this integral has no divergence at x → 1 because the term in squared brackets tends to zero in this limit.
To obtain the last term in the second line we have used 10 10 One can obtain the same result applying the derivative directly into the integral form of δZ 1/ε nc . To that end we change variables h = −βx in the second term of Eq. (4.13). To tame the divergence occurring at x = 1 we add a regulator δ to the argument of the logarithm. Defining G (y, 0) ≡ dG(y, 0)/dy one has where integration by parts has been used to obtain the second line, while to get the third one writes log(δ) = log(1 − δ) − where in the first step the identity a i −b i = (a−b) i−1 j=0 a j b i−1−j is applied, in the second we integrate term by term, and identifying the harmonic number's definition the last expression is obtained. In Eq. (4.16) we switch variables ε = −βx to map the integration to the unit range [0, 1]. If the series has no cusp anomalous dimension G(ε, 0) = 0 and the integral in the second line of Eq. (4.16) vanishes. In this case [dG(−β, u)/du] u=0 = F (−β, 0) and the result in the last line of Eq. (3.14) is recovered.

Exact solution to the RGE
Using the chain rule we can express the µ derivative in terms of a β derivative. Furthermore, we can write the factor log(µ/Q) multiplying Γ A as the difference of inverse β couplings: where in the last line we have pulled out those terms that after integration are not regular if either β µ or β µ 0 tend to zero. All regular contributions are kept in a running kernel which is defined as the difference of two terms, each of them depending on a single renormalization (4.20) Expanding strictly to leading order in large-β 0 one can express the renormalized series as: where we have used 1/β µ − 1/β µ 0 = 2 log(µ/µ 0 ), relation that stems from the solution to the RGE equation for the β coupling in Eq. (2.5). Let us work out each term appearing in Eq. (4.20), starting with the ones involving Γ A . Using the closed form in Eq. (4.14), noting To work out closed expressions for the term involving the non-cusp anomalous dimension we use the expression for γ A as given in Eq. (4.16), taking into account that where to arrive at the second line the β derivative has been computed applying the chain rule. Changing the integration variable ε = −βx in Eqs. (4.16) and (4.23) one finds where to obtain the second equality the following manipulations have been performed to the integral in the second line To obtain the second line we add and subtract G(0, 0) within the square brackets in the first line, collect terms that depend on ε and β separately [ such that G(0, 0) is subtracted in each one ], followed by renaming ε ↔ β in the term that only depends on β (which appears in the third line); to get the last expression we switch the order of integration in the second line integral and flip the sign of the two integration variables ε → −ε, β → −β in the third. Since 0 ≤ −ε ≤ β µ , the two integrands in the last line become equal in the vicinity of the pole that happens at β = −ε, and since the singularities are approached from opposites sides they cancel when the two terms are added up if the integration is performed symmetrically. Equivalently, to obtain the last line of Eq. (4.24) one can use a regulator δ ≥ 0, and keeping track of the changes of variables carried out one obtains the following regularized result: where in the last expression we take the δ → 0 limit, now smooth. Combining the results in Eqs. (4.22) and (4.24) and using the 2-loop cusp anomalous dimension coefficient where the logarithm in the first line appears from the solution to the RGE equation for the strong coupling 1/β Q − 1/β µ = 2 log(Q/µ). The last term in the second line is constant and therefore vanishes when computing ω(Q, µ 1 , µ 2 ). In Appendix B we give an alternative derivation of this result based on series expansions. Introducing Eq. (4.27) and the form shown in Eq. (4.13) for A(µ 0 ) in Eq. (4.21) we obtain a resummed expression for the renormalized series A(µ), with β ≡ β(µ): .
We have displayed all dependence on µ and µ 0 explicitly, and shown that the ambiguous term in the second line is independent of µ. It can be shown that it is also µ 0 independent, and to that end we write the expression between curly brackets in the second line as follows: The twice-subtracted integral resulting from the last line can be solved by Taylor expanding around u = 0, integrating each term and summing up the series: where in the last step we have used Eq. (2.5) to re-write the sum's prefactor. Using Eqs. (4.29), (3.20) and (4.30) in Eq. (4.28) one finds (4.31) The same result can be found using Eqs. (4.29), (3.20) and (4.30) with µ 0 = µ directly in Eq. (4.13). In this closed form ambiguities are manifestly µ 0 independent, but explicitly depend on Λ QCD , hence making the connection to non-perturbative physics manifest. The dependence on µ is fully contained in terms related to the anomalous dimensions. If G(ε, 0) = 0, the first term in the second line vanishes, along with the third andΓ 0,1 A , thus recovering the result in Eq. (3.21).
For the sake of completeness, we present the form of the evolution kernel between two renormalization scales in closed and expanded forms. For the former, defining β i ≡ β(µ i ), we find: where we have intentionally written all dependence on µ i through β i . Expanding G(ε, 0) and integrating each term one arrives at the following perturbative expansion for the evolution kernel: where the first and second terms in square brackets in the first line are the LL and NLL contributions, respectively, while to obtain an N n LL-accurate result one has to add the first n − 1 terms of the sum in the second line.
We finish this section by presenting an expression for the fixed-order coefficients multiplying powers of = log(µ/Q). In this case, due to the cusp anomalous dimension, one has an additional power of at each order as compared to series without cusp: Imposing that the series reproduces both types of anomalous dimensions when acting with a derivative with respect to log µ one finds where we have again equated equal powers of β and in the first line to obtain the relations between coefficients shown below, treating separately the cases j = 0, 1 in the second and third lines, respectively. To obtain the last equality in the second line we have used the expressions forγ i A and a i,i−1 in Eq. (4.16) and (4.6), respectively. Likewise, the expression forΓ i A in Eqs. (4.14) and the result for a i,i−1,1 in the second line has been used to obtain the first equality in the third. To find the last expression in the third line we apply the first relation recursively (j − 1) times until we reach j = 2, such that the previous line can be used (therefore, this last expression is valid for j ≥ 2). One can obtain identical results expanding the logarithms in the last line of Eq. (4.5), complementing the derivation carried out in Eq. (3.25) with a direct computation of G µ i,1 :

Principal value prescription
As already discussed, the inverse Borel transform presented in Eqs. (3.21) and (4.31) is ambiguous (that is, there are poles in the integration path) and one needs a prescription in order to assign a finite value, which in turn defines the so called Borel sum. In this section we present the principal value (PV) prescription, which will be used in the applications discussed throughout the following sections. Let us consider a function f (u), denoting the set of its poles of order m in the real positive axis by u A pole of order m usually implies lower order poles at the same position, so u (m) n stands for all poles of order 1 ≤ i ≤ m at the same position. The asymptotic expansion for f (u) therefore reads: To compute the integral's principal value along the real axis we split the integration in several pieces: In case there are infinitely many poles the last integral in the first line is zero and the first sum in the second line extends to infinity. The result in Eq. (5.2) is independent of the specific values of δ i as long as they verify the condition just stated. Thanks to the pole subtractions, only the last term within the sum in the second line of Eq. (5.2) needs regularization. The PV prescription for an order-i pole, with i a natural number, is defined as 12 where δ > > 0 and C ,un denotes an upper (lower) semi-circumference of radius centered at u n in the (anti-)clockwise direction. In our prescription we want to preserve two important properties of ordinary integration: a) integrating an odd function f (u) with respect to the point u 0 , that is f (u − u 0 ) = −f (u 0 − u), in the interval {u 0 − δ, u 0 + δ} results in zero; b) if the integrand is real, so is the integral. Therefore we will take the average of PV + and PV − as our definition: where in the second line we give the standard prescription to assign an ambiguity to the Borel sum, which turns out to be a circular path enclosing the pole, such that the residue theorem can be used. Therefore, if the integrand is real, so is the ambiguity. Moreover, it can be shown that the numerical value obtained in this way is close to the smallest term of the series (in absolute value), which is normally taken as the intrinsic ambiguity of an 11 Strictly speaking it is enough to require ui + δi ≤ ui+1, but we implement the most restrictive condition to avoid overlapping between consecutive patches. 12 In full QCD poles turn to branch cuts, and the PV prescription has to be adapted accordingly.
asymptotic series [31]. Needless to say, only simple poles contribute to the total ambiguity, which is simply the sum of their residues: such that the divergence 1/ε 2k−1 cancels out when adding the two terms, leaving a finite remainder. We can write a general formula for the PV integral in Eq. (5.3) when the symmetric prescription is used using the modulo operation: with δ 0,mod(i,2) a Kronecker delta symbol which is non-zero only if i is an even number. One can figure out the contribution of each singularity to the asymptotic behavior of the fixed-order coefficients by expanding the pole in powers of u and integrating term by term, which is equivalent to the replacement u n → β 1+n Γ(n + 1): Therefore one has that for larger (smaller) pole's order (u 0 ) the asymptotic behavior becomes worse. Using the asymptotic expansion shown in Eq. (5.1), the contribution to the fixed-order coefficients of a series coming from the poles of its Borel transform takes the following form: The exact coefficients are dominated by the poles closest to the origin of the Borel plane.
We can also work out the case µ = Q taking into account that where we have shifted n → n−i, reversed the order of the sums and expressed the Pochhammer symbol in terms of a binomial to obtain a canonical expansion in powers of β µ and . ℓ ℓ + p p p Applying this result to the contribution coming from a pole of order i as shown in Eq. (5.8), for n 1 we get where we have used that for large n the sum is dominated by small values of j, therefore setting n − j → n (this approximation is more accurate for small values of i). The dependence on µ/Q does not depend on n, and therefore is common to all powers of β µ generated by a given pole. With this result the asymptotic expression for the µ-dependent fixed-order coefficients is obtained: Since the approximation used in Eq. (5.11) is exact if = 0, the result above is in agreement with the second line of Eq. (5.9) for µ = Q. When applied to simple poles one can use the identity (1) n−1 = Γ(n).
For the numerical computations in this article we have developed two independent numeric codes, written in Mathematica [42] and Python [43], that agree within 15 decimal places. While we use built-in native functions in Mathematica, for evaluating special functions and performing numerical integrals in Python, the SciPy module [44], that builds on the NumPy package [45], is used.

Application to the MS and MSR masses
In this section we review the calculation of the relation between the MS and pole quark masses. From this computation we will also obtain the wave-function renormalization for a massive quark in the on-shell scheme, necessary for the SCET to bHQET matching coefficient computation carried out in Sec. 8, as well as the anomalous dimensions for the MS and MSR masses.

Massive quark self-energy and the MS mass
The Feynman diagram shown in Fig. 1 two Dirac structures in the usual way: 13 with 2 F 1 the Gauss hypergeometric function, g 0 and m 0 the bare strong coupling and mass, respectively, and C F = (N 2 c − 1)/(2N c ). The relation between the scale dependent MS mass m(µ) and the scale-independent pole mass m p at one loop with a modified gluon propagator takes the following form where Z MS m (which depends on µ such that the right-hand-side does not) contains only 1/ε n poles and ensures m(µ)/m p is finite. From this result we can identify the function b(ε, u) in Eq. (2.11), with Q = m p , and compute F (ε, u) as defined in Eq. (2.13). To make later formulas as simple as possible we define the function corresponding to the series δ MS (m) ≡ m p − m, which only differs by a sign and has Q = m ≡ m(m): 13 Our definition for the quark self-energy differs from others by a sign.
-25 - From this result one can find a closed expression for the MS mass anomalous dimension γ m , first computed in Ref. [46] (see also Ref. [32]): 14 In the second line we have written the expression in a way that facilitates its expansion in powers of β by computer programs, such that the fixed-order coefficients can be easily obtained. To carry out the expansion of the exponential of a series we use the following recursive relation, which is derived imposing Our result reproduces the full QCD leading flavor structure up to 5 loops [47][48][49][50], collected in the first column of Table 1. Since the pole closest to the origin in the function F MS (β, 0) is located at β = −2.5, the convergence radius of γ m (β) is ∆β = 5/2. Interestingly, since all poles in F MS (β, 0) lay on the positive real axis, γ m has no singularities in the physical region α s > 0. The series is thus not ambiguous and provided that β < 2.5 the partial sum will eventually converge to the exact result.
14 Note that the MS mass anomalous dimension is defined as µ m dm dµ = 2γm, and therefore γm = −γ MS /2.  To obtain the fixed-order relation between the pole and MS masses we can apply Eq. (3.7), and for that we need 15 where in the second line the argument of the exponential has been expanded in powers of u. The fixed-order coefficients obtained using Eq. (6.5) reproduce full QCD leading flavor results up to O(α 4 s ) [53][54][55][56], collected in the first column of Table 2. In Fig. 2(a) we show how the perturbative MS evolution for the bottom quark mass is quickly convergent and already at NLL is very close to the exact value. 16 In Fig. 2(b) we compare the MS and MSR (see Sec. 6.3) evolution for the top quark mass for values of the renormalization scale smaller than m t . Both schemes coincide for µ = m t and, even though the MS evolution is unphysical for those scales, growing out of control and becoming even larger than the pole mass, we observe that m t (m t /2) m pole t . The numerical value of the MSR mass monotonically and smoothly grows as R → 0 to reach m pole t precisely at this limiting value. R-evolution for R < m t is physical and should be used for processes that probe scales smaller than the quark mass itself.
The modified one-loop massive quark wave-function renormalization in the on-shell 15 An equivalent result was found in Refs. [51,52].   scheme Z OS ξ has the form: In this approximation one has Z OS ξ = (1 + u)Z MS m m(µ)/m p , and therefore Z MS m m(µ)/m p and Z OS ξ agree in full QCD at one loop.

Renormalon pole normalization and the RS mass
To study the asymptotic behavior of the series m p − m we look at the poles of its Borel transform. To that end it is enough to analyze m B (u) = F MS (0, u) − F MS (0, 0) /u, a function that can be seen in Fig. 3(a). The subtraction ensures there is no singularity at u = 0, but plays no role for poles with u > 0. We find the following asymptotic expansion: The difference between m B (u) and its asymptotic expansion is shown in Fig. 3(b) (in the range 0 < u < 8 the largest deviation is only 20%). We can apply Eq. (5.7) to the various terms in Eq. (6.8) to obtain the contribution of each pole to the series coefficients a MS n,n−1 defined in Eq. (6.14). This can be seen in Fig. 4(b Table 2. Numeric values for the non-logarithmic coefficients for various short-distance masses a n,n−1 . while the contribution of the first four poles appear in colors other than black. The picture clearly shows how the u = 1/2 pole overly dominates already at low values of n. The black dots represent the non-ambiguous contribution coming from the second term in Eq. (6.13), which does not diverge as n → ∞ and quickly becomes completely negligible. The leading pole is located at u = 1/2, and its residue fixes the norm of the pole-mass renormalon N β 0 to the β 0 -independent value N β 0 1/2 = C F e 5/6 /π 0.976564. In Ref. [35] a sum rule to compute the value of N (n f ) 1/2 was derived using R-evolution. In Fig. 5(a) it is shown how the value of N counterparts 180 MeV, 215 MeV and 250 MeV found in [57], which are a factor of 2.5 larger). In Fig. 4(a)  Equation. (6.9) can be used to write down the relation between the pole and RS masses [58] at leading order in 1/β 0 (see coefficients in the second column of Table 2): 17 where in the second line we show the fixed-order relation between m p and m RS , obtained using Eq. (5.8) with m = 1 and u 0 = 1/2. The series has the same leading renormalon at u = 1/2 as the MS or MSR masses (see Sec. 6.3). Therefore, as can be seen in Fig. 4(b), the RS fixed-order coefficients a RS n (red dots) tend to a MS n (in gray) defined in Eq. (6.14) for large n. The difference m MSR (R 1 ) − m MSR (R 2 ) is free from any renormalon ambiguity and the integration converges for u → ∞ provided that max(R 1 , R 2 ) > Λ QCD . The closed 17 For simplicity we consider only the "non-primed" version of the RS mass.
with Ei(x) the exponential integral function. One can compute the R-anomalous dimension for the RS mass in the way explained in Sec. 6.3, and the result is very simple, β 0 γ RS (β) = 4C F e 5/6 β, with only the one-loop termγ 0 RS non-zero. This implies that Eq. (6.12) is fully compatible with the result shown in Eq. (6.17) truncated at leading-log, such that for the RS mass exact and fixed-order evolution exactly coincide. To relate the RS and MS masses one can take the difference δ MS (m) − δ RS (m): 1 − 2u (6.13) with β m = β(m). The matching relation above is free from the u = 1/2 renormalon, but has higher order singularities. After determining m RS (m) from Eq. (6.13) the RS mass can be computed at any other scale R using Eq. (6.12).

The MSR mass its and R-evolution
The MSR scheme [34,35] is defined from the series relating the pole and MS masses δ MS (m): a MS n,n−1,i log i µ R (6.14) where in the second line we have written an exact form in the large-β 0 limit with β R = β(R).
The fixed-order series us expanded in powers of β R or β µ = β(µ) in the first line, where a MS n,n−1,i are the coefficients that appear in the fixed-order expansion of δ MS . The definition above implies m MSR (m) = m, such that it is convenient to take R = m as a boundary condition to determine the MSR mass at any other scale R. When taking the difference of MSR masses at different values of R, ∆ MSR (R 2 , R 1 ) ≡ m MSR (R 2 ) − m MSR (R 1 ), the inverse-Borel term will not cancel. However, due to the R prefactor in Eq. (6.14), the integrand's pole at u = 1/2 becomes R-independent and thus vanishes when subtracting two series, making the difference less ambiguous than each series taken by itself.
Taking into account that the pole mass is R-independent, one can compute a closed form for the MSR mass R-anomalous dimension using the second line of Eq. (6.14): where the pole at u = 1/2 in the integrand of the first-line cancels, pushing the leading renormalon to u = 3/2, thus making γ R less ambiguous than δ MSR . The corresponding fixed-order coefficients γ n R , collected in the second column of Table 1, can be obtained from to the perturbative expression for δ MSR : reproducing the leading flavor structure of full QCD up to four loops [35]. In Fig 6(a) we compare the exact form for γ R (red dashed line) in Eq. (6.15) with the partial sum of Eq. (6.16) at various orders for a relatively large value of α s . We observe the expected behavior: the partial sum oscillates around the exact value until n 10, and starts diverging after that order. In Fig. 6(b) we show the dependence of the closed form for γ R with α s (gray band) with the partial sum up to γ n R for 5 ≤ n ≤ 10.
One can integrate the R-anomalous RGE equation using γ R as given in Eq. (6.16). Switching the integration variable from R to β through Λ QCD using Eq. (2.5) and reversing the integration order, the β integrals can be carried out analytically. Following this procedure one of course finds ∆ MSR (R 2 , R 1 ) = δ MSR (R 1 ) − δ MSR (R 2 ), with δ MSR defined in Eq. (6.14). This is an important cross check for our derivation. For completeness we provide the perturbative R-evolution running kernel in the large-β 0 limit, simply integrating each term in the fixed-order expansion for γ R : To obtain the result in the second line integration by parts has been used to relate I n (β) ≡ dβ β n−1 exp[−1/(2β)] for consecutive values of n. Iterating this relation until n = 0 each E n (−x) can be written in terms of Ei(x) and a finite sum. 18 Combining this last sum with the one overγ n R one arrives at the final expression. If the sum includes up toγ n R the result is N n LL accurate. In Fig 5(b) we compare the value m MSR t (R) computed through Eq. (6.14) (cyan dashed line) with the R-evolved result obtained with Eq. (6.17) at various orders (black joined dots) for a small value of R.
The pole mass value estimated from the fixed-order expansion of δ MSR (R) together with m MSR (R) obtained with perturbative R-evolution depends at any finite order on R. This matching relation can be expressed in powers of β µ , and as long as µ/R ∼ O(1) there are no large logarithms in the expansion. The numerical value of m p will also depend on µ if the series is truncated. The coefficientsâ MS n,n−1,i can be easily obtained demanding that the series is µ-independent. This provides recursion relation that can be solved exactly: 19 where the first equality is used i times until one reaches i = 0 to obtain the second. This result is equivalent to Eq. (5.10), derived re-expanding β R in powers of β µ and . For the RS mass (or equivalently for the MSR mass if n is large enough) this sum can be carried out exactly, and using Eq. Therefore, at large orders R gets effectively replaced by µ in Eq. (6.11) if we choose to express the series in powers of β µ . This makes clear that when two series with the same leading renormalon are subtracted, the ambiguity will cancel only if both are expressed in terms of α s evaluated at the same renormalization scale µ. At the sight of Eq. (5.12), and given that at large orders perturbative coefficients are dominated by the lowest singularities of the Borel transform, the same conclusion can be drawn from poles located anywhere in the positive real axis and for series carrying an anomalous dimension. In Fig. 7 the pole mass is estimated from its perturbative relation to the MSR mass. The colored dots use N n LL-accurate R-evolution to compute m MSR fix µ and pick three values of R. We observe that for small R the exact value of the pole mass, computed as m t + PV[δ MS (m t )], is properly estimated -within perturbative uncertainties -at low orders (since there are somewhat large logarithms for the choice R = 2 GeV, the partial sum does not exactly hit the exact value for low n). For larger R it takes more orders to adequately predict the pole mass. The explanation for this behavior is simple: since the MSR mass at small values of R is closer to the pole mass (formally both coincide if R = 0) with a few matching corrections one already predicts its value. On the other hand, the large-order behavior depends only on µ. In the right panel we keep R fixed but change the value of µ (green dots are identical in both panels). The three series become flat more or less at the same time (since this depends only on R) but for smaller µ the divergent behavior for n > 8 is more pronounced.

SCET computations
Hadrons produced in e + e − collisions at high energies often adopt di-jet configurations. Jets are sprays of collimated particles traveling very fast in nearly the same direction. In such configurations, final-state particles are either energetic and belong to a jet (collinear) or have small energies, populating the phase-space region in between the jets (soft). Event shapes are kinematic variables constructed from the energy and momentum of the produced particles that measure how di-jet an event looks from geometric considerations. In di-jet configurations one can use SCET to derive factorization theorems that are valid at leading order in the power-expansion parameter of the theory. For thrust one has: In this section we compute the large-β 0 limit of two matrix elements that appear in Eq. (7.1): the hard matching coefficient H Q (common to all jet observables) and the hemisphere jet function J n (relevant for thrust, heavy-jet-mass and C-parameter). The computation corresponds to one-loop diagrams with a modified gluon propagator, and since SCET matrix elements have cusp anomalous dimension, the technology developed in Sec. 4 will be applied.

Hard function
The hard function H Q (Q, µ) is the modulus squared of the matching coefficient C H (Q 2 , µ) between SCET and full QCD di-jet operators. For its computation one needs to determine some matrix element (or Green function) of the vector current operator in both theories.
Since C H does not depend on the specific process used for its computation one usually chooses the simplest possibility, which in this case is the matrix element between vacuum and a pair of on-shell quarks. This choice implies that all diagrams in SCET as well as selfenergy diagrams in QCD are scaleless and vanish. The QCD computation is often referred to as the quark vector form factor. Therefore, to obtain H Q at leading order in the large-β 0 expansion one only needs to compute the diagram shown in Fig. 8 with a modified gluon propagator. The virtual correction δC H (Q) to the bare vector form factor C H = 1 + δC H is then From this result one can obtain the virtual correction to the bare hard matching coefficient defined as H Q (Q) = |C H (Q 2 + i0 + )| 2 ≡ 1 + δH Q (Q) + O(1/β 2 0 ). For e + e − even shapes one has Q 2 > 0 and keeping only terms linear in 1/β 0 the correction δH Q equals twice the real part of δC H : To obtain the cusp and non-cusp anomalous dimensions it is enough to consider C H (Q 2 ).
and taking into account that Γ C Q = −2Γ cusp a closed form for the cusp anomalous dimension is obtained from the equation above: where the second line has been partially expanded to easily compute the fixed-order coefficients recursively. Our result agrees with Ref. [33], and the fixed-order coefficients reproduce the leading flavor structure of full-QCD up to O(α 4 s ) [60][61][62][63], collected in the first column of Table 3. The convergence radius of Γ cusp is set by the distance to the pole closets to the origin of the function G Q (−β, 0), which happens to be at β = −2.5. Furthermore, Γ cusp (β) has no singularities for positive β, and the series is, as expected, free from renormalons. In Fig. 9(a) we compare the exact result for the cusp anomalous dimension to the partial sum of the fixed-order expansion up to 9 loops. Even though we use an unphysically large value for the strong coupling α s = 0.9, the fixed-order expansion converges to the exact value.
To determine the non-cusp hard anomalous dimension γ H (β) the following function is needed: where in the second line we have written a partially expanded expression which is the  starting point of a convenient recursive algorithm to obtain the fixed-order coefficients. This result can be used in Eq. (4.16) together with the analytic form of G Q (−β, 0) provided in Eq. (7.5) to obtain a closed form for γ H (β) ≡ γ C H (β), as well as its fixed-order coefficients. We find full agreement for those when compared to the leading flavor structure of full QCD up to O(α 4 s ) [63][64][65][66], see second column of Table 3. To multiply out two series we use the general result In our case, one of the two series results from the expanding the exponential in Eq. (7.6) using the recursive formula in Eq. (6.5). The convergence radius of the series is again ∆β = 2.5. In Fig. 9(b) we compare for α s = 0.9 the exact and fixed-order results for γ H , finding again that the latter converges to the former already at 6 loops.
The last piece of information necessary to determine  Table 4. To compute δH Q one considers instead G H (0, u) = 2 cos(πu)G Q (0, u), and to determine the fixed-order coefficients the Taylor expansions for cos(πu) and G Q (0, u) must be used in Eq. (7.7).  Table 4. Numeric values for the non-logarithmic coefficients for the various SCET and bHQET coefficients a n,n−1 . u = 1, 2 and an infinite number of simple poles at integer values of u ≥ 3: (7.9) These double poles unambiguously signal the presence of anomalous dimensions with n f dependence at leading order for dimension-2 and -4 operators of the associated OPE. 20 The total ambiguity can be expressed in the following compact form: where in the first line we have explicitly singled out the first two leading contributions, while the second line starts at O(Λ 6 QCD ). The logarithmic enhancement of the two leading terms appears due to the double nature of the poles and is associated to the anomalous dimension just discussed. This entails that Q 2 δ Λ H Q retains some logarithmic dependence on Q in Fig. 10(a) the dependence of δ Λ H Q in Q is shown in red, while the ambiguity when the logarithm is set to zero appears in blue. Despite the logarithmic enhancement, it is negligibly small even for the smallest values of Q for which SCET is applicable.
We compare next the exact result for the hard function H Q (Q, µ) obtained with Eq. (4.31) with its fixed-order expansion. For the latter we compute H Q (Q, µ 0 ) and subsequently run to H Q (Q, µ). A natural choice for the matching scale that avoids explicit large logarithms is µ 0 = Q, but in this way H Q (Q, Q) has analytical continuation π n terms 20 Operators corresponding to simple poles might also carry anomalous dimension, but at leading-order they should not depend on n f . We thank A. Hoang for clarifying this point. Figure 11. Feynman diagrams for the one-loop jet function for massless quarks with off-shell momentum p 2 = s and virtual loop momentum . The cross represents a quark jet-field insertion. Diagram (a) has to be multiplied by a factor of two to account for the symmetric contribution with the gluon line emitted from the right current.
that can be traced as coming from the real part of log n (−µ 2 /Q 2 + i 0 + ) with Q 2 and µ 2 both positive when transitioning from Euclidean to Minkowskian regions, which therefore can also be summed up using RG evolution. In the renormalon formalism these factors come from expanding in powers of u the factor cos(πu) appearing in G H (0, u). To carry out this procedure one performs resummation in C H (Q 2 , µ) by computing first C H (Q 2 , µ 0 ) with µ 0 = −iQ and then running to from µ 0 to µ. From this resummed result for C H (Q 2 , µ) one obtains H Q (Q, µ) simply taking the modulus squared expanded to O(1/β 0 ). We compare these two approaches in Fig. 10(b) for a relatively small value of the center-of-mass energy. The red dots, which include π-resummation, converge to the exact answer much faster than the blue, which do not. Eventually both approaches reproduce the exact result, and start diverging for n > 30. The ambiguity of H Q (Q, µ) is too small to be visible in the plot. Since the hard factor mainly affects the norm of the event-shape distributions, this behavior might explain why the total-hadronic cross section obtained via a direct integration of the differential event-shape cross section, was undershot at low orders in the thrust and C-parameter distributions in Refs. [72,73]. Similarly, Ref. [38] found poor order-by-order convergence in peak differential cross sections for boosted tops unless they were self-normalized.

Jet function
The SCET jet function describes energetic radiation branching from the initiating quark, and depends on the collinear measurement function, which is the same for hemisphere masses, C-parameter and thrust. Furthermore, it is completely inclusive such that the jet function can be computed as the discontinuity of a forward-scattering matrix element: with d = 4 − 2ε, Tr a trace over spin and color, and p 2 = s. We use the gauge-invariant jet field χ n = W † n ξ n , which is the product of a collinear field ξ n and a path-ordered collinear Wilson line W † n . Likewise, χ n,Q is a jet field with total minus label momentanP set to Q with a Dirac delta function. This overly simplifies the computation of the leading 1/β 0 correction as one only needs to modify the gluon propagator of the one-loop diagrams shown in Fig. 11. The vertex diagram 11(a) must be doubled to account for the symmetric contribution, which is omitted along with the diagram in which the gluon is emitted and absorbed from both Wilson since it vanishes in the Feynman gauge. Given that the jet functions for the three event shapes just mentioned are trivially related to one another, in this section we carry out the computation for the invariant mass of a single hemisphere. The discontinuity of J n (s, µ) [ that is, J n (s+i0 + , µ)−J n (s+i0 − , µ) ] yields the momentumspace jet function J n (s, µ), with support for s > 0. It can be converted to the position-space jet functionJ n (y, µ), with y the variable conjugate to s, defined as: While J n (s, µ) is a regular function, J n (y, µ) for massless quarks contains Dirac delta and plus distributions. For its computation we will use the following identity: Both J n (s, µ) and J n (y, µ) have dimensions of an inverse squared mass, while s is a realvalued mass-dimension-two variable. On the contrary,J n (y, µ) is a dimensionless complexvalued regular function depending on the complex variable y, with mass dimension −2. At tree-level one has J tree n (s) = − 1 2π 1 s + i0 + , J tree n (s) = δ(s) ,J tree n (y) = 1 , (7.14) and we define the corrections to the bare jet function in momentum and position space as J n (s, µ) = δ(s) + δJ n (s) andJ n (y) = 1 + δJ n (y), respectively, which are computed next. The vertex contribution shown in Fig. 11(a), and the self-energy correction of Fig. 11(b) contribute to the forward-scattering matrix element as follows Adding the two contributions and using Eq. (7.13) with η = h + ε we obtain 21 Let us compute the anomalous dimensions using the position-space jet function for which we can identify Q 2 = −ie −γ E /y and b(ε, h), to find 22 Taking into account the relation ΓJ = 4Γ cusp we find that βGJ (−β, 0)/β 0 reproduces the result found in Eq. (7.5). To compute the jet non-cusp anomalous dimension we need also where the second equality has partially expanded in powers of β. This expression can be used in Eq. (4.16) together with GJ (−β, 0) = −2G Q (−β, 0) to obtain a closed form for γ J (β), as well as its fixed-order coefficients. We find full agreement with the leading flavor structure of the known full SCET results up to O(α 3 s ) [61], see columns 3 and 4 in Table 3. The soft non-cusp anomalous dimension can be computed from a consistency condition and it reads γ S = −γ H − γ J . In both cases the convergence radius is ∆β = 5/2. In Figs. 9(c) and 9(d) we compare for α s = 0.9 the exact and fixed-order results for the jet and soft noncusp anomalous dimensions, respectively. The soft function has γ 0 S = 0 and therefore the partial sum starting from the first non-vanishing order is shown. We find good convergence in both cases. 21 The (dimensionless) cumulative jet function, defined as ΣJ (sc, µ) ≡ sc 0 ds J(s, µ), has the same cusp and non-cusp anomalous dimension asJn(s, µ). Furthermore, the Borel transform of the series has no poles for u > 0, resulting in a finite convergence radius. Identifying Q 2 = sc one can write down a closed form for ΣJ using GΣ J (0, u) = FJ (u)/Γ(1 + u). 22 The momentum-space jet function for s > behaves as a series without cusp anomalous dimension [ that is, the associated F (ε, u) function is regular at u = 0 ]. Paradoxically, the (non-cusp) anomalous dimension of Jn(s) for s > 0 is precisely −4Γcusp. The series is not ambiguous, and identifying Q 2 = s one can write down a closed form for sJn(s, µ) using Eq. (3.21) and noting that FJ (0, u) = −GΣ J (0, u). To compute the non-logarithmic fixed-order coefficients a J i,i−1,0 of the Fourier-space jet function, which are defined as where the factor 2 −j is necessary to compensate the µ 2 that appears in the argument of the logarithm, we need one last function (7.20) Using its Taylor expansion in Eq. (4.6) we fully reproduce the leading flavor structure of the known coefficients in full SCET up to O(α 3 s ) [74][75][76][77][78], see second column of Table 4. The functionJ B (u) ≡ GJ (0, u) − (GJ ) 0,0 − u(GJ ) 0,1 /u 2 has only two simple poles at u = 1, 2, with the following asymptotic expansion: Therefore the ambiguities in position and momentum space are  Figure 13. SCET Feynman diagrams for the computation of the bHQET to SCET matching coefficient at one-loop for a massive on-shell quark with momentum p 2 = m 2 and virtual loop momentum . We do not show the two symmetric diagrams, as they give identical results, or the contribution in which a soft gluon is interchanged between the two massive collinear quarks, which vanishes, along with all bHQET diagrams.
the jet function is free from the u = 1(2) renormalon, but affected by the u = 2(1) ambiguity. For complex y both real and imaginary parts are affected by the two singularities. In the second line of Eq. (7.22) it is shown that when transforming back to momentum space the leading ambiguity is proportional to δ (s). Therefore, to find an ambiguous series in momentum space one needs to consider a double cumulative (that is, the cumulative version of Σ J ) or a weighted cumulative such as s 0 ds s J n (s, µ). In Fig. 12 we compare the fixed-order expansion for the position-space jet function up to (n+1)-loops with the exact result, splitting in two separate panels the real and imaginary parts, for a real value of the Fourier variable y. The ambiguity is visible for the imaginary part only. We nevertheless observe nice convergence for different values of the matching scale.

bHQET computations
If the quarks primarily produced in the e + e − collision at high energy are massive (e.g. the top quark) a new scale -the quark mass m -becomes relevant in the problem. This mass changes the SCET jet function [ secondary production of heavy quarks through gluon splitting affects also the soft function starting at O(α 2 s ) ]. In the peak of the distribution the jet invariant mass and m are very similar, indicating that quarks are very boosted. In this regime one can match SCET with massive quarks to bHQET, what allows to sum up a new class of large logarithms and treat the top quark decay products in an inclusive way through the following factorization theorem [22,79]: The difference between Eqs. (8.1) and (7.1) is that the former has an additional (universal) hard factor H m that accounts for the matching between SCET and bHQET, and the SCET jet function J n has been replaced by its bHQET counterpart B n (which takes different forms depending on the observable considered, but we focus once more in the hemisphere mass). In this section we compute these two functions at leading order in 1/β 0 .

Hard massive factor
To determine the hard massive function H m one needs to compute the matching coefficient C m between the SCET and bHQET massive di-jet operators. The optimal way of carrying out the computation is using the on-shell scheme and regularizing all infrared divergences within dimensional regularization, making all bHQET diagrams vanish such that only SCET contributions need to be considered. The diagrams to take into account are shown in Fig. 13, where both contributions need to be doubled to account for the symmetric diagrams which are not explicitly shown. Furthermore, the contribution in which a soft gluon is exchanged between the two massive quark lines also vanishes and has been omitted. For diagram 13(b) we use our previous result on the wave-function renormalization Z OS ξ for a massive on-shell quark in Eq. (6.7), since it is the same in SCET and full QCD.
At lowest one has C m = 1. Defining the corrections to the bare matching coefficient as C m = 1 + δC m (m), and using the massive SCET Feynman rules that can be found in Ref. [80], the computation of diagram 13(a) using a modified gluon propagator yields Doubling the result above and including the wave-function renormalization we obtain the modified 1-loop correction necessary to determine the leading 1/β 0 term: Identifying Q = m we obtain the following form for the renormalon master function G(ε, u): For the H m factor one gets an additional factor of two: G Hm = 2G Cm . Equation (8.4) reproduces the cusp anomalous dimension noting that Γ Cm = 2Γ cusp .
To determine the C m matching coefficient, on top of computing the relevant Feynman diagrams and accounting for the wave-function renormalization, one needs to renormalize the SCET current multiplying the bare result with Z C H as computed in Sec. 7.1. At leading order in 1/β 0 this amounts to adding δZ C H and therefore the renormalized C m at this order is obtained simply by dropping the divergent terms, denoted by δZ The second equality is expressed as a Taylor polynomial times the exponential of another series, in a form amenable to be re-expanded by a computer program with the algorithms already explained. This result can be used to obtain a closed form for γ Hm and its fixed-order coefficients, whose series has convergence radius ∆β = 5/2. We find full agreement with the leading flavor structure of the known full bHQET results up to O(α 3 s ) [37,79,81], collected in the fifth column of Table 3. In Fig. 14(a) we compare the exact form of the bHQET non-cusp anomalous dimension for α s = 0.9 with its fixed-order partial sum including up to 23 The SCET current renormalization does not affect the direct determination of the cusp anomalous dimension from Eq. (8.4). γ n Hm . Even for this enormous value of the strong coupling, the expanded result converges to the exact answer.
To obtain the C m matching coefficient at leading order in 1/β 0 we need the following function, which enters the inverse-Borel transform integral: where again we provide a partially expanded expression. With the result above we reproduce the leading flavor structure of the full bHQET coefficients up to two loops [79,81], see column 3 in Table 4. The Borel transform C m B (u) ≡ G Cm (0, u) − (G Cm ) 0,0 − u(G Cm ) 0,1 /u 2 has only simple poles at half-integer values of u, including an O(Λ QCD ) renormalon. Its asymptotic expansion can be written in terms of a finite and an infinite sum The leading ambiguity for H m , which includes a factor of two, is three times as large as the pole mass ambiguity with the same sign: Therefore, the combination H m /m 3 p is free from the leading ambiguity. To illustrate how this renormalon cancellation works in practice we refer to Fig. 15, where the dimensionless quantity (m t /m pole t ) 3 H m (m t , µ) is considered, showing its exact value as a gray dashed line in both panels. In these plots we highlight some important aspects of consistently expanding perturbative series to achieve an effective renormalon cancellation. For the perturbative expansions we use µ H = m t as the matching scale, and therefore for a correct cancellation one needs to expand (m t /m pole t ) 3 in powers of α s (µ H ). In Fig. 15(a) we either keep m p exact expanding in powers of α s (µ H ) only H m (m t , µ H ) (green dots); keep H m (m t , µ) exact expanding only (m t /m pole t ) 3 (blue dots); or expand both (red dots), and only in the last case we see a better-behaved perturbative expansion (keeping both factors exact generates the gray line). If Fig. 15(b) we expand H m and (m t /m pole t ) 3 in powers of α s (µ H ) and α s (µ m ), respectively, with µ m = ξ × µ H and ξ = 4, 1/4, 1 shown as green, blue, and red dots, respectively. As expected, only when ξ = 1 the asymptotic behavior is removed. In both panels we observe that the removal of the leading renormalon ambiguity makes the series approach the exact result after fewer orders have been included in the partial sum. The H m hard factor mainly affects the norm of the distribution. In Ref. [38] this renormalon was not properly accounted for, and this fact, together with the lack of π summation in H Q , might explain the poor convergence found in peak cross sections for boosted top pairs unless the curves were self-normalized.

Jet function
The last matrix element we consider in this article is the bHQET jet function for hemisphere masses. The relevant Feynman diagrams coincide with those of the SCET jet function shown in Fig. 11, but the actual computation uses bHQET Feynman rules in which the heavyquark momentum p = mv + k with v 2 = 1 has a residual component k. We follow the same computational strategy as for the determination of the SCET jet function in Sec. 7.2 based on the following forward-scattering matrix element: with h v an HQET massive quark field,ŝ = 2v · k and W ( †) n the collinear Wilson line already introduced in Sec. 7.2. Its discontinuity along the branch cutŝ > 0 defines the momentumspace jet function B n (ŝ, µ), with support forŝ > 0. Both functions dimensions of an inverse squared mass, although the kinematic variableŝ ≡ (s − m 2 )/m they depend on has dimensions of energy. Since the fixed-order expansion of B n (ŝ, µ) contains Dirac delta and plus distributions [ even though B n (ŝ, µ) is a regular function ], it is convenient to consider its Fourier transformB which due to the prefactor m is dimensionless and depends on x, the variable conjugate tô s with dimensions of an inverse energy. At lowest order one has m B tree n (ŝ) = 1 2π î s + i0 + , mB tree n (ŝ) = δ(ŝ) ,B tree n (y) = 1 . Quantum corrections to the bare jet function, defined as B n (ŝ) = δ(ŝ) + δB n (ŝ) and B n (x) = 1 + δB n (x), are computed next at one-loop with a modified gluon propagator.
The diagram shown in Fig. 11(a) can be expressed as the vertex correction times an HQET propagator: mB a n (ŝ) = 1 2π On the other hand, the self-energy diagram in Fig. 11(b) can be written in terms of the bHQET self-energy M b and a squared HQET propagator Combining the two diagrams with the appropriate factor of two for the vertex correction, and taking the discontinuity using Eq. (7.13) with η = 2(ε + h), one obtains the modified one-loop momentum-space bHQET jet function, which can be converted to position space: .
To obtain the cusp and non-cusp anomalous dimensions for the bHQET jet function we proceed in position-space. To that end, we identify Q = −ie −γ E /x and the 1-loop coefficient b(ε, h), and obtain the following renormalon master function Taking into account the relation ΓB = 2Γ cusp we find that 2βGB(−β, 0)/β 0 reproduces once more the cusp anomalous dimension in Eq. (7.5). To obtain the non-cusp anomalous dimension we need the u derivative of GB(ε, u) at u = 0: where second equality is expressed as a series times the exponential of another series, ideal for a full re-expansion in powers of β. This expression can be used in Eq. (4.16) together with GB(−β, 0) = −G Q (−β, 0) to obtain a closed form for γ B (β), as well as the fixed-order coefficients, with convergence radius ∆β = 5/2. We find full agreement with the leading flavor structure of full bHQET results up to O(α 3 s ) [37,79,81], see last column of Table 3. From the bHQET consistency condition one can deduce the following relation between various non-cusp anomalous dimensions: γ Hm = −2(γ S + γ B ) = 2(γ H + γ J − γ B ), where in the second equality we use the SCET consistency condition to eliminate γ S . This relation is exactly verified by our computations at O(1/β 0 ). In Fig. 14(a) we compare the exact result for γ B with the fixed-order partial sum including up to γ n B . We use α s = 0.9, but nevertheless find that after adding 6 or more terms the perturbative result nicely agrees with the exact computation.
To compute the non-logarithmic fixed-order coefficients a B i,i−1,0 of the Fourier-space bHQET jet function, which we define as we need to evaluate the renormalon master function at β = 0: again provided in a closed form and as the exponential of a series. With this result we correctly reproduce the leading flavor structure of the full theory up to two loops [37,79], see last column of Table 4. The Borel transformB B (u) ≡ GB(0, u)−(GB) 0,0 −u(GB) 0,1 /u 2 has simple poles at u = 1/2 (related to the pole mass ambiguity) and integer values of u ≥ 2.
Its asymptotic expansion is therefore expressed as an isolated term plus an infinite sum: The full ambiguity can be written in terms of Bessel functions of the first kind and generalized hypergeometric functions. The leading ambiguity is, as expected, proportional to Λ QCD and takes the following form in position and momentum space: In position space, except for the factor ix, the ambiguity is exactly twice that of δ MS (m).
If one writes the position-space version of the factorization theorem shown in Eq. (8.1), the bHQET jet function effectively appears in the combinationB n (x, m p , µ) ≡B n (x, µ)e −2ixmp , which is free from the leading renormalon ambiguity. Therefore, in order to have a stable perturbative expansion forB n one needs to express m p in a short-distance scheme and expand consistently in powers of α s (µ) with µ 1/x ŝ. Since the heavy quark mass is no longer a dynamic scale of the effective theory, one should use a low-scale short-distance scheme whose relation to the pole mass is not proportional to the mass itself, such as m MSR . This avoids having renormalon-subtractions that are much larger than the corresponding fixed-order corrections forB n , thus breaking the bHQET power counting. Furthermore, if one chooses R ∼ µ there will be no large logarithms in the renormalon subtraction series. In Fig. 16 we show how the leading renormalon effectively cancels for the case of the top quark, both for the real and imaginary parts in panels 16(a) and 16(b), respectively. Since the real part ofB n does not have the u = 1/2 ambiguity if x is real, we choose a complex number, but for simplicity keep the renormalization scale µ real. When truncating the fixed-order series forB n and δ MSR at any finite order one gets residual dependence on R and µ. In the figures we use n + 1 terms in the fixed-order partial sum, plus N n LL resummation both iñ B n and m MSR t . For this numerical analysis we choose µ 0 = R . If one expands in powers of α s (µ 0 ) onlyB n (green dots) or m p (blue dots) the asymptotic behavior remains in the series, and it takes more orders to approach the exact value (shown as a gray dashed line). The way in which the series diverges in each case seems to be almost exactly opposite. When both functions are expanded consistently (red dots), the series appears much better behaved for large n, and converges to the exact result already at three loops.

Jet masses
As we have seen in Sec. 8.2, in the large-β 0 approximation the leading renormalon of the position-space bHQET jet function is related to that of the pole mass. This relation holds in full QCD as well and therefore one can useB n to define a low-scale short-distance scheme. Using the exponentiation properties ofB n , in Ref. [37] the so-called jet mass, dependent on µ and an infrared scale R analogous to that of the MSR mass, was defined as follows This idea was adapted in Ref. [82] to define a subtraction scheme for the hemisphere-mass soft-function renormalon. In Ref. [38] it was noted that since such subtractions involve a derivative, their asymptotic behavior starts at two loops. Moreover, for the canonical choice µ = R the one-loop soft subtraction vanishes. To solve this shortcoming, a new scheme for the soft function not based on derivatives was introduced in Ref. [38]. In this section we define a new jet mass following this same spirit. We study both schemes in the large-β 0 and derive closed forms for their relation to the pole mass and R-anomalous dimensions.

Derivative scheme
In the large-β 0 one has log(1+δB n ) δB n and the derivative with respect to log(ix) makes the perturbative series associated to m J (R, µ)/(Re γ E ) non-cusp type. Indeed, applying Eq. (9.1) to our result in the second line of Eq. (8.14), the modified one-loop function associated to the jet mass reads From this result, taking Q = R, the following renormalon master function is found which is finite for u → 0. From the equation above we can infer that the µ-anomalous dimension of the jet mass is proportional to the universal cusp anomalous dimension. The evolution in µ for constant R is therefore given by where in the second line the solution at O(1/β 0 ) is provided in terms of β i = β(µ i ). To obtain the fixed-order coefficients of δm J (R, µ) one needs from which one can reproduce the leading flavor structure of the known coefficients for δm J in the full theory up to O(α 2 s ) [37]. We can easily find the asymptotic expansion of the corresponding Borel-transform function m J The leading ambiguity, coming from the u = 1/2 pole, is independent of R and fully coincides with that of δ MS (m) in Eq. (6.10), making δm J (R, µ) − δ MS (m) free from the O(Λ QCD ) renormalon. This renormalon-free difference for R = µ = m can be used as a matching condition between the jet and MS masses. The R-anomalous dimension of the jet mass is defined along the diagonal path µ = R: and in the large-β 0 limit takes the same closed form as Eq. (6.15) with the obvious substitution F MS → e γ E F J m , such that γ J R does not have a u = 1/2 ambiguity. Upon expansion in powers of β R one obtains the fixed-order coefficientsγ n R,J . In Fig. 17(a) we compare the exact value of γ J R (red dashed line), whose ambiguity is too small to be visible on the plot's scale, with the fixed-order partial sum including up to γ n R,J (blue dots) for a relatively large value of the strong coupling α s = 0.25. The divergent behavior does not manifest for this value within the first 17 orders shown in the plot.
The solution to this RGE is the same as that explained already in Sec. 6.3. In fixedorder perturbation theory, starting with the value m J (m, m) = m+ δ MS (m)−δm J (m, m) , one can obtain the jet mass at the scales µ and R using the solutions to the corresponding RG equations. Defining ∆ R J (R 2 , R 1 ) as in Eq. (6.17) with the replacement γ n R → γ n R,J one proceeds as follows: In Fig. 5(b) we compare the exact value of m J (R, µ), shown as a gray dashed line, with the corresponding computation using the fixed-order partial sum for the matching condition including n terms and R-evolution at N n LL. We observe that the perturbative result agrees with the exact value obtained with the PV prescription for 3 ≤ n ≤ 10, and starts diverging for n > 10.

Non-derivative scheme
In this section we consider a low-scale short-distance mass, denoted by m J (R), defined form the bHQET jet function which does not depend explicitly on µ. Based on the prescription introduced in Ref. [38], it is defined as displaying an explicit dependence on the scale R. One of the motivations for the scheme definition in Eq. (9.1) was to have a µ-dependent low-scale short-distance mass. The rationale behind this choice is that, as we have seen, to properly cancel the renormalon of e.g. the bHQET jet function one needs to express δm J (R, µ) in powers of α s (µ), and at first glance it might appear that if the short-distance mass does not depend explicitly on µ this can only be achieved if we take R = µ. However, as shown e.g. in Ref. [38] for soft gap subtractions, or in Ref. [83,84] with the MSR mass, one can still use R = µ as long as δm J (R) is expanded in powers of α s (µ), what implies that the fixed-order coefficients depend logarithmically on µ/R. As long as µ ∼ R these logarithms are not large. From the definition in Eq. (9.9), linearizing again the logarithm at leading order in 1/β 0 and identifying Q = R it is trivial to realize that the master renormalon function for the combination 2δm J /(Re γ E ) coincides with GB shown in Eq. (8.15). From the result in Eq. (8.20) with the replacement ix → 1/(Re γ E ) one infers that the leading ambiguity for δm J (R) exactly matches that of δ MS (m) shown in Eq. (6.10). Therefore, the difference δm J (R) − δ MS (m) is free from the O(Λ QCD ) renormalon and for R = m can be used as the matching condition to relate m J (m) and m.
The corresponding R-anomalous dimension is computed as the R derivative acting on m J (R), and defining h(β) = −2β 0 Γ cusp (−β)/β and b(ε) = d[GB(ε, u) − GB(0, u)]/du u=0 we find the following closed form Once more, the factor (1 − 2u) removes the leading ambiguity from γ J R . To compute the derivative with respect to β of the last term in Eq. (4.13) a regulator has been used such that the logarithm does not diverge at the lower integration limit. After applying the derivative the regulator can be safely set to zero: To obtain the last line we write log(δ) = log(1 + δ) − 1/β 0 −β dε/(1 + ε/β + δ) and set δ = 0 right away. The terms proportional to h (0) are not necessary to make the integrand convergent at the upper integration limit. Furthermore, they cancel one other and therefore do not appear in last line of Eq. (9.10). Equation (9.10) can be integrated analytically in R provided one switches the order of the various double integrals that appear, recovering the difference of jet masses at two values of R. This is an important cross check of our computation. In Fig. 17(b) we perform a comparison of the fixed-order partial sum and exact results for γ J R . The description and observations are identical to those already presented for Fig. 17(a).
Defining ∆ R J (R 2 , R 1 ) in a way analogous to ∆ R J one can obtain the jet mass in fixedorder perturbation theory at any scale R simply using R-evolution and the matching relation: m J (R) = m + δ MS (m) − δm J (m) + ∆ R J (R, m). In Fig. 5(b) we compare the exact value of m J t (R), shown as a red dashed line, with the fixed-order partial sum, in green. Explanations and remarks are identical to those made for the derivative-scheme jet mass.
either scheme. This might have no effect in practical applications since the value of α s rarely becomes so large and at most four perturbative orders are known in full QCD. In any case, this worse perturbative behavior of γ R seems to affect the stability of perturbative R-evolution when used to determine the MSR and jet masses, as can be seen in Fig. 5(b): if one evolves the various short-distance masses to a very small scale R, the asymptotic behavior for the MSR mass sets in after 7 orders, while jet masses remain stable until the 12th loop. On the other hand, from Fig. 18(a) it is hard to conclude that either of them is better suited to estimate the pole mass: starting from the second loop all schemes go hand in hand. In particular, the MSR (red) and derivative-scheme jet (green) masses are almost on top of each other, while the non-derivative jet mass is closer to the exact value (gray dashed line with a fade gray band) at one loop. Likewise, one can see in Fig. 18(b) (with the same color convention) how the three short-scale low-scale masses successfully remove the leading renormalon in the real part of the position-space bHQET jet function (the imaginary part is not shown as it does not add any substantial additional information), which uses the same values as in Fig. 16. Even though the MSR mass appears to reach the exact value earlier, after the second loop all schemes look almost identical.

Estimate of bHQET jet function at three-loops
Even though the bHQET jet function has been computed to two-loop accuracy, its cusp and non-cusp anomalous dimensions are known up to O(α 3 s ), which determine the logarithmic coefficients at this order. SinceB n depends on a single scale only, the missing nonlogarithmic term is a real number. In the previous section we have determined the leading flavor structure at three-loops, α 3 s n 2 while in Ref. [37] it was shown that the position-space jet function obeys non-abelian exponentiation [85,86]. In particular this encompasses that in the abelian limit (obtained taking {n , C A } → 0) the function log(m pBn ) is one-loop exact, unambiguously fixing the C 3 F color structure. The pieces which are left to determine are thus C 2 F n , C F C A n , C F C 2 A and C 2 F C A , although in practice one can simply set C A = 3 and C F = 4/3 and estimate the linear and n -independent terms. The exponentiation theorem and the fact that log[B n (x, µ)]−2ixm p is free from the leading u = 1/2 renormalon can be exploited to estimate the non-logarithmic term. To that end we use a modified version of Eq. (4.21) in Ref. [35] applied to the jet mass in the non-derivative scheme defined in the previous section. We start by writing the corresponding perturbative series as which is formally independent of λ. The coefficientsb i ≡b i (λ = 1) correspond to the non-logarithmic fixed-order terms of log[B n (x, µ)]. The λ-dependent coefficientsb i (λ) are computed asb withĉ i,0 ≡b i . From theb i (λ) coefficients one computes λ-dependent S k coefficients using e.g. Eq. (A.18) of Ref. [35]. In our case only S λ 1 and S λ 2 can be determined and accordingly the upper limit on the first sum must be set to k = 2. In order to study the convergence of our estimates we will consider setting the upper limit also to k = 1, which is useful to check if the known two-loop result is reproduced, within uncertainties, if only S 1 is used as input. All in all, the m-loop estimate of the coefficientb n (λ) takes the following form: where the g coefficients can be computed with the recursive relation provided in Eq. (A.14) of Ref. [35]. The benefit of using this formula -rather than the more conventional (4.20) therein -is that it incorporates information beyond leading-renormalon dominance. In particular, assuming that m ≤ 4, all coefficients with n ≤ m are reproduced exactly. From b m n (λ) one can obtain an estimate for the canonicalb m n (λ = 1) ≡b m n coefficients simply inverting Eq. (10.2) [ that is, derivingĉ m i,j (λ) coefficients from the recursive expression in Eq. (10.2) and using the first relation with the replacement log(λ) → − log(λ) ]. Since the sum over S λ k in Eq. (10.3) is convergent, truncating it at k = m (rather than the exact value n − 1) inflicts a small uncertainty that can be quantified directly onb j n varying the parameter λ in the range (1/2, 2).
Let us present our numerical analysis. In first place we fix the active number of flavors to the physically relevant cases for charm, bottom and top n = 3, 4, 5 and directly estimatê b 3 (n ) using 1-and 2-loop information. We find that the 2-loop results are contained in the 1-loop uncertainty bands. As an additional check we estimateb 2 (n ) using 1-loop information and find that the known exact values are contained within the estimated uncertainties. Our results are summarized in Table 5. Next we try to estimate the coefficients of the various powers of n , which we denote by 24 b i (n ) =  To that end we determine the value ofb 3 (n ) in the range −10 ≤ n ≤ 6, and given thatb 3,2 has been computed in Sec. 8.2, least-square minimization is used to determine the remaininĝ b 3,0 andb 3,1 . We have checked that varying the n fit-range boundaries by one (upper limit) or several (lower limit) units changes the estimated values by an amount much smaller than the uncertainty caused by λ variation. Therefore even if included in our final analysis, we will not discuss this source of error any further. Considering the individualb 3 (n ) errors (obtained by varying λ) as uncorrelated greatly underestimates the uncertainties of the coefficients. This is expected since there is of course a strong correlation between the individual uncertainties. One can take the extreme situation in which the errors are considered as 100% correlated by taking the same value of λ for all n values, determining thenb 3,0 andb 3,1 as a function of λ. By varying λ in the range specified above one can estimate the uncertainties of each coefficient and the correlation between them. We find a very strong negative correlation, very close to 100%, which is responsible for anb 3 (n ) uncertainty in the range 3 ≤ n ≤ 5 smaller than that of the direct determination. This is caused by our rough assumption of a 100% correlation in λ. To correct this deficiency we keep the individual uncertainties (to which we add in quadrature those coming from the linear regression) ofb 3,0 andb 3,1 but rescale the correlation between them by a factor that makes the direct and indirect uncertainties match for the physically relevant values of n , see Fig. 19(a). This rescaling also makes the 1-and 2-loop estimates ofb 3,0 andb 3,1 compatible, as shown in Fig. 19(b). Our numerical results are presented in Table 6, where we also shown that using 1-loop input the two-loop n -independent coefficient is correctly predicted within its uncertainty.
Carrying out the same exercise for the thrust and C-parameter soft function, under the assumption of u = 1/2 renormalon dominance one obtains less precise results, particularly orderb 3,0b3,1 correlationb 3,2b2,0b2,1  Table 6. Numerical values, uncertainties and correlation for the sub-leading flavor structures of the position-space bHQET jet function's exponent coefficients at 3-loops using at 1-and 2-loop input. For completeness we also shown the known values forb 3,2 andb 2,1 , and the estimated two-loop coefficientb 2,0 if 1-loop input is used.
so for the former. Since the leading flavor coefficient is not known, we do not attempt to determine the individual n f coefficients. Following the notation in Refs. [73,88] we find The thrust result for n f = 5 is however not compatible with the determination carried out in Ref. [77], s τ, n =5 3 = −10000 ± 2100, which in turn is based on the numerical analysis of Ref. [89] which uses NNLO binned distributions generated with EERAD3 [90].

Conclusions
In this article we have adapted the methodology to compute the leading term in the 1/β 0 expansion for series with cusp anomalous dimension. Such series are ubiquitous in effective field theories for massive and massless jets such as SCET and bHQET. We have found closed expressions for the renormalized series, its multiplicative renormalization factor as well as their cusp and non-cusp anomalous dimensions. Furthermore, we have solved the corresponding RGE equations in a closed form in this limit, and shown that the ambiguity does not depend on the renormalization scale. We have provided an optimal algorithm for the principal value prescription used to compute the inverse Borel integral when there are poles in the integration path.
As a first application, the series that relates the pole and MS masses has been revisited, recovering a known exact from for the MS mass anomalous dimension. We have derived, for the first time, a closed expressions for the MSR mass R-anomalous dimension γ R in the large-β 0 limit, and provided an analytic solution for the corresponding RG equation at any perturbative order. We have investigated how the fixed-order expansion of the series defining this scheme can be used to estimate a numerical value for the pole mass, finding that for smaller R the series quickly converges to the exact value. On the other hand, the asymptotic behavior of the series only depends on the renormalization scale µ at which the strong coupling is evaluated. The RS mass has also been worked out in this limit, and we have demonstrated that the difference of two such masses for different values of R obeys a perturbative R-evolution equation. Finally, it has been shown that the renormalon normalization constant N 1/2 as computed with the sum rule derived in Ref. [35] tends to the large-β 0 result, which can be obtained from our analytic expressions.
The hard and jet soft functions for massless (in SCET) and massive (in bHQET) quarks have been computed in the large-β 0 with a modified gluon propagator, applying the technology developed in this article. Our determination of the universal cusp anomalous dimensions agrees with a previous computation, and our results for the non-cusp anomalous dimension and matrix elements agree with the leading flavor structure of known results computed in the full theory. Our results for SCET reveal that the leading ambiguities happen at u = 1, so they are proportional to Λ 2 QCD . We find that the Borel transform of the SCET hard function has double poles, associated to operators with non-trivial anomalous dimension in the OPE which depend on n f at leading order. We also find that its fixedorder expansion converges faster if π n terms are summed up choosing a complex matching scale for the Wilson coefficient C H . The SCET jet function, on the other hand, has only two simple poles. Our computation of the hard and jet functions in bHQET shows that in both cases the leading renormalon is of the most severe type, that is, proportional to Λ QCD . We have seen that the fixed-order expansion of the renormalon-free quantity H m /m 3 p shows nice convergence if the pole mass is expressed in terms of the MS mass. On the other hand, the fixed-order expansion for the combinationB n (x, µ)e −im 2 x appearing in the bHQET factorization theorem, which is also free from the leading renormalon, behaves well if the pole mass is expressed in terms of the MSR mass (or other R-dependent low-scale shortdistance masses). We have estimated the unknown three-loop non-logarithmic coefficient ofB n in full QCD with 20% accuracy using renormalon dominance.
The ambiguities found in the EFT computations can have two physical meanings, depending if Λ QCD is interpreted as the scale of non-perturbative (hadronic) physics Λ had or as the soft scale µ s of SCET or bHQET. For the hard function H Q , since it comes from matching two operators and its definition does not involve a matrix element (even though, in practice, for its computation matrix elements are taken) it appears natural to interpret Λ QCD as the soft scale, and since µ s /µ H ∝ λ 2 with λ the SCET power-counting parameter, ambiguities are likely related to sub-leading operators in the EFT power expansion. The jet function is indeed a matrix element and in this case Λ QCD could signal either powersuppressed jet functions or effects from collinear hadronization corrections. In either case, given the relation µ s /µ J ∝ λ, these are power suppressed contributions. For the H m hard massive factor and the bHQET jet function, the leading u = 1/2 ambiguity is obviously related to the well-understood pole mass renormalon. Higher-order ambiguities can be interpreted in the same way as just done for H Q and J n , respectively.
From the position-space bHQET hemisphere jet function one can define the so called jet mass, which involves a derivative of the logarithm ofB n . We have provided an alternative scheme for a µ-independent jet mass in which no derivative is taken. We have worked out closed expressions for the relation of these two schemes to the pole mass, showing they have the same leading ambiguity as the corresponding MS series, making possible an unambiguous matching to m. Closed expressions for their R-anomalous dimensions have also been derived, and it has been shown that the two schemes work similarly to the MSR mass when it comes to removing the renormalon ofB n or estimating a numerical value for the pole mass in fixed-order perturbation theory.
Our results can be used to derive an exact prediction at leading order in 1/β 0 for the singular structure of massless event-shape distributions in the di-jet limit, as well as massive event shapes in the peak region, once the large-β 0 soft function is known. Such computation is more involved than those carried out in this article as it cannot be cast as the discontinuity of a forward-scattering matrix element. Instead, one has to deal with two-body real-radiation phase-space integrals (when "cutting" the quark bubble), which are similar in complexity to a two-loop integration, with the additional complication of modified gluon propagators. This computation is therefore left to future work. Our results also suggest that studying the effect of π summation in the SCET hard function for a full QCD event-shape distribution is warranted, and we will do so in the near future. Likewise, it appears important to figure out a practical implementation of the bHQET hard factor renormalon subtraction within the corresponding factorization theorem, and how this is linked to reparametrization invariance [91], a transformation that connects different orders in the 1/m expansion. where in the first step the definition of binomial is used, in the second we shift n → n − 1, pull out the largest factor in (i − 1)! = (i − 1)(i − 2)!, and combine it with the rest of factorials as a new binomial. To obtain the first equality in the second line we use the binomial expansion for (n + 1) m−1 . For i = 1 the last equality automatically vanishes due to the (i − 1) factor multiplying the sum. Since for i ≥ 2 one has J i−1,0 = 0 it is implied that J i≥2,1 = 0; likewise, given that for i ≥ 32 both J i−2 = J i−1 = 0 it immediately follows that J i≥3,2 = 0. Proceeding in this fashion one concludes that J i≥m+1,m = 0 for any non-negative m, and this proofs I i,m = 0 for 1 ≤ m ≤ i − 1.
In the second part of this appendix we provide a proof for Eq. (3.5), making use again of recursion relations. Noting that K 1 = −1, assuming i ≥ 2 and shifting n → n + 1 in Eq. (3.5) one gets where in the second equality we have used the binomial expansion for (n + 1) i−1 (assuming n j = 1 for n = j = 0); to get the third equality the definition of J i,j in Eq. (A.3) is used.
Our previous result for J i,j (that is, only if j = i − 1 it is non-zero) is used to obtain the second equality in the second line: only the larger value of j in the sum is kept and the definition for J i,i−1 is used. To obtain the one-to-last equality we write the binomial in terms of factorials and pull out the largest factor in both (i − 1)! and n!; writing the various factorials back as a binomial one arrives at the last equality. Applying the recursion relation i − 1 times we obtain the result in Eq. (3.5).
Next we show a result involving the sum of binomials, necessary to obtain the final form in Eq. (4.3). We start by noting that the harmonic numbers can be written as the following integral where in the first step we have used a n − b n = (a − b) n−1 i=0 a i b n−1−i with a = 1. To find the relation we seek one switches variables x → 1 − u in the integration above where we have used the binomial expansion to obtain the second equality, and noticed that the i = 0 contribution in the sum cancels against the term outside; after integrating the remaining terms the final expression is obtained.
In this final part of the appendix we show how to convert a sum over harmonic numbers into an integral, that is, we provide a proof for Eq.
where in the first step we switch variables ε = −β(1 − x); to obtain the first expression in the second line we perform the binomial expansion of (1 − x) i−1 ; next we write the logarithm in terms of a derivative and integrate each term in the sum. After carrying out the derivative we write the binomial in terms of factorials and use (i − 1)! = i!/i and j 2 (j − 1)! = j × j! to obtain yet another binomial. Finally, using the result in (A.9) we conclude our proof for Eq. (4.9).