Massive Event-Shape Distributions at N$^2$LL

In a recent paper we have shown how to optimally compute the differential and cumulative cross sections for massive event-shapes at $\mathcal{O}(\alpha_s)$ in full QCD. In the present article we complete our study by obtaining resummed expressions for non-recoil-sensitive observables to N$^2$LL + $\mathcal{O}(\alpha_s)$ precision. Our results can be used for thrust, heavy jet mass and C-parameter distributions in any massive scheme, and are easily generalized to angularities and other event shapes. We show that the so-called E- and P-schemes coincide in the collinear limit, and compute the missing pieces to achieve this level of accuracy: the \mbox{P-scheme} massive jet function in Soft-Collinear Effective Theory (SCET) and boosted Heavy Quark Effective Theory (bHQET). The resummed expression is subsequently matched into fixed-order QCD to extend its validity towards the tail and far-tail of the distribution. The computation of the jet function cannot be cast as the discontinuity of a forward-scattering matrix element, and involves phase space integrals in $d=4-2\varepsilon$ dimensions. We show how to analytically solve the renormalization group equation for the P-scheme SCET jet function, which is significantly more complicated than its 2-jettiness counterpart, and derive rapidly-convergent expansions in various kinematic regimes. Finally, we perform a numerical study to pin down when mass effects become more relevant.


Introduction
Since the late 70s, a class of observables called event shapes has been used to test and determine fundamental properties of QCD (for a review see [1,2]), most notably to measure the strong coupling. As the name suggests, these observables contain information about the geometric momentum distribution of the final-state particle momenta. The (historic) main fields of application are e + e − collisions and deep inelastic scattering (DIS), but today there also exist adaptations developed specifically for pp colliders. In high-energy experiments most of the time it is sufficient to use the approximation that all particles in the final state are massless. If one is interested in high-precision calculations or in cases where the quark mass is a dominant effect, this approximation is no longer valid. While the theoretical computation of event-shape distributions for massless quarks at e + e − colliders has been pushed to unprecedented precision in recent years, including fixed-order results to O(α 3 s ) [3][4][5][6][7][8] and resummation at next-to-next-to-leading-logarithm (N 2 LL) [9][10][11] and next-to-next-to-next-to-leading-logarithm (N 3 LL) [12][13][14][15][16], computations for massive quarks remain at a lower precision. Such high level of resummation for massless quarks has been achieved in most cases using Soft-Collinear Effective Theory (SCET) [17][18][19][20][21] (or in the equivalent formalism by Collins, Soper and Sterman [22][23][24][25][26]), in which the most singular terms of the distribution are written in a factorized form, and summation of large logarithms is carried out with standard renormalization-group evolution. On the other hand, using the coherent branching formalism [27], the resummation can be automatized up to N 2 LL with numeric codes [28,29]. The results have been extended to the case of oriented event shapes in Ref. [30], and can be used to extract the strong coupling with high precision by comparing the theoretical expressions to LEP data, see e.g. [10,13,[31][32][33][34][35].
A first step towards having a firmer theoretical control over massive event shapes was taken in Ref. [36]. In that article, the fixed-order differential and cumulative cross sections were computed at O(α s ) for any event-shape. Since quark masses screen collinear divergences, there are only soft singularities that translate at threshold into two types of singular terms: a Dirac delta and a plus function (at higher orders other singular functions might appear). In Ref. [36] the coefficients of the delta and plus functions were computed analytically, providing a closed 1-dimensional integral form for the former, and showing that the latter is universal. Furthermore, the delta function coefficient was provided in an analytic form for most event shapes. An algorithm was devised to compute the nonsingular terms through a single integral (which can be carried out analytically in a few cases), a method much faster and accurate than binning the distribution in conjunction with a Monte Carlo integrator. In the present paper, we complement those developments by adding resummation of large Sudakov logarithms at N 2 LL. This kind of resummation for boosted quarks has been worked out only for the hemisphere-mass (doubly differential) distribution in Refs. [37,38], which can be easily marginalized into heavy jet mass and 2-jettiness (and with extra little work, also into a variable called C-Jettiness in Ref. [36]). In both cases, and only through recent computations that will be reviewed in the next paragraph, N 3 LL precision has been achieved. In this article we take the necessary steps to bring the accuracy to N 2 LL + O(α s ) for SCET-I type observables in the massive E-and P-schemes.
Although SCET was first developed as an effective theory for massless particles, it was quickly generalized to include massive quarks in Ref. [39]. Resummation is achieved in SCET by writing the most singular terms of the cross section as the product or convolution of various pieces: the hard matrix element, which is the square of the matching coefficient between SCET and QCD; the jet function, describing radiation collinear to the initiating partons; and the soft function, accounting for wide-angle soft radiation. Quark masses are infrared modes of the effective theory and therefore do not show up in the hard matching coefficient, but they do contribute to the jet (starting at one-loop) and soft (starting at two loops) functions. In the limit in which the difference between the jet and the heavyquark massess is much smaller than either of those, a new physical scale emerges, requiring additional resummation of large logarithms. This is achieved by matching SCET into boosted heavy-quark effective theory (bHQET). For jettiness all the ingredients necessary to build an N 3 LL-accurate cross section are known: bHQET jet function [40] SCET massive jet function [41], and bHQET matching coefficient [42], as well as the secondary production contributions [43,44] (all up to two-loop level). The computations carried out in this article make possible N 2 LL resummation both in SCET and bHQET for other classes of event shapes. Masses can also modify the endpoint of the distribution, which is their only effect at LL and NLL. This threshold modification depends solely on the scheme used to define the massive event shape.
In the approach followed in [45] to build the differential cross section, kinematic power corrections are taken into account in fixed-order perturbation theory by matching the SCET cross section onto full QCD. This setup was successfully applied to phenomenological analyses of massless event shapes, and is crucial to obtain reliable predictions in the tail and far-tail of the distribution. In the case of massive event shapes the situation is a bit more complex for two reasons: a) the partonic threshold is different in full QCD and the EFT, and b) the EFT prediction does not completely reproduce the singular terms at threshold. The reason is simple to understand: since the mass is an infrared mode, it is power-counted together with other infrared scales such as the quark virtuality, in such a way that powersuppressed terms include kinematic as well as mass corrections. Since the main focus of this article are event shapes with mass-independent thresholds, issue a) is not relevant. At leading power, SCET and bHQET only predict the leadingm behavior of the singular terms in QCD. 1 It is however possible to modify the hard and jet functions to fully account for the singular terms in the factorization theorems, such that power corrections are not distributions and behave well close to threshold. We will show how to apply this prescription to event shapes in the E-and P-schemes.
Non-vanishing quark masses imply that the mass sensitivity can be tuned using different schemes for the event-shape definition. To the best of our knowledge, this possibility has only been studied in the context of fixed-order perturbation theory in Ref. [36], and in the present work this analysis is continued by considering resummation of large logarithms. Furthermore, we complete the list of ingredients necessary for a full N 2 LL + O(α s ) compu-tation, which are applied to the case study of P-scheme thrust, making our work a useful reference for upcoming studies. Our result is essential to consistently include bottom-quark mass effects in analyses that aim to extract α s from fits to LEP data, and to clarify the top quark mass interpretation problem.
This article is organized as follows: in Sec. 2 we discuss kinematics in the dijet limit, which is the core for factorizing cross sections; in Sec. 3 we review the concept of massive schemes, explore how they affect the mass sensitivity of cross sections, and derive the implications they have for the SCET power counting; the factorization theorems in SCET and bHQET are presented in Sec. 4, in which also large-logarithm resummation is introduced; we carry out the computation of the SCET massive jet function in Sec. 5, and use this result in Sec. 6 to write down the fixed-order expression of the cross section in the dijet limit; in Sec. 7 we compute the bHQET jet function, and derive analytic results for the running of the non-distributional pieces of the SCET jet function in Sec. 8, presenting some useful expansions that can be implemented in different regions of the spectrum; kinematic, massive and non-perturbative power corrections are discussed in Sec. 9, while some numerical investigations are carried out in Sec. 10; finally, our conclusions are contained in Sec. 11. Some technical aspects of the computations are relegated to Appendices A and B.

