Boosted Top Quarks in the Peak Region with N$^3$LL Resummation

We present results for the 2-jettiness differential distribution for boosted top quark pairs produced in $e^+e^-$ collisions in the peak region accounting for QCD large-logarithm resummation at next-to-next-to-next-to-leading logarithmic (N$^3$LL) order and fixed-order corrections to matrix elements at next-to-next-to-leading order (NNLO) calculated in the framework of soft-collinear effective theory and boosted heavy quark effective theory. Electroweak and finite-width effects are included at leading order. We study the perturbative convergence of the cross section in the pole and MSR mass schemes, with and without soft gap subtractions. We find that there is a partial cancellation between the pole mass and soft function renormalons. When renormalon subtractions concerning the top mass and the soft function are implemented, the perturbative uncertainties are, however, systematically smaller and an improvement in the stability of the peak position is observed. We find that the top MSR mass may be determined with perturbative uncertainties well below $100$\,MeV from the peak position of the 2-jettiness distribution. This result has important applications for Monte Carlo top quark mass calibrations.


I. INTRODUCTION
The top quark mass m t is one of the most important parameters of the Standard Model (SM). In conjunction with the Higgs boson mass, it is an essential input for studies of the stability of the SM electroweak vacuum [1][2][3][4][5][6], and it plays an important role in precision electroweak fits [7]. The most precise determinations of the top mass to-date come from so-called "direct measurements", that are based on the kinematic reconstruction of the final-state top quark decay products and the comparison of the resulting kinematic distributions with parton-shower Monte Carlo (MC) simulations. The current world average for direct measurements reads m MC t = 172.76 ± 0.30 GeV [8] and projections for the HL-LHC indicate that uncertainties as small as 200 MeV can be reached for individual measurements [9].
The interpretation of these measurements are, however, (as reviewed below) impacted by an additional ambiguity from the lack of understanding of the field theoretic meaning of the top mass parameter encoded in the MC event generators [9][10][11][12]. This ambiguity is not yet precisely quantified and should be considered at the GeV level, i.e. it is comparable to the uncertainties quoted by the experimental analyses [13]. Carrying out firstprinciple theoretical predictions of kinematic distributions that exhibit high sensitivity to the top mass is a challenging program in the light of the disparate energy scales that enter the top production and the measure-ments on the decay products. These lead to large logarithms of ratios of these scales that require resummation. Furthermore, due to nonperturbative corrections it is necessary to take into account the hadronic nature of the final state. While these effects can be accounted for via MC event generators, it comes at the cost of the limited perturbative [14,15] and conceptual precision of the MC description of perturbative and nonperturbative effects. This is the origin of the interpretation problem of the top mass parameter employed in the MC generators (see [13] for further discussion).
In order to achieve a precise top quark mass determination we therefore need an observable that has the required kinematic sensitivity and is theoretically tractable, such that (a) it can be reliably calculated in perturbation theory within a specific short-distance top quark mass scheme, and (b) nonperturbative effects can be consistently quantified from a field theory perspective. Such calculations can be carried out in the framework of effective field theories (EFTs) that are systematically improvable in their power counting expansion, as well as in perturbation theory where resummation of large logarithms now reaches next-to-next-to-next-to leading logarithmic (N 3 LL) accuracy for a number of applications [16][17][18][19]. Moreover, the framework of EFTs offers ways to rigorously describe and quantify nonperturbative effects due to hadronization [20][21][22][23].
Therefore, for such an observable, EFT based calculations may naturally incorporate all the distinct features of parton-shower MC simulations while retaining a sys-tematic connection to field theory and thus theoretical control. Hence, EFTs offer promising prospects for precision collider physics and for developing diagnostic tools for improving MC simulations [23,24]. Such a framework has recently been applied to the broad program of top quark mass measurements, including a factorized description of a hadron-level differential top jet mass spectrum for jets initiated by boosted top quarks in the peak region at a future e + e − collider [25,26], boosted top jets with soft drop grooming at the LHC [15], a calibration of the MC top quark mass parameter [27] (based on the work of Refs. [25,26]), and parton-level studies of the correlation of the MC top quark mass parameter with the parton shower evolution cutoff [15]. The sensitivity of event-shape variable definitions with respect to quark mass effects has also been recently studied in Refs. [28,29] in the context of fixed-order and resummed perturbation theory.
In this work we continue this effort by extending the perturbative calculations of boosted tops in the peak region at an e + e − collider from the N 2 LL resummation of logarithms and O(α s ) fixed-order matrix elements used in [27], to N 3 LL resummation with NNLO [ O(α 2 s ) ] fixedorder matrix elements.

A. Status of top mass measurements
Recent direct top quark mass measurements have yielded the results m MC t = 172.26±0.61 GeV (CMS) [30], m MC t = 172.69±0.48 GeV (ATLAS) [31] at the LHC, and m MC t = 174.34 ± 0.64 GeV [32] at the Tevatron. The top mass superscript MC signifies that the direct measurements extract the top quark mass parameter coded in the MC generators used for the analyses. The lack of understanding of a precise field-theoretic definition of m MC t (and hence its relation to short-distance masses defined in the context of quantum filed theory, which is the preferred input parameter for high-precision theoretical predictions) results in an additional conceptual uncertainty in how m MC t should be related to the pole mass m pole t or a short-distance mass such as the MS mass m (6) t (µ) or the MSR mass m MSR t (R) [33][34][35]. This ambiguity is not included in the individual quoted experimental uncertainties but should be considered to be (at least) comparable [13]. Analyses that can shed light on a more complete and quantitative understanding of these issues, as well as more first-principle aspects of MC event generators, are underway [15,27,36]. Other recent related studies include analyses of the theoretical limitations concerning the modeling of the dynamics in the top quark production and decay [37], finite lifetime [38], hadronization effects [39] and observable infrared sensitivity [40].
The so-called top quark pole-mass measurements based on the total cross section [41], for which precise theoretical predictions expressed in terms of the pole mass renormalization scheme have been employed, have generally yielded results with larger uncertainties than the direct measurements: m pole t = 172.9 +2.5 −2.6 GeV (ATLAS) [42] and m pole t = 172.7 +2. 4 −2.7 GeV (CMS) [43]. The relatively large errors result from the uncertainty in the normalization of the inclusive cross section (both in the measurement and the theoretical prediction) and its relatively weak dependence on m t . Employing differential cross sections such as leptonic distributions, and using matched NLO+PS (parton shower) MC simulations results in an enhanced top mass sensitivity [44]. Such differential measurements have been included in the world average m pole t = 172.4 ± 0.7 GeV [8] for the pole mass measurements which is in good agreement with the corresponding world average of direct measurements mentioned above. Recent precision measurements, also using matched NLO+PS MC simulations for the theoretical predictions, accounting for tt+jet final states [45] (m pole t = 171.1 +1.2 −1.0 GeV [46]) and the reconstructed top-antitop invariant mass distribution (m pole t = 170.5 ± 0.8 GeV from a simultaneous α s fit [47]) have comparable uncertainties but have, however, resulted in significantly lower m t values, posing some tension that may partly arise from missing theoretical input in predictions for the corresponding differential cross sections [13,48,49].
It is well known that the pole mass additionally suffers from a conceptual O(Λ QCD ) renormalon ambiguity which in dedicated analyses was estimated to amount to 110 MeV in Ref. [50] and to 250 MeV in Ref. [35]. 1 This ambiguity does not represent an uncertainty due to perturbative truncation or limited information, but signifies the principle conceptual imprecision in assigning a definite value to m pole t . The ambiguity is not related to any physical effect, but inherent to the unphysical nature of the pole mass renormalization condition. It can therefore be avoided by expressing cross sections obtained in perturbation theory in terms of a short-distance mass scheme at an appropriate renormalization scale, which can also improve the overall convergence of the perturbative series at the first few orders. As far as the spread of the above-mentioned recent pole mass measurements is concerned, the pole mass renormalon ambiguity likely plays no role, because it is smaller than the quoted uncertainties of these measurements. For the direct top quark mass measurements, the pole mass renormalon problem has been discussed intensely in the context of the frequently used approach of identifying m MC t and m pole t . However, in Ref. [15] it was shown analytically for the simple case of the 2-jettiness distribution in e + e − → tt+X, that the quark mass parameter associated to a NLL-precise parton shower is in general not the pole mass, but a low-scale short distance mass that depends on the value of the shower cut and may differ from the pole mass by an amount larger than the pole mass ambiguity. 2 So the problem of how to properly interpret the MC top quark mass m MC t in terms of a well-defined and ambiguity-free field theory mass is not related to the pole mass renormalon ambiguity, but to the limited theoretical precision of the state-of-the-art MC event generators and to ignorance concerning MC systematics.
Motivated by the interpretation issues of the direct top mass measurements and the still large uncertainties in the pole mass measurements from inclusive and differential cross sections, a number of alternative methods to measure m t have been proposed, which are based on differential cross sections with respect to variables constructed from top decay products exhibiting strong kinematic top mass sensitivity. The observables these analyses are based on include the M T 2 variable and variants of it [51,52], the lepton b-jet invariant mass [53], the shape of b-jet and B-meson energy distributions [54], and the J/ψ and lepton invariant masses [55,56]. Conceptually, these observables are governed by parton shower dynamics as well as various nonperturbative effects, in a way analogous to the direct reconstruction method (albeit with differing systematics). This is because they are also based on the idea of assessing the kinematics of decaying (colored) top quark particles through simulations obtained from parton-shower MCs. Their reliance on these MCs further makes their potential extension to higher logarithmic precision nontrivial, as the theoretical precision of parton showers is quite observable-dependent. Additionally, the nonperturbative effects of hadronization are accounted for through a multi-parameter MC hadronization model. Here, since a systematic way to quantify the intrinsic MC hadronization uncertainty does not yet exist, uncertainties are typically being estimated by comparing predictions based on different models.

