A formalism for the resummation of non-factorizable observables in SCET

In the framework of soft-collinear effective theory (SCET), we show how to formulate the resummation for a broad family of final-state, global observables in high-energy collisions in a general way that is suitable for a numerical calculation. Contrary to the standard SCET approach, this results in a method that does not require an observable-specific factorization theorem. We present a complete formulation at next-to-next-to-leading logarithmic order for $e^+e^-$ observables, and show how to systematically extend the framework to higher orders. This work paves the way to automated resummation in SCET for several physical observables within a single framework.


5.
General expressions for the fully differential transfer function 16

Introduction
Perturbation theory is one of the most widely used techniques to make predictions for interacting quantum field theories. In fixed order (FO) perturbation theory, one expands the physical observables in powers of the coupling constants of the theory, where leading order (LO) predictions describe a given process to lowest order, while higher order corrections are suppressed by additional powers of α i = g 2 i /4π. The dominant corrections are in most cases due to the strong interaction, where the expansion is in terms of α s . This approach is well suited for sufficiently inclusive observables. On the other hand, for observables that depend on several scales, it is well known that the perturbative expressions contain logarithmic dependence on the ratio of these scales r, with up to two powers of such logarithms for each power of the coupling constant (α n s ln 2n r). In this case the FO perturbative expressions can become unreliable and a different expansion is required to make accurate predictions. Resummed perturbation theory is a rearrangement of the perturbative series such that the dominant logarithmic terms are resummed to all orders in the coupling constant α s .
In this article we consider the family of global, recursively infrared and collinear (rIRC) safe [1] observables v whose leading order perturbative series features double logarithms in the limit v → 0. For such observables, the perturbative series of the cumulative distribution with σ B being the Born cross section, can be shown to exponentiate at the leading logarithmic order such that one can write Σ(v) = exp [Lg 0 (α s L) + g 1 (α s L) + α s g 2 (α s L) + . . .] , (1.2) where L = ln 1 v .
(1. 3) In the limit L ∼ 1/α s , such that α s L ∼ 1, each term in a Taylor expansion of the functions g i (α s L) is of the same size, such that arbitrarily large powers of α s are required to determine these functions. However, the contribution of the function g n is suppressed by one power of α s compared to the function g n−1 . The goal of resummation is to compute these functions, and the N k LL order is defined by the inclusion of the functions g m (α s L) with m ≤ k in the prediction (e.g. at NLL one needs g 0 (α s L) and g 1 (α s L)).
Two main approaches to calculate the resummed expressions exist. The branching formalism [1,7,12,35,36] uses the factorization properties of squared amplitudes in QCD to describe the radiation dynamics at all perturbative orders. This approach is widely used also in parton shower algorithms, which are typically accurate to (N)LL for specific observables [35,[37][38][39]. However, in contrast to a parton shower, in a resummed calculation one does not need to simulate physically sensible events which exhibit momentum conservation. This makes the branching formalism usable in principle to any logarithmic accuracy. The main idea relies on the computation of the squared amplitudes with an arbitrary number of soft or collinear emissions, but requiring them to be correct only to a given logarithmic accuracy. With these simplifications one can formulate a Monte-Carlo (MC) algorithm to simulate the radiation above a certain resolution scale q 0 , while the emissions below q 0 are treated analytically. For sufficiently simple observables one can rewrite the branching formalism in terms of differential equations, which can be solved in closed form.
An alternative approach is based on the derivation of a factorization theorem for the observable under consideration, in which the cross section is decomposed in terms of ingredients that are sensitive to a given infrared scale. Logarithms are then resummed through evolution equations of the various ingredients. Collins, Soper and Sterman [40] pioneered this technique with their seminal calculation of the transverse momentum distribution of the vector boson in W and Z production. Subsequently, this approach was formulated [13] in soft-collinear effective theory (SCET) [41][42][43][44]. In SCET, the interactions between particles can be factorized at the level of the Lagrangian [43] and, in the case of observables with factorizing measurement functions [45,46], one can separate the cross section into contributions from hard, collinear and soft degrees of freedom. In dimensional regularization, each element of the factorization theorem then only depends on a single scale through its ratio to the renormalization scale [41,42]. Therefore, the renormalization (RG) group evolution can be used to resum the logarithmic dependence of each of the ingredients in the factorization theorem, and the resummation is achieved via the computation of anomalous dimensions of SCET operators.
Using the thrust observable as a case study, in [47] it was shown how these two formulations of resummation can be combined into a hybrid method. This approach proceeds by defining a simpler version of thrust, defined in order to obey a particularly simple factorization theorem that can be handled using the standard SCET formalism. In the resulting factorisation theorem, the different ingredients are combined through a multiplication, rather than a more complicated convolution as in the case of the full thrust observable, hence leading to a much simpler structure of the RG evolution.
Given a resummed result for the the cumulative distribution of the simplified thrust observable Σ max (τ ), one could then obtain the actual thrust cumulative distribution Σ(τ ) as where the transfer function F(τ ) can be calculated numerically using an efficient MC algorithm.
In [47] it was argued that this procedure is rather general, and can be applied to a wide variety of observables, without giving the details of this generalization. The purpose of this paper is to derive this formulation for a generic recursive infrared and collinear (rIRC) safe observable v and to give all the details required for a numerical implementation up to NNLL accuracy. In particular, we derive the generalization of Eq. (1.4) namely give a particularly general definition of the simplified observable Σ max (v), and provide a detailed derivation of the transfer function F(v). The paper is organized as follows: We start in Section 2 with some comments about the counting of the logarithms in a resummed expressions beyond NLL accuracy. After defining the kinematics of the process and detailing the notation used throughout the paper in Section 3, we introduce the basic idea behind the numerical resummation in Section 4. In that section, we briefly touch on all aspects of the derivation: we introduce a generic and fully differential factorization approach derived from SCET, from which any observable can be calculated, define the simple observable and how it can be resummed using SCET and finish by an SCET definition of the transfer function. In Section 5 we then derive a general expression for a fully differential transfer function, from which the transfer function for any observable can be calculated. While in Section 5 the derivation is purposely kept rather general, we give explicit and detailed results for the transfer function in Section 6 at both NLL and NNLL accuracy. In Section 7 we discuss various ways one can simplify the observable dependence in the transfer function. We present our Conclusions in Section 8.

Remarks on the counting of logarithms
In Eq. (1.2) we presented the general form of the resummed expression for a cumulative distribution of a general observable. To obtain N k LL accuracy one needs to include the functions g n (α s L) with n ≤ k, which ensure a reliable theory description in the regime L ∼ 1/α s . Within this counting, NLL accuracy is defined by including the functions g 0 (α s L) and g 1 (α s L), which are required to obtain a prediction that is accurate up to O(α s ) corrections. Beyond NLL accuracy, one includes the functions g n (α s L) with n ≥ 2, which lead to progressively suppressed perturbative corrections.
This allows for different ways of consistently defining the resummed cross section beyond NLL, such that they differ from one another by terms beyond the considered order. Consider as an example the case of a NNLL resummation. The first scheme is to keep the full exponential form of Eq. (1.2), which we repeat here for convenience: (2.1) A second counting scheme can be defined by expanding the perturbative corrections introduced by the g 2 term in Eq. (2.1) in powers of α s , but keeping the scaling L ∼ 1/α s intact. This allows one to write It is easy to verify that the difference between eq. (2.2) and eq. (2.1) amounts to terms which contribute only to N 3 LL.
Finally, one can define a third scheme, which is a hybrid between the two schemes defined so far. In this scheme one writes whereg 2 (x) = g 2 (x). Therefore, this scheme keeps part of the NNLL corrections in the exponent, and expands out the remainder. Once again the difference between eq. (2.4) and the previous two schemes is of order N 3 LL. All of these counting schemes give equivalent predictions if α s L ∼ 1, since they only differ by terms that are of higher order than those considered. Generally, most observables do not naturally satisfy the fully exponentiated form of eq. (2.1), and some NNLL contributions (and beyond) are usually expanded out in a fashion similar to scheme (2.4). However, the exact definition of the functiong 2 (α s L) is intrinsically ambiguous, in that one can always reshuffle NNLL terms from the exponent to the expanded remainder and vice versa. More precisely, different approaches to resummation lead to different definitions of the exponentiatedg 2 (α s L) function, while yielding the same NNLL result for the resummed cross section. For instance, standard SCET resummations keep all constant terms and their own running-coupling corrections outside of the exponent, while some keep part of the constants (specifically those coming from the virtual corrections) exponentiated. In the approach of refs. [7,9,12], conversely, only the universal contributionsg 2 (α s L) are exponentiated, while the observable-dependent NNLL corrections are expanded out at fixed α s L.
This freedom can be exploited in the formulation of numerical approaches to resummation. This can easily be understood from the way numerical resummation approaches work by relating two different observables that have the same LL structure. Given two such observables Σ and Σ , whose logarithmic resummation is defined by the functions g i and g i with g 0 = g 0 , the two resummations can be related as This means that the difference between the g 2 and g 2 functions, rather than being exponentiated, can be evaluated at fixed order while keeping α s L fixed. This will considerably simplify the numerical algorithm. This feature has been exploited to formulate approaches that push the automatic NLL resummation of CAESAR [1] beyond NLL accuracy, in particular in the algorithms adopted in the programs ARES [7] or RadISH [29,30]. The above conclusions will be used in this paper to formulate a resummation method in the framework of SCET that can be implemented in a numerical algorithm. Note that the same considerations apply to any logarithmic accuracy beyond NLL, for which different resummation approaches predict resummed cross sections that differ by corrections beyond the considered order. It follows that a comparison between resummed predictions must be performed up to the nominal logarithmic order, while subleading differences may well be present.

Kinematics and notation
In this section we introduce a convenient notation to formulate the resummation in a way that is applicable to broad classes of observables. We consider e + e − collisions producing a Z/γ * which subsequently decays into hadrons, i.e. e + e − → q n qn + X. We begin by introducing the usual light-cone parametrization of four-momenta with respect to a lightlike direction n µ = (1, n) with n 2 = 0. This allows one to write wheren µ = (1, − n), such thatn 2 = 0, n ·n = 2. The momentum k t is a space-like momentum satisfying n · k t =n · k t = 0. We also define k 2 ⊥ = −k 2 t . A common choice is to align n with the z-axis, resulting in n = (1, 0, 0, 1), but the notation is generic.
Next, we introduce the notation that allows us to describe the fully differential energy distribution of the considered process. We consider a state X containing M massless particles with four-momenta k i = (E i , k i ). The total energy ω X and 3-momentum k X at a given solid angle Ω in the final state are given by where the solid angle of each particle is determined by its 3-momentum Ω i = Ω( k i ). One can now define a functional integration measure by discretization. We divide Ω into infinitesimally small bins {Ω k }, and define the set of discrete variables {ω k , k k } as the integrals of ω(Ω) and k(Ω) over the bins {Ω k }, Now we define an integration measure over the energies and 3-momenta at a given solid angle Ω k as Restricting the particles inside of each solid angle to be on-shell (which is justified since we have infinitesimally small solid angles Ω k ) one can write With this notation, the phase space for a state X N , containing N final state particle plus the two Born particles q n and qn is given by where the symmetry factor S(N ) takes into account identical particles in the final state (for example, for N gluons we have S(N ) = 1/N !).
Since a given energy density is identified by its momenta, as discussed above, the sum of two energy distributions is simply given by the the combined momenta of the two individual distributions One simple consequence of this result is that the total momentum is conserved where P measures the total 4-momentum in a given energy distribution. Any observable is defined by the way it acts on a set of final state particles. We consider observables whose value only depends on the momenta of the particles, and not on any additional quantum numbers. Such observables can be defined in terms of a fully differential energy distribution, following the discussion of [45,46].
We introduce the phase space Φ B of the underlying Born configuration, that in our case is defined by the two back-to-back quarks p andp prior to any emission. It is constructed in terms of the energy distribution through the total momentum P µ [ω (N ) ] = Q(1, 0) and the thrust axis n. The cumulative distribution for a given observable can be constructed from the projection of the cross section that is fully differential in the energy distribution of the event, onto the value of the observable. We express Σ(v) as 1 and define Here Φ B [ω (N ) ] denotes the direction of the initial (Born) quarks that can be entirely reconstructed from the final-state emissions k i and from the final (i.e. after all radiation has been emitted) quark momenta q n , qn. The integration measure Dω (N ) was defined in Eq. (3.6), and the observable for any fixed N can again be written in terms of the momenta k i , q n and qn (3.11) The fully differential cross section is given by the square of the matrix element as δσ δω (N ) ≡

2F
|M (q n , qn; k 1 , . . . , k N )| 2 , (3.12) where F denotes the flux factor, and |M (q n , qn; k 1 , . . . , k N )| 2 the squared matrix elements to all orders in perturbation theory with a fixed number N of emissions. We conclude the notation section by giving our definition of the running coupling constant that will be used throughout the paper. The dependence of α s on the renormalization scale is given in general by where (3.14)

Basic idea of numerical resummation for rIRC safe observables
To resum a cumulative distribution of a given observable in SCET, one typically proceeds in various steps. First, one needs to identify the degrees of freedom of the effective theory, which allows one to reproduce the singular behavior of the observable under consideration. Typically, this requires a combination of collinear modes along the directions of the hard particles of the underlying Born process, and soft modes, which can provide interactions between different collinear directions. The scaling of the collinear and soft modes depends on the observable. For example, for thrust (τ = 1 − T ) the scaling is for soft and collinear modes, respectively. These modes are then used to define the effective SCET Lagrangian [41][42][43][44]. Resummation is then normally achieved by factorizing the cumulative distribution into a hard matching coefficient, together with soft and collinear (jet) functions, and then solving RG equations for the various components of the factorization theorem. The details of this factorization theorem as well as that of the RG equations depend again on the observable. As stated in the introduction, in this paper we approach the resummation from a different angle. Within the framework of SCET, one proceeds as follows: 1. Given an observable V (q n , qn, k 1 , . . . , k N ), build a simplified version of the observable V max (q n , qn, k 1 , . . . , k N ), which shares the same leading logarithmic structure with V .
2. Identify the degrees of freedom for the simplified observable V max , and build the corresponding SCET Lagrangian. Use this Lagrangian to obtain a factorization formula for the fully differential energy distribution used in Eq. (3.10).
3. Decompose the resummed cumulative distribution as where Σ max is the cumulative distribution for V max , while F is a transfer function that relates the latter to the full observable.
4. The resummation is achieved by performing the resummation for Σ max (Φ B ; v) analitically using the usual SCET technology, and then computing the corresponding transfer function F(Φ B ; v) using the above Lagrangian.
As we will show, the factorization theorem for the simple observable in SCET, and hence its resummation, is straightforward and can be obtained in a general manner. The feature that makes this approach powerful is that the transfer function can be computed numerically for wide classes of observables. As a consequence, this technique can be used to resum observables without the need for an observable-dependent factorisation theorem, as it is usually done in SCET.
In the following sections we will go through the derivation of each of the above steps in some more detail.

Factorization of the fully differential energy distribution
The starting point for our formulation is the expression for the cumulative distribution given in Eq. (3.10). As was shown in [46], the fully differential energy distribution δσ/δω can be written in terms of the energy momentum tensor of SCET. Following the notation of [46], we define where X denotes a general state. The energy flow operator E 0 (Ω) can be defined in terms of the energy-momentum tensor as [45,48] E 0 (Ω) = lim Here, u ≡ u(Ω) is the unit three-vector pointing in the direction identified by Ω. Therefore, E 0 (Ω) measures the total energy arriving over time at infinity in the direction Ω. An explicit proof of Eq. (4.5) for scalars and Dirac fermions can be found in Ref. [45].
As it was shown in [46], given that the energy momentum tensor is linear in the Lagrangian of the theory, it can be written as where E 0 n ,s (Ω) are defined analogously to Eq. (4.5), but using the energy-momentum tensor obtained from the Lagrangian L ,s only. This then implies that the fully differential energy distribution can be written in the factorized form [45,46] Here |C(Φ B )| 2 denotes the matching coefficient describing the short distance fluctuations in the full theory that are not included in SCET, which depends on the underlying Born configuration of the process under consideration, but is independent of the definition of the observable. The terms δσ F /δω F denote the fully differential cross section as computed from the part of the SCET Lagrangian describing sector F ∈ {s, n,n}. In each sector, one can write In the soft sector, one finds where the virtual corrections V s and the real emission matrix element are computed using the Feynman rules of the given sector. The only difference with Eq. (3.12) is the absence of the differential Born cross section, which is contained in the matching coefficient |C(Φ B )| 2 .
In the collinear one instead has An important difference between the collinear and the soft sector is that in the latter there is still a single quark q n contributing to the energy distribution (ω n = ω[q n ]), even in the absence of any radiation. Combining Eq. (4.7) with Eq. (3.10) one can write where the sum of energy densities is obtained by taking the union of the momenta in each sector, as defined in Eq. (3.7). Note that Eq. (4.11) does not imply a factorization formula for Σ(Φ B ; v) in the commonly used sense, since the observable V [ω n + ωn + ω s ] does not factorize in general. Only observables which do not combine momenta from different sectors in a non-trivial way, such that they can be written as satisfy a commonly used factorization theorem. One then obtains (4.14) Thus, for factorizable observables one reproduces the factorized result that was the starting point in [47] when discussing the thrust distribution at NLL accuracy. For example, an additive observable, such as thrust satisfies Another important consequence of Eq. (4.11) is that, in general, the delta function constraining the total energy density to the sum of the ones in each sector introduces a kinematic cross-talk between the soft and collinear sectors. This can be understood by noting that Thus, the quark that is initiating each collinear sector recoils against the soft radiation in the corresponding hemisphere. This is important for observables that are sensitive to the recoil of the Born quarks, such as the jet broadening.

Definition of the simple observable
The next essential step is the definition of the simple observable V max . The only constraint on such an observable is that it needs to have the same leading logarithmic structure as the full observable V . Many choices are of course possible, but it is convenient to adopt a simple definition that is generic and can be easily handled analytically. A particularly simple and systematic procedure is to define a global and rIRC safe observable in each sector individually, and then to compute the final value of the simple observable as the maximum of the value in each sector. In other words, we define where the definition of V max [ω F ] in each sector F is given below.

The soft sector
Since the observable needs to be defined in resummed perturbation theory, one needs to define it for an arbitrary number of final state soft particles. We decompose the squared amplitude |M s (k 1 , . . . , k N )| 2 for producing a set of particles with momenta k 1 to k N into soft webs, that in the following will be referred to as correlated clusters 2 Each cluster is recursively defined as the portion of the squared amplitude that can not be written in terms of products of lower-multiplicity clusters. The 1-particle correlated clusterM 2 s (k) describes the square of the eikonal amplitude for emitting a single gluon from the Wilson line. For the emission of two gluons, the correlated clusterM 2 s (k 1 , k 2 ) is defined as the difference of the 2-gluon squared eikonal amplitude |M s (k 1 , k 2 )| 2 and the product of two 1-gluon squared eikonal amplitudes determined in the previous step. For the emission of a quark anti-quark pair, the correlated cluster is given by the full qq eikonal squared amplitude, since the 1-particle eikonal squared amplitude with only a quark does not exist. The extension to higher multiplicities proceeds in a straightforward fashion. Each correlated cluster admits a perturbative expansion arising from the virtual corrections, and we writeM Using this structure, one can define the simple observable V max in the soft sector through its action on products of correlated clusters One therefore only needs a definition ofṼ s (k), acting on the total momentum of each correlated cluster. A convenient choice is to consider the following generalisation of the moments of energyenergy correlation [1] 22) and to define the simple observable using the following procedure: 1. Associate each momentum k i (for correlated clusters with more than one momentum, k i corresponds to the total momentum of the correlated cluster) with the Wilson lines it is closest in rapidity. This effectively divides the event into two hemispheres.
2. For each hemisphere, compute the contribution of each correlated cluster to eq. (4.22) in the soft-collinear limit. If the correlated cluster contains more than one particle, i.e. it has a non-vanishing invariant mass, we evaluate eq. (4.22) in the massless approximation, i.e. we change the energy of the cluster k 0 such that k 2 = 0. This amounts to considering the energy-energy correlator between each of the clusters and the quark moving in the same hemisphere (after all soft radiation has occurred), that can be recast asṼ where k ⊥ and η are the transverse momentum and pseudo-rapidity of the considered correlated cluster with respect to the above quark, labelled by = {n,n}. 3 The final momentum of the two quarksq n ,qn after the emission of the soft radiation (note that these differ from the final state quarks q n and qn by the recoil associated with collinear emissions) is simply obtained by the Born momenta as where p µ andp µ are the Born momenta prior to any emissions. The directions of the two momentaq n andqn then define the direction of the Wilson lines used in the (recoil-free) definition of the SCET Lagrangian.
The parameters a and b must be chosen in order to match the soft-collinear scaling of the full observable V , which guarantees that the simple observable has the same leading logarithmic structure as the full observable.
3. Finally, define the simple observable v s ≡ V max [ω s ] as the largest of all of theṼ s (k) values.
The above definition corresponds to a toy observable that will lead to important simplifications both in its analytical resummation and in the numerical evaluation of the transfer function.