Dijet Kinematics
In this section we study the kinematics of dijet events: two narrow, nearly back-to-back jets, plus additional soft radiation. In this situation the value of the event-shape is not far from its minimal value, which we call collectively e min , and the SCET power-counting rules apply. Particles in such events can be either soft, n-collinear orn-collinear, with momenta whose light-cone coordinates p µ = (p + , p − , p ⊥ ) scale like p µ s ∼ Q(λ 2 , λ 2 , λ 2 ), p µ n ∼ Q(λ 2 , 1, λ) and p μ n ∼ Q(1, λ 2 , λ) respectively, with λ the SCET power-counting parameter. Since we are interested in the primary production of heavy quarks and p 2 = m 2 , for consistency one has thatm ∼ λ and soft particles are either massless or havem s λ 2 . The SCET scaling holds for momenta defined in the P-and E-schemes as well. In this limit it is possible to write the event-shape measurement e as a sum of contributions from collinear (in both directions) and soft particles [46,47]: e = e n + en + e s , (2.1) where e denotes the event-shape measurement in the dijet limit at leading power, and we consider only soft perturbative particles. For SCET-I type observables, the three terms scale like λ 2 and are equally important. Moreover, e n , en and e s can be written as a sum over single-particle contributions. If the event shape is already defined by a single sum of finalstate particle momenta (as it is the case for thrust or angularities) this statement is trivial. When the event shape correlates momenta of final-state particles [ e.g. the definition involves a double sum (C-parameter) or there is a single sum squared (jet masses) ] the situation is more complicated. In the latter case one has to show explicitly that in this the dijet limit the leading contribution to the event shape can be written in the form of Eq. (2.1), as was done e.g. in Ref. [14] for C-parameter. Let us work out this decomposition for hemisphere masses, which are defined as the square of the total four-momentum flowing into one of the hemispheres, being those delimited by the plane normal to the thrust axis. Taking the z axis in the thrust direction, the mass of the plus hemisphere takes the following form in light-cone coordinates: where we have used that the component of the total hemisphere momenta normal to the thrust axis is identically zero. Let us assume the negative direction of the z axis pointing towards the plus hemisphere such that it does not containn-collinear particles (with p z > 0). In the dijet limit, particles i and j can be either soft of collinear, but if both are soft the corresponding contribution is ρ s,s ∝ O(λ 4 ) and therefore power suppressed. Next we consider that i is soft and j is collinear, such that the leading contribution comes from the p − j term: where we have used that up two power corrections in the dijet limit p − i = 2E i and the total available energy Q is carried by collinear particles only and equally divided into each hemisphere, such that the total minus momentum flowing into the plus hemisphere is Q up to power corrections. The set of soft particles that belong to the plus hemisphere is denoted by s + . With an identical computation we get for the collinear-collinear contribution with c + the set of n-collinear particles. In the dijet limit we have Qρ + = P + s + + P + n , which is equal to the total plus momentum entering the minus hemisphere. An identical reasoning leads to Qρ − = P − s − + P − n . Since we have not made any assumption on the mass of the particles, our result is valid for massless and massive quarks.

Massive Schemes
In this section we review the generalization of event-shape measurements for massive particles that go under the name of "massive schemes" (which should not be mistaken with the quantum field theory mass schemes such as pole, PS, 1S, MSR, MS, . . . ). These schemes were introduced in an article by Salam and Wicke [48] to study the effects of hadron masses on hadronization power corrections, which were further studied in [47]. Both studies consider light quarks only, such that massive schemes have no effect on partonic computations, but change the size of non-perturbative power corrections. For massive quarks, however, switching schemes can dramatically change the cross section, in particular its sensitivity to the quark mass, which is obviously of high interest either in cases where very accurate computations demand including quark mass corrections, or when they are the leading effect. A propaedeutic study of mass schemes for heavy quarks in event-shape cross sections was carried out in Ref. [36], where fixed-order results where computed and subsequently analyzed. A more complete study demands resummation of large logarithms in the peak and tail of the distribution, and we aim to fill in this gap here.
Massless particles travel at the speed of light and their four-momentum satisfies p 2 = 0, which when translated into energy and momentum in a given frame implies E p = | p |. Therefore on can interchange E p ↔ | p | in the event-shape definition with no visible effect. If particles become massive, p 2 = m 2 and one has E p > | p |, meaning in turn that replacing the momentum magnitude by the energy (or vice versa) changes the value of the event shape. In this context, various mass schemes where defined to quantify this possibility, which reproduce the original "massless" definition in the limit m → 0: 1. E-scheme: As indicated by its name, one replaces momenta by energies with the substitution The scalar product takes the following form: such that p 2 E = 0 even for massive particles. One advantage of this prescription is that hadronization corrections become universal. The variables called angularities [49] were originally defined in this scheme.
2. P-scheme: Again the name suggests that energy gets replaced by momenta as , and it happens that most of the classical event shapes were originally defined in this scheme: thrust [50], C-parameter [51,52] and broadening [53]. The scalar product now reads p P ·q P = | p || q | − p · q, which again implies p 2 P = 0 for massless and massive particles.
3. M-scheme: The name "massive scheme" is used for event shapes that in their original definition where neither in the P-nor in the E-scheme, such as heavy jet mass [54][55][56]. Their definition contains both energy and momentum and are the most sensitive to quark masses, in particular because in this scheme the usual relation p 2 = m 2 is satisfied.
It is important to realize that four-momenta as defined in the E-and P-schemes are framedependent, and that event shapes are usually defined in terms of magnitudes measured in the center-of-mass frame. The usual light-cone decomposition applies in either scheme S with S = E, P . The specific definition of the event shapes just introduced can be found e.g. in the original papers and will not be repeated here unless necessary. They are also summarized, including a discussion on massive schemes, in Ref. [36].

Mass sensitivity
When studying the sensitivity of event shapes at parton level 2 the first possible contribution for e + e − annihilation comes from the production of a heavy quark-antiquark pair without 12m 2 (1 −m 2 )m 2 P-and E-schemes 0 0 0 0 Table 1. Threshold position for various event shapes in the case of primary production of a stable quark-antiquark pair in different massive schemes. We use β = √ 1 − 4m 2 , the velocity of the quarks at threshold in natural units. additional radiation. In this case, the thrust axis is parallel to the three-momenta of the quarks, which makes trivial to calculate the threshold for two particles in the final state with equal mass m. Moreover, this simple computation sets the lower threshold even if additional gluons and (massless) quarks are radiated. The results in Tab. 1 show that for events in which a massive stable quark-antiquark pair is produced (primary production) only the M-scheme is sensitive to the quark mass while P-and E-schemes are not. In most of the events there will be some extra radiation present which will modify the former dijet into two fatter jets or an even more isotropic momentum distribution. For the observables we study, such processes will mainly contribute for event-shape values away from threshold value adding subleading mass sensitivity (i.e. suppressed by a factor of α s ) even in the Pand E-schemes, but will not substantially change the leading sensitivity of the M-scheme definition since it comes from the tree-level peak position. From this we can conclude that the M-scheme is preferred if the aim is a mass-sensitive observable (e.g. for quark mass determinations), but in case that one seeks a mass-insensitive observable, the P-and E-schemes are a better choice. 3

Mass schemes in the collinear limit
Collinear particles in the n direction satisfy E p = (p + + p − )/2 = p − /2 + O(λ 2 ) and also such that the E-scheme (Q = i E i ) and P-scheme (Q p ≡ i | p i |) normalizations are the same at leading power. Let us compute the fourmomenta of massive n-collinear particles in the E-and P-schemes. Since we have seen that E p /| p | = 1 + O(λ 2 ), the "large" (or label) components p − and p ⊥ are the same in any scheme. Let us then focus in the small p + momenta, which in the massive scheme takes the following form (for n-collinear particles the z component of momenta is negative): where we have used the so-called perpendicular energy E ⊥ ≡ | p ⊥ | 2 + m 2 . The computation of p + in the P-scheme is very similar since one only has to use | p | instead of

4)
3 If the massive partons enter the final state via gluon splitting in a massive quark-antiquark pair (that is, through secondary production) the sensitivity to the quark mass will again be subleading (now suppressed by a factor of α 2 s ).
where we notice that in the second step p + P takes the same analytic form as in Eq. (3.3) with the replacements E p → | p | and E ⊥ → | p ⊥ |, in the third step we use | p | = E p + O(λ 2 ), and in the last step we replace | p ⊥ | 2 = p + p − − m 2 and E p = p − /2 + O(λ 2 ). It is important to notice that in general p + P = p + because m 2 /p + O(1). Moreover, the mass that appears in Eq. (3.4) (and in other any event-shape measurement function) comes from kinematic on-shell considerations and therefore corresponds to the pole scheme. Finally, let us compute the E-scheme p + component, and for that we only need to make the replacement where in the last step we again have used E p /| p | = 1 + O(λ 2 ). Since for collinear particles (in any direction) p µ P = p µ E at leading power, we can safely conclude that at this order the collinear measurements for all event shapes take the same form in the P-or E-schemes, but is in general different from the M-scheme. In what follows we will work out the collinear measurement for a few event shapes.

Thrust
The original thrust definition is already in the P-scheme and reads witht the thrust axis. For collinear particles we have that |t · p i | = −p z and therefore up to power corrections we have where we already indicate that the collinear measurement is the same in the E-scheme. In Ref. [57] an M-scheme generalization of thrust, dubbed 2-jettiness, was introduced such that its collinear limit is the total plus momentum flowing into the plus hemisphere Qτ J c = i∈+ p + i . Since the measurement is completely inclusive, the computation of the jet function can be carried out as the imaginary part of a forward-scattering matrix element. This is not the case for the E-and P-schemes if quark masses are non-vanishing.

Hemisphere Jet Masses
We already worked out the collinear measurement for heavy jet mass in Eq. (2.4), and getting the P-scheme measurement is equally simple since Eq. (2.2) still applies with minimal modifications: where we again use that the total perpendicular momentum vanishes and use the identity . With this result we can trivially obtain the collinear measurement using that p − P,j = 2E j + O(λ 2 ): that matches the P-scheme thrust result. The total perpendicular momentum does not vanish in the E-scheme, since there is not such thing as E-scheme three-momentum conservation (in the same way, P-scheme energy is not conserved either). However, in the dijet limit the perpendicular components are already O(λ) and therefore p ⊥ , thence power suppressed, such that the result in Eq. (3.10) is also valid for the E-scheme.