B. Top mass determination using effective field theories
In Refs. [25,26] a framework of EFTs was developed to describe boosted top quarks in the peak region at a future e + e − collider. Using soft-collinear effective theory (SCET) [57][58][59][60][61] a factorization theorem for the double-differential hemisphere-mass cross section in e + e − → tt + X was derived in the boosted top quark limit, with center of mass energy E cm = Q m t . The two invariant masses M t and Mt are defined using all particles in each of the two hemispheres that are determined by the event's thrust axis described below. The peak refers to the region where the M t -Mt double differential distribution exhibits the hemisphere mass top and antitop resonances. One key feature of this factorization formula is that the most important hadronization effects are parametrized via a convolution with a top-mass-and Q-independent nonperturbative shape function, which is field theoretically defined from a vacuum matrix element of Wilson-lines. The universality of the factorization formula states that, at least in principle, data for massless dijet events (obtained even at past experiments such at LEP) could be used to fix the nonperturbative shape function and to make the analysis independent of estimates of nonperturbative corrections obtained from MC event generators. The other key feature is that the dependence on the renormalization scheme of the top quark mass is fully controlled so that one can make and test predictions in any scheme to the extent that higher order perturbative corrections are incorporated. In the tail of the distribution an operator expansion can be applied such that the nonperturbative corrections are dominated by a single parameter.
The dijet hemisphere mass cross section exhibits a clear peak at the top and antitop resonances which are directly sensitive to the value of the top quark mass. However, the resonance location in M t and Mt is not directly at the top quark mass due to radiative effects related to ultra-collinear and large-angle soft radiation, as well as hadronization and finite width effects. The ultra-collinear radiation refers to radiation that is soft in the (anti)top quark rest frame, but becomes collinear due to the (anti)top quark boost. In the factorization theorem this ultra-collinear radiation is described by boosted versions of heavy quark effective theory (HQET) [62][63][64][65][66] called boosted HQET or simply bHQET. 3 The typical scales of the ultra-collinear radiation in the peak region range between the top width Γ t and ∼ 10 GeV, which also quantifies the typical 'off-shellness' ∼ (q 2 − m 2 t )/m t of the decaying top quarks. The top mass dependence of the location and shape of the observable peak, along with the top quark scheme dependence, is specified by ultracollinear radiation effects which can be calculated perturbatively. Furthermore, the leading-order electroweak effects come from the top width which is fully encoded in Breit-Wigner-modified top quark propagators. The large-angle soft radiation is only sensitive to the collinear top quark color flow (which is fully taken over by the top decay products in the boosted limit) and describes soft momentum exchange between the two hemispheres. It is governed by scales below the top quark width, and its nonperturbative contributions constitute the shape function. Due to the boost of the top quarks, the effect of large-angle soft radiation on the peak locations is enhanced by a factor Q/m t . 3 Here the letter "b" in bHQET stands for the fact that the two HQET theories for the top and antitop quarks are boosted in opposite directions.
The above-mentioned framework was employed in Ref. [27], where a calibration of the m MC t parameter in Pythia 8.205 was carried out using the N 2 LL + O(α s ) prediction for the 2-jettiness event shape in e + e − → tt + X collisions defined as where the sum runs over all produced particles i in the event andn t is the thrust axis that maximizes the sum in Eq. (1). Since a sum over all final-state momenta is involved in Eq. (1), we restrict ourselves to hadronically decaying top and anti-top jets, which is accomplished by simply including the corresponding branching fraction in the Born cross section σ 0 in the factorization formula discussed below. In the limit τ 2 − 2m 2 t /Q 2 1, i.e. when focusing on the peak region, the event shape corresponds (up to power corrections) to the sum of the hemisphere invariant masses, such that and thus has the same kinematic sensitivity to the top mass as the double differential cross section considered in Refs. [25,26]. Being a global event shape, the τ 2 differential cross section is furthermore free from non-global logarithms [67,68]. The results of the calibration carried out in Ref. [27] showed that m MC t in Pythia 8.205 cannot be simply identified with the pole mass, but is numerically close to the MSR mass m MSR t (1 GeV). Compatible numerical results were obtained in the analysis of jets from top quarks produced in pp collisions with soft drop grooming in Ref. [36] based on theoretical calculations at NLL order.
In this work we improve the calculation of the 2jettiness event shape cross section in e + e − → tt + X by including NNLO fixed-order matrix elements with N 3 LL resummation [ referred to as N 3 LL + O(α 2 s ) ]. By now, all the ingredients that enter the factorization formula in Eq. (3) below are individually known to the accuracy needed to achieve this precision for the 2-jettiness cross section, in particular the two-loop heavy quark jet function [69], full two-loop thrust soft function [70,71], and the N 3 LL result for the Wilson coefficient for the matching at the top mass scale [72]. We also include the recently calculated analytic 4-loop result for the cusp anomalous dimension [73], which despite having a tiny numerical impact, is an important formal ingredient to obtain N 3 LL accuracy. We consistently combine these ingredients to obtain the N 3 LL + O(α 2 s ) boosted top cross section fully analytically, and study the convergence of the resummed perturbation theory. We exclusively consider the e + e − 2-jettiness cross section in the boosted and bHQET limits where power corrections in m t /Q and in the top quarks' off-shellness ∼ (q 2 − m 2 t )/m t are neglected. These corrections, even though there are formally power-suppressed, can be non-negligible for phe-nomenological analyses in the peak region. With inclusion of these power corrections, the results obtained in this article will serve to improve analyses such as the top quark mass calibration carried out in Ref. [27], which we leave to future work.
The outline of the paper is as follows: We first present the factorization formula in Sec. II describing its key features and elements. We describe in Sec. III the implementation of the cross section in terms of short distance mass schemes and a renormalon-free soft function in Sec. IV. The setup of observable-dependent renormalization scales is discussed in Sec. VI. Finally, in Sec. VII we combine all the pieces to calculate the resummed cross section and study its perturbative convergence when the renormalons in the pole mass and the soft function are either subtracted or left unsubtracted. Using an analysis of the 2-jettiness peak location we draw conclusions on the perturbative uncertainties of a determination of the MSR mass and the pole mass. In the appendices we review and state additional details for the various ingredients that are needed in this analysis. We conclude in Sec. VIII.