The collinear sector
As will be discussed later, the collinear sector only contributes non-trivially to NNLL accuracy and beyond. Given the discussion in Section 2, this implies that in the contribution of the collinear sector to the transfer function one needs to consider only a fixed number of emissions at a given logarithmic order, with r emissions needed to obtain N r+1 LL accuracy. Thus, to obtain NNLL accuracy a single collinear emission suffices, while at N 3 LL accuracy at most two collinear emissions are required. For this reason, a decomposition of the collinear squared amplitudes in terms of correlated clusters is not necessary, and we can rather define the simple observable in terms of its action on the full collinear squared amplitude. In each of the collinear sectors, say n, the simple observable V max [ω n ] is then simply given by eq. (4.22) in the small θ ij approximation, namely where the sum runs over all particles belonging to the collinear sector n. For a single emission, a useful parametrisation of v n ≡Ṽ n (k n ) is given in eq. (B.10). We note that this simplified observable in the collinear sector only depends on the momenta in the final state, but is independent of the axis used to define the collinear Lagrangian. The corresponding jet function is therefore independent of any recoil against the soft sector.

The SCET Lagrangian and resummation of the simple observable
From the above definition of the simple observable, it is now clear that the scaling of the soft and collinear modes in SCET can be obtained in terms of the a and b parameters as The SCET Lagrangian is uniquely defined by the modes given above The measurement function (4.17) allows for a simple factorization formula in the traditional sense. Using 28) and the fact that the simple observable is recoil insensitive [49], one finds where the expression in each sector is given by the obvious result Eq. (4.29) can be handled with standard RG techniques as it is usually done in SCET, as discussed in more detail in Appendix A. From now on we restrict ourselves to observables whose infrared dynamics is described by the degrees of freedom of this SCET Lagrangian. It is well known that more involved observables may require additional modes. An example is given by joint resummations [27,50,51], or by some cases of single resummations (such as negative-b angularities defined w.r.t. the thrust axis 4 ) where some extra care must be taken when constructing the effective Lagrangian. We believe that the formulation presented in this article can be adapted to such cases, although we will not discuss this further in the present article.