C-parameter
In Ref. [14] it was shown how the C-parameter measurement splits into the sum of soft and collinear contributions in the dijet limit. The proof relied on the particles being massless, so it cannot be taken for granted that it will work when quarks have a non-zero mass. Here we carry out a similar proof valid for massive particles as well. C-parameter is defined already in the P-scheme as where we simply use sin 2 (θ ij ) = 1 − cos 2 (θ ij ) = [1 + cos(θ ij )][1 − cos(θ ij )] and the definition of the euclidean scalar product to get to the final form. Next one can express the result in terms of P-scheme light-cone coordinates using Eq. (3.2) as follows , with a similar result in the E-scheme. Arguments analogous to those used in Ref. [14] apply, and we focus on the collinear measurement only. First consider that both i and j are n-collinear such that the SCET scaling implies where once again we use that the total collinear perpendicular momenta flowing into the plus hemisphere is zero up to power corrections. One gets an analogous result for C P nn , while if i is n-collinear and j isn-collinear we get (we already include a factor of 2 to account for the case in which i isn-collinear and j is n-collinear) (3.14) Summing up the contributions in Eqs. (3.13) and (3.14) we obtain the collinear measurement in the P-and E-schemes: 15) identical to that of thrust or the hemisphere masses up to a factor of 6. The computation for the E-scheme is identical, relies on arguments already exposed, and therefore will not be repeated. For the M-scheme C-jettiness variable introduced in Ref. [58] it was shown in Ref. [59] that the collinear measurement takes the simple form QC J c = 6 i∈+ p + i . Therefore, for the reduced C-parameter variable C ≡ C/6 the collinear measurement for the three event shapes we consider coincide in every massive scheme.

Factorization Theorems
For simplicity we consider the well-known case of thrust (in either massive scheme), which can be easily modified to obtain the corresponding factorized results for C-parameter or heavy jet mass. After having shown that in the three schemes considered Eq. (2.1) holds, the derivation of the factorization theorem is obtained following the steps outlined in Ref. [46]. The leading hadronization corrections (which are soft) can also be factorized as an extra convolution with the so-called shape function, and even though they are included in our numerical analysis, for the sake of conciseness we ignore them in this section. Likewise, we include kinematic and mass power corrections in our final analysis, but postpone their discussion until Sec. 9

SCET
The value of the quark mass can have different hierarchies with respect to the (EFT) hard, jet, and soft scales. In Refs. [43,44] it was extensively discussed how to setup a consistent variable-flavor number scheme for final-state jets accounting for primarily and secondarily produced massive quarks. Four scenarios can be defined for the cases in which m > µ H (scenario I, which is of no interest for primary quarks since there is no energy to produce them), µ H > m > µ J (scenario II, relevant for very boosted heavy quarks, and better described in bHQET), µ J > m > µ S (scenario III) and m < µ S (scenario IV). Each of them has a different factorization theorem and renormalization group evolution setup. Even though the heavy quark mass is a fixed parameter, the jet and soft scales depend on the event-shape value and therefore they change along the spectrum, such that several scenarios might occur in a given distribution. For simplicity, we assume the quark mass is always smaller than the soft scale, such that we stay in scenario IV even in the peak of the distribution. In this way, we avoid having to deal with integrating out the heavy quark mass and the partonic factorization formula reads 4 with σ 0 the Born or point-like (massless) cross section, H and S τ the hard and soft functions, respectively, and the thrust jet function, which is the convolution of two single-hemisphere jet functions. The M-scheme hemisphere jet function has support for s > s min = m 2 , what sets the integration limits in Eq. (4.2). Accordingly, the thrust jet function has support for s > 2s min , implying that the minimal value for 2-jettiness is τ J min = 2m 2 . We shall present the computation of the M-and P-scheme SCET jet function in Sec. 5. The definition of the soft function in terms of Wilson lines can be found e.g. in Ref. [46] and the corresponding expression for the jet function will be given in Sec. 5. The factorization formula takes a simpler form in Fourier space In Eqs. (4.1) and (4.3) all matrix elements appear evaluated at the same renormalization scale µ. In order to minimize large logarithms that appear in each of them one should use RGE equations to evaluate them at their respective natural scales, denoted by µ H ∼ Q, µ J ∼ Q √ τ and µ S ∼ Qτ , such that for small τ there is a strict hierarchy among those: µ H > µ J > µ S and the SCET scaling parameter takes the value λ ∼ √ τ . The form and solution of the renormalization group equations is also simpler in position space. Using those and changing variables to y = x/p one arrives at where µ i denotes collectively all renormalization scales (including the common µ) and we use the following compact notation withω andk the exponential running kernels defined in terms of integrals over the SCET and QCD anomalous dimensions as follows .
Here γ can refer to cusp or non-cusp anomalous dimensions, and their dependence on α is in the form of perturbative series that define their respective coefficients The integrals in Eq. (4.7) can be solved analytically in terms of the anomalous-dimension coefficients if an expansion in α s is carried out. Their explicit form up to N 3 LL can be found e.g. in Ref. [31]. General expressions valid for arbitrarily high order can also be derived and will be given elsewhere.
The jet function of a massive quark contains terms which are distributions, and hence easy to Fourier transform, plus others which are regular functions, and to the best of our knowledge it seems impossible to find an analytic expression in position space for them. Up to one loop, the momentum-space hemisphere jet function can be decomposed in the following form: where the massive corrections, either with distributions J m or fully non-distributional J nd , are µ independent dimensionless functions with support for positive values of their (dimensionless) arguments. The µ dependence of J dist is entirely determined from the jet and QCD anomalous dimensions, does not depend on the quark mass, and therefore can be fully accounted for in the massless jet function of Eq. (4.9). The only piece that needs an explicit computation in the E-and P-schemes is J nd , since the rest can be obtained using consistency conditions and results obtained in Refs. [38] and [36].
The integral in Eq. (4.5) can be easily solved for all terms involving only distributions, and generic formulas can be found for instance in Ref. [12]. For the non-distributional piece of the jet function we carry out resummation in momentum space, and at one loop it is multiplied by the hard and soft functions at tree-level only. Therefore, using Eq. (4.4) in (4.5) and carrying out the y integration, the non-distributional part of the 1-loop partonic cross section reads The lower limit of integration in the first line has been moved to 2s min since below that value the jet function has no support. In the E-and P-schemes s min = 0 so we have not lost any generality. To get to the second line we have switched variables in the integral to s = xm 2 + 2s min , and to obtain the second expression for I np (ω, y) we switch variables to x = zy. For the partonic cumulative distribution one gets instead To make the function I(ω, y) smooth in the no resummation limit, achieved whenω → 0, one can integrate by parts to obtain with J np the derivative of the J np function. This form is particularly useful if the integration has to be carried out numerically, making it more convergent and defining its analytic continuation to values 0 <ω < 1. Further integration by parts can be implemented to define the integral for even larger values ofω. If a closed analytical form is found, this procedure is unnecessary. Although the discussion in this section has been carried out assuming the pole mass for the heavy quark, it is straightforward to convert the result to a short-distance scheme. In the scenarios in which SCET applies, the MS scheme is perfectly adequate. In scenario II it is more convenient to employ low-scale short-distance schemes such as the MSR mass [60,61].