II. FACTORIZATION THEOREM IN THE PEAK REGION
In Ref. [25] two factorized expressions for the differential cross section were derived that are valid in the peak and tail regions of the double hemisphere mass, or equivalently 2-jettiness, distribution. In the tail region the fluctuations in the top mass can be large, such that M 2 t,t − m 2 t ∼ m 2 t , and a factorization formula based on SCET with massive particles can be derived. On the other hand, in the peak region, the off-shellness of the decaying top is constrained such that M 2 t,t − m 2 t m 2 t , necessitating an additional factorization and a resummation accounting for the top width as an additional relevant scale that is carried out in the bHQET framework. As already mentioned, in this work we focus on the factorization in the peak region accounting exclusively for the bHQET contributions, leaving off-shellness power corrections described in the SCET factorization (that become essential in the tail region) and the inclusion of m t /Q SCET power corrections to future work. The factorization formula in the peak region is given by where we have the boost parameter given by and have defined the off-shellness variableŝ τ as 4 Eq. (3) involves various perturbative ingredients, including an evolved matching function H (5,6) evol , jet and soft functions J (5) B,τ2 andŜ (5) τ2 , and evolution kernels U (5) B and U (5) S . It also includes a non-perturbative shape function F , whose independence from other parameters is a prediction of the factorization theorem. These ingredients will be discussed in detail in subsections below. Eq. (3) applies in the peak region whereŝ τ ∼ Γ m t , with Γ 2Γ t being the effective width of the distribution broadened by radiation, hadronization as well as finite width effects. In this region, where the τ 2 distribution exhibits a resonance, the 2-jettiness variable is (up to power corrections) directly related to the sum of the squared hemisphere masses defined with respect to the thrust axis, as given in Eq. (2). It is therefore convenient to define the inclusive jet mass variable M J which inherits some of the features of a reconstructed top invariant mass, albeit being based on a hemisphere top jet. The M J -distribution peaks close to m t , but is in addition affected (with respect to peak position as well as the width of the observed peak resonance) by large-angle soft radiation exchanged between the two hemispheres. The widening of the peak due to top-decay width and soft QCD effects, however, does not affect the kinematic sensitivity of the observable, and normalizing the M J distribution enables the uncertainties in the M J peak location to be taken as a direct measure for the uncertainties in the associated top mass determination. The factorization formula separates perturbative contributions from the hard local interactions involving the scales Q and m t (encoded in the hard factor H (5,6) evol ), dynamical effects associated to large-angle soft radiation (accounted for in the soft function S (5) τ2 ), and dynamical effects due to ultra-collinear radiation (contained in the jet function J B,τ2 incorporates the leading-order effects due to the top quark width Γ t and also carries, because of its peaked structure, the main top quark mass sensitivity of the τ 2 distribution. This allows testing at high precision the impact of either using the pole mass scheme m pole t or a suitable shortdistance mass m sd t . This is indicated by the argument is the perturbative series for the difference between the pole and the adopted short-distance masses. The choice δm = 0 implies the use of the pole mass scheme. We note that there is also top quark mass (scheme) dependence in the hard factor H (5,6) evol , indicated by the argument m t , which, however, only affects the normalization of the τ 2 distribution and is very subdominant compared to the main sensitivity to the top quark mass. The different character of the top mass dependence in the hard and the bHQET jet functions is discussed in more detail below.
For the τ 2 distribution in the peak region for boosted top quarks we have the hierarchy Q m t {ŝ τ , Γ}. This leads to large logarithms, which are resummed via renormalization group (RG) equations for the corresponding perturbative matrix elements. This resummation is implemented through the evolution factors U The choice of µ is arbitrary and the factorized prediction is (strictly) invariant under changes of this µ. In contrast the dependence on the initial scales µ i only cancels out order-by-order in resummed perturbation theory. The scale µ is typically set equal to one of the renormalization scales µ i (i = H, m, S, J) such that one of the renormalization evolution factors disappears.
The superscripts "(5)" and "(5, 6)" indicate the number of active dynamical flavors relevant for the momentum scales of the respective quantum effects, where "(5)" and "(6)" indicate scales below and above the top mass, respectively.
Due to the simple inclusive character of the 2-jettiness (or the M J ) distribution, the leading-order nonperturbative effects arise from the low-energy dynamics of the large-angle soft radiation and are encoded in the (hadronization) shape function F , which is convolved with the perturbative soft function S (5) τ2 . Physically, the shape function incorporates the leading effects of hadronization and controls the amount of nonperturbative radiation being exchanged between the two hemispheres. Even though the details of the shape function form must be determined from experimental data, the way how the shape function appears in the factorized cross section represents a very strong theoretical constraint on hadronization.
In the following we briefly review the physical aspects of all functions appearing in the factorization theorem one by one. For a detailed discussion on how the bHQET factorization theorem of Eq. (3) connects to the corresponding formula in the tail region, and on possible alternative versions to organize the renormalization group evolution, we refer to Ref. [74].

A. The Hard Function
The factorized cross section involves a two-step matching from QCD to SCET, and then from SCET to bHQET, at the scales µ H ∼ Q and µ m ∼ m t , respectively. Therefore, the difference between 6 or 5 flavors is related to the top quark being a dynamical degree of freedom or not. The resulting hard matching coefficients in Eq. (3), together with their renormalization group (RG) evolution kernels, are collectively written as The term U H Q evolves the SCET hard function H Q [75][76][77][78][79][80] from µ H ∼ Q to µ m ∼ m t and resums large logarithms of Q/m t in the cross section. U v is responsible for the evolution of the bHQET current between µ m ∼ m t and µ, which we assume is smaller than µ m , and only depends on the top quark boost factor defined in Eq. (4).
The hard matching at the top quark mass scale µ m , which encodes off-shell top quark quantum fluctations that arise in the heavy quark limit, is given by H m [26,72]. Since the matching is performed at the top quark mass, one can express H m in terms of α s with either 5 or 6 active flavors. In the numerical analysis below we choose 6. The effect of this freedom in the scheme choice is, however, tiny and numerically irrelevant. We now discuss the parameter = Q/m t [ see Eq. (4) ] appearing in the mass mode matching factor H (6) m and the bHQET current evolution kernel U (5) v . In the peak region, for the bHQET factorization treatment of the top quark dynamics, the momentum of the nearly onshell top quarks is parameterized as is the 4-velocity of the energetic (anti)top quarks and k µ is a small residual momentum accounting for the fluctuations caused by the low-energy radiation [ in the (anti)top quark rest frame ], such that k µ m t and one can expand the dynamical effects to using light-cone coordinates defined relative to the thrust axis n µ = (1, n t ) (which we take along the top direction), such that p µ = (n · p,n · p, p ⊥ ), wheren µ = (1, − n t ) is an auxiliary vector satisfying n·n = 2 and n 2 =n 2 = 0. The parameter appearing in the definition of the reference velocities is related to the top quark boost and defined in Eq. (4), such that the on-shell (anti)top 4-velocity would approach v µ tt in the boosted limit Q m t in the absence of any radiation. The mass mode matching factor H (6) m and the bHQET current evolution kernel U depend on the reference velocities v µ t and v μ t and thus on through the scalar products v t,t ·k appearing in the bHQET Feynman diagrams.
The choice of v µ t,t is ambiguous with respect to higherorder power corrections of O(k µ /m t ) [ ∼ O(Γ t /m t ) in the (anti)top rest frame ] and entails a symmetry of bHQET factorization with respect to changes of v µ t,t and thus of , which is one aspect of a more general class of symmetry transformations called reparametrization invariance [81] that connects different orders in the bHQET 1/m t expansion.
As a consequence, there is a power-suppressed freedom in the choice of the boost parameter , which in the factorization theorem of Eq.
(3) appears in the function H (5,6) evol and the momentum argument of the jet-function evolution factor U (5) B . The m t parameter appearing in is, however, not associated to any particular top mass renormalization scheme, but for the final evaluation it just has to be numerically chosen close to the invariant mass of the top quark in the resonance region. It is therefore not mandatory to re-expand the boost parameter when expressing the top quark mass in a short-distance renormalization scheme. The other consequence is, that variations → + δ with δ ∼ QΓ t /m 2 t in the peak region account for uncertainties due to the truncation of power corrections in the context of the leading power bHQET factorization. We conclude by noting that actually setting the parameter m t in equal to the pole mass artificially reintroduces the pole mass renormalon (albeit with formally power suppressed numerical effects) which cannot be cured by a re-expansion in terms of a short-distance mass, because there are no associated factorially divergent perturbative corrections in other parts of the factorization theorem. It is therefore not advisable to identify the parameter m t in as the pole mass.
The top mass parameter m t that appears in the argument of H (6) m in Eq. (8) is identical to the mass used in the threshold decoupling relation of the strong coupling when transitioning between 6 and 5 flavors. This m t can be expressed in either the pole or MS scheme (with the same mass scheme employed in the α (6) s ↔ α (5) s decoupling relation). Since H m only involves logarithms of the ratio µ m /m t (in addition to logs of ) it is also insensitive to the O(Γ t ) fluctuations at the top threshold that are suppressed by O(Γ t /m t ). The difference between the choice of pole or MS masses is accordingly found to be also numerically subleading.

B. The Ultra-Collinear Sector
Upon integrating out the off-shell (small) components of the top quark field in bHQET, the leading-order dynamics of the remaining ultra-collinear fluctuations along the top and anti-top quark directions is captured by the 2-jettiness bHQET jet function, defined as Here "ultra" distinguishes the collinear modes in the peak regionŝ τ ∼ > 2Γ t [ which are soft in the (anti)top quark rest frame ] from higher virtuality collinear modes appearing in the tail regionŝ τ ∼ m t , where bHQET off-shellness power corrections become large and SCET provides the adequate description for collinear radiation. Here J B (ŝ t , δm t , Γ t , µ B ) is the familiar bHQET jet function that appears in di-hemisphere mass and jet mass distributions [15,26,69]. At leading order, the finite top width effects of the bHQET jet function J where the factor 2Γ t arises from accounting for the widths from both top and antitop quarks. This implies that the measurement on tt final states is fully inclusive in the decay products as well as any radiation from them. The consistency of this inclusive treatment is ensured by considering the boosted limit Q m t where the decay products from the top and anti-top quarks are collimated back-to-back along the direction of the thrust axis in distinct hemispheres. Both the boosted top decay products as well as gluon and light-quark radiation encoded in J B have typical angles ∼ 2m t /Q relative to these axes. 6 The variableŝ t appearing in the (single) hemisphere jet functions J where k is the total residual momentum of the collimated system (after removing the contribution from the top quark mass). This captures the invariant mass of the decaying top quark system together with its ultra-collinear radiation, up to terms of O(ŝ 2 t /m t ) [25]. In the peak region, the ultracollinear fluctuations have virtuality of has support only for non-negativeŝ τ and equals δ(ŝ τ )/m t at tree-level. It is this threshold behavior which causes the strong top mass sensitivity of the τ 2 distribution in the peak region. We emphasize, however, that while the partonic threshold is atŝ τ = 0, the observable peak position exhibited by the entire factorization theorem of Eq. (3) is determined coherently from the effects of the ultra-collinear and large-angle soft radiation together with the Breit-Wigner smearing. The well-known pole mass renormalon problem arises from higher order perturbative corrections in J (5) B,τ2 in the pole mass scheme, and is encoded in the size of the coefficients of plus distributions inŝ τ . The pole mass renormalon can be remedied by using, instead of the pole mass m pole t , a suitable short-distance mass scheme m sd t in the definition ofŝ τ in Eq. (5). Because the bHQET jet functions are defined strictly at leading order in the 1/m t expansion, switching the top mass renormalization scheme requires that one also accounts for the mass scheme correction strictly to leading order in the 1/m t expansion. This leads to the generic form where the B ij are constant coefficients. This yields the 2-jettiness bHQET jet function for stable top quarks in an arbitrary (short-distance) mass scheme m sd t , where it is still strictly mandatory to expand the dependence on δm consistently in powers of the strong coupling α (5) s (µ B ) such that the pole mass renormalon cancels order by order. Setting δm = 0, one recovers the corresponding result in the pole mass scheme, which was calculated up to O(α 2 s ) in Refs. [26,69].
x + is the standard plus function distribution with a vanishing integral over the range x ∈ [0, 1] for k ≥ 0 and From the expression in Eq. (13) we can also clearly see that the bHQET power counting requires that the perburbative series for the mass scheme correction δm(R) obeys the scaling δm ∼ŝ τ m t . This shows from a power-counting point of view why low-scale shortdistance masses have to be employed and the MS mass, which has δm ∼ m t , is forbidden. Such low-scale short-distance masses always involve an infrared subtraction scale R, which is necessary to eliminate the large O(Λ QCD ) renormalon corrections appearing in the pole mass scheme [33]. To avoid upsetting the bHQET jet function power counting it is important that the series coefficients defining the low-scale short-distance mass in Eq. (12) have the property δm i (R) ∝ R, where the infrared scale R is parametrically close to the typical ultracollinear scale, i.e. R ∼ µ B ∼ŝ τ . In our analysis we em-ploy the MSR mass m MSR t (R) [33][34][35]83] that satisfies this requirement as explained in more detail in Sec. III.
We finally note that the overall factor (1/m t ) 2 appearing in the generic expression for the 2-jettiness bHQET jet function J (5) B,τ2 in Eq. (13) arises from Eq. (10) since each J B has a 1/m t factor [26]. The renomalization scheme dependence of this factor is power suppressed in the leading order bHQET expansion and therefore in principle not specified at the level of the factorization theorem of Eq. (3).