Definition of the transfer function
The second ingredient for the resummation is the transfer function F, which relates the resummation of the simple observable V max to that of the desired observable V . Most global rIRC safe observables have the property that the cumulative distribution is exponentially suppressed [e.g. Σ(v) ∼ exp(−α s ln 2 v)] in the limit v → 0. Such observables satisfy the basic relation to the cumulative distribution of their simple observables [47] (4.31) Using Eqs. (4.11) and (4.29) one can write a general factorization theorem for the transfer function 33) and the individual fully differential transfer functions are given by .

(4.34)
Note that the form of the general factorization theorem given in Eq. (4.32) does not depend on the observable, and in particular it does not depend on whether the observable factorizes or not.
A couple of considerations are in order. In the transfer functions for each separate sector, the virtual contributions to the numerator and denominator cancel in the ratio, such that only real emissions contribute. As explained in [47], and shown in the following sections, Eq. (4.34) can be computed numerically provided that the UV divergences in the SCET real radiation are handled by means of a regulator different from dimensional regularization, such as a cutoff or an analytic regulator. In this case the transfer function is finite in 4 dimensions and can be evaluated efficiently.

General expressions for the fully differential transfer function
In this section we derive the fully differential transfer function at NLL and NNLL accuracy. In our derivation, we only use the general factorization theorem given in Eq. (4.7) that is valid for any observable that depends only on the kinematics of the final state, but without assuming that the measurement function of the desired observable itself factorizes. This therefore extends the results of [47], in which similar results were derived for the thrust distribution, to generic observables. Note that in order to keep this section as simple to understand as possible, we write the results in terms of general squared amplitudes of the effective theory, and we defer giving the results using the explicit expressions for them obtained in SCET to Section 6.
In order to obtain an expression that is suitable for a numerical integration, one needs to address the cancellation of both IRC and UV divergences. We discuss them separately in the following.