bHQET
If the heavy quark mass is large enough or if the jet is very narrow one enters scenario II, in which the jet and heavy quark masses are close to each other, corresponding to very boosted quarks. In this kinematic situation, a new physical scale emerges µ B ∼ Qτ /m ∼ µ S /m, such that there is a new hierarchy between scales: µ H > m > µ B > µ S . 5 A practical way to see how this becomes manifest is looking at the structure of the one-loop jet function in Eq. (4.9). Since the non-distributional terms are power suppressed when s → s min , it is enough to focus on the terms with distributions, which generically read where the coefficients A, B and C (with or without subindex m) depend neither in µ nor in m. In scenarios III and IV one has s m 2 and therefore the choice µ 2 ∼ s makes sure there are no large logarithms in neither term (the massless limit is smooth since J m + J nd → 0 when m → 0 and no new class of large logarithm emerges). On the other hand, if s m 2 the choice µ 2 ∼ s cannot prevent the logarithms in J m from becoming large. 6 The heavy quark carrying momentum p = mv + k with v 2 = 1 gets integrated out as a dynamical degree of freedom giving raise to heavy-quark effective theory [62][63][64][65]. The remaining degrees of freedom are referred to as ultracollinear and carry residual momentum k. In the heavy quark rest frame [ in which v µ = (1, 1, 0 ⊥ ) ] they are soft k µ = ∆(1, 1, 1), with ∆ m a low-energy scale, 7 being able to interact with each other and with color sources representing the integrated-out heavy quarks. In the center-of-mass frame these momenta get boosted and a hierarchy is generated among their light-cone components 8 where we have included momenta q s which is soft in the center-of-mass frame. The two boosted copies of HQET are matched onto SCET in order to account for global soft radiation, such that the heavy quark and ultracollinear particles can interact with soft degrees of freedom. The typical off-shellness of ultracollinear particles is softer than for collinear degrees of freedom which are part of SCET. Using this framework, it is possible to derive a factorization theorem for the partonic cross section [37] which effectively separates physics at the different involved scales.
In practice one can still identify µB with the jet scale (µB = µJ ), and we will do so in what follows. 6 The bHQET limit should not be confused with the threshold limit, yet another interesting physical situation in which radiation other than the heavy quark is soft as compared with the quark mass, while there is no hierarchy between m and Q. 7 For an unstable top quark, this scale is of the order of its width ∆ ∼ Γ, but for a stable bottom quark it can be identified with ∆ ∼ Q 2 τ /m ∼ (Q 2 τJ − 2m 2 )/m for thrust and 2-jettiness, respectively. 8 If only two back-to-back quarks are produced, their velocity equals β = √ 1 − 4m 2 1 − 2m 2 , and therefore the boost factor reads γ = 1/(2m). When boosting momenta in light-cone coordinates the plus/minus components get multiplied/divided by γ( where the hard and soft functions are the same as in the SCET factorization theorem, but there is an additional matching coefficient H m between SCET and bHQET. The jet function B τ (ŝ) is different from J τ in SCET, has support forŝ > 0, its mass dependence is only through a global 1/m factor and contains only distributions. It is also the convolution of two hemisphere bHQET jet functions B n : whose operator definition shall be given in section 7. Since H and H m are the same for all event shapes and S τ does not depend on the quark mass (it sees only light degrees of freedom), the anomalous dimension of the B n function is the same in any massive scheme.
In turn this implies that all terms in the jet function except for the Dirac delta are fixed by consistency and hence are the same in the M-, P-and E-schemes. Knowing H, H m and S τ at one loop, the delta function coefficient in B n can be obtained taking them → 0 limit of the result quoted in Ref. [36] for the full QCD prediction for the threshold delta function coefficient. In this sense, our computation in section 7 will be just a sanity check. The H m matching coefficient and bHQET jet function satisfy , such that both factorization theorems smoothly join. For the various arguments exposed in this paragraph, this also implies that the coefficients B m and C m in Eq. (4.13) do not depend on the massive scheme, and they are known from the 2-jettiness computation of Ref. [38]. Therefore, different massive schemes can differ only in A m and the J nd , but knowing the one-loop hard and soft functions, A m can be obtained again from the massless limit of the full-QCD threshold result. Carrying out resummation in bHQET is identical to massless SCET. Since we will limit our numerical analysis to situations in which it is sufficient to use SCET with masses, we will not give further details on how to solve the corresponding RGE equations, which can be found elsewhere. Moreover, we will not provide a detailed discussion on how to switch to a short-distance mass scheme in this setup.

SCET Jet Function Computation
The jet function accounts for the dynamics of collinear particles within the hemisphere. Since the collinear measurement function in the P-and E-schemes is not the total plus momentum, it cannot be computed as the discontinuity of a forward scattering amplitude, as was done in [37,38]. Instead, one has to use the definition given in Ref. [46], which after a small modification to match the form of our factorization theorem and minimal manipulations can be cast into the form.
with χ n,Q the jet field with total minus momentum equal to Q, d = 4 − 2ε the spacetime dimension in dimensional regularization, − = Q and ⊥ = 0 due to label momentum conservation, and the trace is taken over spin and color indices. The n-collinear eventshape operatorê n acting on some final state |X pulls out the value of the event shape for the e n (X), the contribution from n-collinear particles to the value e of the event shape: e n |X = e n (X) |X . To simplify the expression of in the second line of Eq. (5.1) we insert the identity I = X |X X| after the delta function and shift the field χ n,Q (x) to x = 0 employing the momentum operator. Using the label operators for the large components of the momenta, the sum over X can be carried out and we obtain the following convenient expression For practical computations one inserts a complete set of states after δ(P − Q) where we exclude the vacuum from the sum because it does not contribute to the jet function. Each term in the sum over n can include several contributions, accounting for various particle species (heavy/light quarks and gluons), and the sum over polarization affects all particles in the final state. The perturbative expansion of the jet function in powers of α s is obtained by adding more particles to the sum as well as more virtual (loop) contributions to the matrix elements that appear after inserting the identity, which in compact form can be written as For the computation of the P-scheme hemisphere jet function one does not need any regularization beyond taking the space-time dimension from 4 to d = 4 − 2ε. In the following we carry out the computation of the jet function using Eq. (5.2) for both 2-jettiness and P-scheme thrust. Although the result is already known for the former, it is instructive to repeat its computation to highlight the difference between the two approaches. In a way, the computation that uses Eq. (5.2) can be obtained applying Cutkosky rules to directly compute the imaginary part of the forward scattering matrix element. It is important to remember that χ n in Eq. (5.2) is composed of bare SCET (quark and gluon) fields, and that it is convenient to carry out our computations using perturbation theory "around" those (that is, we will not use the so-called renormalized perturbation theory). For the jettiness computation through the discontinuity of the forward matrix element, this entails that the wave-function renormalization factor computed with the diagrams shown in Fig. 1 (the soft-gluon contribution vanishes), never appears directly. The mass in Eq. (5.5) should be understood as the pole scheme. When using Eq. (5.2) one needs to account for Z 1/2 ξ since this factor is precisely the overlap between the quantum (bare) collinear field ξ n and the physical collinear state |q n : 0| ξ n |q n ( p, s) = Z 1/2 ξ u s ( p ), with u a particle spinor in the collinear limit, and s, p the spin and 3-momentum of the on-shell collinear quark. On the other hand, when using Eq. (5.2) self-energy diagrams on external legs are not included, since their effect is already accounted for in the Z 1/2 ξ factor, and it is in this way that one has a one-to-one correspondence with the computation through the imaginary part of the forward matrix element. The computation at leading order is simple enough that can be carried out for the two massive schemes simultaneously. The corresponding tree-level diagram is shown in Fig. 2, where the double line represents a heavy quark and the dashed line marks which particles are on-shell. To compute the phase-space integration it is convenient to use the following parametrization which implies that p + has to be expressed in terms of the minus and perpendicular components through the on-shell condition p + = (m 2 + | p ⊥ | 2 )/p − , and since the mass appears through on-shell kinematic relations it corresponds always to the pole scheme. In the second line we have carried out the angular integrals for the perpendicular momentum, assuming that matrix elements depend only on its magnitude. We then obtain (5.7) where we have used that the trace of the polarization sum equals 2p − and have integrated all delta functions except the one with the measurement. The color trace cancels the 1/N c prefactor, and the on-shell condition implies p + = m 2 /Q, such that for the 2-particle collinear measurement we get which correspond to e min . To include the wave-function renormalization at O(α s ) one only needs to multiply this result by Z ξ . The contribution from virtual gluons can be carried out for the two massive schemes simultaneously since the phase-space integration is identical to the tree-level computation. There are two diagrams contributing, as shown in Fig. 3, which yield the same result, so we will compute only one of them which will be multiplied by a factor of 2. Pulling out a collinear gluon field from the Wilson line and using the Feynman rules for massive collinear quarks we obtain the following integral for the leftmost diagram:

Virtual radiation
whereμ 2 = µ 2 e γ E /(4π), the factor of 2 comes from the productnn and the Casimir C F from the color trace with two Gell-Mann matrices. We are left with two master integrals I 1 and I 2 that can be solved using Feynman parameters for the former and with a combination of Feynman and Georgi parameters for the latter Adding those two results we find a closed expression for J virt 1 : Interestingly, this result is zero in the massless limit (which has to be taken before expanding in ε). Therefore, using dimensional regularization, only real-radiation diagrams contribute. The m appearing in Eq. (5.12) is, strictly speaking bare, but since we limit our computation to O(α s ) we can safely take it as the pole mass as the difference between these two is a higher order correction. Implementing this result to the jet function computation and integrating the real momentum results in adding a factor of δ(e n − e min ). Multiplying by 2, expanding in ε and adding the wave-function renormalization, which is obviously a virtual contribution, we obtain which is the final result of this section. Equation (5.13) should be valid also for SCET-II type observables since the 1-particle phase space is not yet afflicted by rapidity divergences. Since the phase-space integrals with two particles do not fully collapse with the Dirac delta functions, the real radiation contributions differ depending on the collinear measurement. The diagrams that contribute at O(α s ) are shown in Fig. 4, where we have omitted the term in which both gluons are radiated from the Wilson line since it vanishes. Diagrams (a) and (b) give identical contributions and therefore we will compute one of them which will be multiplied by a factor of 2. For all the real contributions label-momentum conservation implies q ⊥ = − p ⊥ and q − = Q − p − , which together with the Heaviside function in Eq. (5.6) sets the integration limits for p − between 0 and Q. For the first diagram, after integrating the gluon momenta with the delta functions and carrying out the angular perpendicular integration one gets

Real radiation
Since this diagram involves gluons radiated from Wilson lines, one expects 1/ε 2 poles, implying also harder integrals. On the other hand, since the corresponding Feynman rule is simpler, the result is also shorter. For the second (symmetric) diagram, which does not need a factor of two, the result reads Since Feynman rules are more cumbersome for gluons radiated from a massive quark, the result is lengthier. On the other hand, since these correspond to boosted QCD processes single 1/ε poles are expected to appear, meaning that no special treatment for the integral is necessary.

P-scheme Thrust
To solve the integrals in the previous section, the thrust measurement must be expressed in terms of p − and | p ⊥ |. Since the gluon is massless and the quark has mass m we have which is mass independent. With this result we can use the measurement delta function to Switching variables to p − = Qx we find the following 1-dimensional integral for J a : J real a,P (s, µ) = The complication arises because the integral diverges as 1/ε 2 and contains distributions. The divergence comes from the (1−x) −1−ε factor, but the subtraction around x = 1 behaves as 1/s, which combined with the s −ε prefactor implies a new divergence and invalidates the subtraction. This pathological behavior is usual in two-loop computations involving double integrals, and the standard way to solving it is using sector decomposition [66].
To do so, one needs to get rid of distributions by considering the cumulative jet function, which converts Eq. (5.18) into a double integral. In Appendix A we show how to use this general method to solve the integral, and follow in this section an easier, albeit less general, procedure. Before that, we solve the integral for the massless case, which is not affected by the problem just described, and is valid for 2-jettiness and thrust: In the second line we have expanded the result around ε = 0 to obtain distributions using the identity For the m > 0 case we can transform the integral using hypergeometric function identities. Since these special functions will appear also in Sec. (8), we remind here its integral definition (the hypergeometric function is symmetric with respect to its first two arguments): where the second line is the so called Euler transformation. 9 The integral in Eq. (5.18) is already into this canonical form and therefore we can write where in the second step we have used Euler's identity and in the third we write the hypergeometric function back as an integral. The integration in the last term can be easily expanded in ε using Eq. (5.20), and definings ≡ s/m 2 we have 10 with f 1 (s) a function involving a dilogarithm. This result can be reexpanded in ε together with the prefactors −1−2ε , responsible for the appearance of distributions, finding theñ This property can be easily shown as follows: switching variables x → 1 − x and rearranging terms one where the second term is obtained using the symmetry a ↔ b on the first equality. Using the first relation followed from the second one arrives to the second line of Eq. (5.21). 10 Alternatively one can use the Mathematica package HypExp [67] to expand directly the hypergeometric function.
Therefore one only needs f 1 (0), which takes a simple form Putting all partial results together and expanding in ε we get the following expression The result is divergent for s → m 2 , although at that kinematic point there is no physical phenomenon that implies a singularity. We therefore expect that the singularity will cancel when adding together all real-radiation diagrams.
For the cut self-energy diagram in Fig. 4 (c), performing the same change of variable as in Eq. (5.18) we arrive at To see how the 1/ε divergence occurs, we compute first the massless limit of J real b , for which we get At the light of this result one can realize that switching variables to x = ys/m 2 exposes the divergence, factoring it out front the integral: In the one-to-last step we have used Eq. (5.20) to partially expand in ε and in the last step we use Eq. (3.26) of Ref. [36]. The Dirac delta function δ(s) sets the upper integration limit to infinity and s = 0 in the integrand, becoming the integral so simple that no further expansion in ε is necessary for this term. For the contribution proportional to the plus distribution we can set ε = 0 right away, and solve the integral with standard methods. In the last step we have consistently expanded in ε the full result. The expression is again divergent as s → m 2 , but as anticipated, the full real-radiation contribution is regular in this limit: For completeness, we also provide the real-radiation contribution for the massless case, which coincides with the full jet function. Adding the tree-level result we recover the known result

Jettiness
Let us express the jettiness measurement in terms of minus and perpendicular components: which can be used to solve the measurement delta function for the magnitude of the perpendicular momentum The argument of this delta function can be zero only if sp − > Q m 2 , what sets the lower limit of integration. Therefore, changing variables to p − = Q(1 − x) we obtain for the diagram in which the gluons are radiated from the Wilson line and the quark particle the following result where in the second line we have switched variables to x = y(1 − m 2 /s). Performing the same change of variables in the diagram in which the gluons are radiated from both quark lines one gets where we have carried out the same manipulations as in Eq. (5.35). Adding the two results we obtain the total contribution for 2-jettiness:

Final Result for the Jet Function
Adding together the contributions from real and virtual corrections one obtains the complete jet function. The divergences are now entirely of UV origin and can be renormalized multiplicatively (with a convolution). Since in either massive scheme these are the same as for massless quarks, the renormalization factor is the same, along with the anomalous dimension. Therefore we quote directly the result for the renormalized jet functions [ α s ≡ α s (µ) ]: with A J = 2π 2 /3 + 1 and A P = 4π 2 /3 − 3. We shall see that J nd (x → ∞) = −J m (x) for both schemes, and show this behavior graphically in Fig. 5(a). This is expected since it corresponds to taking the massless limit, and therefore mass corrections should vanish such that the jet function becomes equal to the (renormalized) massless result of Eq. (5.32).
Since J nd contains distributions in this limit, it is advantageous to work with the cumulative jet function which is a regular function. Likewise, one can define the cumulative functions for J nd and J m , which are also shown in Fig. 5(b). The result can be obtained easily and involves polylogarithms: As expected, in the m → 0 limit both results reduce to the (same) massless cumulative jet function To take the derivative one needs to recall that the jet function has support only for positive s, such that it is effectively proportional to an (implicit) Heaviside function θ(s). Using the following relations: one readily arrives at Eq. (5.32). For s > 0 one can expand around m = 0 to find the following compact series with similar results for the cumulative jet functions. Since individual pieces of the P-scheme thrust jet function have divergences at s = m 2 it is convenient to compute the expansion of J P nd (x) around x = 1, which can be cast as (1 − x) i 9 + 5i (i + 1)(i + 2)(i + 3) .
one arrives at the fixed-order prediction for the partonic singular terms of the P-scheme thrust differential cross section: 11 In the same way, one can get a similar expression for the cumulative distribution Σ SCET P , which is among other things useful to take the m → 0 limit. The differential cross section has a similar structure in full QCD, although it is different for vector and axial-vector currents as discussed in Ref. [36], and for P-scheme thrust takes the following form 12 where C = V, A labels the type of current and with R C 0 the tree-level massive R-ratio. Analytic results for A C and B plus C can be found in Ref. [36] and we quote here the universal value for the latter: with β = √ 1 − 4m 2 , and where the first and second line of the expression in big parentheses correspond to the vector and axial-vector currents, respectively. One recovers the SCET result for small masses, A C (m → 0) = A SCET (m) and B C plus (m → 0) = B SCET plus (m), and also lim τ →0 The partonic fixed-order bHQET cross section is identical to the SCET one dropping F SCET NS . 12 In the threshold limit one gets the same result as in full QCD dropping F C NS (τ,m), which is a power correction. with α ∼ O(1). Since for thrust F C NS is only known numerically, in Fig. 6 we show a comparison of QCD and SCET results for the NLO corrections scaling the mass as indicated in Eq. (6.5). Excellent numerical agreement is found as τ → 0 for various values of α between 1.2 and 15. We show only the vector current as for small values of τ it is indistinguishable from the axial-vector one.

bHQET Jet Function Computation
The computation of the bHQET jet function is significantly simpler that for the SCET counterpart since in this EFT the mass is no longer a dynamical scale and we are left with tadpole-like integrals. As an immediate consequence of that, much as it happened for the massless SCET jet function, all virtual graphs are automatically zero in dimensional regularization since they are scaleless (this includes the wave-function renormalization factor). We are then left with the tree-level, which is common for both massive schemes, and realradiation diagrams. The collinear event-shape measurements are the same in SCET and bHQET, although the contribution of massive particles needs to be power expanded, such that using p = mv + k we obtain for thrust and 2-jettiness the following results where to get to the last expression in the second lines we use the on-shell condition for heavy quarks v· k = 0, to be discussed later in this section. The field-theoretical definition of the bHQET jet function can be obtained from the expression given in Eq. (5.2), and taking into account that applying the bHQET power counting to the minus component of momenta one obtains δ(P − Q) → δ(k − + q − ) the jet functions can be written as (7.2) Here K is an operator that pulls out the residual momenta of the heavy quarks and the (full) momenta of ultra-collinear particles. We have also used v ⊥ = 0 and W v + has the same functional form as W n replacing the collinear gluons by ultra-collinear fields: A n → A + : The bHQET phase-space integration involving a heavy quark gets also simplified, and using again p = mv + k one has that where in the second line we have used Eq. (4.14) and in the third we power-count away the k 2 in the delta function argument along with all terms but Q inside the Heaviside function. The on-shell condition for heavy quarks written in light-cone coordinates implies v − k + +k − v + = 0, in agreement with the argument of the delta function. When comparing to Eq. (5.6) we observe that the p − in the denominator got replaced by Q and that the p − integration is not limited to positive values only. The phase-space integration for ultracollinear particles stays the same as in SCET. Feynman diagrams look exactly the same in SCET and bHQET, with the replacement p → k for the heavy quark momenta. Let us compute the tree-level contribution as given in Fig. 2, which is analogous to the corresponding SCET calculation: Tr u s (p)u s (p) = δ(ŝ) , (7.5) where we have used that the trace of the polarization sum equals 4m and have integrated all delta functions except the one with the measurement. The on-shell condition k − = 0 makes both (shifted) measurements coincide at tree-level, see Eq. (7.1). There are some generic features to be learned from this diagram: since there is no Dirac structure in the diagram, the trace of the polarization sum will be always 4m at any loop order, and since there is always one heavy quark which brings an inverse power of Q through its phase space one has the following combination: Tr u s (p)u s (p) = 1 , (7.6) which eliminates the spurious dependence on m and Q, ultraviolet scales that should not appear in EFT computations. To make this non-dependence explicit at higher orders one can rescale the minus component of ultracollinear real particles as q − i = (Q/m) i , as we shall do in the rest of the section.
We turn our attention now to real-radiation contributions, for which we can simplify the heavy-quark propagator using v · k = 0. We start with diagram (a) of Fig. 4, that after applying the bHQET Feynman rules becomes The on-shell condition on the ultra-collinear gluon momenta implies in light-cone coordinates: Let us work out the measurements for thrust and 2-jettiness where we have used Eq. (7.1), the on-shell condition for heavy quarks and ultra-collinear massless gluons, and the fact that label momentum conservation implies k − = −q − . With this result it is very simple to solve the measurement delta function in terms of the perpendicular gluon momenta and we will use these results to compute the jet functions in the next two sub-sections.

Thrust
We start with the diagram in which the gluon is radiated from the Wilson line. Switching variables to − =ŝx we arrive at With an identical change of variables we arrive at the following result for diagram (b): Adding both diagrams with the appropriate factors we obtain the final expression for the P-scheme hemisphere jet function:

Jettiness
The Dirac delta function in Eq. (7.10) implies that there is a solution for | q ⊥ | only if − <ŝ, which can be implemented through a Heaviside function and bounds the upper integration limit for − . With the change of variables implemented in the previous section the integration limits are mapped to the interval (0, 1) and we get the following result for diagram (a): that, as expected, differs from the expression in Eq. (7.11) only in the non-divergent term of the delta-function coefficient. Similarly, we obtain for diagram (b) again almost identical to the corresponding P-scheme computation. Adding twice the first diagram plus the second we recover the known result for the 2-jettiness bHQET jet function: Both schemes have the same divergent structure and hence their anomalous dimension, as expected, is the same. Furthermore, the difference between the respective delta coefficients is the same as the same difference for the SCET jet functions. This result was also expected since both theories should smoothly match in the bHQET limit.

RG Evolution of the SCET Jet Function
In this section we solve the renormalization group equation for the non-distributional part of the jet function for thrust and 2-jettiness. This amounts to finding an analytic expression for the function I np defined in the last line of Eq. (4.10). Even though the result for I J np has been already worked out in Ref. [38], we present here the main steps to find the solution as they are illustrative. Using the rightmost integral expression of the bottom line in Eq. (4.10) we find brings the second term also into to a canonical form that we can easily integrate, finding where in the last step we have used the integral representation of the 3 F 2 function: a 1 , a 2 , b 1 , tz) , (8.4) with a 1 = a 2 = a 3 = 1, b 1 = 1 −ω and b 2 = 2. After adding the result for the first term we find an expression slightly simpler than that quoted in Ref. [38], although fully equivalent: which has a smoothω → 0 limit. For a numerical implementation, one can use standard routines to evaluate 2 F 1 hypergeometric functions in programming languages such as Mathematica, Fortran, Python, or C++. For 3 F 2 there are built-in routines in Mathematica and Python, while for other languages one can use a numerical integration over 2 F 1 as shown in Eq. (8.4). For P-scheme thrust we can write the logarithm as a derivative to bring all terms into a canonical form: For y > 1 each of the terms in the integral diverges when z = 1/y. We can regularize the divergence adding a small imaginary part y → y + i . This makes each integral complex, although the sum is real when → 0. To express our result in terms of a minimal set of hypergeometric functions, we use the following identity: Furthermore, one can use an additional identity to make the final result manifestly real also for the case y > 1 After solving all integrals, recursively applying Eq. (8.7) and transforming the hypergeometric functions using Eq. (8.8), one arrives at the second important result of this article: 13 with H a the harmonic number, which for non-integer values of a can be expressed in terms of the digamma function: H a = ψ (0) (1 + a) + γ E . Equation (8.10) has been cast in a way in which the no-resummation limitω → 0 is smooth. The singularities that appear in individual terms of J P nd for x = 1 manifest themselves now as a double pole in I P np at y = 1, which however is fictitious, as the result is indeed smooth at this value. To solve this problem in numerical implementations we provide in Sec. 8.2 an expansion of this result around y = 1 at arbitrarily high order. The result in Eq. (8.10) is adequate for a numerical implementation since the derivative with respect to ε can be taken numerically through finite differences. It can be also performed analytically, using Eq. (8.8) in 2 F 1 (1, 1+ε, 2+ω, 1−y) and the following identity The result as given in this equation is very convenient for a numerical implementation, since one only needs to evaluate two hypergeometric functions (which might be numerically expensive) using the following approximations: with a value of ε which can be safely taken as small as 10 −6 .
to arrive at the equivalent expression: The (1 − y) −3−ω factor and both hypergeometric functions are complex for y > 1 but the combination is real. To have all terms explicitly real for y > 1 one can use the following relation 1 +ω (2 +ω) 2 3 F 2 (2 +ω, 2 +ω, 2 +ω, 3 +ω, 3 +ω, 1 − y) , (8.13) which does not rely on numerical derivatives, is manifestly real for all positive values of y but is numerically unstable if y → 0. This poses no problem in practice, since for y < 1 one can simply switch to Eq. (8.12). To derive the result in Eq. (8.13) we proceed as follows: where in the first step we use Eq. (5.22), in the second we apply the chain rule on derivatives, and in the third line we use again Eq. (5.22) on the first term. Using the identity d da 2 F 1 (a, b, a + 1, z) = bz (a + 1) 2 3 F 2 (a + 1, a + 1, b + 1, a + 2, a + 2, z) , (8.15) in the two terms of Eq. (8.14) we arrive at the result quoted in Eq. (8.13). In Appendix B we present an alternative (although more complicated) expression for I P np which does not involve numerical derivatives and with every term manifestly real for y > 1. We use this result as an additional cross check of our analytic derivations. In any case, we shall see that for numerical implementations one never needs to use expressions involving hypergeometric functions.