C. The Soft Sector and Nonperturbative Effects
The thrust partonic soft functionŜ (5) τ2 accounts for the effects of large-angle soft radiation with respect to the thrust axis [16,17,70,71,84,85]. In the context of inclusive observables, such as the global event shape 2-jettiness, it accounts for the so-called ultrasoft modes, which have scale fluctuations of the order of µ S ∼ŝ τ / , parametrically smaller than the typical scale µ B for the ultra-collinear modes. Thus µ S represents the smallest perturbative scale relevant for the τ 2 distribution.
For global event shapes such as τ 2 or jet-based observables without jet grooming, the leading nonperturbative effects from hadronization can be described via convolution with a shape function F (k − 2∆), where F (k) has support for k ≥ 0, peaks at k ∼ Λ QCD [21,86] and is by definition normalized to unity, ∞ 0 dkF (k) = 1. Here, ∆ is a model parameter which accounts for the minimum hadronic energy deposit in each hemisphere (hence the factor of 2), referred to as the 'gap'. In the tail region where ∼ µ S Λ QCD , one can expand for large and the most important nonperturbative effect is encoded in the first moment Ω 1 of the shape function, where Ω 1 can also be expressed as a vacuum matrix element of soft Wilson lines [20,22], and the ellipses represent higher order terms in the expansion. The concrete form of the shape function we use for the generic numerical examination carried out in Sec. VII is provided in App. G.
As we move further into the peak region, µ S decreases and eventually approaches the nonperturbative scale Λ QCD . Thus in the peak region µ S ∼ 1 GeV > Λ QCD and the effects of the shape function have to be accounted for exactly in terms of the convolution of Eq. (14). Even though this in principle implies that an infinite amount of information could be required to fix the analytic form of the shape function, the normalization and the requirement that all moments of F exist, together with the fact that factorization strictly demands convolution, allow us to reliably constrain the form of the shape function in terms of a few parameters by means of an expansion in optimally designed basis functions [87], see e.g. [27] and [88]. In practice, determining the first moment of the shape function fixes the bulk of the information encoded in it for the whole τ 2 spectrum.
The factorization into partonic soft and nonperturbative shape functions displayed in Eq. (14) depends on the regularization and renormalization schemes that are employed for the computatation of the partonic soft function. The argument 'δ = 0' shown in Eq. (14) stands for the standard MS scheme. As was shown in Ref. [21], these prescriptions entail that the partonic soft function has an O(Λ QCD ) renormalon which affects the partonic threshold at = 0 for the perturbative ultrasoft radiation [21,89] in a way very similar to how the O(Λ QCD ) pole mass renormalon affects the thresholdŝ τ for the partonic ultra-collinear radiation. Both features are physically disentangled by the fact that the factorization theorem in Eq. (3) predicts that the ultrasoft effects on the hemisphere jet masses (and thus τ 2 ) are enhanced by a factor compared to the ultra-collinear effects. This also entails that for a top quark mass determination, different c.m. energies Q and simultaneous fits including parameters of the shape function must be considered to lift the degeneracy between m t and hadronization effects. This is in close analogy to the α s determinations from e + e − event-shape data carried out in Refs. [19,90,91], where nonperturbative effects were not fixed from hadronization corrections in MC event generators, but from a simultanous fit using data obtained for different c.m. energies Q.
The difference between the O(Λ QCD ) renormalons affecting the partonic soft and bHQET jet functions is that the former cancels inside Eq. (15) with the nonperturbative matrix element Ω 1 , while the latter is only an artifact of the pole mass and nonexistent when employing a short-distance mass scheme. One can remove the partonic soft function renormalon in an analagous way by re-expressing the first moment Ω 1 in a new scheme that includes a perturbative subtraction: is a perturbative series constructed such that it has exactly the same renormalon as Ω 1 and the MS renormalized partonic soft functionŜ (5) τ2 ( ,δ = 0, µ S ). Eq. (16) implies that This entails the introduction of the scale R s Λ QCD (in analogy to the scale R for the MSR mass dicussed below), which effectively represents an infrared cut for the partonic soft function which is then free from the O(Λ QCD ) renormalon ambiguity. This results in ∆(R s ), and hence Ω 1 (R s ), being renormalon free, and that instead of the MS renormalized partonic soft functionŜ (5) τ2 ( ,δ = 0, µ S ), one employs the 'gap subtracted' partonic soft function [21]: In analogy to the partonic 2-jettiness bHQET jet function it is strictly mandatory to expand the dependence on δ consistently in powers of the strong coupling α (5) s (µ S ) such that the soft function O(Λ QCD ) renormalon consistently cancels order by order. Furthermore, to avoid upsetting the soft function power counting and to avoid the appearance of large logarithms in the subtraction it is mandatory that R s is parametrically close to the typical soft scale, i.e. R s ∼ µ S ∼ŝ τ / . Settingδ = 0 one recovers the MS renormalized partonic soft function. Hence, the hadron-level soft function becomes We discuss the precise definition of the scheme that defines ∆(R s ) in Sec. IV.
At this point we note that there is a partial cancellation between the O(Λ QCD ) renormalons in the partonic soft and bHQET jet functions as the corresponding ambiguities and the associated diverging behavior of the perturbative series are equally severe but have an opposite sign [15]. The extent of the cancellation is Q-dependent through the factor that enters the factorization convolution between the partonic soft and bHQET jet functions in Eq. (3). As a consequence, the impact of the soft function renormalon in the 2-jettiness distribution increases with Q, while the impact of the jet function renormalon does not. This means that the overall effect of both O(Λ QCD ) renormalons may be hidden for a certain range of Q values when the pole mass is used in the jet function and no gap subtraction is carried out for the partonic soft function. As we show in Sec. VII this indeed happens within the range of Q values relevant top mass determinations. However, since simultaneous fits for different c.m. energies Q are mandatory to independently determine the top mass and the parameters of the shape function without degeneracy, the impact of the individual O(Λ QCD ) renormalons cannot be avoided, so that the perturbative uncertainties in analyses accounting for both types of renormalon subtractions are systematically smaller than those when renormalon subtractions are not implemented.

III. TOP MASS SCHEMES
In our analysis we employ two renormalization schemes for the top quark mass: the pole mass m pole t and the scale-dependent MSR mass m MSR,(5) t (R) [33][34][35]. The pole mass scheme has -from a technical point of viewthe simplest implementation because we can set the residual mass term δm appearing in Eqs. (3) and (13) to zero and all entries of m t discussed before to m pole t . However, employing the pole mass scheme entails an O(Λ QCD ) renormalon ambiguity, which visibly destabilizes the order-by-order behavior of the cross section in the peak region as we show explicitly in Sec. VII. This can be systematically avoided by using a suitable shortdistance mass scheme. The MSR mass m MSR,(5) t (R) is a short-distance scheme that is derived from the MS mass and represents an extension of the MS mass concept: While the MS mass m (6) t (µ) is suitable for scales µ ≥ m t , the MSR mass is appropriate for scales R ≤ m t .
For the MS mass we will use the notation m We give the details of the MSR scheme definition and the numerical values of the coefficients in App. D. This means that the MSR mass is a scheme derived from the MS mass, but where the virtual off-shell fluctuations in the on-shell self-energy from scales beyond R (which includes virtual top quark effects) are integrated out. The MSR mass is therefore designed for top mass dependent observables sensitive to soft QCD dynamics.
The scale R can be interpreted as the resolution scale below which virtual self-energy and real ultra-collinear radiation is treated as unresolved, so that only self-energy contributions above R are absorbed into the mass. This means that m from scales above R. This interpretation entails that in the limit R → 0, where all virtual self-energy and real ultra-collinear radiation is treated as resolved and all virtual self-energy contributions are absorbed in the mass, we approach the pole mass, which is precisely expressed in Eq. (21). The renormalon ambiguity of the pole mass can thus be seen to be associated with the problem that the limit R → 0 involves crossing the Landau pole of the strong coupling which a priori cannot be carried out in an unambiguous way.
Here we use the interpretation in [34], where the top MSR mass m  t . The matching relation is given in Eq. (D3). For scales R < m t the MSR mass evolves with the R-evolution equation where γ R n 's are obtained from the coefficients a i 's in Eq. (21) using the procedure outlined in App. F.
As explained in Sec. II B, the consistent use of a shortdistance mass in Eq. (3) entails that the MSR mass scale R satisfies the parametric relation R ∼ µ B ∼ŝ τ , which means that R depends on τ 2 and m MSR,(5) t adopts the status of a dynamical scale-dependent 'mass coupling' in complete analogy to the well-known concept of the scale-and flavor-number-dependent strong coupling α s . This dynamical treatment of the top quark MSR scheme resums important large logarithms via the R-evolution equation Eq. (22). The reader should note that the RHS of the R-evolution equation is linear in R, which differs from the common logarithmic renormalization group equations. This linear evolution is an essential aspect of properly treating the physical mass effects that govern the resonance/close-to-mass-shell dynamics of heavy colored particles.
In our numerical analysis we use m t ) as the standard reference mass value, which we quote as our main input and from which we then calculate the MSR or MS masses at the respective scales needed within the factorization formula in Eq. (3). For the flavor-number-dependent-strong coupling α (5,6) s (µ) we always use matching at 4-loops [92][93][94] and running at 4loops [95,96], where the flavor matching is carried out at µ = m (6) t . We finally note that the analytic properties of the pole mass O(Λ QCD ) renormalon in terms of knowledge on the large-order behavior of the perturbation series obtained in the pole mass scheme are by now very well understood, see e.g. Refs. [34,50] for the case where all quarks except for the top are assumed massless and Ref. [35] where finite bottom and charm masses are included systemat-ically. 8 One of the most interesting observations in this context is that the large-order asymptotic behavior (for some unknown reason) universally sets in at O(α s ) and is already well saturated at O(α 2 s ), which is one of the reasons why the pole mass renormalon has received significant attention in the literature. Another useful (but also confusing) consequence of this fact is that is it very easy to devise different types of low-scale short-distance masses from either physical [25,69,[97][98][99][100][101] or conceptual considerations [33][34][35]102]. All are -as long as the resolution scale R is assigned appropriate valuessimilarly effective in minimizing mass-related QCD corrections and stabilizing the perturbation series already at O(α 1,2 s ). The MSR mass m MSR,(5) t (R) provides a unifying concept to connect all low-scale short-distance masses with the MS mass via its renormalization group equation given in Eq. (22).