Numerical treatment of the UV divergences
As already discussed in ref. [47], UV divergences appear both in the real and in the virtual corrections in a SCET calculation. In SCET I theories, where soft and collinear degrees of freedom have different virtualities, both of these UV divergences are regulated by dimensional regularization, which then in turn leads to the anomalous dimensions of the soft and jet functions. On the other hand, in SCET II theories, where soft and collinear degrees of freedom have the same virtuality, an additional regulator must be introduced to regularize the rapidity UV divergences that arise from the real radiation. Possible choices involve different analytic regulators [52][53][54][55][56].
Since in the numerical approach to resummation one needs to be able to perform real phase space integration numerically in 4 dimensions (after IRC divergences have been properly subtracted), the UV divergences in SCET I theories need to be regulated with a regulator other than dimensional regularisation. This guarantees that real and virtual UV divergences are regulated independently, and thus can be handled separately. The exponential regulator introduced in [56] regulates the UV rapidity divergence by shifting the separation of the soft Wilson lines that define the soft operator, and the collinear fields that define the collinear operator, by a small complex amount τ ∼ 1/Λ in both the x ± light-cone directions. In momentum space, this gives rise to an exponential suppression of the form where k 0 denotes the total energy that flows across the real phase-space cut (the total energy of the real radiation). Effectively, this amounts to regulating the UV divergences at a scale ∼ Λ. Such divergences will then appear as logarithms of Λ. Since the rapidity divergences cancel between the soft and collinear sectors, the dependence on ln Λ vanishes once these sectors are combined, but it does change the RG equations of the soft and collinear sector individually.
To isolate the logarithmic dependence on Λ which characterizes these extra UV divergences (and neglect any power suppressed dependence), one needs to expand the integrals obtained in SCET around the limit Λ → ∞. In a numerical implementation such an expansion is more difficult to do. One could simply choose a very large numerical value of Λ, such that the power suppressed effects would be numerically small, but this can result in an inefficient numerical algorithm. A better approach is to take formally the limit Λ → ∞ of the exponential suppression factor before performing the integration, such that it only gives rise to the logarithmic terms upon integration. For the soft integrals this amounts to making the asymptotic replacement 6 This converts the regulator into a cutoff on the larger light-cone component of each real momentum. One restricts the momentum k + in the integration region k + > k − and or k − for k − > k + . For the collinear integrals, the purely collinear contributions are already regulated in the UV by the phase space limits, and therefore the exponential regulator does not give rise to any logarithmic contributions and it can be entirely dropped. On the other hand, for the 0-bin subtraction in the collinear sector in the n direction, one proceeds as in the soft sector and uses the constraint Θ soft [n · k, Λ] for each emission. The above prescription amounts to defining the 0-bin subtraction integrals by replacing the collinear squared amplitude by In practical calculations this reproduces precisely the cutoff prescription discussed in [47], while retaining the nice properties of the exponential regulator in the analytic formulation. Since the dependence on Λ cancels between the soft and collinear sectors, one can choose any value of the UV scale for the practical numerical calculation. As discussed in ref. [47], a particularly convenient choice is to set it to its characteristic scale of the collinear sector, i.e. Λ = Q. This choice makes the collinear sector single logarithmic, while all double logarithms are included in the soft sector of the SCET Lagrangian. In the following we adopt this scheme in the derivation of the expressions for the transfer function.
Since the collinear sector is single logarithmic, the fully differential collinear transfer function at NLL reduces to a δ function, giving such that it only depends on the energy density introduced by the single recoiling quark of the underlying Born configuration. Conversely, the soft transfer function at NLL requires an infinite number of 1-particle soft clusters at tree-level. At NNLL, the collinear transfer function requires a single collinear emission at tree level, while for the soft transfer function one has to add a single insertion of the 1-loop contribution to the soft 1-particle cluster, and of the soft 2-particle cluster at tree-level. This generalizes in a systematic way to any logarithmic order. The transfer function at a general order N k LL requires a certain number of insertions of n-particle soft clusters with n = n particles computed to n loop order [see Eq. (4. 19)]. Given the NLL transfer functions, described above, the general expression for k > 1 can be written recursively as The insertions needed the computation of the corrections δF N k LL F (ω F , v) are summarized in Table 1.

Numerical treatment of the IRC divergences
We start by recalling that, in the definition of the transfer function (4.34), Σ max (Φ B ; v) is chosen to have the same LL structure as Σ(Φ B ; v). An important consequence of this fact soft sector (clusters) collinear sector n loops + n particles times n loops + n particles times is that the pure LL contribution completely cancels in the ratio and the transfer function differs from one by NLL corrections. A second important observation is that, as it will be shown shortly, the contribution from very unresolved emissions entirely cancels in the ratio when computing the transfer function for a given observable V < v. Therefore, F only receives a relevant contribution from emissions (or correlated clusters in the case of the soft sector) withṼ F ∼ v. A crucial consequence of this fact is that the phase space integration associated with each emission loses one logarithmic power, which in turn implies that the computation of the transfer function to N k LL accuracy requires the knowledge of the numerator and denominator of eq. (4.34) only to order N k−1 LL, that is . (5.7) In order to make the cancellation of the IRC divergences explicit in the transfer function, we introduce a parametrically small quantity δ 1 which needs to be independent of the observable v in the v → 0 limit. This allows us to recast Eq. (5.7) as .

(5.8)
In this way each of the two factors in the above equation is IRC finite, and the dependence on the regulator δ cancels in the product. It is important to notice that the Σ max F factor in eq. (5.8) is to be computed with the rapidity regulator defined in the previous section. Although this regularization procedure can be also adopted in the resummation of the simple observable in eq. (4.31), it is not always necessary (for instance in the case of SCET I observables) in the analytic calculation. However, it is crucial that the UV divergences of the real radiation are regulated as in Section 5.1 in the computation of the transfer function (5.8). In the following two subsections, we work out the results at NLL and NNLL in detail.