Expansion around s = 0
For numerical implementation purposes, it might be convenient to obtain an analytic expansion of I np (ω, y) around y = 0. One can do so by using the known expansions for the hypergeometric functions, e.g. (8.16) but in order to have a relation valid at arbitrarily high orders it is simpler to use the expansion of J nd around s = 0 on the leftmost expression of the bottom line in Eq. (4.10) and integrate analytically term by term. It turns out that one can sum up the corresponding series using Eq. (8.16) to recover an expression analytically equivalent to Eq. (8.12). The master integrals that we will need are where the bottom line can be obtained from the top one acting with a derivative with respect to i. We then arrive at where we have used the Pochhammer symbol (a) n = Γ(a + n)/Γ(a) since it is convenient for an optimized numerical implementation. Both series converge well for |y| < 1, and therefore apply mainly in the peak of the distribution. For 2-jettiness the series can be easily summed up using Eq. (8.16) and the series expansion for the 3 F 2 hypergeometric function: For P-scheme thrust one can convert the term involving harmonic numbers into the derivative of ratios of gamma functions to use Eq. (8.16) and recover the result we already obtained with a direct integration. The numerical implementation of Ref. [45] (which dealt with 2-jettiness) did not use this expansion and the evaluation of the non-distributional jet function running was the most severe performance bottleneck for the code.