IV. SOFT GAP SUBTRACTION
The O(Λ QCD ) renormalon of the soft functionŜ has large-order properties very similar to the O(Λ QCD ) pole mass renormalon in the bHQET jet function. It is known to differ from the latter just due to a different normalization [21] and the sign difference already mentioned at the end of Sec. II. However, compared to the pole mass renormalon, the soft function O(Λ QCD ) renormalon is harder to pinpoint quantitatively at low orders due to the soft function anomalous dimension and hadron mass effects [22,103]. Therefore, the soft function renormalon does not seem to exhibit the same universality as the pole mass renormalon (even though it is phenomenologically equally relevant). Apart from e + e − event shapes [104], it is still largely unknown whether or in which way the O(Λ QCD ) renormalon inŜ (5) τ2 appears in a universal manner in soft functions relevant for collider observables which are linearly sensitive to large-angle soft radiation.
It has so far been the practice to define the renormalon subtraction seriesδ(R s ) in Eq. (17) directly from the soft function, see Refs. [21,84] for such approaches. In our analysis we employ a modified version of the prescription suggested in Ref. [84]. It is based on the Fourier transform of the renormalon-subtracted partonic soft functioñ where the gap subtraction seriesδ(R s ) is factored into the exponential factor shown in the second line. It is therefore possible to define an expression forδ(R s ) which cancels the O(Λ QCD ) renormalon of the soft function by the condition where s ij are coefficients of the fixed-order series expansion of log[S (5) τ2 (y, µ)] shown below explicitly in Eq. (28). When the gap subtraction seriesδ(R s ) is used in the factorization theorem, it is crucial that the renormalization scale of the strong coupling α s (R s ) is re-expressed in terms of α s (µ S ), the coupling used in the series for the soft function, to ensure a systematic order-by-order cancellation of the renormalon. This is detailed in App. E.
This definition of the gap subtraction can be contrasted with the one used in Ref. [84] where the subtraction series was instead related to a derivative of the soft function logarithm: In this definition, the scale of the strong coupling is µ S by construction, and the gap subtraction inherits a nontrivial anomalous dimension in µ S from the soft function. In App. E we describe a set of generic gap subtraction schemes that include Eqs. (24) and (25) as special cases.
While both definitions in Eqs. (24) and (25) are perfectly viable subtraction schemes (i.e. equally effective at asymptotic large orders), the seriesδ Ref. [84] in Eq. (25) is numerically zero at O(α s ) for µ S = R s because the oneloop non-cusp anomalous dimension vanishes, γ Sτ 0 = 0. This necessitates choosing R s strictly below µ S in the peak region to reduce the size of the O(α s ) correction. This can, however, be problematic when considering O(α 2 s ) corrections because in the peak region the soft scale µ S ∼ŝ τ / is already parametrically smaller than the top quark width such that setting R s > µ S can lead to instabilities. On the other hand, for the gap subtraction definition in Eq. (24), we have s 10 = 15.053 at one loop which, in addition to being non-zero, is also numerically sizable allowing for the implementation of an effective gap subtraction for the more natural setting R s = µ S . We will therefore adopt the gap subtraction scheme defined in Eq. (24) in our analysis.
The R-evolution of the gap parameter corresponding to the scheme defined byδ in Eq. (24) is given by where the anomalous dimension coefficientsγ R n are derived from the fixed-order coefficients s ij in Eq. (24) following the steps laid out in App. F. In our numerical analysis we take the first moment value for Ω (2 GeV) obtained in Ref. [91] as the input and determine from it Ω 1 (2 GeV) in the gap scheme defined in Eq. (24) as well as Ω 1 for the case of no gap subtraction. For the no-gap case we use ∆ = 0.1 GeV as the input value of the gap parameter in Eq. (15) which fixes the analytic form of the shape function F (k). From this we calculate ∆(R s = 2 GeV) in the scheme of Eq. (24) using Eq. (16). This input then unambiguously fixes the value of ∆(R s ) at any scale R s using the R-evolution equation Eq. (26), thus determining the form of the gapsubtracted shape function F (k − 2∆(R s )) in Eq. (20). For the R-evolution in Eq. (26) we employ 2-loop precision [70,71,84]. For F (k) we adopt the parametrization given in Ref. [87], see App. G for explicit expressions. We note that this approach implies that the form of the shape function F (k − 2∆(R s )) entering the factorization theorem Eq. (3) depends dynamically on the value of the physical 2-jettiness variable τ 2 . In the peak region this dependence is, however, quite weak because the scale µ S saturates, see Sec. VI.

V. N 3 LL RESUMMED CROSS SECTION
For the numerical evaluation of the factorization formula in Eq. (3), we find it convenient to work in the position (Fourier) space where the convolutions involving the bHQET jet and and soft functions and their respective renormalization evolution factors become simple products. (The results in this section could equally well be expressed using the Laplace transform.)

A. Stable-top cross section without renormalon subtractions
We first discuss the stable top quark cross section for the case of having neither soft function gap nor mass subtractions. The Fourier transforms of the stable top 2-jettiness bHQET jet and soft functions are defined via where for brevity we have dropped the zero arguments δm = 0,δ = 0 and Γ t = 0. Thus the position space jet and soft functions have the form: such that Eq. (3) can be written as × H (5,6) evol (Q, m t , , µ; µ H , µ m )e K (5) where due the convolution in Eq.
Here µ is the common final renormalization scale of all the RG evolution factors. Taking the inverse Fourier transform back to distribution space we find where for later convenience we have defined the following function of the derivative operator ∂ Ω : It acts on the function of Ω shown in Eq. (31), and the outcome is evaluated at the following µ-independent evolution kernel between the bHQET jet and soft scales: S (µ, µ S ) + ω We note that since the jet scale is always above the soft scale, one hasω (5) < 0. In Eq. (32) the arguments cusp non-cusp matching β[αs] γR δ LL 1  -tree  1  --NLL  2  1  tree  2  1  -N 2 LL  3  2  1  3  2  1  N 3 LL  4  3  2  4  3  2  NLL  2  1  1  2  1  1  N 2 LL  3  2  2  3  2  2   TABLE I. Loop corrections required for specified orders. In the last two columns γR and δ refer to either soft or MSRmass subtractions.
of the position-space jet and soft functions are understood to replace the corresponding logarithms shown in Eqs. (28). For sake of brevity we have suppressed the arguments µ B and µ S that appear in the running coupling, as shown in Eq. (28). The meaning of the superscript '(0)' onσ will be clarified below in Eq. (39). We also note that the dependence of the result in Eq. (32) onŝ τ is defined in terms of rational power plus-distributions which have support forŝ τ ≥ 0 in the case of stable top quarks [26]. To account for the fixed-order corrections contained in the product of the functions H Finally, since we carry out resummation at the level of the differential cross section, when implementing the cross section at N k LL+O(α k−1 s ) accuracy for k ≥ 1 (referred to as 'unprimed' orders), we explicitly incorporate the O(α k s ) plus-function boundary condition [105] in order to correctly sum up logarithms that are counted as N k LL in the exponent of the cumulative distribution. This amounts to including the single logarithmic terms appearing at O(α k s ) in the jet and soft functions. For the 'primed' orders N k LL , or equivalently, N k LL + O(α k s ) accuracy, this is not necessary as the O(α k s ) fixed-order matching already includes this single logarithmic term. The loop-order of the theoretical ingredients for the primed and unprimed orders are summarized in Tab. I. We refer to Ref. [91] for further details on primed and unprimed orders.