General expression for the transfer function at NLL
According to Eq. (5.7), the computation of the transfer function to NLL accuracy requires the knowledge of δσ F /δω F and Σ max F (v) to LL accuracy. As discussed above, the collinear sector gives rise to only single logarithmic behavior such that where we have used the fact that V max is zero for a state that only includes a single quark. From this one immediately obtains F NLL n,n (ω n,n ) = δ (ω n,n − ω[q n,n ]) . (5.10) Thus, to NLL accuracy, only the transfer function F s (ω s , v) is required and one finds To compute the soft transfer function to NLL accuracy, we start from the general expression (see eq. (4.8) and (4.9)) The virtual corrections V s in the soft sector of SCET are given by a scaleless integral, and therefore contain both UV and IR divergences, while with NLL accuracy (as discussed above) we can approximate This allows us to write where V LL s denotes the virtual corrections at LL accuracy. The explicit expression for this soft matrix element will be given in Sec. 6.1. The precise form of the virtual corrections is not relevant, since it will drop out of the transfer function to all orders. The running coupling in each 1-particle correlated cluster needs to be evaluated at the transverse momentum, which is the only Lorentz invariant and reparametrization invariant (RPI) scale available in the soft sector of SCET. Σ max,LL s in eq. (5.8) can be written as whereṼ s (k) is defined in Eq. (4.23). We have used that the N = 0 term in the sum is just equal to unity, and that to LL accuracy all final state particles are gluons, given the symmetry factor S(N ) = 1/N !.
In the combination of Eq. (5.14) and Eq. (5.15) the virtual correction drops out, and one obtains where we have used the shorthand notation v i ≡Ṽ s (k i ), and ∆ LL s (v, δv) has been defined by Eq. (5.8) and can be written as As a next step, we split the numerator of Eq. (5.16) into an unresolved subset (where each k satisfiesṼ s (k) < δv) and a resolved subset (where each k satisfiesṼ s (k) > δv) Now we use that this fully differential transfer function will be integrated against global rIRC safe observables. They have the property that unresolved emissions [withṼ s (k) < δv] can be neglected up to power corrections in δ in the measurement function of the observable [1,7] Θ(Ṽ s (k i ) < δv) V (q n , qn; k 1 , . . . , k N ) where h is a positive constant. Thus, upon integration against such an observable, and getting rid of power corrections in δ, the unresolved terms cancel between the numerator and denominator and one can write This gives the final equation In Sec. 6.1 we will use the SCET expression for the matrix elementM 2 s,0 (k) and a particularly suitable parametrization of the phase space integration to derive an explicit expression for this result that can be implemented in a straightforward numerical algorithm.

General expression for the transfer function at NNLL
The NNLL transfer function necessitates the calculation of δσ F /δω F and Σ max F (v) to NLL accuracy. This requires higher order corrections to the correlated clusters used for the soft transfer function, as well as the tree level single-emission amplitude for the collinear transfer function. We will discuss these two cases in turn.

The soft transfer function
As discussed in Section 2, in order to compute the cumulative distribution beyond NLL, it is sufficient to expand the NNLL correction perturbatively for fixed α s ln v. In the context of the soft transfer function, this amounts to correcting the NLL result by adding either a single one-loop correction to the 1-particle correlated cluster, or adding a single 2-particle correlated cluster at tree level. This allows one to write the numerator of the fully differential transfer function as where we have expanded the virtual contribution beyond LL accuracy perturbatively s + . . . just as for the real emission matrix elements. As for the LL case, the precise form of the virtuals is not important, as these will drop out of the final result.
In Eq. (5.22) there are IRC divergences in all of the three terms of the r.h.s., and one needs to work out their cancellation. We start by considering the denominator Σ max in the definition of the transfer function of eq. (5.8), and we will return to the general cancellation (which is required for the numerator of the transfer function) later. Using Eq. (4.30) one finds Σ max,NLL  In the square bracket of Eq. (5.23), the IRC divergences cancel between the various contributions. We start by canceling the divergence between the double real and real-virtual terms. For this, we use the well known result where the . . . represent terms beyond NNLL, and K is given by the ratio of the 2-loop to the 1-loop cusp anomalous dimension To obtain Eq. (5.24) it is crucial that the coupling in M 2 s,0 (k) is evaluated at µ = k ⊥ , which, as already discussed, is the only possible choice of scales in SCET for this quantity. While this result is typically derived in a setting outside of SCET [12,35,57], we have explicitly checked that it is obtained using the Feynman rules of the SCET soft Lagrangian.
Combining Eqs. (5.23) and (5.24), we obtain This gives the denominator of the general expression for the transfer function of Eq. (5.8), as well as the Sudakov factor at NLL accuracy Note that the first, second and third terms in the square bracket contain N , N + 1 and N + 2 momenta, respectively, and the action of the observable therefore is different in each of these terms, which are again separately IRC divergent.
The IRC divergence of V (1) s ultimately cancels in the ratio to Σ max in Eq. (5.22), so we ignore it for the time being. To handle the cancellation of the IRC divergence between the second and third term, we introduce a simple subtraction term in the square bracket, where the projection P + (k a , k b ) is defined through its action on the observable V (q n , qn; {k i }, k a , k b ) upon integration over the phase space of k a and k b as where {k i } represent a set of other emissions, and k ab is the massless momentum that is constructed from the transverse momentum and pseudo-rapidity of the vector k a + k b .
With this subtraction term, Eq. (5.28) can be recast as where we have used Eq. (5.24) to combine the integrated subtraction withM 2 s,1 (k a ). Combining Eqs. (5.26) and (5.31), the virtual correction drops out again, and one obtains after a few lines of algebra the following soft transfer function and we have taken the limit δ → 0 in the term that is now regulated through the P + projection. As a final step, we replace the ∆ NLL s with the expression given in Eq. (5.27). The final result can be conveniently expressed as and, as above, we have taken the limit δ → 0 in the term regulated by the P 0 projection. This is the final result for the fully differential soft transfer function at NNLL written in terms of a general expression for the correlated clusters. In Section 6.1 we will use the explicit SCET results together with a suitable phase space parametrization to simplify these results further.

The collinear transfer function
Since, as already discussed, collinear emissions only start contributing at NNLL order to the transfer function, it is sufficient to include the effect of one single emission at this order. One can therefore write where the subscript n denotes the direction of the collinear sector. The Σ max n is given by Using Eq. (5.8), and following the same steps as in the soft case, one finds where P 0 (k n ) was defined in Eq. (5.34). We stress that the collinear sector must undergo a zero-bin subtraction to avoid double counting with the soft sector. This will be directly performed in Section 6.2.3, where the explicit parametrisation of the collinear transfer function is given.

Explicit results for transfer function
In this section we formulate the somewhat abstract results obtained in the previous section in a concrete way that is suitable for a numerical evaluation. 7 We rely on the explicit parametrization of the phase space and SCET matrix elements reported in Appendices B and C. We start by making some general remarks. The considerations of the previous section have resulted in equations where we have neglected squared amplitudes which are of subleading logarithmic order. However, each of the transfer functions given in Eqs. The first source of subleading corrections is due to the presence of the coupling constant in each correlated cluster. One obviously needs to evaluate the running coupling at the desired resummation accuracy, i.e. only include β 0 at NLL, add β 1 at NNLL and so on. The second source is due to the rapidity boundaries, which in general are such that each emission develops a subleading logarithmic contribution upon phase space integration. The third source is arising from the fact that the transfer functions are integrated against the observable's measurement function. In general, this integral will again give rise to subleading logarithmic and power suppressed effects.
As we will show in the following, it is possible to address each of the above points in such a way that no subleading logarithms are contained in the final transfer function. While this sounds very different from what is done in the usual SCET resummation, it is in fact quite similar. An SCET resummation relies on fixed order computations of anomalous dimensions (extracted from the divergent structure of the result) as well as finite contributions. In these calculations, regulators and observable constraints are expanded around their limiting value. To a given order of resummation, one deliberately picks the anomalous dimensions and finite pieces at different orders from the fixed order calculation. In a numerical procedure as employed here, one then needs to find a way to reproduce these results and to ensure that no subleading logarithms are produced upon phase space integration.

Explicit expressions at NLL
The final expression of the transfer function in Eq. (5.11) is given in terms of Eq. (5.21), with the Sudakov defined in Eq. (5.17). One therefore needs to parametrise the phase space integral over the single particle soft correlated cluster. Using the expressions reported in Appendices B and C, we can write where C is the colour factor associated with leg . Note that in this expression it should be understood that the variables v i , χ i and φ i are always defined with respect to the leg relative to the emission k i , and we only suppress this dependence here to simplify the notation. To NLL accuracy, Eq. (6.1) can be simplified by performing the following Taylor expansion 2C a α s (k i,⊥ ) π and the running coupling at LL accuracy is simply given by One can easily see that the integration over the second term leads to a subleading logarithmic contribution, which can be dropped to the order we are working. It will be included in the next section when going to NNLL. Using this, one immediately obtains for the Sudakov factor in Eq. (5.17) that can be computed analytically yielding Combining these results into Eq. (5.11), one finds In the last line we have made explicit that the functional form of F NLL (Φ B ; v) is determined entirely in terms of P LL s, and the observable V .