Expansion around s = m 2
The results obtained in Eqs. (8.10) and (8.12) are not useful for a numerical implementation in the vicinity of y = 1. When y is sufficiently close to unity one can switch to a series expansion to arbitrary high power using the change of variables z → 1 − z in the rightmost expression at the bottom of Eq. (4.10) and the following expansion: Terms have been combined such that the coefficient of each power in y has a well-defined z → 0 limit and therefore we can integrate coefficient by coefficient. In practice one can integrate each piece assuming a non-integer value of i and subsequently convert the gamma functions that would become divergent ifω = 0 using the identity As expected, there are large cancellations among the various terms for a given power of y, but when adding all contributions one gets the following nicely convergent series: −(ω + 6)(5ω + 9)] + [20i − 2 iω(i(3ω + 7) + 7ω + 9) + 4(5ω + 9)] where again special care has been taken to write the expression in a manner in which one can setω = 0 without any worries. The series converges well for |1 − y| < 1, and therefore, combined with the expansion worked out in the previous section, for P-scheme thrust one can use expansions if y < 2.

Expansion around s = ∞
Since the numerical evaluation of hypergeometric functions is slow, it is convenient to figure out another series expression (in this case, of asymptotic type) around s = ∞, which is of course tantamount to m = 0. This limit is very relevant, since it can be applied in the tail of the distribution and almost everywhere if the heavy quark mass is much smaller than the center-of-mass energy, as is the case for bottom quarks at LEP. Such asymptotic expansion was not known by the time in which Refs. [31,32,45] were published, and was significantly affecting the performance of the respective numerical codes. Even though one could, in principle, use known results for the asymptotic expansions of 2 F 1 and 3 F 2 hypergeometric functions, it is in practice simpler and more efficient to compute the series directly on its integral form. This is complicated since, as we shall see, the expansions involve powers of log(y), and so one cannot simply expand the integrand and integrate term by term, as we did in Secs. 8.1 and 8.2. We found out that the Mellin-Barnes representation with 0 < c < ν, is optimal to achieve our goal [68]. 14 After applying Eq. (8.24), the asymptotic expansion around X 1 is obtained integrating by residues the poles that appear on the real axis for t > ν (the poles for t ≤ 0 correspond to the expansion X 1). We work out this expansion for thrust and 2-jettiness in the rest of this section, but before that we note that the asymptotic expansion is well convergent if 1/y < 1, which for P-scheme thrust means that in numerical evaluations one can always use one of the three expansions discussed in this section and never needs to evaluate hypergeometric functions with dedicated routines. For jettiness the same statement is almost true, except in a small vicinity of y = 1 in which, to the best of our knowledge, no expansion can be used.

2-Jettiness
We start from the integral form given in Eq. (8.1). The only complication in this case is that we have to deal with a logarithm, which does not have the form in Eq. (8.24). However, it can be brought to the standard form using Eq. (8.2) 14 This representation can also be used to solve the RG equation exactly. Applying a Mellin transformation to the first line of Eq. (8.6) in the y-variable, solving the z-integral and transforming back one gets a closed (and rather short) analytic expression for I P nd in terms of MeijerG functions, which are not very convenient for a direct numerical evaluation, but can be related to hypergeometric functions.
with 0 < d < 1. Since the denominator of the first term in Eq. (8.1) is squared, when applying the Mellin-Barnes representation (8.24) the first pole appears at t = 2. This is accompanied by an extra power of y, such that we can nicely map the poles of first term into those of the second by shifting the integration variable t → t + 1 in the former. After integrating over z we obtain with 0 < c < 1. The integrand has a triple pole at t = 1 and double poles at natural values of t larger than 1. We compute the triple pole by itself and treat the rest generically using With this result we obtain the following asymptotic expansion: using again the Pochhammer symbol. We have written each coefficient in a form such that theω → 0 limit, relevant in the far tail of the distribution, is smooth.