B. Renormalon subtractions
We now describe how the renormalon subtractions with respect to the top quark mass and soft function gap are to be included starting from the unsubtracted stable top cross section in Eq. (31). First we recall that the δm dependence in the bHQET jet function J The two terms shown in the second line of Eq. (34) represent those to be accounted for in Eq. (3) since powersuppressed contributions in the peak region must be systematically dropped for consistency. For the term δm(R), which contains the pole mass renormalon ambiguity, this is particularly important to achieve the order-by-order cancellation of the pole mass renormalon. Thus, the stable-top bHQET jet function using Eq. (13) in position space can be expressed up to NNLO as where we have dropped the zero argument Γ t = 0 in the jet function and the argument R in the mass subtraction δm(R) for simplicity. We have kept only terms at most quadratic in δm(R) so that the pole mass renormalon can be consistently canceled to O(α 2 s ). Likewise, the soft gap subtraction can be incorporated using Eq. (23) such that S (5) τ2 ( x,δ, µ S ) = e −2i δ xS(5) τ2 ( x, µ S ) where we dropped the argument R s in the gap subtractionδ(R s ). Including the subtraction terms in Eqs. (36) and (37) and strictly expanding to O(α 2 s ) we arrive at the following expression for the renormalon-subtracted cross section for stable top quarks: Here we have where the total subtraction series δ tot has the form and dσ (0) /dτ 2 is given in Eq. (32). It is essential to consistently drop terms of O(α 3 s ) and higher (concerning the fixed-order corrections in the hard, jet and soft functions as well as the renormalon subtractions) in the product in Eq. (39). Note that the derivatives in the sum over n in Eq. (38) are evaluated at Ω =ω + n. For simplicity we will display the functionsŝ MSR τ , δ m ,δ and δ tot without their arguments R and R s below.

C. Including the top width
From Eq. (11) we see that the cross section for unstable top quarks involves an additional Breit-Wigner convolution such that 9 dσ where as a result of the stable top cross section, the integration is bounded below byŝ = 0. However,ŝ MSR τ (orŝ τ in the pole mass scheme) for the unstable top 2jettiness cross section can be negative as well. To incorporate this convolution in Eq. (38) we first note that the cross section is proportional to 1/(ŝ )ω (5) +1 , as can be seen from Eq. (32) using Eq. (33). This factor can be brought to the right of the ∂ Ω derivatives using the identity where we remind the reader thatŝ is the integration variable in Eq. (41) whereasŝ τ (orŝ MSR τ in the MSR mass scheme) is related to the τ 2 measurement as defined in Eq. (5). The analogous relation also holds for the ∂ Ω derivatives associated to the soft function. Hence, we need to evaluate the following convolution: which has a smooth Ω → 0 limit, and where we have defined In the limit Γ t → 0 one smoothly recovers the stable-top results. For Γ t = 0 we find that I = 0 whenŝ τ < 0. Using these expressions in Eqs. (38) and (41) we now arrive at the final expression for the unstable-top cross section with renormalon subtractions: Finally, the hadron-level cross section is obtained from the partonic cross section via convolution with the nonperturbative shape function:

VI. PROFILE FUNCTIONS
To properly sum large logarithms we use τ 2 -dependent renormalization scales µ i (τ 2 ), R(τ 2 ) and R s (τ 2 ), called profile functions [87,91]. They have canonical scaling in the resummation regions, and freeze at a perturbative scale in the resonance region to avoid the breakdown of perturbation theory for anomalous dimensions. In the far tail region, they become equal to the hard scale to reproduce the fixed-order perturbative expansion with a common scale µ. They are expressed in terms of 7 parameters which can be varied to estimate perturbative uncertainties. Following Refs. [27,106] we employ a natural generalization of the profile functions devised for massless event shapes in [19], to which they reduce in the massless limit.
The strategy to estimate perturbative uncertainties involves varying all the profile functions up and down by at most a factor of 2 and 1/2, respectively, including a shape dependent variation in the jet scale, as well as varying the value at which the soft scale freezes in the nonperturbative region. We show bands for the latter two variations in Fig. 1, and indicate the factor of 2 variations by arrows. We scan over these profile variations by generating a sample of 500 profiles were all their parameters are varied simultaneously with random choices within their allowed ranges (see Ref. [91] for details on this general approach). The concrete form of the profile functions and how their parameters are varied are given below. The total uncertainty is determined by the envelope of the resulting cross sections.
For the hard and mass matching scales we use Our τ 2 -dependent profile functions are implemented through the soft scale µ S (τ 2 ), which in the peak region is parametrized with the following piecewise function where τ min 2 = 2m 2 t withm t ≡ m t /Q. We refer to the three corresponding intervals as the "non-perturbative", "bHQET", and "SCET-resummation" regions. The function ζ[f a (τ 2 ), f b (τ 2 ), t a , t b , τ 2 ] smoothly connects any two linear functions f a,b (τ 2 ) that end/begin at the points t a,b . This is achieved by means of two quadratic polynomials of τ 2 smoothly joined at τ 2 = (t a + t b )/2, where the explicit formula can be found in Ref. [19]. The parameter n s has the default value 0.5 and is varied in the range |n s − 0.5| ≤ 0.025. Its effect in the peak region is relatively mild. The default slope in the SCET-resummation region is set by r s = 2 and using the default value 0 for the variable e slope . Slope variations are implemented by varying e slope in the interval [1/1.13 − 1, 1.13 − 1]. Note that the rescaling factor to the left of r s approaches 1 in the massless limit. In a similar way, the parameter affecting the flat non-perturbative region is e np . Its default value is 1 and it is varied in the interval [1/ √ 2, √ 2] with µ 0 = 3 GeV. The values of the transition points t 0,1 (m t , Q) depend on the mass and center-of-mass en-ergy and have the form with |d 0,1 | ≤ 0.05 and |n 2 − 0.25| ≤ 0.025. Their default values are d 0,1 = 0.05 and n 2 = 0.25. For the energies and masses considered in this article one has t 0 < t 1 < t 2 . For the jet scale profile function µ J (τ 2 ) we first definẽ µ J (τ 2 ) = √ e H µ S (τ 2 )/m t and t s = n s +m t , and then use the piecewise function Here the jet-function parameterẽ J is defined with a rescaling factor with variations |e J | ≤ 1.5 used for assessing uncertainties and the default value e J = 0. For m t = 0 we recover e J =ẽ J . We set the mass and soft-function renormalon subtraction scales to their respective canonical values: R s (τ 2 ) = µ S (τ 2 ) and R(τ 2 ) = µ J (τ 2 ).