Explicit expressions at NNLL
In this section we derive explicit expressions for the NNLL corrections to the transfer function. They can be classified in 3 broad terms. First, one needs to re-compute the NLL transfer function including the higher order terms to the various expansions that were made in Eq. (6.2). One therefore writes As previously discussed, to obtain NNLL accuracy it is sufficient to describe a single emission with the NLL splitting function, while using LL splitting functions for all other emissions. Therefore, one could expand P NLL s, about P LL s, in Eq. (6.9), and retain only the linear term at this accuracy, i.e.
where the dots represent subleading logarithmic corrections. We leave the details of this operation to Appendix D.3, and only present the expression for P NLL s, below. We stress that this approximation is by no means mandatory, and one can directly evaluate the NLL transfer function with the NLL splitting function. 8 As a second NNLL term, we need the higher order corrections to the soft n-particle correlated clusters, given by in Eq. (5.33). As discussed there, this involves either one or two soft emissions on top of the ones required for the NLL transfer function. Finally, one needs the collinear transfer function, given in Eq. (5.37). As before, we will denote these two corrections by with F ∈ {s, n,n}. Here we explicitly state the dependence on the observable, as this will become important later. Since one only needs to include one correction at a time, the NNLL corrections to the soft transfer function are always combined with the NLL result for the collinear transfer function (which is just a delta function), while the NNLL correction to the collinear transfer function is combined with the NLL soft transfer function derived the previous section. We now discuss these three pieces in turn.

s,
In this subsection we just provide the expression for the NLL splitting function P NLL s, (v, χ a ), described above. This is simply defined by keeping the NLL term in the expansion Eq. (6.2), that is To use this result in Eq. (6.8), we need the integral of Eq. (6.12) over χ a and φ a . One finds . (6.14)

s (v): The soft contribution
We recall once more that in order to obtain a NNLL result for the cumulative distribution, the NNLL correction to the soft transfer function is to be combined with the NLL collinear transfer function, which at this order is just a delta function. Starting from Eq. (5.33) and keeping only the term involving a single extra emission we find, using eq. (6.2) 2C a For the term involving two extra emissions, we use the phase space parametrization given in Appendix B and the soft 2-particle correlated cluster given in Appendix C. We can write whereM 2 s,0 (k a , k b ) is given in Eq. (C.10), and we defined κ = k ⊥,b /k ⊥,a . As done for the single emission contribution above, we now Taylor expand the arguments of the couplings and the rapidity bounds in order to retain only NNLL terms, hence obtaining s,0 (k a , k b ) + . . . , (6.18) where . . . refers to higher logarithmic terms. Using eq. (5.33) leads to where the second line is as in Eq. (6.17), and as before k ab is the massless momentum that is constructed from the transverse momentum and pseudo-rapidity of the vector k a + k b . The total NNLL result originating from the higher order soft correlated clusters is then The NNLL collinear correction is obtained by combining the fully differential NNLL collinear transfer function with the NLL soft transfer function, obtained in Section 6.1. Using the collinear phase space parametrization and collinear matrix elements given in Appendices B and C, respectively, we obtain As for the soft transfer function, we have expanded the integrand retaining only the most singular terms. Thus, we find To obtain the final collinear contribution to the NNLL transfer function, we first need to combine this with the fully differential soft transfer function at NLL obtained in Section 6.1. This gives, for each collinear sector, × Θ(V (q n , qn; k 1 , . . . , k N , k a ) < v) − Θ(V (q n , qn; k 1 , . . . , k N ) < v) Θ(V max (k a ) < v) .

(6.24)
Finally, we need to perform the 0-bin subtraction to avoid the double counting with the soft sector. This amounts to subtracting from eq. (6.24) its purely soft approximation, defined as where (6.26) and V max (k a ) in eq. (6.25) is to be considered in the soft limit, i.e. V max (k a ) =Ṽ s (k a ), defined in eq. (4.23). The NNLL collinear correction is simply given by the difference of eqs. (6.24) and (6.25), namely Note that the same definition of the observable is used in both (6.24) and (6.25). However, as we will discuss in the next section, the whole observable in the zero-bin contribution has to be expanded about its soft limit.

Dealing with the observable
The expressions for the transfer function derived in the previous section use the generic rIRC safe observable V (q n , qn; k 1 , . . . , k N ). In this section we discuss how to evaluate the observable on a set of emissions k i corresponding to a given sector of the SCET Lagrangian.
One possible choice, convenient for a flexible MC implementation, would be to use the full definition of the observable, which often is know in the form of an algorithm. For example, a jet observable requires a sequential clustering algorithm, and even a simple observable such as thrust involves a minimization procedure to find the thrust axis. An issue with using the full definition is that it often relies on the fact that the final state satisfies momentum conservation, which is explicitly invalidated when constructing the kinematics of the momenta of the fully differential transfer function. This is because phase space limits were explicitly taken in the soft and collinear limits. Such a problem can be overcome by reconstructing the final state kinematics after the generation of the radiation, following a procedure similar in spirit to what is currently used in some parton shower generators [58,59]. However, a drawback of this approach is that the corresponding resummed prediction will include uncontrolled subleading-logarithmic and subleading-power corrections that are ultimately eliminated with the matching to fixed order at the relevant perturbative accuracy.
A different recipe, commonly adopted in resummations, is to expand the observable around the relevant kinematical limits. As in a common SCET factorisation theorem, this implies that the soft and each collinear sector is evaluated using the observable expanded in the corresponding soft or collinear limits. Performing the explicit expansion of the observables has the advantage of dropping undesired subleading logarithmic corrections, as well as all power corrections, which are retained when using the full observable dependence. This is important in a purely resummed calculation. Furthermore, in the context of the method presented in this article, it allows for the numerical extraction of specific ingredients to other resummation approaches. In the following we will briefly discuss the implications of the latter choice in the calculation of the transfer functions defined in the previous section.

The soft sector
The evaluation of the transfer functions involving the soft sector, namely eqs. (6.8), (6.16), and (6.19) can be carried out by taking the soft limit of the observable V , namely V (q n , qn; k 1 , . . . , k N ) → V s (q n , qn; k 1 , . . . , k N ) . (7.1) While this choice automatically eliminates any source of subleading-power corrections in the resummed prediction, for generic observables it may still lead to the inclusion of subleading logarithmic effects in the result. As discussed in Section 2, this is a common feature in any resummation approach.
Nevertheless, it is interesting to work out a simple recipe to neglect sources of subleading logarithmic corrections due to the approximation of the observable. This can be easily done by noticing that, at NLL, the relevant kinematics is described by an ensemble of soft radiation (E i Q) that is strongly separated in angle and collinear to either of the legs in the hard process. Since η i ∼ χ i ln 1/v, for v → 0 the collinear limit implies that all χ i are of order 1, while the strong separation ordering implies that all χ i 's are different from each other. In other words, the soft-collinear approximation for the observable V sc is obtained by expanding the full observable about the limit V sc (q n , qn; k 1 , . . . , k N ) : The corresponding approximation can be adopted in the calculation of the NLL transfer function F NLL (6.8) and (6.10), that can be recast as 3) The first line provides the correct NLL result, while the second line contributes at NNLL and beyond. The same approximation can be used in evaluating the correction δF NLL defined in eq. (6.10), as well as the single-soft correction δF NNLL s,1 at NNLL (6.16). The second line of eq. (7.3) can be once again expanded systematically beyond NLL. At NNLL, the soft sector is described once again by an ensemble of soft and collinear radiation strongly separated in angle, and by a single soft emission that can either be close in angle to any of the above, or be emitted with a large angle (away from the collinear limit). The two limits can be promptly translated into two approximations for the observable V , namely: • Two of the existing emissions k a and k b are soft and collinear, but have similar rapidities V sr (q n , qn; k 1 , . . . , k N ; k a , k b ) : This configuration is irrelevant for event shapes, but it contributes for observables that explicitly depend on the relative rapidity between emissions. A typical case is when jet algorithms are involved in the definition of the observable [9].
• One of the existing emissions k a is soft and strongly separated in angle from all other emissions, but with large angle V wa (q n , qn; k 1 , . . . , k N ; a) : Accordingly, at NNLL, eq. (7.3) can be expanded as Finally, the double-soft NNLL correction (6.19) involves two correlated emissions that are close in rapidity, and therefore it should be naturally evaluated using the approximation V sr of the observable described above.
We stress that the approximations described in this section are useful to neglect subleading logarithmic corrections but they are, strictly speaking, unnecessary. The kinematic expansion described here can be systematically extended to higher logarithmic orders by progressively relaxing the strong angular separation and the collinearity constraint for an increasing number of emissions. The downside of this approach is that, while systematically extendable to higher orders, it requires working out the necessary limits for the observable at each new logarithmic order.
The restrictions on the observable that we just discussed can be simply understood in SCET terms as follows. In traditional SCET resummation the NLL result is entirely given in terms of the anomalous dimensions. Observable dependence in the anomalous dimension only arises at 1-loop, and therefore only depends on a single emission. This picture is consistent with taking the strong angular separation limit. Taking the collinear limit of the observable guarantees that no finite contributions are generated, which corresponds to setting the initial conditions to one in a NLL resummation. Thus, if one is concerned with generating subleading logarithmic corrections in addition to the desired NLL terms, the above mentioned limit of the observable should be taken.

