A numerical formulation of resummation in effective field theory

In this article we show how the resummation of infrared and collinear logarithms within Soft-Collinear Effective Theory (SCET) can be formulated in a way that makes it suitable for a Monte-Carlo implementation. This is done by applying the techniques developed for automated resummation using the branching formalism, which have resulted in the general resummation approach CAESAR/ARES. This work builds a connection between the two resummation approaches, and paves the way to automated resummation in SCET. As a case study we consider the resummation of the thrust distribution in electron-positron collisions at next-to-leading logarithm (NLL). However, the results presented here are easily generalizable to more complicated observables and processes as well as to higher orders in the logarithmic accuracy.


Introduction
A well known fact of perturbation theory is the presence of logarithmic terms, sensitive to the ratio of scales in the problem, whose power grows with the perturbative order. For most processes of interest at high energy colliders, two powers of such logarithms (L) arise for every power of the strong coupling constant, and the numerical size of these logarithms can be of order L ∼ 1/α s . This makes fixed order (FO) perturbation theory for such processes ill behaved, requiring a rearrangement of the perturbative expansion, in which these large logarithms are resummed to all orders in perturbation theory. Instead of simply counting the powers of the strong coupling constant, where N n LO refers to a calculation satisfying σ exact /σ N n LO − 1 ∼ α n+1 s , one instead performs a logarithmic resummation, in which N n LL implies that ln[σ exact /σ N n LL ] ∼ α m s L m−n for any m ≥ n. As long as the ratio of scales in the logarithm is large (L 1), this reorganization of perturbation theory provides a sensible expansion.
In this paper we will study the 2-jet cross section at lepton colliders (equivalent considerations apply to the 0-jet cross sections at hadron colliders), that we denote as where σ B denotes the Born cross section and V (Φ n ) an observable that goes to zero in the Born kinematics. The radiation phase space Φ n has the property that in the limit v → 0 each strongly interacting particle is either infinitely soft or it is collinear to either of the two Born legs (either the directions of the 2-jets in lepton colliders or the beam directions for hadron colliders). The resummed expression commonly takes the form Resummation to order N k LL requires knowledge of the functions g m (α s L) with m ≤ k + 1 and the coefficients C (m) with m ≤ k − 1. So for example, for NLL resummation one needs g 1 (α s L) and g 2 (α s L), but only the leading order coefficient C (0) = 1.
There are two main approaches to calculate the resummed expressions. The first is based on deriving a factorization theorem for the specific cross section under consideration, and then using evolution equations to resum the logarithms of the various ingredients of the factorization theorem. This approach was started by Collins, Soper and Sterman [1], where the transverse momentum distribution of the vector boson in W and Z production was studied. The development of soft-collinear effective theory (SCET) [2][3][4][5] has formulated this approach in the framework of effective field theories. In SCET, the interactions between particles are already factorized at the level of the Lagrangian [4]. As long as the measurement function can be shown to factorize [6,7], the relevant interactions can be separated as arising from collinear and soft particles, which do not interact with one another directly. Once a factorization theorem for a given observable has been derived, the important property of SCET is that in dimensional regularization each element of the factorization theorem only depends on a single scale through the ratio to the renormalization scale [2,3]. This means that all logarithmic terms can be obtained from the renormalization (RG) group evolution of each element in the factorization theorem, resulting in analytic results for resummed cross sections.
An alternate approach to resummation is based on the branching formalism [8,9], which is built on the factorization properties of the QCD squared amplitudes. It relies on slicing the radiation phase space by means of a resolution scale q 0 and then describing the radiation of particles above and below that scale separately. At lowest order this approach is the same as what is used in a parton shower, but by not requiring to produce fully exclusive final states which satisfy momentum conservation, one can systematically improve the branching formalism to any logarithmic accuracy desired. The logarithmic dependence on the resolution scale q 0 cancels for infrared and collinear (IRC) safe observables, leaving only a power suppressed dependence on q 0 , which allows one to take the limit for q 0 → 0. The key difference between the branching formalism and the resummation from factorization is that in the factorization approach the resummation is obtained by solving a set of differential equations (RG equations in the case of SCET), while in the branching formalism resummation is usually performed through an all order calculation in perturbation theory using a Monte-Carlo (MC) algorithm. For sufficiently simple observables one can rewrite the branching formalism in terms of differential equations, reproducing the results from the factorization approach.
Each of these two approaches has its advantages and disadvantages. The factorization approach has the advantage that to perform resummation only requires to obtain the desired differential RG equations that, in particular using SCET, can be obtained by computing well defined loop diagrams. On the other hand, the formalism only works for observables for which a factorization theorem is known. While these factorization theorems can easily be obtained for the simplest observables, for more complicated observables the factorization formulae are not known. In fact, the most difficult part to obtain resummation in SCET is often the derivation of the appropriate factorization theorem. Since the branching formalism does not rely on a specific factorization theorem, it can be applied to a rather wide class of observables. In fact, one can show that for any continuously global, recursively infrared and collinear (rIRC) safe observable [10], resummed results can be obtained using an appropriate MC algorithm. The downside of the branching formalism is that going to higher logarithmic accuracy requires a careful definition of the necessary squared amplitudes and phase space constraints, making extensions to arbitrarily high accuracy less systematic. Besides the formulation of resummation, the two approaches differ in a number of important aspects. For instance, while in the branching formalism the uncertainties associated with higher-order terms are estimated by varying the renormalization (and factorization in the hadronic case) and resummation scales, in SCET each of the ingredients of the factorization theorem has a characteristic scale that can be separately varied (commonly done through profile functions [11]) to estimate the perturbative uncertainty. Another important difference is in the way the two methods are sensitive to non-perturbative effects. While in the branching formalism this sensitivity comes from the scale at which the running coupling is evaluated, in a factorized approach one has an operator definition of the factorization ingredients in the non-perturbative regime.
As was shown in refs [10,12,13], one can reduce the requirements of the accuracy of the ingredients in the MC algorithm by one order if one has analytic knowledge of the resummation of any other observable that has the same double logarithmic structure as the observable one is interested in resumming. In other words, if one wants to resum Σ(v) for a given observable v to order N n LL, and one knows the resummmation analytically for Σ(v s ), where v s is a different (simpler) observable that has the same double logarithmic structure as Σ(v), then the ratio Σ(v)/Σ(v s ) can be computed using a MC that only requires ingredients at N n−1 LL accuracy. This was used in [10] to construct a completely generic method capable of computing any IRC safe observable to NLL accuracy. The simplified observable was chosen such that the branching formalism could be solved analytically. This approach was later reformulated to NNLL for generic rIRC safe observables in [13,14] and to even higher orders for specific observables [15].
By combining the two approaches with one another one keeps the advantages of each, while removing the main obstacle. If one can find for any observable v a simplified observable v s with a simple factorization formula, one can use SCET to obtain the analytical resummation of this simplified observable, while using the branching formalism to relate this analytical result to the resummed result for the more complicated observable, for which a factorization formula might be difficult or impossible to obtain. This is in spirit very similar to the CAESAR/ARES approach, but instead of finding a simplified observable for which the branching formalism can be solved analytically, one chooses the simplified observable such that a factorization theorem is easily derived and its resummation performed in SCET. In fact, there is a large body of observables for which high logarithmic accuracy is known in SCET (see references above). By combining this with the branching formalism allows one to obtain results for all related observables to the same level of accuracy.
It is this combination of SCET results with the branching formalism that we will address in this paper. A major part of this discussion will explain how to deal with UV divergent phase space regions that are crucial in the SCET approach (since it is the UV divergences from these regions that give rise to the anomalous dimensions leading to the RG equations in SCET), but can not be present in a MC approach which has to integrate over physical regions of phase space. We will explain this combination using the thrust distribution (using τ = 1 − T ) as an explicit example. Although the ingredients for a N 3 LL resummation are currently known [11,16,[40][41][42] (with the sole exception of the fourloop cusp anomalous dimension), in this paper we limit ourselves to NLL for the sake of simplicity. However our results are easily generalizable to more complicated observables of interest as well as to the computation of higher-order corrections. We leave this to a forthcoming publication [43]. This paper is organized as follows: We review the resummation using the branching formalism in Section 2 and the SCET approach to resummation in Section 3. The main part of the paper is contained in Section 4, where we discuss how to combine the two approaches to obtain a numerical approach to resummation in SCET. Conclusions and an outlook to future work is presented in Section 5.

Review of QCD resummation in the CAESAR/ARES framework
In this section we briefly review how resummation is carried out in the approach of refs. [10,13]. We begin by discussing the basic setup of the formalism for a general observable v and to arbitrary order in the resummation, but then restrict ourselves to the specific case of the thrust distribution when deriving the NLL result in more detail. We consider observables that vanish in the 2-jet limit, and when considering the thrust distribution we use τ ≡ 1−T .
At Born level, the final state consists of two back-to-back particles along the thrust axis with center of mass energy Q. Beyond Born level, further radiation is present and the final state consists in general of n secondary emissions, k 1 , . . . , k n , and of the primary quark and antiquark which recoil against these additional emissions. We denote the value of the observable by V (Φ B ; k 1 , . . . , k n ), where the Born phase space Φ B contains the dependence on the two Born momenta.
In order to single out the dependence on the Born phase space, we write and we will work with the expression Σ(Φ B ; v) for most of this paper. This means that the Φ B phase space integral needs to be performed at the end, and the final result needs to be divided by the Born cross section σ B . In the infrared and collinear limit, Σ(Φ B ; v) receives contributions from either virtual, or soft and/or collinear real corrections. In the following, we denote by V(Φ B ) the quark form factor at all orders (see e.g. [44]). Therefore we can write where M is the matrix element for n real emissions (the case with n = 0 reduces to the Born matrix element), and [dk i ] denotes the phase space for the emission k i . Each matrix element in Eq. (2.2) receives higher-order virtual corrections, while the number of real emissions is fixed by the index of the sum. The Θ function represents the measurement function for the observable under consideration.

The simplified observable and the transfer function
The general strategy of the CAESAR/ARES approach is to write the cross section for a rIRC safe observable v into the cross section of a simpler observable v s which has the same logarithmic structure as v at lowest order 1 , and a transfer function that accounts for the difference between the two observables v and v s . The latter is formulated in such a way that it can be evaluated efficiently using Monte Carlo methods.
As we will discuss in a little while, a good choice for such a simple observable is v s = v max , where the simple observable is defined by taking a suitably defined maximum of the observable calculated for independent emissions. We use the following notation from now on. A detailed definition of this observable will follow shortly. Using some trivial manipulations, one can write where we introduced the small positive parameter δ 1 that is independent of the observable's value v. An important comment at this stage is in order. Eq. (2.4) is strictly valid only for observables that admit a resummed cross section of the type (1.2). This is not the case for observables which can vanish also in the presence of resolved emissions due to kinematic cancellations (for instance for the transverse momentum of a color singlet in pp collisions). In this case Eq. (2.4) takes the form of a convolution between Σ max and F. In this article we limit ourselves to observables that behave like Eq. (1.2) in the v → 0 limit, and leave the study of the above observables for a future publication.
The product of ratios gives the relation between the cross section of the desired observable Σ(Φ B ; v) and the cross section of the simplified observable Σ max (Φ B ; v), and it accounts for the exact behavior of the observable in the presence of radiation. For this reason it is normally referred to as multiple-emissions functions. For the sake of brevity, we have dubbed it transfer function An important property of this transfer function is that it is IRC and UV finite, and as long as Σ max (Φ B ; v) has the same LL structure as Σ(Φ B ; v) (as we are assuming), its contribution starts at NLL (as will be shown later). The small parameter δ was introduced in order to allow for the transfer function to be easily calculable in a MC framework, as we will discuss below. The basic idea is that in the second ratio, Σ(Φ B ; v)/Σ max (Φ B ; δv), the denominator removes the unresolved emissions with V < δv, such that this ratio is IRC finite. The resulting dependence on the resolution parameter δ is cancelled against the first ratio, which can be calculated analytically once Σ max is known.
The first step in computing a resummed expression for Σ(Φ B ; v) is to build an explicit logarithmic counting for the squared matrix elements. To achieve this, one introduces the n-particle correlated (nPC) squared matrix elements |M (k 1 , . . . , k n )| 2 , which are defined recursively as These denote the parts of the squared amplitudes with n real emissions that can not be obtained by multiplying together squared amplitudes with less than n real emissions, and therefore represent the fully correlated part. With these definitions, the renormalized squared amplitude for n real emissions can be decomposed as Each of the correlated squared amplitudes admits a perturbative expansion where the index j denotes the order of virtual corrections to the squared amplitude with n real emissions. The rIRC safety of the observables considered here guarantees a hierarchy between the different blocks in the decomposition (2.7), in the sense that correlated blocks with n particles generally start contributing at one logarithmic order higher than correlated blocks with n − 1 particles [10,13].
Having introduced the nPC decomposition of the squared matrix elements allows one to precisely define the simplified observable V max . It is defined to be the maximum value of the full observable V calculated on the sum of momenta in each correlated block. In equations, this becomes (2.9) where V max (Φ B ; k 1 , . . . , k n ) is defined through its action on the n-particle correlated blocks as whereṼ (Φ B ; q) denotes the functional form of the observable on a single emission. It is important that the content of each correlated block is treated inclusively when computing the relativeṼ so that the above definition is collinear safe and can be extended at all orders. It is obvious that, in general, the cross section Σ max has no physical meaning, but it simply defines one ingredient for our resummation approach.
It is worth stressing that one has some freedom in choosing the form ofṼ . Conceptually, the simplest choice is to setṼ = V evaluated on the inclusive content of each block. In general, however, the only important feature is that it shares the same leading logarithms with the observable we are ultimately interested in. It is therefore very useful to define the simple observable such that the corresponding Σ max can be used for a whole class of more complicated observables. This can be achieved, for instance, by using the soft-collinear approximation of the full observableṼ = V sc instead of its full form V . This indeed, besides simplifying further the computation of Σ max , guarantees that this ingredient can be directly used for the resummation of all observables that share the same soft-collinear limit for a single emission, which defines a much broader class than the first definition given above. For the sake of simplicity, however, we avoid this technical complication in the rest of this article, and refer the interested reader to refs. [10,13] for more details.
We now continue with the derivation of the master formula. Using Eq. (2.9) together with Eqs. (2.7) and (2.10), Σ max (Φ B ; v) can be written in terms of the nPC (j) blocks. Since the observable acts separately on each nPC block, the expression in Eq. (2.9) can easily be shown to exponentiate. To perform a N k LL resummation for global, rIRC-safe observables, one needs to include nPC (j) blocks with n + j ≤ k + 1, but additional simplifications might be made on each nPC (j) block. For example, to LL accuracy, only the 1PC (0) block is required, and one only needs to keep the soft-collinear limit of it. Thus, one obtains where the 1/n! prefactor accounts for n identical gluons in the final state. The LL radiator function R LL (Φ B ; v) is the combination of the virtual and real contribution which at this order simply reads The definition of V max ensures the exponentiation at higher orders as well, such that one can always write where R(Φ B ; v) is called the radiator function, and Σ 0 max denote constant terms. Σ 0 max differs from one starting at NNLL order. Using the expression for Σ max (Φ B ; v) obtained just above, the first ratio in Eq. (2.5) can easily be computed, and we give the explicit expression when deriving the results at NLL accuracy below.
To compute the second ratio in Eq. (2.5) we need to carefully define the notion of resolved (unresolved) momenta, by demanding that the value V max evaluated on the set of momenta is above (below) a resolution scale δv (2.14) A key point now is to notice that for δ → 0 one can neglect the unresolved real emissions in the observable measurement function as they are much softer and/or more collinear than any other resolved emission in the final state. rIRC safety then guarantees that where p is a positive and real constant, independent of v. This allows one to split the total momentum q of each nPC block in Σ(Φ B ; v) into a resolved and unresolved component (depending on whether the value ofṼ (Φ B ; q) is greater or less than δv). 2 Thus, for each nPC block we use and so on. Up to power corrections in the small parameter δ, this allows us to separate Σ(Φ B ; v) into a resolved component (where all momenta are resolved) and an unresolved component, One therefore finds for the second ratio in Eq. (2.5) (2.18) 2 A second comment is in order. Once again, for observables that can vanish because of kinematic cancellations (a primary example being the transverse momentum of a color singlet in pp collisions), our choice of the simple observable can have issues when the above cancellations occur. A more appropriate, and more general, prescription would be to use δV (q1) as a resolution scale, where q1 is the total four momentum of the hardest correlated block. In this case Eq. (2.4) takes the form of a convolution as discussed in Refs. [12,39]. We will however proceed with the initial choice in the rest of this article for the sake of simplicity, since all of the other considerations made here are fully general. Note that the above discussion holds to any logarithmic accuracy. To go to a given order in the resummation of Σ(v) or Σ max (v) one needs to rewrite the full matrix element in terms of the nPC (j) blocks, and only keep the blocks that are relevant at the desired logarithmic order.
For the two ratios required in the transfer function Eq. (2.5), the argument of the numerator and denominator scale with the observable v. This implies that to compute the ratio to a given logarithmic accuracy, one needs the numerator and denominator at one logarithmic order lower [10,13]. To understand this fact better, let us consider the first ratio in Eq. (2.5) as an example. At NLL order, we can write where we have dropped all terms contributing beyond NLL. One can clearly see that the result depends only on the LL function g 1 (α s L v ), such that each term is only required to LL accuracy. Furthermore, one can perform additional kinematical expansions to simplify the expressions of the nPC (j) blocks, and we decompose each block nPC (j) by singling out its most singular (hence leading) term [nPC (j) ] sc , that is obtained by taking the soft and collinear limit of all emissions, i.e.
In summary, the ingredients required to a given order in logarithmic counting are summarized in Table 1. In the next section we perform the calculation at NLL for the thrust event shape.

Resumming the thrust distribution to NLL accuracy
In this section we compute all ingredients necessary to obtain Σ(Φ B ; τ ) for the thrust distribution to NLL accuracy, using Eq. (2.4). The thrust distribution is an additive observable, which satisfies The first ingredient is the NLL resummation of the simplified observable. To NLL accuracy one obtains [using the obvious extension of the LL result given in Eq. (2.11)] where one keeps the full kinematical dependence in the tree level contribution of |M (k i )| 2 , but only the soft-collinear limit of the one-loop contribution to |M (k i )| 2 and of One can evaluate the integrals involving |M sc (k a , k b )| 2 in dimensional regularization, and neglecting NNLL corrections one finds For the computation of the transfer function one only needs to keep the 1PC (0) block in its soft-collinear limit. Therefore, the first ratio in F NLL (τ ) can be written as where R LL (Φ B ; v) was given in Eq. (2.12).
To compute the second ratio of the transfer function to NLL accuracy, one uses Eq. (2.18), which leads to (up to power corrections in δ) Trading the 1/n! in Eq. (2.28) with an ordering in v i , this allows us to write the transfer function in the form where to obtain the second identity we have used this strongly resembles the standard parton shower evolution. It is therefore solved using the usual parton shower algorithm: 1. Start with i = 0 and τ 0 = τ 2. Increase i by one 3. Generate τ i randomly according to 3 4. If τ i < δτ exit the algorithm, otherwise go back to step 2 If the sum over all generated τ i is less than τ , accept the event, otherwise reject it. The value of F(τ ) is equal to the fraction of the accepted events.

Neglecting subleading effects: the CAESAR formula
The expression for the transfer function obtained in the previous section contains effects beyond NLL. It is often useful to be able to neglect all subleading effects and hence have a pure NLL answer. This can be done through a set of simplifications that we briefly summarize below. We stress that the operations performed in the present section are not, strictly speaking, necessary, but they can simplify considerably the numerical evaluation of the transfer function, and allow for an analytic solution in some cases. There are two important sources of subleading corrections in the treatment presented in the previous section. First, since the relevant squared amplitudes in the transfer function 1PC (0) are taken in the soft-collinear limit, it is natural to also approximate the observable V in the same limit. It is convenient to parametrize the emissions' momenta as where κ t,i is a space-like four-vector κ t,i = (0, k t,i ), orthogonal to the two reference momenta p 1 and p 2 that are aligned with the thrust axis n T Finally, we introduce the emission's rapidity η i with respect to the thrust axis, which is given by where the boundary for η i is obtained by imposing z ( ) i < 1 for any leg = 1, 2. Using the additivity of thrust one finds The above expression for the observable can be used in the evaluation of both the Sudakov radiator and the transfer function, neglecting terms beyond NLL order. Starting again from Σ NLL max , we can evaluate the integral (2.24) by parametrizing the phase space [dq] in terms of the transverse momentum q t and rapidity η of the emission q in the centre-of-mass frame of the emitting dipole. For the soft-collinear contribution one finds where the sum runs over the two Born emitters. The remaining hard-collinear contribution can be recast as [10,13] [dk]|M Using the above integrals, one can express where L = ln 1 τ and the expressions of g i have been long known [9] and are summarized in Appendix A. The derivative of the LL radiator, required in the transfer function, is then given by The second source of subleading corrections in the formulation of Section 2.2 has to do with the phase space bounds of the resolved radiation. In particular, we see from Eq. (2.37) that |η where in the last step we used Eq. (2.36). At NLL the upper rapidity bound can be approximated as which is then common to all resolved emissions. In our notation, this operation amounts to Taylor expanding the functions R LL (Φ B ; v i ) in the resolved radiation as where all terms in the r.h.s. beyond the first one are logarithmically subleading (each extra derivative suppresses the contribution by one logarithmic order). Similarly, the first ratio in the transfer function can be expanded about τ , in order to retain only the actual NLL terms necessary to cancel the δ dependence of the resolved radiation With these simplifications, we can recast the transfer function as which can be written as Eq. (2.45) is purely NLL, and does not contain any correction of subleading logarithmic nature, but still has the same general form as Eq. (2.30). The algorithm to compute it simplifies considerably: 1. Start with i = 0 and τ 0 = τ 2. Increase i by one 3. Generate τ i randomly according to 4. If τ i < δτ exit the algorithm, otherwise go back to step 2 If the sum over all generated τ i is less than τ , accept the event, otherwise reject it. The value of F(v) is equal to the fraction of the accepted events. The form of the transfer function can be manipulated further in order to make its numerical evaluation more efficient by getting rid of the Θ function in Eq. (2.46), as shown in refs. [10,13].
There is a second advantage of using Eq. (2.46) rather than Eq. (2.30). We notice that the starting equation (2.30) involves the function R (and hence the running coupling) evaluated at scales Qτ i that can get as small as Qδτ . When δ → 0 the above scale hits the Landau pole of the theory, which requires a prescription to deal with the non-perturbative region (e.g. a cutoff or a non-perturbative model) if this equation is implemented in a Monte Carlo method. On the other hand, the final equation (2.46) does not have this issue since we expanded the arguments of the couplings about Qτ Λ QCD , hence avoiding the Landau pole as long as the observable τ is sufficiently large.
For an additive observable such as thrust considered here, further manipulations are possible to obtain an analytic solution which reads , (2.47) which leads to the following NLL formula for the thrust cumulative distribution (2.48)

Review of resummation in SCET
SCET [2][3][4][5] is an effective field theory of QCD constructed to capture the long distance physics arising from soft and collinear radiation. To describe these different types of long distance effects requires two separate types of fields in the effective theory: soft and collinear fields. All short distance physics is integrated out of the theory, and contributes only via short distance matching coefficients.
Given that SCET has several degrees of freedom and exhibits a rich gauge structure, a detailed derivation of it is beyond the scope of this work and we refer the reader to the original literature [2][3][4][5] for details. One important feature, however, is that by defining the collinear fields in an appropriate way [5], the SCET Lagrangian can be written in a way that at leading power the collinear and soft degrees of freedom can be completely separated, giving where the soft Lagrangian is identical to the full QCD Lagrangian. In the following we are going to use the collinear fermionic Lagrangian which can be written as Here ξ n denotes a collinear fermion field after a so-called BPS [5] field redefinition, and the derivatives D n are covariant with respect to collinear gauge transformations and therefore only include collinear gluons fields. Note that a single collinear fermion field ξ n can be made invariant under collinear gauge transformations by combining it with a collinear Wilson line W n to define Operators in SCET are typically constructed out these gauge invariant fields. The Feynman rules that are obtained from the SCET Lagrangian are given in Fig. 1. The starting point for resummation in SCET is the derivation of a factorization theorem that expresses the cross section as a combination of contributions arising from three different singular sectors: hard, soft and collinear. Although such a type of separation is already performed at the level of the SCET Lagrangian [4], the observable under consideration mixes the various soft and collinear modes in its definition. Therefore, in order to derive a factorization theorem one must decompose the observable into soft and collinear contributions [6,7], which can be treated separately and then combined to give the final value of the observable. This allows for a separation between the phase space of the soft and collinear sectors, and therefore makes the factorization manifest. A clear complication arises for complex observables, for which the separation between soft and collinear modes in the measurement function can become quite cumbersome.
We express a generic factorization theorem as The hard function H n 1 n 2 only depends on the directions n i , but is independent of the observable. The jet functions describe the evolution of the radiation collinear to the directions n i , and the soft function describes the soft interaction between the two jet functions. The precise definition of the jet and soft functions, as well as of the convolution ⊗ in Eq. (3.4), depend on the definition of the observable whose value is required to be less than v in the integrated cross section Σ(v). In this paper we will need two types of observables. The first is the thrust observable we intend to resum, which is an additive observable for which the total value of the observable is the sum over the contributions from each particle. The factorization formula for such an additive observable takes the form [45][46][47] The two collinear directions n 1 and n 2 are back to back along the thrust axis t such that n 1 = n, n 2 =n with n = (1,t),n = (1, −t) and n·n = 2. Suppressing the dependence of the hard and soft function on the directions n andn, we write Σ(τ ) = H(µ) dτ n J n (τ n ; µ) dτn Jn(τn; µ) dτ s S(τ s ; µ) Θ[τ > τ n + τn + τ s ] . (3.6) We will also need an expression for the simple observable used to define Σ max in the previous section. This is defined by first grouping the various collinear and soft emissions separately into clusters in an infrared and collinear safe manner, computing the observable from each cluster and taking the maximum value of those. Such an observable factorizes in a multiplicative way, such that no convolutions are required We start by discussing the resummation for the additive observables described by the factorization formula (3.5), and then we comment in more detail on the definition of Σ max in SCET.
The soft and jet functions that appear in the above factorization theorems have the operator definition [46,[48][49][50][51] where J n (τ n , l + ; µ) / n αβ and Y n (x) denotes a soft Wilson line along the n direction. V soft , V n and Vn denote the expression of either thrust V or the simple observable V max as function of the final state momenta in the soft and collinear approximations, respectively. For notational simplicity, from now on we will omit the trace operation as well as the 1/N c prefactor in the color average of the above expressions, which will be understood in the rest of this article.

Resummation via Renormalization group equations
Once a factorization theorem has been obtained, one can use the renormalization group equations to resum the logarithmic dependence in the various contributions to the factorized cross sections. For this to work, it is crucial that each contribution depends kinematically on only a single scale µ F . This ensures that the logarithmic dependence in each contribution is directly tied to the dependence on the renormalization scale, since it can only occur in the form ln(µ/µ F ). It immediately follows that each contribution is free from logarithmic dependence if one chooses µ = µ F (the initial condition), and that the logarithms can be resummed using the RG equations. Before we discuss this in more detail, we take a short digression and discuss a feature of SCET that will be important later. In SCET, both the physical phase space and the observable's measurement function are expanded out according to the scaling of soft and collinear modes, since it ensures that each ingredient in the factorization formula depends on only a single scale 4 . Written in terms of the invariants y qg = s qg /Q 2 and yq g = sq g /Q 2 , the matrix element squared of the real radiation behaves as 1/(y qg yq g ), such that divergences arise both in the IR (y → 0) or UV (y → ∞) limit. To understand the consequences of this, we investigate the phase space boundary of a single emission, which are given in full QCD as QCD : where y ij denotes both y qg and yq g . They are shown graphically in Fig 2a). Clearly, neither of the two Mandelstam variables can exceed the physical bound set by the total energy in the event Q 2 , and therefore the phase space integration over each variable is bounded from above. The phase space boundary of the soft function in SCET is obtained by expanding the full QCD phase space boundary about the limit y qg , yq g 1. This gives which is shown graphically in Fig 2b). This implies that the larger of the two Mandelstam variables y qg or yq g is unbounded from above, leading to a UV divergence. The first collinear limit is obtained by taking the limit y qg yq g ∼ 1 (the second is the same under the replacement y qg ↔ yq g ). This gives (3.12) The collinear regions are shown by the hatched region in Fig 2 c) and d). In this case both variables are bounded from above, just as in the case of the full theory. However, adding the soft and collinear regions naively, leads to a double counting of the soft-collinear region [53], which is handled in SCET by subtracting a 0-bin region from the collinear integrals, which is nothing but the soft limit of the collinear integral. The soft limit of the first collinear phase space region (with the obvious replacement to the obtain the soft limit of the second collinear phase space region) is given by such that the integral over yq g is again unbounded from above, leading to a UV divergence. Diagrammatically, the 0-bin regions are summarized by the gray region in Fig. 2 c) and d).
While UV divergences are present in SCET as just discussed, each of the terms in the factorization formula Eq. (3.5) is IRC finite. Thus all divergences are of UV origin and are removed by renormalization. The renormalization of the UV divergences leads to renormalization group equations (RGE) for each component. As already discussed, each of the ingredients of the factorization theorem has its own characteristic scales that we denote by µ H , µ J , and µ S for the hard, jet, and soft functions respectively. At these scales no logarithms are present to any order in perturbation theory. For the thrust observable considered in this work, the scales are [16,47] (3.14) The resummation in SCET is then performed by evolving the hard, soft and jet functions from their characteristic scales to a common renormalization scale µ. The evolution is simply obtained by solving the corresponding RGE. The hard function is always multiplicatively renormalized, giving the following evolution equation whereJ andS denote the Laplace transform of the jet and soft functions, u is the Laplace variable conjugate to τ , and u 0 = e −γ E . Since the cross section Σ(v) is independent of the renormalization scale, the anomalous dimensions of the various pieces satisfy the consistency condition An analogous condition, trivially satistfied, holds for the terms in the anomalous dimension proportional to Γ cusp . We write the solution to the RGEs in Eqs. (3.15), (3.16) and (3.17) as where, as discussed above, all logarithms arise from the RG Kernels U F (. . . ; µ, µ F ). This leads to the final resummed Σ(τ ), which takes the form For an observable that is multiplicatively renormalized, such as Σ max (τ ), one finds the simpler expression The boundary conditions F (. . . ; µ F ), as well as the anomalous dimensions Γ cusp and γ F (for F = H, J, S) have a perturbative expansion The logarithmic accuracy is determined by the perturbative order with which the anomalous dimensions and boundary conditions are determined. For example, to achieve LL accuracy, one only needs Γ SCET and resummation based on factorization theorems in general is extremely powerful. Since higher logarithmic accuracy is achieved simply by computing anomalous dimensions and boundary conditions at higher perturbative accuracy, progress in our ability to perform fixed order calculations directly leads to higher logarithmic resummation, and some of the highest logarithmic accuracy has been achieved for several observables using this approach. The main drawback is that only observables for which a factorization theorem is known can be resummed using this approach. Deriving such a factorization theorem is often quite complicated, and for many observables it is not known.

NLL resummation for thrust
In this section we give the result for the thrust distribution at NLL accuracy, repeating the example of Section 2.  We will perform all required calculations at 1-loop order, but will include the 2-loop cusp anomalous dimension when giving the final result. We parametrize the generic momentum q as and define In the following we also use We start with the computation of the soft function. The diagrams contributing to the one-loop corrections are reported in Figure 3. The virtual correction is given by a scaleless integral, and hence vanishes in dimensional regularization, allowing us to set IR = UV = . Conversely, the real correction (plus its conjugate) is obtained by cutting the gluon propagator. Using the Feynman rules given in Fig 1 this gives After renormalizing the coupling in the MS scheme (α s (4π) → α s (µ)e γ E ) we take the Laplace transform of the result and expand it in α s (µ) to obtaiñ where u 0 = e −γ E . One can see that this soft function does not contain any logarithms at the characteristic scale which corresponds to µ S = Qτ in thrust space. Next, we compute the jet function along the n direction (analogous considerations apply to Jn), whose one-loop corrections are given by the diagrams of Figure 4. Virtual corrections are again given by scaleless integrals, so the only non-vanishing contribution is obtained by cutting through the loop in the diagrams of Figure 4. The sum of the diagrams (a) and (c) can be obtained by using the SCET Feynman rules reported in Fig. 1. We obtain (using l − Q) where we have defined x = k − /Q. The calculation of the remaining two diagrams ((b) and (d) in Figure 4) can be simplified further by noticing that their sum is related to the QCD wave function [3] as follows where the projectors P n and Pn read The result therefore reads Since the integrals above include a contribution where the momentum k µ becomes soft (and these effects have already been included in the soft function), this soft contribution needs to be subtracted. In SCET this procedure is called zero-bin subtraction [53], but in this case is given by scaleless integrals and hence vanishes. Combining Eqs. (3.31) and (3.33) (after the usual MS renormalization), performing the Laplace transform and expanding the result in α s (µ) one finds (3.34) One can see that the jet function does not contain any logarithmically enhanced terms at the characteristic scale which corresponds to µ S = Q √ τ in thrust space. The 1/ divergences are of UV origin, and in Laplace space can be renormalized with a multiplicative renormalization constant as follows where we defined and Z ψ is the wave-function renormalization By imposing the RG invariance for the bare soft and jet functions one can obtain the RGE of the renormalized ones  (3.34) in thrust space, which contain plus distributions. In this case the resulting RGEs take the form reported in Eqs. (3.16) and (3.17).
Using only these one-loop results, the solution to the previous RGE reads (where (3.40) where at NLL the initial conditions read Now Eq. (3.40) can be inverted to thrust space. One can decide to set the scales as in Eqs. (3.30), (3.35) and perform the Laplace transform or, rather, to first perform the inverse Laplace transform with symbolic µ S and µ J and then set the scales to µ S = τ s Q, µ J = √ τ n Q directly in thrust space. The difference between the two procedures is subleading, therefore we choose the latter which yields , (3.43) where µ S = τ s Q, µ J = √ τ n Q and Combining all results together, setting the common renormalization scale to µ = µ H = Q, such that the hard function contains no logarithmically enhanced terms and to NLL order can be set to unity, and including the 2-loop cusp anomalous dimension, one obtains , (3.46) where the expressions for the anomalous dimensions are reported in Appendix A. After evaluating the integrals in the exponent, and neglecting terms beyond NLL, one finds , (3.47) where the functions g i are reported in Appendix A. One can easily show that the above equation is equivalent to the QCD result of Eq (2.48) by writing which is equal to R LL (Φ B ; v) given in (2.29) Before moving on, we report the result for Σ max (τ ), which enters as an ingredient of the decomposition that will be used in Section 4. The simple observable used to define Σ max (τ ) is such that its UV divergences can be renormalized in a multiplicative way in thrust space, that is, the corresponding factorization theorem is multiplicative (see Eq. (3.7)). To the order we are working, the resulting soft and jet functions are trivially obtained from the Laplace space results reported above by simply evaluating them directly in thrust space, i.e. S max (τ ; µ) =S(u = u 0 /τ ; µ) J max (τ ; µ) =J(u = u 0 /τ ; µ) . (3.50) This gives At higher orders the initial conditions in Laplace space are different than they are in thrust space, such that the Eq. (4.23) is no longer exactly correct. To obtain the correct expression requires to perform the calculation of S max and J max directly in thrust space according to the factorization theorem (3.7).

Neglecting subleading logarithmic effects
The exact definition of the logarithmic order in resummation is somewhat convention dependent, and different prescriptions can be found in the literature. The prescription given in the previous section in Eq. (3.47) includes in fact various subleading logarithmic terms. For example, the cusp anomalous dimensions at 2-loop order is only required for the contribution in the first line of Eqs. (3.16) and (3.17), while in the second line it is enough to include the cusp anomalous dimension at 1-loop order. This implies that, instead of using the full expression for η ≡ 2η j + η s in the term e −γ E η /Γ(1 + η), one can perform a Taylor expansion of this result. For example, to NNLL accuracy one has Also, in general one finds differences depending on how the RG equations are solved. As already mentioned, performing resummation to a given order in Laplace space and then inverting the Laplace transform, gives results that differ beyond the order one is working compared to solving the RG equations directly in thrust space. A second example is that resumming the thrust distribution dσ/dτ (by setting the scales to the characteristic scales of the distribution) and then computing Σ(τ ) by integrating over 0 < τ < τ yields results that again differ at higher logarithmic order from those obtained by directly resumming the distribution Σ(τ ). For a detailed discussion of differences in logarithmic counting, see [56].
This existence of different conventions needs to be kept in mind in the next section when comparing the results obtained from an automated SCET resummation with the analytical results. In particular, a consistent comparison between different approaches can be only performed up to formally subleading terms.

Automated resummation in SCET
The starting equation for the automated resummation in Section 2 was the separation of the desired cross section Σ(v) into the product of the simplified cross section Σ max (v) and the transfer function F(v) given in Eq. (2.5). The resummation of the simplified observable was computed analytically, while the transfer function could be obtained numerically. In this section we derive a similar result, but where all ingredients are defined within SCET.
To simplify the discussion, we consider here a factorizable observable (such as thrust) and perform a similar decomposition at the level of the individual soft and jet functions. The SCET factorization theorem (3.5) for the thrust event shape that can be recast as (note that we drop the Φ B dependence from now on) where we expressed the soft and jet functions as (with F = S, J n , Jn) Next, we define .
The goal is to compute each of the transfer functions through a MC algorithm defined uniquely in terms of either soft or collinear fields, in a way that is similar to Section 2. We will show in Section 4.2 that in the framework of SCET one can compute each of the transfer functions F J (τ n , τ, µ) and F S (τ s , τ, µ) through a separate MC. This ensures that all observable dependence is restricted to the numerical MC algorithm.
The computation of Eqs. (4.4) via MC methods requires that each can be obtained in 4 dimensions by recursively computing real emissions. This relies on two important facts: First, the transfer function has to be determined entirely through the real radiation, and second, each contribution needs to be finite in 4 dimensions. The first fact is trivially satisfied, since in the ratios Σ F (τ )/Σ max F (δτ ) the purely virtual corrections cancel exactly. The second requirement deserves some closer investigation.
The IRC divergences cancel quite trivially in the ratio Σ F (τ )/Σ max F (δτ ), since the numerator and denominator include the same unresolved real radiation (for rIRC safe observables). However, as we discussed in Section 3.1 and contrary to full QCD, in the standard formulation of SCET real radiation is UV divergent. The resulting UV divergences of the real radiation appear both in the soft and in the jet functions and they cancel entirely only in their combination to give the physical cross section. The existence of the above divergences is a feature of the effective theory formulation in which the UV bounds of the theory are completely integrated out into Wilson coefficients. This guarantees that each of the soft and jet functions only depends on a single characteristic scale, which allows for the resummation of the dominant logarithmic terms via RG equations. In the usual formulation of SCET the UV divergences from the real radiation are regulated using dimensional regularization, and they constitute a crucial contribution to the anomalous dimensions which resummation is based on. However, the presence of the additional UV divergences prohibits a MC formulation of the problem that requires the phase-space integrals of the real radiation to be computable in 4 dimensions.
We solve this problem by introducing an explicit UV regulator for real phase space integrals into SCET. In this formulation of SCET, the UV divergences from virtual diagrams are regulated in dimensional regularization, just as before, while those from the real radiation are regulated with an alternative regulator, which can be chosen to be either a physical cutoff or an analytic regulator. This will give rise to a different RG structure in SCET, resulting in different logarithmic structures for the soft and jet functions individually. However, when soft and jet functions are combined into physical observables, one reproduces the same result as in the standard SCET formulation. By introducing such a regulator, we make sure that the UV divergences in the real radiation are now regulated in 4 dimensions, hence allowing for a formulation of the resummation through a MC algorithm.
In Section 4.1 we discuss the standard resummation in SCET in the presence of this new UV regulator for real-emission phase space integrals. We perform an explicit computation of the relevant soft and jet functions at one loop, and we show how the resummation can be performed through RG evolution. In Section 4.2 we show how to formulate a MC solution to the corresponding RG equations. We briefly comment on the extension to other observables in the conclusions, while the detailed generalization will be treated in a future publication. An alternative interpretation of the results presented in this section in the context of SCET is reported in Appendix C, where we comment on the structure of the theory when a IRC regulator δτ is included.

SCET with a UV regulator for real radiation
In this section we perform the calculation of the one-loop soft and jet functions by using an additional UV regulator for the phase-space integrals of the real radiation. This can be compared directly with Section 3.2, where the same calculations were performed without the additional UV regulator. The infrared and collinear divergences, as well as the UV divergences of the virtual corrections, are regularized by conventional dimensional regularization as before. One has some freedom in choosing the form of the UV regulator. In what follows we employ a cutoff Λ on the light-cone components of the emissions' momenta which is assumed to be larger than any other scale in the problem, which implies Λ ≥ Q . This mimics what happens in the full theory where the upper bound is set by the centerof-mass energy of the reaction. As a cross check, we have performed the calculations shown below using the exponential regulator proposed in Ref. [57], and found analogous results. As we will see, introducing this new regulator moves UV divergences between the soft and the jet functions, but of course does not affect the result after soft and jet functions have been combined into the total cross section. Since the UV divergences determine the RG equations, and therefore the logarithmic structure, this also implies that logarithmic contributions are moved between the soft and the jet functions. That is of course not a problem, since the separation into the logarithms of contributions from the various ingredients of the factorization theorem is to some extent arbitrary. Even in standard SCET one can move contributions between the different ingredients by changing the choice of the common renormalization scale µ.

The soft and jet functions at one loop
Consider the soft function of the factorization theorem given in Eq. (3.5) or Eq. (3.7). As before, the virtual contribution (plus its conjugate) is given by This integral is scaleless and therefore vanishes, hence setting UV = IR . The real contribution to the soft function is obtained by cutting the gluon propagator and imposing that the contribution to thrust from the real emission is smaller than τ s . This gives (remember that we impose k + , k − < Λ) After renormalization, we take the Laplace transform and expand in α s (µ). We obtaiñ From the above expression one can see that the soft function does not contain any logarithmically enhanced terms at the characteristic scales two separate evolution equations for each subprocess, the first of which describes the evolution in the dimensional regularization scale µ and the second one describes the dependence on the UV cutoff Λ. This is in spirit similar to what happens in SCET II problems [52] where a rapidity regulator is introduced to regularize the additional UV divergence of the real radiation [53,[57][58][59][60][61]. In fact, the same conclusions that follow would apply in that case. One finds for the soft and jet functions 5 where the lower bound of the Λ RGE arises from the fact that the corresponding anomalous dimension vanishes at this scale, which occurs at (4.16) From the system of equations (4.15) one can easily see that the Λ dependence cancels in the combination of the soft and two jet functions at a given µ.
One can now solve the RG equations by evolving the jet and soft functions simultaneously from their characteristic scales in µ and Λ to a common scale. Since the order of taking the derivatives with respect to ln µ and ln Λ commutes, this evolution is independent of the path chosen in the 2-dimensional µ − Λ plane. We can therefore write the evolution kernels relating the soft and jet functions at the characteristic scales µ F and Λ F to the common scales µ and Λ as From Eq. (4.20) we observe that the evolution of the jet function now starts at NLL. All double logarithms are entirely contained in the soft function, contrary to the case of standard SCET, where both the soft and jet function contained double logarithmic terms when evolved to the hard scale. However, as shown in Appendix B once the soft and jet functions are combined into a physical cross section, the logarithmic terms in the two formulations of SCET agree to all orders in perturbation theory, and one reproduces again the result given in Eq. (3.46). The logarithmic structure of the soft and jet radiation in this formulation of SCET reproduce exactly the physical structure in full QCD, allowing us to establish a one-to-one correspondence between the two formulations. Notably, this makes it possible to formulate the resummation in SCET via a Monte-Carlo approach, as it will be described in the next section.

Monte-Carlo resummation of the transfer functions at NLL
Using the decomposition in Eqs. (4.1) and (4.4), we first consider the soft transfer function. We start from the definition of the thrust soft function as the following expectation value where |k denotes a generic state with a fixed number of soft particles, e.g. |k = |k 1 for a single real emission, |k = |k 1 , k 2 for two real emissions and so on. In Eq. (4.22), V soft (p n , pn) denotes the expression of thrust in the soft limit, and p n (pn) denotes the sum of the momenta of the soft particles in the n (n) hemisphere. Our goal is to use Eq. (4.22) such that the soft transfer function F S (τ, τ s , µ) required in Eq. (4.5) can be computed via a MC algorithm. The first ingredient to evaluate Eq. (4.3) for the soft function is Σ max S . To the order we are working, the result for Σ max both for soft (jet) function is trivially obtained from the Laplace space results reported in the previous section by simply evaluating the soft (jet) function with the Λ regulator directly in thrust space, i.e. S max (τ ; µ, Λ) =S(u = u 0 /τ ; µ, Λ) J max (τ ; µ, Λ) =J(u = u 0 /τ ; µ, Λ) . (4.23) At higher orders the initial conditions in Laplace space are different than they are in thrust space, such that the Eq. (4.23) is no longer exactly correct. To obtain the correct expression requires to perform the calculation of S max and J max directly in thrust space according to the factorization theorem (3.7) with the additional UV regulator Λ. Eq. (4.23) leads to . (4.25) The first ingredient is the ratio Σ max S (δτ, µ)/Σ max S (τ, µ), which is obtained using Eq. (4.24). Since the leading logarithms cancel in the ratio of the Σ max S (δτ ) and Σ max S (τ ), one needs the resummed expression only to LL, which is given by Eq. (4.24) where only the one-loop cusp anomalous dimension is considered.
The second ratio Σ S (τ s )/Σ max S (δτ ) can be now computed numerically, as it is both IRC and UV finite. Indeed, the UV finiteness is guaranteed by the presence of the cutoff Λ in the real radiation, while the IRC finiteness is due to the rIRC safety of the observable that ensures that the radiation below the resolution scale δτ cancels out completely in the ratio. To achieve this, we introduce the decomposition for the squared amplitude in Eq. (4.22), similar to what was done in Section 2 We recall that the reason for the above decomposition is that squared amplitudesM S with n correlated real emissions start contributing at N n−1 LL to the evolution of the soft function for all rIRC safe observables. Using SCET with a UV regulator for the real emissions, as discussed in Section 4.1, these are now in clear correspondence with the QCD counterparts discussed in Section 2. Just as before, each of the squared amplitudes in the r.h.s. of Eq. (4.26) admits a perturbative expansion in powers of α s due to virtual corrections S (k 1 , . . . , k n ) . (4.27) The notation nPC S in Eq. (4.26) denotes the soft n-particle correlated blocks. In order to compute the transfer function to NLL accuracy, we only require the 1PC (0) S block, or in other words the squared amplitude |M S (k 1 )| 2 at tree level.
Putting all this together one finds (4.28) which can be evaluated with the same MC algorithm described in Section 2.2. Next we consider the jet function. As for the soft function, the expression for Σ max J (τ ) is immediately obtained from the results of Section 4.1.1 [see Eq. (4.23)] (4.29) The computation of the jet transfer function is trivial at NLL order. As we have seen in Section 4.1.1, the jet function is only single logarithmic once the additional UV regulator Λ has been introduced, since the only kinematic region of phase space giving rise to large logarithms is of hard-collinear origin. Given that the collinear sensitivity is the same in Σ J (v) and Σ max J (v), the resulting logarithmic dependence due to the phase space bounds cancels in their ratio, and the only logarithmic sensitivity in the jet transfer function comes from the running coupling constant. This implies that each additional emission is suppressed by an additional power of α s , such that only a finite number of emissions need to be taken into account at a given order N k LL. In particular, to NLL accuracy, the jet transfer function does not contribute for the reasons stated above, and one has the trivial result We can now combine the result for the two jet functions just computed with the NLL soft function as in Eq. (4.1) obtaining where we have performed the trivial integrations over τ n and τn. We note that the prefactor given by the product Σ max = Σ max S Σ max Jn Σ max Jn can be directly computed in the standard SCET without the need for the Λ regulator, whose dependence will completely cancel in the product of the three terms.
Using the same steps as in Section 2, one can neglect terms that only contribute to order NNLL and higher, such that one can write the above result in a way that allows for a simpler MC implementation. To this end, we define where we have evaluated the scale of the running coupling constant at k t = √ k + k − . This is the only available choice in the soft function differential in the two light-cone components, since it is the only possible scale which is invariant under a rescaling of the directions of the Wilson lines. F NLL S (τ, τ, Q) becomes which can be solved with the following MC procedure:  Although the extension to the general case is beyond the scope of this article, we do want to mention that it is possible to apply the above method to a more complicated observable than thrust. In general, if one is able to find an SCET Lagrangian for the simple observable and define Σmax which by definition contains the same LL as the full observable v, then the resummation for v can be obtained by means of a transfer function that is defined in terms of the fields of the same Lagrangian, and can be computed via Monte Carlo methods.

Conclusions and Outlook
In this work we have shown how to formulate a numerical approach to resummation in SCET using the example of NLL resummation of the thrust distribution. This was achieved by combining the automated CAESAR/ARES approach to resummation with the factorization of the long distance degrees of freedom in SCET.
In SCET, resummation is obtained by first factorizing the required cross section such that a process independent hard function (H) multiplies the convolution over jet (J) and soft (S) functions The jet and soft functions describe the long distance physics of the process and therefore contain all observable dependence. Each of the factorization ingredients depend on only a single scale, and logarithms can be resummed by solving RG equations for each of the factorization ingredients separately. This general approach makes resummation relatively straightforward once the appropriate factorization formula has been obtained, and simply requires the computation of anomalous dimensions at a given order in perturbation theory.
In the numerical approach introduced in this paper, we identify a simplified observable, which has the same leading logarithmic structure as the thrust distribution, and for which a factorization theorem can be built in SCET. This simplified (max) observable is constructed such that it has a very simple multiplicative factorization theorem, which is just the product of the same hard function (H) multiplied by jet (J max ) and soft (S max ) functions Due to this simple multiplicative form of the factorization theorem, resummation is achieved in a straightforward manner by solving multiplicative renormalization group equations. The ratio between the full and simplified jet and soft functions defines a transfer function where F = J, S. The main result of this work was to show how to compute this transfer function by performing the phase space integration over real emission diagrams to all orders in perturbation theory. To NLL accuracy, this was shown to result in a rather simple expression, which can be numerically implemented into a straightforward MC algorithm. In order to compute the phase space integrals numerically, we needed to ensure that they are finite in 4 dimensions. This is not the case in regular SCET, where the multipole expansion of the phase space limits of the soft function (and the 0-bin of the jet functions) leads to UV divergences in the real integration. In order to overcome this, we introduced an additional regulator to control the UV divergences in the soft real phase space integrations. While this modifies the UV structure of the theory, and requires to perform the RG evolution in two different variables, we showed that the results obtained in SCET with and without this extra regulator are in fact equivalent.
While we have focused for simplicity only on the NLL resummation of the thrust distribution, our results are very general and are readily extended to higher orders in resummation accuracy and to more general observables as long as one can find a simple observable that has the same LL as the full observable, and it is factorizable in SCET. Using the general definition of Σ max given in this article, this is an almost trivial task for most observables. In the approach discussed in this paper, the degrees of freedom required, and hence the effective Lagrangian, are determined by the simple observable which can be resummed analytically through a factorization theorem. The required transfer function, that relates the simple observable to the desired observable, can then be computed numerically using the Feynman rules of the above Lagrangian. Moreover, owing to the fact that the UV limit is now separately regularized, SCET II problems can be formulated exactly on the same footing as SCET I ones.
Higher-logarithmic accuracy can be obtained in a relatively straightforward manner by keeping subleading terms in the expansions performed in this work. Furthermore, the numerical approach to resummation in SCET applies even to observables for which a factorization theorem is not known. This opens the door to a systematic resummation for a wide class of observables by combining the analytical power of SCET with numerical MC integrations, which can be automated in an algorithmic way. The details of the general formulation are discussed in a forthcoming paper [43].
where we have used that Λ J = Λ H , such that one does not need any Λ evolution for the jet function, and we have defined The two anomalous dimensions γ S and γ J are different from the usual SCET ones starting from their NLO expression. However, they satisfy γ S + 2γ J = γ S + 2γ J . The integration over Λ in the Λ evolution kernelsŨ (Λ) F can be performed analytically by changing the order of integration. For this, we write This gives for theŨ  where we have used in the last line µ 2 J = µ S µ H . This combined evolution factor therefore contains the complete evolution due to the cusp anomalous dimension, which is usually split between the soft and the jet evolution kernels, as well as the non-cusp part of the soft evolution. The evolution factor of the jet function only contains the non-cusp part of the collinear evolution which reads which shows that the physical combination of the evolution factors is identical to the standard SCET one at all orders.

C. RGE of the thrust soft function with a IRC resolution scale
In this appendix we wish to comment more on the logarithmic structure of the decomposition (4.5) for the resummed cross section. It is instructive to consider the soft function as a case study, although the same conclusions apply to the two jet functions. We express the soft function as in Eq. (4.3), namely . (C.1) In this paper we have shown how Σ max S is computed analytically, while the transfer function F S is obtained via MC methods. The latter task requires the introduction of an explicit UV regulator in the theory, which led to the results in Section 4.2. The goal of this appendix is to study how the decomposition (C.1) modifies the soft function's RGE if it were used in the standard SCET formulation, i.e. without the Λ cutoff. This is a useful exercise to understand from a different viewpoint the method proposed in this article.

(C.2)
The full soft function for thrust in momentum space fulfills the non-local RGE [see Eqs. where the precise observable dependence (in particular its additive nature) is reflected in the second term in the right-hand-side of the above equation, that essentially shows how an extra real emission modifies the existing value of the observable (i.e. in an additive way). Now we study how the RGE is modified by the introduction of the IR resolution scale δτ s . We first consider the contribution Σ max S (δτ, µ), which fulfills the following local RG equation in momentum space [see Eqs. (3.19)  which can be easily solved by setting the initial conditions at µ = δτ Q.
Next we consider the remaining ratio Σ S (τ s , µ)/Σ max S (δτ, µ) which, by means of Eqs. Thus, the decomposition of (C.1) allowed us to separate the initial evolution equation into a local piece (C.4) and a non-local piece (C.8). The precise definition of the non-local piece depends on the form of the observable, and the result given here holds for an additive observable. The decomposition of Eq. (C.1) allows to separate the RGE into a local piece, which is independent of the definition of the observable and can therefore easily be solved analytically, and a purely non-local piece, which contains all the observable dependence. While we have discussed how to compute the non-local piece via a MC algorithm, for an additive observable it can easily be calculated analytically, as we now show. We take the Laplace transform L of the ratio defined as Σ S (τ s , µ) Σ max S (δτ, µ) By writing the logarithm in the exponential function as ln δτ Q µ = ln τsQ µ + ln δτ τs we obtain S(τ s ; µ) = exp µ µ 0 dµ µ 4Γ cusp (α s (µ )) ln τ s Q µ 1 τ s Q e −γ E η(µ,µ 0 ) Γ(η(µ, µ 0 )) , (C. 16) which reproduces the NLL result for the thrust soft function with µ 0 = τ s Q.