P-scheme thrust
Applying the Mellin-Barnes representation in Eq. (8.24) to the first line of Eq. (8.6) and integrating over z we arrive at an expression that involves different powers of y with poles shifted accordingly. Therefore, using the same strategy as in the previous section, we can shift the integration variable by one or two units such that poles and powers of y in each term exactly match. This is very important, since the expansion in 1/y must be carried out consistently given the large cancellations that take place among the various terms due to the divergence at x = 1 of individual terms in J P nd (exactly as it happened for the expansion around s = m 2 ). After some work we arrive at the following expression: We have already implemented a few simplifications because we assume the result is real, and therefore discarded the imaginary parts that would arise from (−y) −t . We have checked that indeed this is the case as long as one expands strictly in y without mixing any powers. Harmonic numbers are caused by the term in J P nd proportional to log(z). The integrand has now double and triple poles, located at natural values of t, the latter arising precisely because of the harmonic numbers. There are no poles arising from H −t−ω because of the corresponding gamma function in the denominator which has poles at the same values, making the ratio regular. To compute the residues of the poles we need, on top of Eq. (8.27), the following expansion which can be obtained from the relation between harmonic numbers and the digamma function and a bit of algebra. Using these results we arrive at the following expression, in which again special care has been taken to make theω → 0 limit smooth:

Kinematic, Mass and Hadronization Power Corrections
The resummed SCET cross section can be matched to full QCD such that its validity is extended beyond the peak and tail into the far tail. The usual procedure is to add in fixedorder those terms not included in the factorization theorem. For massless quarks these are usually denoted as non-singular contributions, since singular terms (that is, delta or plus functions) are fully accounted for in SCET. For massive quarks, terms not contained in the factorization theorem can be singular, and therefore these will be referred to as non-SCET.
In the far tail one sets all renormalization scales equal due to the large cancellations that take place between SCET and non-SCET terms around τ ∼ 1/3 that would be spoiled by resummation. To ensure this cancellation when including hadronization effects one usually convolves the added terms with the same shape function. As already explained in the introduction and discussed in further detail in Sec. 5.3, the fixed-order QCD prediction contains terms which are singular as τ approaches 0. The reason is that when including quark masses there are two kind of power corrections to the leadingorder EFT prediction: kinematic and massive. Both are power-counted equally in SCET and therefore, when considering the singular cross section one necessarily neglects higher powers of m. Mass power corrections can be singular, while kinematic power corrections are genuinely non-singular. It is in general desirable to absorb mass power corrections into the EFT description (although this is not strictly speaking necessary), since resummation turns log n (τ )/τ into an integrable singularity. This prescription has been adopted already for 2-jettiness in Refs. [31,32,45], and here we succinctly explain how this is implemented.
One can modify the SCET matrix elements to absorb all singular pieces. Since at treelevel all matrix elements are either 1 or Dirac delta functions, one only needs to multiply the tree-level SCET results by R C 0 to fully account for massive power corrections at this order. The O(α s ) massive corrections can be implemented modifying the hard function, that we write as such that it includes the tree-level mass modification, 15 and the jet function (at this order one cannot have mass corrections from soft dynamics). For the latter one only needs to modify J m , the mass correction with distributions, which we write as Corrections coming from B C plus are easy to implement, since they only come from the jet function (the hard function does not contain any distribution). It is convenient to define Implementing these modifications into the SCET factorization theorem with fixed scales one arrives at Since one cannot compute separately h C m and j C m , we make an ansatz and split it according to a parameter ξ This parameter reflects our lack of knowledge on the structure of mass power corrections. To estimate the associated uncertainty we vary it between 0 and 1, such that for the extreme values the correction is fully contained in the hard or jet functions. Once all singular 15 Therefore one has to rescale the non-distributional jet function J nd → J nd /R C 0 (m) as well.
corrections have been absorbed into the SCET matrix elements we can incorporate the truly non-singular corrections as an additive term where dσ C SCET /dτ refers to the mass-corrected (resummed) partonic SCET cross section, that is, to Eq. (4.1) with the substitutions H → H C and J m → J C m and with resummation kernels implemented. A similar strategy can be carried out to add power corrections to the bHQET cross section, but these will be discussed elsewhere.
So far we have dealt with partonic cross sections. Although for infrared-and collinearsafe observables partonic predictions are already a good description of the full result, for a precision analysis hadronization cannot be ignored. Here we will be concerned with the dominant effect of hadronization, that comes from soft dynamics and is already contained in the leading-power SCET factorization theorem. It is well known that for Qτ Λ QCD the main effect of soft hadronization is shifting the distribution to the right τ → τ − Ω 1 /Q (due to the operator product expansion or OPE), with Ω 1 a non-perturbative parameter that can be defined in terms of field-theory matrix elements. In the peak of the distribution, hadronization is more complex and must be taken into account by convolving the partonic result with a hadronic shape function F (p): The shape function has support for p > 0 and is normalized ∞ 0 dpF (p) = 1. It is strongly peaked at p ∼ Λ QCD and has an exponential tail extending towards infinity that ensures any of its moments is well defined. This behavior enforces the OPE and one has ∞ 0 dp pF (p) = Ω 1 . As discussed in Ref. [69], Ω 1 is afflicted by an u = 1/2 renormalon that can be removed with appropriate subtractions defined on the partonic soft function [70]. Since at the order we are working these effects are not yet relevant they will not be discussed any longer.
So far we have presented all our results in the pole scheme for the heavy quark mass m ≡ m p . Expressing our cross section in the MS scheme is trivial at this order, since for Pand E-scheme thrust there is no mass dependence at lowest order except in R C 0 . Therefore one simply has to substitute m → m(µ) everywhere the mass appears: jet function, mass power corrections and fixed-order kinematic power corrections, and add an α s correction from the conversion of R C 0 to the MS scheme. We associate the MS mas renormalization scale to µ J since the jet function is the main responsible for mass effects at the order we are working.

Numerical analysis
We have implemented our result for the cross section as given in Eqs. (9.6) and (9.7) in two independent numerical codes, which use Mathematica [71] and Python [72], respectively, that agree with each other within more than 8 significant digits. For the evaluation of dilogarithms Li 2 , hypergeometric 2 F 1 , 3 F 2 and polygamma functions ψ (n) , as well as for interpolations and numerical integration in Python we use the scipy module [73], that builds on the numpy package [74], which is also used for numerical constants such as π or γ E . In Mathematica we simply use built-in native functions.
While for the partonic SCET cross section all ingredients are analytic, the partonic non-SCET cross section is only known numerically through the results of Ref. [36]. The algorithm used in that article allows to determine the fixed-order cross section with highprecision, and in practice numerical errors are negligibly small. Our strategy to parametrize the fixed-order cross section is based on a combination of fit functions and interpolations. In a first step we make the curve less divergent at threshold by subtracting the known singular structures. This leaves integrable log-type singularities that cannot be described with an interpolation. To make the curve smoother as m → 0 we also subtract the SCET non-distributional contribution. We split this subtracted cross section into two regions that meet at τ = τ lim = 0.0016661, and use a fit function below τ lim and an interpolation above, constructed in a way such that the curve is smooth at the junction. The fit function is the sum of a term linear in log(τ ) multiplied by a degree-7 polynomial in τ plus a second polynomial of the same degree. The logarithm contains the expected behavior of nonsingular terms. The coefficients of these two polynomials are functions of the reduced mass, and each one of them is parametrized with a fit function ofm, which again consists on the sum of a 10 th -degree polynomial inm, and log(m) times another polynomial of identical degree. For the interpolation we take an evenly spaced grid with 0.033 bin-size for τ < 0.026 and a finer grid below. While the values of τ in the grid do not depend on the mass, the height of each node does, and we use fit functions ofm to parametrize this dependence with the same functional form used for τ < τ min : a logarithm of the mass times a polynomial of degree 10, plus another polynomial of the same order. To code this parametrization in Python in a way which is flexible and efficient we use object-oriented programming.
The convolution of the (now fully analytic) partonic cross sections with the shape function is performed numerically. To ensure that resummation is properly implemented in the peak and tail of the distribution, being smoothly switched away in the multi-tail, we employ the profile functions introduced in Ref. [14]. It is reasonable to think that the presence of a non-zero quark mass should modify the profile functions, but since we consider here physical situations in which the mass is still small, we stick to mass-independent parametrizations. More sophisticated profile parametrizations depending on the value of m were employed e.g. in Ref. [45]. All plots and analyses carried out in the rest of this section use the MS scheme for the heavy quark mass. Furthermore, we do not implement gap subtractions since they are not very relevant when matrix elements are used at the one-loop level. Unless stated otherwise, we take n = 4 massless quarks and a massive bottom with m b (m b ) = 4.2 GeV. For the strong coupling we use 4-loop running with the boundary condition α We start our numerical discussion by analyzing the size of each term in Fig. 7, which shows differential cross sections for Q = 20 GeV and 40 GeV for vector and axial-vector currents. We use only the default parameters for the profiles and set the parameter ξ defined in Eq. (9.5) to its canonical value 0.5. We split the distribution as follows (to alleviate notation, in the remainder of this section we drop the superscript C that indicates the current): where each term contains hadronization power corrections computed as a convolution with the same shape function. In the first equality we split the full cross section (shown as a black solid line) in SCET and non-SCET contributions, shown in magenta and green, respectively. The SCET cross section can be further divided into the sum of singular dσ sing m /dτ (shown as a red solid line) and non-distributional dσ nd m /dτ (in solid blue) contributions. In our setup, the singular cross section is defined as the contribution from terms in the SCET factorization theorem which are singular at threshold if no resummation is implemented. At N 2 LL, these  correspond to the distributions that arise from the hard function, the 1-loop soft function, and the J m=0 , J m pieces of the one-loop jet function, with the modifications discussed in Sec. 9 to absorb the relevant mass corrections, integrated against the resummation kernels. The non-distributional terms (shown in blue) are defined as the resummed contribution from J nd defined in Eq. (4.10). We observe that while the singular contribution is positive, the non-distributional is negative, and they significantly cancel each other when added together. The singular distribution can be cast as the sum of the massless approximation dσ sing m=0 /dτ (shown as a dashed gray line) and singular massive corrections dσ sing m /dτ (cyan solid line). The massless approximation is quite close to the SCET cross section (specially for the vector current), as expected, since the P-scheme decreases the sensitivity to the quark mass, and the singular massive corrections are very similar to the non-distributional term up to a global sign. We define the SCET massive corrections dσ SCET m /dτ (pink solid line) as the sum of the singular massive corrections and the non-distributional terms, which turns out to be rather small, specially for larger values of the center-of-mass energy. The non-SCET cross section has been defined in Eq. (9.6) and contains non-distributional kinematic corrections coming from the QCD fixed-order cross section. Interestingly, once we absorb all singular terms into the SCET factorization theorem, the non-SCET corrections are absolutely negligible everywhere except in the far tail. Within the setup defined in Sec. 9 only the non-distributional cross section and the massless approximation is the same for vector and axial-vector currents. We study next the convergence of the (resummed) perturbative series for the differential cross section. To that end, we generate bands randomly modifying our profile functions via a flat scan on their parameters, varying them within the ranges specified in Ref. [14], with the exception of the non-singular scale, for which we use the following continuous variation with −1 ≤ n s ≤ 1. In our scan we also randomly vary ξ between 0 and 1. In Fig. 8 we show the resulting perturbative bands at LL (green), NLL (blue) and N 2 LL + O(α s ) (red) for the vector and axial-vector currents, at two center-of-mass energies: Q = 40 GeV and 80 GeV. Our curves are not self-normalized, but we nevertheless observe and excellent convergence in all cases (even at low energies) in the tail of the distribution, where higher-order bands are nicely contained in lower-order ones. In the peak we see a big jump between LL and the two highest orders, and the convergence is not as good as in the tail, what might  Figure 10. Difference between the vector and axial-vector differential cross sections normalized to the average of the two currents. We show results for Q = 20 GeV (red), Q = 40 GeV (blue) and Q = 80 GeV (green).
indicate that the parameters affecting mainly the peak should be varied in wider ranges. A careful inspection of the error bands reveals that the relative uncertainties for LL and NLL are nearly identical in the whole range, and both monotonically increase as τ grows: at Q = 40 [80] GeV they change from 36 [45]% at τ = 0.07 to 84 [80]% at τ = 0.28. On the other hand, at N 2 LL the relative uncertainty is completely flat between 0.07 ≤ τ ≤ 0.3, and smaller than the two lower orders: 36 [30]% for Q = 40 [80] GeV. We observe the same relative uncertainties for vector and axial-vector currents.
In Fig. 9 we compare the 2-jettiness (blue) and thrust (red) cross sections for massive quarks produced through the vector current at various center of mass energies, as indicated in the caption of the plot. As a reference, we also show in green the massless cross section. We observe that the 2-jettiness cross section has a negative deep which becomes more pronounced at low energies. It is produced by large logarithms which could be summed up by matching SCET to bHQET, as discussed in Sec. 4.2. While the massless cross section is always quite similar to massive P-scheme thrust, the 2-jettiness distribution gets quite different at low energies, with a higher peak shifted to the right. We will study this behavior in further detail later in this section. As energies become larger, the three cross sections become similar to one another, but P-scheme thrust is always closer to the massless result. In fig. 10 we plot twice the difference between the vector and axial-vector currents, normalized to the sum of the currents. To make the figure clearer, we use a logarithmic scale on the y axis. We observe, as expected, that for larger energies the difference becomes smaller, since both currents approach the (current-independent) massless result.
In our last analysis we study the dependence of the peak position and peak height with the heavy quark mass. Since the peak position retains some dependence on Q from soft hadronization, we fix the value of the center-of-mass energy to 40 GeV, such that we can make sure the peak moves only due to changes in the mass. In this case we compare the results for thrust and 2-jettiness, since the former is relatively mass insensitive while the latter has been designed to measure the top quark mass in future linear colliders, see e.g. Ref. [37]. We restrict the values of the bottom quark below m b = 14 GeV to make sure we can still apply SCET and scenario II, which should be described using bHQET, is unimportant. The results of our study are summarized in Fig. 11, where one  Figure 11. Peak position (a) and peak height (b) for 2-jettiness (blue) and P-scheme thrust (red) massive cross section. Results correspond to default profiles, vector current and a center-of-mass energy of 40 GeV, and with m b ≡ m b (m b ). We vary the bottom mass between 0 and 14 GeV, such that SCET still applies.
can clearly observe a flat behavior for P-scheme and an obvious quadratic dependence for 2-jettiness. For the latter this is nothing but expected, since the peak position is shifted by In fact, if we perform a fit to the 2-jettiness peak position we find τ max 0.0255 + 1.75m 2 , which follows almost exactly the blue line in Fig. 11(a) and is in fair agreement with our expectations [the small disagreement is expected since the peak position should be computed with m b (µ J ) and not with m b (m b )]. The dependence of the peak height on the bottom mass is also much larger in jettiness than P-scheme thrust.

Conclusions
When considering heavy quarks in the context of event shapes, depending on the scheme used in their definition the mass sensitivity of the cross section can vary significantly. This sensitivity manifests itself already at lowest order by setting the threshold to a non-zero value, and of course increases as the mass grows. This shifted threshold comes solely from the jet function.
While in a recent paper we discussed how to compute these distributions in fixedorder at NLO, in this article we have shown how to analytically compute the differential and cumulative cross sections in the E-and P-schemes at N 2 LL + O(α s ) accuracy in SCET and bHQET. To achieve this goal, we have calculated the missing pieces, namely the NLO jet function in those two effective field theories. We have shown that in the collinear limit the heavy quark momenta expressed in the E-and P-schemes coincide, but are different from the original (massive) definition. This entails that for any event shape, the jet function will be identical in the former two schemes, but in the case of thrust, heavy jet mass and C-parameter the measurement function is no longer completely inclusive, meaning that one needs to compute the jet function with cut diagrams, integrating over phase space rather than loop momenta. We provide an optimized and compact form for the jet function definition in each EFT, written in terms of quantum and kinematic operators, that facilitates their computation. For the bHQET jet function we explain how to rescale the integrated light-particle momenta such that the heavy quark mass drops out before any integration is carried out. In the computation of the P-scheme SCET jet function one needs to use either sector decomposition or hypergeometric function identities to properly expand in ε and extract the distributions that appear in this limit.
Shifting its argument to relocate the threshold back to zero, the 1-loop SCET massive jet functions can be written as the sum of the massless jet function plus mass corrections. The latter can be further divided into terms with distributions (or singular) and terms without distributions. While carrying out resummation for the former is already well known, the terms with regular functions need to be treated in a case-by-case basis. In this article we show how to analytically RG evolve the non-distributional terms for 2-jettiness and P-scheme thrust, and derive rapidly-convergent expansions that can be carried out around the heavy-quark (p 2 = 0), "threshold" (p 2 = m 2 ) [ only for P-and E-schemes ] and massless (p 2 = ∞) limits, which nicely overlap with one another such that in numerical implementations there is no need to explicitly evaluate hypergeometric functions at all (for 2-jettiness there is a small region around p 2 = m 2 for which one cannot use expansions). Our expansions can be carried out up to any order such that the result is also arbitrarily precise. This is much faster than a direct evaluation of 2 F 1 and 3 F 2 functions, which were the bottleneck of the analysis carried out in Ref. [45].
We show how to absorb into the SCET factorization theorem those mass-suppressed singular terms that appear in fixed-order corrections by a suitable redefinition of the hard and jet functions. After this procedure is carried out, we complete our resummed expression with purely kinematic corrections (which are now entirely non-singular), which become relevant in the far tail. Hadronization power corrections can be incorporated in the usual way by convolving with a shape function, and with this complete description we have performed some numerical investigations. We have shown that there are strong cancellations taking place between the two types of massive corrections to the SCET factorization theorem (with or without distributions) everywhere except in the peak, and that the remaining nonsingular corrections are immaterial everywhere except in the far tail. The cancellations are stronger at larger energies, where also vector and axial-vector currents yield similar results. We have demonstrated that the P-scheme thrust cross section is much closer to the massless prediction than for 2-jettiness by comparing cross sections as well as investigating the peak position and height as a function of the heavy quark mass. We have also observed a nice convergence of the cross section when adding perturbative orders.