VII. NUMERICAL ANALYSIS
In this section we present a numerical analysis of the bHQET N 3 LL + O(α 2 s ) 2-jettiness peak region cross section based on the factorization formula of Eq. (3). We remind the reader that this factorization theorem is based on the bHQET limit and does not account for subleading terms related to higher powers ofŝ τ (kinematic power corrections) and m t /Q (mass power corrections). As was shown in Refs. [27,82,106], these corrections are formally power-suppressed in the peak region, but still important numerically for a realistic phenomenological analysis concerning the top mass dependence of the 2-jettiness peakregion line shape. In the following we therefore carry out a generic numerical analysis pointing out important features of the bHQET 2-jettiness peak region cross section at this order related to the convergence of the perturbative series, as well as the impact of the jet and soft function renormalons and the improvement related to their subtractions. A phenomenological analysis aiming for a systematic study of other sources of theoretical uncertainties is postponed to future work.
For our analysis we have two independent codes to implement all cross sections, one in fortran [107] and one in C++ using SCETlib [108]. Numerical integrations are carried out using either quadpak [109] or the gsl library [110].
Before getting to results some comments on our parametric input are in order. For the subsequent discussion we use α (5) s (m Z ) = 0.118 with m Z = 91.1876 GeV as input for the strong coupling. As input for the top quark mass we take the standard MS mass m (6) t ≡ m for the MSR top mass at 2 GeV, which one can interpret as the (renormalon-free) kinematic mass that governs the top-mass dependence of the peak position. The same top mass value is used in the boost parameter defined in Eq. (4) when the MSR mass scheme is used. For comparison we will also discuss results for the cross section in the pole mass scheme and without gap subtraction. At this point one has to recall that the pole mass is, due to its renormalon ambiguity, an order-dependent concept, where the size of the fixed-order corrections in its relation to a short-distance mass at a given order depends on the renormalization scale of the short-distance mass. Thus to achieve comparable theoretical predictions employing the MSR and the pole mass schemes (i.e. with peak positions that are compatible), it is essential to apply fixed-order conversion from the MSR to the pole mass at the renormalization scale of the MSR mass that is employed in the peak region of the distribution (which is the part of the distribution that carries the highest top quark mass sensitivity). Furthermore the order of conversion has to match the fixed-order input used for the theoretical calculation [13,35]. In the peak region of the distribution the appropriate renormalization scale is just the bHQET jet function scale µ J , around 10 GeV, see Sec. VI. Thus, the proper way to determine the pole mass is to use the fixed-order conversion from m MSR t (10 GeV) as was pointed out in Ref. [27]. Therefore, to determine the pole mass for the following analysis we first determine m MSR t (10 GeV) using Eq. (51) and R-evolution and then apply fixed-order conversion to the pole scheme at the scale R = 10 GeV. As we are employing O(α 2 s ) fixedorder corrections to the bHQET jet function at the highest N 3 LL order of our analysis, this fixed order conversion must be carried out at 2-loops. Using this procedure, we find This input value for the top quark pole mass ensures that at N 3 LL the peak position in the pole mass scheme is compatible with that obtained in the MSR mass scheme.
Lastly, in order to fix the form of the nonperturbative model function F (k), see Eq. (14) and App. G, we have to specify values for the first moment Ω 1 . In analogy to the top quark mass there are also different schemes for Ω 1 related to the definition of the gap subtraction parameter. Thus the corresponding values for Ω 1 also have to be fixed with some care. Here we aim to adopt values for Ω 1 consistent with the e + e − thrust analysis of Ref. [91], using the fact that the same universal soft function that is given in Eq. (14) also appears in the thrust factorization theorem for massless quark production. The value of Ω 1 determined in Ref. [91], Ω Ref.
[84] 1 (R s = 2 GeV) = 0.323 ± 0.045 GeV, was based on the gap subtraction scheme suggested in Ref. [84] [explained in detail in and below Eq. (25)]. For our analysis it needs to be converted to the gap scheme of Eq. (24) adopted in this work, referred to as Ω 1 (R s ), as well as to the unsubtracted MS gap scheme Ω 1 (which still contains the soft function renormalon). For these conversions we must use the 2-loop fixed-order formulas in Eqs. (E5) and (G4), respectively, since at the highest N 3 LL order of our analysis we employ the soft function at O(α 2 s ). This gives where we convert at the scale R s = 2 GeV, the typical value of the soft function profile in the peak region. Fixing ∆ = 0.1 GeV for the unsubtracted gap [ see Eq. (15) ] and using Eq. (16) we obtain ∆(R s = 2 GeV) = 0.563 GeV for the gap term entering Eq. (20). With this choice of parameters we can use the same shape function parameters as employed in Ref. [91], which corresponds to taking λ = 0.349 GeV and c 2 = 0.05 in Eq. (G1). We then use Eq. (26) with O(α 2 s ) running to determine ∆(R s ) at the τ 2 -dependent R s values required by the gap subtraction profile function.
In Fig. 2 the 2-jettiness differential cross section in the MSR mass scheme with gap-subtraction is shown for Q = 700 GeV as a function of the inclusive jet mass variable M J [ see Eq. (6) ] at all primed and unprimed orders up to N 3 LL for the default set of profile functions (see Sec. VI). Here the primed orders include all contributions of the unprimed with the addition of the fixed-order matrix elements at one higher order in α s . Our results at NNLL and N 3 LL are new and have not been analyzed before in the literature.
In Fig. 2 all curves are normalized to the Born-level massless total cross section σ 0 . The behavior of the curves at the different orders therefore reflects the effects of the perturbative corrections to the shape as well as to the normalization of the cross section. We see that, apart from sizable normalization corrections which happen to be positive for all subsequent orders, the convergence with respect to the peak location and the shape is excellent. From the difference in the curves obtained from primed and unprimed orders we can also see that Perturbative convergence of the 2-jettiness cross section, with all curves normalized with the Born cross section σ0. The peak location shows excellent convergence, while the normalization corrections are significant and exhibit slower convergence.
the effects on the normalization from higher order corrections in the renormalization group equations (and thus the resummation of large logarithmic terms) are smaller than those in the fixed-order matrix elements.
In order to analyse the effects of the higher-order corrections to the distribution shape and its orderdependent perturbative uncertainty, it is useful to normalize the curves from the different orders to a common M J interval. In Fig. 3 the 2-jettiness differential cross sections at Q = 700 GeV (upper panels) and Q = 2000 GeV (lower panels) are shown in the MSR mass scheme with gap subtractions using default profile functions. The results are normalized to the M J interval displayed in the respective panels at NLL (green dotted line), NNLL (blue dash-dotted line) and N 3 LL (red solid line). Primed orders are not displayed to avoid cluttering. We also display uncertainty bands with the corresponding colors at each of these three orders. These bands are derived by determining the upper and lower value of the distributions (for each M J value) obtained by considering 500 profile functions generated randomly within the profile function parameter ranges given in Sec. VI. To generate the bands, each cross section from a given profile is normalized to the displayed M J range. The central curves exhibit excellent perturbative convergence for the shape. The width of each band illustrates the size of the perturbative uncertainty, which nicely decreases with increasing order. For better visibility the error bands and lines are displayed once more in the lower parts of each plot showing the fractional deviation from the central N 3 LL curve. At Q = 700 GeV for M J ≥ 171 GeV the relative uncertainty in the peak region is ±(4-10)% at NNLL and ±(3-7)% at N 3 LL. In contrast, at Q = 2000 GeV for M J ≥ 175 GeV the relative uncertainty in the peak region is ±(3-8)% at NNLL and ±(1-5)% at N 3 LL.
In Fig. 4 we show for comparison the analogous results for cross sections in which the pole scheme for the top quark mass is employed without gap subtractions. We will return and discuss this figure in more detail below. It is instructive to first examine the importance and interplay of the O(Λ QCD ) renormalons contained in the perturbative fixed-order series of the bHQET jet and par- s ) (solid red) and consistently using gap subtractions. Thus, renormalon subtractions associated with the pole mass renormalon are not included, while the soft function renormalon is still removed systematically. Since the renormalization group evolution, which predominantly effects the normalization, does not contain any renormalon effects, we adopt the highest order N 3 LL anomalous dimensions for all renormalization-group resummation factors, so that the behavior of the three curves is focused on the pole mass renormalon in the unsubtracted jet function. The curves clearly exhibit the well-known pole mass renormalon problem which causes the peak position to systematically shift towards smaller jet masses with increasing order. At the level of the twoloop bHQET jet function itself, this behavior was discussed in [69]. In a fit to data this behavior would correspond to a pole mass value that systematically increases with the perturbative order. This is the known behavior of the perturbative series for the pole mass in terms of a short distance mass [13,35,50]. Furthermore, the curves show some instabilities in its shape, in particular in the form of the distribution at and above the peak.
In order to illustrate the impact of the soft function renormalon we display in the middle panel of Fig. 5 the 2-jettiness cross section without gap subtractions, again for Q = 700 GeV using the default profile functions and consistently expanding all fixed-order matrix elements entering the factorization theorem, but this time using the MSR top quark mass scheme (and the standard MS mass in the hard function). Here subtractions associated to the soft function renormalon are not included, while the pole mass renormalon is removed systematically. As in the upper panel, we adopt the highest order N 3 LL anomalous dimensions for all renormalization-group resummation factors, so that the three curves focus on the behavior due to the soft function renormalon. We see that the soft function renormalon causes the peak position to systematically shift towards larger jet masses with increasing order. In a fit to data this behavior would correspond to an MSR mass value that systematically decreases with the perturbative order. Furthermore, the curves at O(α s ) and O(α 2 s ) show considerable shape instabilities in the region below the peak, where the cross section can even become negative. We note that the impact of the soft function renormalon increases with the c.m. energy Q. This dependence on Q arises from the boost factor = Q/m t appearing in the convolution integral shown in Eq. (3), which is also manifest in Eq. (40).
Finally, in the lower panel of Fig. 5 we have displayed the corresponding three curves once again, but systematically accounting for the subtractions associated to both the pole mass and soft function renormalons, by using the MSR top quark mass scheme and gap subtractions, respectively. We now observe very good convergence of the peak position and, furthermore, no instabilities in the shape of the distribution are visible.
The upper and middle panels of Fig. 5 also nicely illus- trate the presence of a partial cancellation of the jet and soft function renormalon effects since they have opposite signs. In the combined order-by-order cross sections both corrections thus partially cancel when the pole mass scheme is employed and no gap subtraction is applied for the soft function. Even though this partical cancella-tion arises between two independent physical effects and should therefore be considered as accidental from a principle point of view, it does undeniably take place in the physical regions of c.m. energies where high-precision extractions of the top mass can be carried out. One may therefore ask the question whether this cancellation may in principle allow for a pole mass determination where the impact of the pole mass renormalon could be tamed or even avoided altogether. At this point we would like to remind the reader that for a top mass determination from data (or MC pseudo data) simultaneous fits of the peak region 2-jettiness distribution for several Q values are needed to disentangle the dependence on the top quark mass and the shape function parameters. So there is a strong degeneracy concerning the dependence on the top quark mass and the shape function parameters, and in particular its first moment Ω 1 . Given that there are strong cancellations between corrections affecting both of these dependences, it can be expected that they degrade the overall precision of such an analysis. Furthermore, since the amount of mutual cancellation between the pole mass and soft function renormalons is Q-dependent, it is expected that such fits for theoretical predictions without any renormalon subtractions will exhibit larger theoretical uncertainties compared to those where the pole mass and soft function renormalons are separately and independently subtracted. Such an extensive analysis is, however, beyond the sope of this work.
Furthermore, it should also be pointed out that the shape function appearing in Eq. (3) is universal and appears in the same form also in the factorization theorem for the e + e − thrust distribution below the top pair threshold, where precise information on its parameters can be extracted from available e + e − data [91] at significantly smaller values of Q. The thrust distribution for massless quark production is sensitive to the same soft function renormalon, but does not have any top mass dependence. Thus if information on the renormalon-free shape function parameters obtained from e + e − thrust data are systematically accounted for, it is unavoidable that renormalon effects must be properly handled for top quark mass determinations from the 2-jettiness distribution.
Finally, it is instructive to also have a closer look at the 2-jettiness cross section without any renormalon subtraction. In Fig. 4    ther the pole mass nor the soft function renormalons are subtracted, the setup used for all curves and uncertainty bands in Fig. 4 is precisely the same as the one used for Fig. 3. We see that the perturbative behavior concerning the convergence and the perturbative uncertainties is also good even without any renormalon subtraction. This underlines the partial cancellation of the jet and soft function renormalons. However, a closer inspection shows that the perturbative uncertainty bands are narrower when the subtraction of all renormalons is taken care of systematically. This is visible in the fractional deviation plots, where the N 3 LL renormalon-subtracted predictions in Fig. 3 exhibit an average uncertainty of ±3.8% at Q = 700 GeV for M J ≥ 171 GeV compared to ±5.5% for the predictions without any renormalon subtraction in Fig. 4. In contrast, for Q = 2000 GeV Fig. 3 has an average uncertainty of ±2.4% for M J ≥ 175 GeV compared to ±2.9% for the predictions without any renormalon subtraction in Fig. 4. An interesting aspect of our definition for the jet mass variable M J [ defined using 2-jettiness in Eq. (6) ], is that it is normalized in a way such that it can be seen as a direct measure for the top quark mass. Therefore, the behavior of the peak position for the M J distribution allows us to draw conclusions on the size of the perturbative uncertainties of a top mass determination from the peak position. To avoid outliers we discard the two highest and two lowest points in the scan so as to better represent the bulk of the points. In Fig. 6 we show the peak positions of the curves at Q = 700 GeV and Q = 2000 GeV for the default profile functions and their perturbative uncertainties, estimated from the 500 random profile functions. Results are shown in the MSR mass scheme with gap subtraction (MSR, red) and in the pole mass scheme without gap subractions (Pole, blue) at NLL, NNLL and N 3 LL. The central values, shown by dots, correspond to the default profile scales, so the perturbative uncertainties are asymmetric. The corresponding numbers are also given in Tab. II. The results show that the perturbative uncertainty is systematically smaller when the top quark mass and soft function renormalons are subtracted. For Q = 700 GeV, where we have the highest top quark mass sensitivity, using renormalon subtractions leads to an uncertainty in the peak location of around ±85 MeV at N 3 LL order. Without renormalon subtractions the uncertainty at this order increases to around ±150 MeV, which is almost factor of two larger. For Q = 2000 GeV the uncertainties are larger for the analysis without renormalon subtraction as well (around ±450 MeV with renormalon subtraction compared to around ±650 MeV without renormalon subtractions, both at N 3 LL order). Here the difference is less pronounced because the overall top mass quark sensitivity decreases for larger Q values and the overall uncertainties increase. The analogous behavior is also visible at lower orders. Our results indicate that the MSR mass may be extracted with an uncertainty of well below 100 MeV, while the pole mass uncertainty is at the level of 150 MeV. Interestingly, this is about the size of the top quark pole mass renormalon ambiguity of 166 MeV that was estimated recently in Ref. [35] for the case of massless charm and bottom quarks (which is the approximation we use in our analysis). Note that in an earlier analysis in Ref. [50] the top quark pole mass renormalon ambiguity was estimated as the smaller value of 67 MeV (also for massless charm and bottom quarks).
The results we have obtained in this simple analysis of the peak positions do -taken by themselves -not contradict the view that the top quark pole mass can be extracted from the 2-jettiness cross section with perturbative uncertainties below the pole mass renormalon ambiguity, but they also show that at least at N 3 LL or-der the precision is not (yet) sufficient to achieve that goal and that higher-order corrections beyond this order would be mandatory to get there. On the other hand, the results also support the view that, even though the 2-jettiness cross section exhibits a cancellation between the pole mass and soft function renormalons, the pole mass can still not be extracted with a precision below its renormalon ambiguity. In any case, using renormalon subtractions, and in particular the MSR mass scheme, will yield substantially higher precision and smaller perturbative uncertainties.
At this point we would like to again mention, that in mass determinations from data (or MC pseudo data) simultaneous fits of the peak region 2-jettiness distribution for several Q values are needed to disentangle the dependence on the top quark mass and the shape function parameters, and that the whole distribution in the peak region (rather than just the peak position) would enter such fits. As mentioned before, however, this kind of study requires that also off-shell and m t /Q powersuppressed contributions are included, as their effects can be non-negligible depending on how the cross section is normalized.
The dominant such QCD corrections to the factorization theorem in the bHQET region come from two sources, mass power corrections appearing as higher order terms in Eq. (2), and corrections to the perturbative singular structures. (Additional non-singular kinematic power corrections are very small at one-loop and hence irrelevant.) The former are universal at any order in α s and shift the distribution to the right by O(m 4 t /Q 4 ) but are trivial to incorporate. The latter are known analytically to O(α s ) is QCD [28]: at tree-level one gets a modification of the coefficient of the delta function, while at O(α s ) also the plus distribution coefficient is affected. One can include these mass corrections by a suitable modification of the hard and jet functions (see e.g. Ref. [29] for more details). These mass power corrections decrease the cross sections in Fig. 2 by 5% (beyond NLL) everywhere except for the region to the left of the peak where the effects are smaller. However if the cross section is normalized the effect of these power corrections drops below a percent becoming negligible in all relevant regions. It is reasonable to believe that these power corrections will be of similar form and size once O(α 2 s ) are added. A complete analysis that accounts for these effects will require the inclusion of the O(α 2 s ) correction to the primary massive quark SCET jet function that was computed recently in Ref. [74] and shall be carried out in future work.