The collinear sector and the zero-bin subtraction
The collinear transfer function at NNLL (6.24) involves a set of soft emissions and one collinear emission. Therefore, the observable must be evaluated in the corresponding kinematic limit.
As for the soft transfer function, it is possible to make minor approximations in order to neglect subleading logarithmic corrections from the final result. Following the discussion of the previous subsection, in the computation of the observable the soft emissions can be simply approximated by their collinear limit and by assuming that they are strongly separated in angle, hence neglecting corrections of N 3 LL order and higher.
The collinear transfer function (6.24) at NNLL is then approximated by where V hc is defined on a set of soft-collinear emission which are strongly separated in angle, and a single (hard) collinear emission k a . Finally, similar considerations can be made for the corresponding zero-bin subtraction (6.25), that at NNLL can be evaluated using the V sc approximation of the observable, i.e.
where the replacement also implies that one should take the soft limit of V max in eq. (6.25).
The final NNLL collinear correction is given by the difference of the above two equations.

Conclusions
In this work we outlined a general formalism to resum rIRC safe observables in soft collinear effective theory. This formalism combines two quite different approaches to resummation with one another and results in a very general and powerful tool that allows one to resum in a relatively straightforward way most observables of interest to higher logarithmic order. The first approach used is the standard SCET technique, which commonly relies on the derivation of a factorization theorem for a given observable, and then uses RG equations to resum logarithms at all orders in each of the components of the factorization theorem. The second approach is the branching formalism, which uses the factorization properties of squared amplitudes in QCD to model the soft and collinear radiation at all-orders, often simulated numerically using a MC algorithm. By combining these two approaches, our new formalism exploits the strengths of both methods, namely the generic numerical solution provided by the branching formalism with the ability to easily and systematically go to higher orders of the SCET approach.
The main idea behind this method is inspired by numerical resummation techniques [1,7,12] and consists of deriving the resummation for a simplified version of the considered observable, which shares the same LL structure with the full one. For a suitably chosen simplified observable, this resummation is much simpler than for the full observable, and can be always handled analytically. Moreover, the same simplified observable can be adopted in many cases. One can then relate the resummation of the full observable to that of the simplified one multiplicatively via a transfer function, which can be computed numerically. While this approach was already introduced in [47] for the thrust distribution, in this paper we discussed how it can be generalized to any rIRC safe observable. The main insight that allowed this generalization was the introduction of a fully differential energy distribution, which allowed us to formulate the transfer function for any observable.
We described all details of this method, and provided all equations necessary for an explicit implementation of the calculation of the transfer function up to NNLL accuracy. By collecting the results obtained in the article, we arrive to the following master formula for the NNLL transfer function for a generic rIRC safe observable V As discussed in Sections 6.2 and 7, it can be sometimes convenient to neglect any source of subleading logarithmic corrections that may originate from evaluating eq. (8.1) directly. While this step is, strictly speaking, not necessary, we find it useful to provide a simple recipe to neglect all sources of N 3 LL contamination from eq. (8.1). To this end, we can further decompose the above result as follows While we did not give any explicit results to N 3 LL and beyond, the derivation of our results using the language of effective field theory was presented in a way that should make the systematic extension to higher logarithmic accuracies obvious. In future work, we will present applications of this approach to a range of experimentally relevant observables and provide generic computer code that can be used by non-experts to obtain the resummed prediction for a given rIRC safe observable.

A. Resummation for Σ max (v)
In this appendix we discuss in some more detail the resummation for the simple observable V max defined in Section 4.2.
Given the simple form of the measurement function (4.17), one is able to find an analytic result for this observable at all logarithmic orders. From the definition, it is easy to see that the characteristic scale of the soft sector is µ s = Q v 1/a , while each collinear sector has µ J = Q v 1/(a+b ) . By observing that the observable definition does not mix soft and collinear modes, and in particular that there is no kinematic cross talk between the different modes, we obtain the simple multiplicative factorisation theorem The jet and soft functions are defined as where Y n denotes a soft Wilson line along the n direction, while χ n ≡ W † n ξ n , with ξ n being a collinear fermion field after the BPS [44] field redefinition, and W n being a collinear Wilson line. By comparing eq. (A.1) and eq. (4.29) we obtain the obvious relations Resummation proceeds as usual by renormalizing each of the above building blocks up to a given perturbative order. Given the multiplicative structure of the factorisation theorem, after renormalisation each of them satisfies an RGE of the type with F = {s, n,n}. Its solution yields where the anomalous dimensions and the boundary conditions Σ max F (µ F ) have the following perturbative expansions The logarithmic accuracy is defined in terms of the perturbative order of the anomalous dimensions and boundary conditions, as summarized in Table 2. For example, to achieve LL accuracy, one only needs Γ (1) cusp , while for NLL accuracy Γ An important property of the simple observable defined in Section 4.2 is that the same definition can be used for the resummation of any rIRC safe observable V , with the only observable dependence being encoded in the a and b parameters. A consequence of this property is that the anomalous dimensions and initial conditions (A.6) depend on the observable under consideration only via the same parameters, which allows for a generic calculation of Σ max at once.

B. Kinematics and phase space parametrization
In this appendix we report the parametrizations used for the phase space in each sector of the SCET Lagrangian with the UV regulated by a cutoff at Q, as described in [47]. The only phase space we need to compute for the NNLL transfer function is a single soft or collinear momentum, and a correlated pair of soft momenta.

Single soft emission
The most natural phase space representation in SCET is in terms of the light-cone components and the transverse momentum, giving where the rapidity is defined in terms of the light-cone components as In the above equation the rapidity limit comes from eq. (5.2) where we have explicitly used Λ = Q as discussed in Section 5.1. The unit 4-vector n defining the light cone decomposition of the momentum k is in principle arbitrary. In the following we will separate the phase space into two regions associated with either of the legs according to the definition of the simple observable given in Section (4.2). For two legs this corresponds to separating into two hemispheres, giving Now one can change the phase space variables from k ⊥, to v ≡Ṽ s (k), withṼ s (k) given in Eq. (4.23). We introduce the rapidity fraction (with support from 0 to 1) as This allows us to write where η and φ are measured with respect to leg and Single collinear emission Starting again from the standard SCET parametrization in terms of the light-cone components and the transverse momentum, we parametrize the collinear momentum as where now unit 4-vector n is defined by the direction of the collinear sector It is useful to express the simple observable in the collinear limit (4.25) as a function of k ⊥ and z n , where k ⊥ is defined with respect to the collinear direction n. For a single collinear emission we obtainṼ By changing the phase space variables from k ⊥ to the value of the observable v n ≡Ṽ n (k), one can write where k n,⊥ ≡ k n,⊥ [v n , z n ] = Qv 1/(a+bn) n z bn/(a+bn) n (1 − z n ) bn/(a+bn) . (B.12) Note that the limit z n < 1 comes about automatically from the constraintn·k n < Q in the collinear sector of SCET. Conversely, when performing the zero-bin subtraction, the above limit comes from the extra constraint due to the UV regulator in eq. (5.4), using Λ = Q as discussed there. This implies that the same phase space can be used both for the collinear contribution and its corresponding zero-bin subtraction.

correlated soft emissions
The phase space of 2 correlated soft emission has 6 independent variables. They can be parametrized in many different ways. The simplest parametrization is to use the usual phase space given in eq. (B.6) for one (say k a ) of the two soft emissions, and parametrize the remaining phase space in terms of the relative rapidity between the two emissions ∆η, the ratio of the two transverse momenta 13) and the azimuthal angle (φ b ) of the second emission k b . We find At NNLL accuracy, we are interested in integrating over phase space regions where the two emissions k a and k b are in the same hemisphere. All configurations in which the two emissions propagate into opposite hemisphere start contributing at N 3 LL, and hence we purposely omit them in our parametrisation.
We therefore approximate Eq. (B.14) by The argument of the running coupling constant again is the natural scale of the collinear sector. At NNLL, the transverse momentum k ⊥ (B.12) in the collinear squared amplitude can be approximated with the scaling Eq. (4.26), i.e. Q v 1/(a+bn) , with the difference between the two scales being of N 3 LL order.

correlated soft emissions
The tree level expression for the 2-body amplitude M s (Φ B ; k a , k b ) in the soft sector of SCET was computed in [60][61][62]. We write where the three different color structures C 2 F , C F C A and C F n f T F are given by (C.8) When calculating the 2-particle soft correlated cluster the Abelian contribution cancels and one obtainsM where the factor of 1/2 takes into account the symmetry factor for having identical gluon fields.

D. Numerical algorithms
In this section we discuss how to generate the transfer function in a simple numerical algorithm. The algorithms given below are meant as a guideline for a numerical implementation of this formalism. Furthermore, as discussed in the main text, several choices can be made at NNLL, and below we simply limit ourselves to the main formulae discussed in the text. We stress that more efficient implementations can be adopted, like the ones proposed in refs. [1,7], but their discussion is beyond the scope of this article.
In the algorithms defined in the remainder of this section we adopt the prescription outlined in Section 7, according to which the observable in each sector of the SCET Lagrangian is expanded about its NLL approximation (7.2), denoted by V sc . We stress that this choice does not affect the generality of the algorithms given below.

D.1 Results at NLL
The starting point of the numerical algorithm is Eq. (6.8). Using the simple expression of the Sudakov factor given in Eq. (6.6), one easily finds Furthermore, we use .

(D.2)
This allows us to write Eq. (6.8) as where we have eliminated the 1/N ! symmetry factor by ordering the emissions according to v i , and we have used that v 0 ≡ v.
Thus, the emissions with momenta k i (v i , χ i , φ i ) have a distribution that is very similar to that of a parton shower algorithm. They involve an evolution variable (v i ) that is monotonically decreasing with each emission, together with a splitting function P LL s, (v; χ i ) and associated Sudakov factor ∆ LL s . The emissions can therefore be generated by the following algorithm: Multiply the event weight w by P LL s, (v, χ i )/ P LL s, (v) × n legs ; Determine k i = k(v i , χ i , φ i ) and add to the list of emissions; end Return the list of momenta {k i } and associate weight w; From the emissions generated with Algorithm 1 one then immediately computes the NLL transfer function by taking the average weight of all emissions with the constraint on the observable V sc (q n , qn; {k i }) < v. This is results in the algorithm if V sc (q n , qn; {k i }) < v then Increase W by w; Increase W Sq by w 2 ; end end Compute F NLL ± ∆F NLL from the average value of W and its standard deviation;

D.2 Results at NNLL
As discussed in Section 6.2, there are 3 broad classes of terms contributing to the transfer function at NNLL. The first is F NLL [P NLL As discussed in Section 6.2, keeping the NLL terms in the expansions made in Eq. (6.2) results in the NLL splitting function given in Eq. (6.12). While the LL splitting function given in Eq. (6.3) only depended on the values of v and χ a , the NLL splitting function depends on the value of v as well as all three emission variables v a , χ a and φ a . This means that after performing the analytical integration over χ i and φ i , the function still depends on v a . The Sudakov factor is then given by In terms of these expressions, the algorithm can be written as The transfer function is then computed according to Algorithm 2, where the emissions are generated as just described. The NNLL soft correction derived in Section 6.2.2 is further divided into two contributions (see Eq. (6.20)). The correction involving a single emission reads δF NNLL s,1 × Θ(V sc (q n , qn; k 1 , . . . , k N , k a ) < v) − Θ(V sc (q n , qn; k 1 , . . . , k N ) < v) Θ(Ṽ s (k a ) < v) , Choose v a ∈ [0, v], χ a ∈ [0, 1], φ a ∈ [0, 2π] randomly from a flat distribution; Construct the momentum k a from the three emission variables; Multiply the event weight w by P LL s, (v; χ a ) K α LL s (κ a,⊥ ) 2π × n legs ; Compute Θ ≡ Θ(V sc (q n , qn; {k i }, k a ) < v) − Θ(V sc (q n , qn; {k i }) < v) Θ(Ṽ s (k a ) < v); if Θ = 0 then Add w Θ to W ; Add w 2 to W Sq ; // Θ = ±1 end end Compute δF NLL ± ∆δF NLL from average value of W and its standard deviation; The second term, involving two extra emissions, reads δF NNLL S,2 × Θ(V sr (q n , qn; k 1 , . . . , k N , k a , k b ) < v) − Θ(V sc (q n , qn; k 1 , . . . , k N , k ab ) < v) , where k ab is the massless momentum that is constructed from from the transverse momentum and rapidity of the vector k a + k b , and the squared matrix element is reported in Eq. (6.17). As for the previous one, this correction can be computed with the following algorithm: Algorithm 5: Computing the NNLL correction δF NNLL S,2 [V ](Φ B ; v) Set weight W = 1, W Sq = 1; for i = 1 . . . N do Generate a set of soft-collinear emissions {k i } with weight w using Algorithm 1; Choose the leg randomly from a flat distribution; Choose v a ∈ [0, v], χ a ∈ [0, 1], φ a , φ b ∈ [0, 2π] randomly from a flat distribution; Choose κ ∈ [0, 1] uniformly and ∆η ∈ (−∞, ∞); Construct the momenta k a , k b and k ab from the above emission variables; Build the event weight w according to eq. (6.18) ; Compute Θ ≡ Θ(V sr (q n , qn; {k i }, k a , k b ) < v) − Θ(V sc (q n , qn; {k i }, k ab ) < v); if Θ = 0 then Add w Θ to W ; Add w 2 to W Sq ; // Θ = ±1 end end Compute δF NLL ± ∆δF NLL from average value of W and its standard deviation; × Θ(V hc (q n , qn; k 1 , . . . , k N , k a ) < v) − Θ(V sc (q n , qn; k 1 , . . . , k N ) < v) Θ(Ṽ (k a ) < v) . Once again, the algorithm is a simple adaptation of Algorithm 1: an analogous expansion for the real contribution, and after some manipulations we obtain Thus, the correction arising from the higher terms in the expansion of the splitting function (6.9) is given by × Θ(V sc (q n , qn; k 1 , . . . , k N , k a ) < v) − Θ(Ṽ s (k a ) < v)Θ(V sc (q n , qn; k 1 , . . . , k N ) < v) .
This can be implemented using the same variant of Algorithm 1 used in the computation of δF NNLL (v; v a , χ a , φ a ) × n legs ; Compute Θ ≡ Θ(V sc (q n , qn; {k i } , k a ) < v) − Θ(Ṽ s (k a ) < v)Θ(V sc (q n , qn; {k i }) < v); if Θ = 0 then Add w Θ to W ; Add w 2 to W Sq ; // Θ = ±1 end end Compute δF NLL ± ∆δF NLL from average value of W and its standard deviation;