VIII. CONCLUSIONS
In this article we have presented results for the 2-jettiness differential distribution for boosted tops produced in e + e − collisions in the peak region, accounting for the resummation of large QCD logarithms at next-to-next-to-next-to-leading logarithmic (N 3 LL) order and fixed-order corrections to the hard, soft and jet function matrix elements at next-to-next-to-leading order [ O(α 2 s ) ], calculated in the framework of soft-collinear effective theory and boosted heavy quark effective theory. We have systematically removed the O(Λ QCD ) renormalons contained in the soft and jet functions, by using a gap subtraction as well as the MSR mass, and have provided a numerical analysis indicating that the perturbative uncertainties of a determination of the top quark MSR mass from the N 3 LL + O(α 2 s ) prediction at a c.m. energy of Q = 700 GeV are well below the level of 100 MeV. For future reference all theoretical formulae have been given explicitly in several appendices.
An interesting aspect of the 2-jettiness distribution is that the soft-and jet-function renormalons partially cancel each other for center-of-mass energies above 700 GeV where the boosted top quark approximation is valid and precise top mass determinations can be carried out. This cancellation arises because the soft-function and polemass renormalons enter with different signs. While these two renormalons represent two physically independent infrared sensitivities, the cancellation allows for rather stable and convergent predictions in the pole mass scheme, if at the same time also the soft function renormalon is left unsubtracted. However, the resulting perturbative uncertainties are still systematically larger compared to the predictions where both renormalons are independently removed.
The analysis done here based on boosted heavy quark effective theory neglects subleading collinear off-shell corrections, which have been determined recently at O(α 2 s ) in Ref. [74] and shall be accounted for in future work. The factorization theorem presented in Eq. (3) is expressed in momentum space where the jet and soft functions, along with the evolution factors, are distributions involving series in plus and delta functions. We find it convenient to combine the ingredients in position space where the convolutions become simple products. Below our notation and definitions with variable mass dimension follow Ref. [26]. For earlier work on the associated resummation formulae see Refs. [112][113][114]. For a function F(q, µ) that depends on a momentum-space variable q with mass dimensions j F and the renormalization scale µ, the Fourier transform in position space is defined as where x has mass dimensions −j F . The position-space anomalous dimension permits writing the corresponding RG relation as a local equality: such that RG-evolved position-space soft and jet functions can be expressed as a regular product: F(q, µ) = dq U F (q − q , µ, µ 0 ) F(q , µ 0 ) (A3) = dx 2π e iqxŨ F (x, µ, µ 0 )F(x, µ 0 ) .

Generating fixed order terms
We now describe a helpful algorithm that allows one to generate the fixed-order expansion of the positionspace matrix elements and factorization functions discussed here from the non-logarithmic coefficients and their anomalous dimensions in Eq. (B6). This applies to the matching coefficient H (6) Q , the bHQET jet and the soft matrix elements, and can be used to reproduce the results stated above in Eqs. (C1), (C4) and (C5). Note that the algorithm described below must be generalized in an obvious way for H (6) m to obtain the results in Eqs. (C2) and (C3), since its running results from the anomalous dimensions of the bHQET current and the SCET matching coefficient (with different number of dynamical flavors in the running coupling) and includes an additional rapidity logarithm.
Consider a position-space matrix element or a factorization function F (µ, Q) having the generic form where we have A = 1/m 2 t for the thrust position-space bHQET jet function and A = 1 in the other cases, see Eq. (28). The constant terms a F m0 serve as independent data, whereas other coefficients can then be determined by anomalous dimensions and the beta function. Thus the a F m0 serve as boundary condition data for the RG differential equations in Eq. (A4). The logarithmic terms a F mn for n ≥ 1 then can be expressed as where the three terms result from the running of the coupling, the non-cusp and cusp pieces of the anomalous dimension of the given function. The coefficients in Eq. (C7) can be obtained via the following recursion relations: with m > 1 in order to have a sensible upper limit. The starting values of the three series (with m ≥ 1) are given by Here, the integer j F corresponds to the dimension of the momentum-space variable q as it appears naturally in logarithms, 10 as shown in Eq. (A4). The dependence on j F is factorized as in Eq. (C7) and only enters the a F mn [β, Γ F ] coefficients through the boundary condition in Eq. (C9). For the factorization functions in Eq. (A4) this is simply set to 1. The constant terms of the SCET Therefore the final expression used for our analysis is with d ij = 2 j i−1 k=j k d k(j−1) β (5) i−k−1 .
Finally, we can relate the leading power correction Ω 1 (R s ) between the two subtraction schemes: