The Sudakov radiator for jet observables and the soft physical coupling

We present a procedure to calculate the Sudakov radiator for a generic recursive infrared and collinear (rIRC) safe observable in two-scale problems. We give closed formulae for the radiator at next-to-next-to-leading-logarithmic (NNLL) accuracy, which completes the general NNLL resummation for this class of observables in the {\tt ARES} method for processes with two emitters at the Born level. As a byproduct, we define a physical coupling in the soft limit, and we provide an explicit expression for its relation to the $\overline{\rm MS}$ coupling up to ${\cal O}(\alpha_s^3)$. This physical coupling constitutes one of the ingredients for a NNLL accurate parton shower algorithm. As an application we obtain analytic NNLL results, of which several are new, for all angularities $\tau_x$ defined with respect to both the thrust axis and the winner-take-all axis, and for the moments of energy-energy correlation $FC_x$ in $e^+e^-$ annihilation. For the latter observables we find that, for some values of $x$, an accurate prediction of the peak of the differential distribution requires a simultaneous resummation of the logarithmic terms originating from the two-jet limit and at the Sudakov shoulder.


Introduction
Distributions in event shapes and jet resolution parameters, collectively jet observables, are among the most studied QCD observables. Since they are continuous measures of the hadronic energymomentum flow in jet events at colliders, they constitute a powerful probe of the dynamics of strong interactions, from high scales where fixed-order perturbative calculations can be applied, down to low scales where the yet unexplained phenomenon of hadronisation plays a decisive role.
Jet observables play a major role in measurements of the QCD coupling α s , and in testing non-perturbative hadronisation models (see e.g. Ref. [1] and references therein). The study of jet observables also led to important advances in the understanding of all-order properties of QCD radiation, which lead to the discovery of the so-called non-global logarithms [2][3][4]. Distributions in jet observables can be computed at fixed order in QCD perturbation theory. Such calculations have reached next-to-next-leading order (NNLO) accuracy for a number of relevant QCD processes. In particular, for e + e − annihilation, NNLO corrections to three-jet production have been computed in refs. [5][6][7][8][9].
While fixed-order calculations provide a reliable tool to describe jet observables in the region where their values are large, the bulk of data lies in a region where multiple soft-collinear emissions give rise to large logarithms of the jet observable at all orders in perturbation theory. To be precise, given a generic jet observable, let us consider its cumulative distribution Σ(v), the fraction of events such that the observable's value is less than v. This quantity exhibits logarithmic terms as large as α n s L 2n , where L = − ln v and n is the order in QCD perturbation theory. Resumming those large logarithms means reorganising ln Σ in such a way that it can be written as α s g 1 (α s L) + g 2 (α s L) + α s g 3 (α s L) + . . . , where g 1 (α s L) resums the so-called leading logarithmic (LL) contributions, α n s L n+1 , g 2 (α s L) the NLL ones, α n s L n , g 3 (α s L) the NNLL ones, α n s L n−1 , and so on.
Next-to-leading logarithmic (NLL) resummations, that include all terms O(α n s L n ) in the logarithm of cumulative distributions, have been available for many years for specific observables [10][11][12][13][14][15][16]. Nowadays, NLL resummation for jet observables that have the properties of recursive infrared and collinear (rIRC) safety and continuous globalness [18,19] is a solved problem. The general solution is based on a semi-numerical approach developed for e + e − event-shapes and jet rates in Ref. [20], and later extended to any suitable jet observable in any QCD hard process [18,19]. The method is implemented in the computer program CAESAR [18], that also verifies whether a given observable is rIRC safe and continuously global. This led to a first systematic study of event shapes in hadronic dijet production at NLL accuracy matched to next-to-leading order (NLO) results at hadron colliders [21,22].
NLL predictions have a sizeable theoretical uncertainty. Given the precision of current experiments, theoretical accuracy for resummations should aim at NNLL, and in some cases beyond. Most NNLL resummations are observable specific, and rely on the properties of the observable to achieve resummation through factorisation theorems that separate different kinematical configurations (e.g. hard, soft, collinear), and appropriate renormalisation group equations based on the fact that physical distributions do not depend on the unphysical scales that need to be introduced to achieve such separation. Such approaches made it possible to obtain full next-to-next-to-leading logarithmic (NNLL) predictions for a number of global e + e − event shapes such as thrust 1 − T [23][24][25], heavy jet mass ρ H [26], jet broadenings B T , B W [27], C-parameter [28], energy-energy correlation (EEC) [29][30][31], heavy hemisphere groomed mass [32], and angularities [33]. Among the above examples, for 1 − T , ρ H , C-parameter, and EEC, all N 3 LL corrections are also known, except for the four-loop cusp anomalous dimension, that has been computed numerically more recently [34]. Jet observables have been resummed at NNLL accuracy also in deep inelastic scattering [35][36][37]. For hadronic collisions, full NNLL resummations are available for processes where a colour singlet is produced at Born level, specifically for a boson's transverse momentum [38,39] and φ * [40], the beam thrust [41,42], transverse thrust [43], and the leading jet's transverse momentum [44][45][46][47], and for heavy quark pair's transverse momentum [48,49]. For an arbitrary number of legs, a NNLL accurate resummation is available for the N -jettiness variable [50]. Very recently, resummations for the boson's transverse momenta and φ * have been pushed to N 3 LL accuracy [51][52][53].
Despite these remarkable results, most jet observables are beyond the scope of factorisation theorems. This is especially true for those observables which cannot be expressed in terms of simple analytic functions of momenta, e.g. event shapes like the thrust major, or the two-jet rate in the Durham algorithm. For rIRC safe observables, it is possible to achieve NNLL accuracy by means of the semi-numerical method ARES (Automated Resummer for Event Shapes), developed for e + e − event shapes in Ref. [54], and later extended to jet rates [55]. These publications focused on NNLL corrections induced by resolved real radiation. There, the cancellation of infrared singularities between unresolved real and virtual corrections to the Born process was parametrised in an observable dependent Sudakov form factor called "radiator", that was extracted from existing calculations. This made it possible to study only observables that scale like powers of the transverse momentum or the invariant mass of the jet, which constitute a vast set of the phenomenologically relevant observables. This led to the first resummations for complicated observables such as the thrust major and the two-jet rate with various jet algorithms in e + e − [54,55].
In this paper we complete the last analytic ingredient necessary to have a fully general formula for the resummation of rIRC safe jet observables at NNLL in processes with two hard legs at the Born level. This makes it possible to handle all known observables of this type in a single framework, and hence paves the way to systematic phenomenological applications. We formulate the Sudakov radiator at all orders for a generic rIRC jet observable in QCD, and we explicitly compute it at NNLL accuracy. Our calculation of the Sudakov radiator must be then supplemented with the finite contributions coming from real radiation that are computed as in refs. [54,55].
The paper is organised as follows. In section 2 we give the formulation of a general procedure for the resummation of jet observables with the ARES formalism. There we define the main object of our paper, the Sudakov radiator, which we compute at NNLL accuracy in section 3. The computation of the radiator leads to a definition of the physical coupling for soft radiation at higher orders. This generalises the scheme of Ref. [61], and constitutes an ingredient for future NNLL accurate parton shower algorithms. In that same section, we present a new formulation of the NNLL correction δF correl introduced in Ref. [54], that we redefined in order to make sure that the Sudakov radiator can be computed analytically for an arbitrary rIRC safe jet observable. In section 4, we apply our method to angularities and moments of energy-energy correlation in e + e − annihilation, and outline briefly the main features of their phenomenology at present and future colliders. Section 5 contains our conclusions.

Resummation in the ARES formalism
In this section we summarise the structure of the NNLL resummed cross section for a generic rIRC observable in e + e − annihilation. We first set up the relevant notation and kinematics, and then we move on to derive the resummed cross section up to NNLL accuracy.

Kinematics and notation
Let V ({p}, k 1 , . . . , k n ) be a generic continuously global, rIRC safe final-state observable, a function of all final-state momenta. 1 Here {p} denotes {p 1 ,p 2 }, which are the quark-antiquark pair initiating the process after all radiation has been emitted, and k 1 , . . . , k n are the momenta of the additional radiation. At Born level V ({p}) = 0. In order to parametrise the radiation momenta k i , we introduce the following Sudakov decomposition where p µ 1 , p µ 2 are two light-like reference vectors, and κ µ i is a space-like four-vector and its magnitude is denoted by k ti ≡ −κ 2 i . The on-shell condition k 2 i = 0 implies that z The choice of the reference momenta p 1 and p 2 is arbitrary, and determines the mapping between the Sudakov variables {z (1) i , z (2) i , κ i } and the actual final-state momenta {p 1 ,p 2 , k 1 , . . . , k n }. A particular choice is therefore motivated by computational convenience. In our case, we observe that the product of the squared amplitude and phase space for an emission k collinear to eitherp 1 orp 2 is proportional to wherek We choose the reference vectors p 1 and p 2 such that The actual inputs of final-state observables are hadron momenta. However, it is well known that using parton momenta gives distributions that differ from the measured ones by corrections suppressed by powers of the typical hard scale of the process, in this case the centre-of-mass energy of the e + e − collision.
up to corrections that vanish as a power of k t in the limit k t → 0. One possible solution, adopted in Ref. [18], is to choose the two reference vectors p 1 and p 2 in Eq. (2.1) to be the momenta of the emitter of parton k i and the corresponding spectator. The precise definition of emitter and spectator requires specifying an ordering of insertion of the emissions in the event. A natural way of reconstructing the kinematics is to insert the emissions as follows: • For each leg , insert the emissions into the initial event starting from the one at smaller η ( ) . The reference momentum p in Eq. (2.1) represents the emitter and takes the transverse recoil. The longitudinal recoil is shared between the emitter and a spectator momentum (i.e. the remaining reference vector in Eq. (2.1)). However, for the purpose of the analytical method presented in this paper, one can safely neglect the longitudinal recoil of the spectator (which is proportional to k 2 ti ) that would otherwise give rise to regular terms in the cross section. For instance, for an emission k emitted off leg 1, and parametrised by Eq. (2.1), we have The resulting momenta p 1 and p 2 , after the emission, will be massless up to O(k 2 ti ) corrections.
• Update the reference momenta, and proceed with the insertion of emissions at progressively larger η This procedure guarantees the validity of Eq. (2.3), and implies that the reference vectors p 1 and p 2 are different for each emission [54]. Before proceeding, we stress that the ordering chosen for the insertion of the emissions in the event has nothing to with the way the triplets {z (1) i , z (2) i , κ i } are produced. For instance, when computing the resummation via a Monte Carlo algorithm, it is convenient to generate this triplet according to an ordering in the observable's value (see, e.g. [18,54]). This ordering is however unrelated to the ordering with which the emissions are inserted, which follows the angular ordering arguments outlined above. This is essential to specify a correct recoil scheme that preserves the simple factorised form of QCD matrix elements.

General structure of NNLL resummation
Any rIRC safe observable V ({p}, k 1 , . . . , k n ), can be parametrised in the following way for a single soft and collinear emission k collinear to each leg Here k ( ) t , η ( ) , φ ( ) are the transverse momentum, rapidity and azimuth of k with respect to the emitter p as defined in the previous section, whereas a, b , d are constants. 2 The scale Q represents a typical hard scale of the process, in our case the centre-of-mass energy of the e + e − collision.
Our aim is to resum large logarithms in the cumulative distribution Σ(v), the fraction of events for which V ({p}, k 1 , . . . , k n ) < v, in the region v 1. This is given by where H(Q) includes all virtual corrections to the Born process (normalised to the total cross section σ) and M 2 (k 1 , ..., k n ) is the amplitude squared for n real emissions. The Lorentz-invariant phase-space in d = 4 − 2 dimensions is denoted by [dk] and defined as Given the Sudakov decomposition of any four-momentum k of (squared) invariant mass k 2 = m 2 as in Eq. (2.1), the measure d d k can be expressed as The variable y is the rapidity of k (with respect to some reference light-like directions p 1 and p 2 ), and for real emissions, i.e. k 2 = 0, it is bounded by 3 It is immediate to link y to the rapidity of a massless emission k with respect to a given leg . In fact, η (1) = y for y > 0, and η (2) = −y for y < 0. This implies that the phase space [dk] in Eq. (2.9) can be written as follows So far, all momenta have to be considered in d-dimensions, because dimensional regularisation is needed to regulate the IRC divergences present both in H(Q) and in the real radiation. The virtual corrections H(Q) can be expressed as [56,57] The virtual corrections are parametrised in three objects. First, w(m 2 , k 2 t + m 2 ; ) denotes the soft web [58,59] of total momentum k and squared invariant mass m 2 , in d = 4 − 2 dimensions. The web is obtained by considering the Feynman graphs with two eikonal lines that cannot be further decomposed into subgraphs by cutting each eikonal line once. For example, at lowest order in perturbation theory we find where C = C F in the case of quarks, and C = C A in the case of emitting gluons. Note that the web does not depend on the rapidity y, due to the properties of eikonal Feynman rules. Remarkably, at each order in perturbation theory, the web has a finite limit for → 0, which we will simply denote by w(m 2 , m 2 + k 2 t ). Second, the function γ (α s (k, )) coincides, up to an overall sign change, with the coefficient of the δ(1 − x) term of the regularised splitting functions P qq (x) and P gg (x), according to whether leg is a quark or a gluon, respectively. The strong coupling α s (k, ) is defined as the solution of the d-dimensional renormalisation group equation: where β (d=4) is the beta function in four dimensions, given by the following expansion In our representation of H(Q), the upper integration bound for the k t -integral of the web is set by the centre-of-mass energy Q, and the upper bound for the rapidity integral to |y| < ln(Q/ k 2 t + m 2 ). Finally, the overall quantity C(α s (Q)) is a multiplicative constant that is obtained by matching Eq. (2.13) at each order in perturbation theory to the quark or gluon form factor computed in the MS scheme.
In order to proceed, we need to define a procedure to cancel the IRC singularities in Eq. (2.13) against those in the real emissions. This can be done by introducing a resolution parameter that is engineered in such a way to divide the real radiation into a resolved set and an unresolved one. The idea behind this procedure is to handle the unresolved part of the radiation analytically, and hence cancel the divergences against the virtual corrections. The cancellation is performed in a manner that the resolved contributions could subsequently be computed numerically in d = 4 dimensions. The resolution parameter is defined through its action on the soft and/or hard-collinear contributions to the squared amplitudes, as outlined below.
We start by considering soft radiation. The soft squared amplitudes for n emissions, denoted hereafter as M 2 s (k 1 , . . . , k n ), can be iteratively reorganised as follows The quantitiesM 2 s (k 1 , . . . , k n ) represent the correlated portion of the n-emission soft amplitude squared, together with its virtual corrections. 4 This is strongly suppressed unless all emissions k 1 , . . . , k n are close in angle. We refer to the latter as soft correlated blocks and they play a dominant role in constructing the webs, which are the building objects of the Sudakov radiator to be defined below. Each correlated block admits a perturbative expansion in α s due to virtual corrections, hencẽ For instance, at tree level, the squared matrix element for the emission of a single soft gluon is given byM (2.20) In the following the coupling will be always renormalised in the MS scheme, i.e. we replace where β 0 is the first coefficient of the beta function in four dimensions of Eq. (2.16), given by The tree-level correlated block with two emissionsM 2 s,0 (k 1 , k 2 ) is reported in Appendix A, and will be useful later. This decomposition is particularly convenient to define a logarithmic counting. Each correlated blockM 2 s (k 1 , . . . , k n ) will contribute to ln Σ(v) at most with a factor α n s ln n+1 (v), with n powers of ln(v) coming from the soft singularities and an extra power from the only collinear singularity.
The definition of the resolution parameter proceeds as follows. First, we define a clustering algorithm that combines together the momenta of all particles emitted according to each correlated blockM 2 s (k 1 , . . . , k n ). For example consider the simple case of two emissions, the clustering is assigned as follows The property of rIRC safety [18] implies that all particles in a cluster are both close in angle and have commensurate transverse momenta. This allows one to evaluate the QCD running couplings of each cluster at the transverse momenta of the corresponding emissions. This procedure allows us to absorb all logarithms of µ R /k ti into the running of the coupling. As a consequence of this procedure, for rIRC safe observables, using the decomposition (2.17), every correlated blockM 2 s (k 1 , . . . , k n ) (when combined with the corresponding virtual corrections) will contribute to ln Σ(v) with terms of order α m s ln m+2−n (v) for m ≥ n. This allows us to build a logarithmic counting at the level of the squared amplitude, which defines which contributions must be considered at a given logarithmic order.
In order to proceed with the calculation of Σ(v), we then choose a resolution parameter δ 1 such that all clusters of total momentum k clust. satisfying are labelled as unresolved. This choice guarantees that one is able to compute analytically the contribution of unresolved emissions for an arbitrary rIRC safe observable. For this class of observables, the unresolved clusters can be neglected from the Θ function in Eq. (2.8) since they do not contribute to the observable V up to corrections suppressed by powers of δ p v, with p being a positive parameter. The above definition of the resolution parameter allows us to exponentiate the contribution of unresolved soft blocks. From the decomposition of Eq. (2.17), it is straightforward to connect the correlated blocksM 2 s (k 1 , . . . , k n ) to the webs introduced in Eq. (2.13). In fact where the factor S(n) represents the multiplicity coefficient for each soft final state (quarks or gluons). For instance, for n identical gluons, S(n) = 1/n!. Therefore, the contribution of an arbitrary number of soft clusters (and no hard-collinear clusters) gives rise to the following exponential factor Eq. (2.26) can be promptly combined with the virtual corrections in (2.13) to give where we defined the soft radiator R s as and we took the four-dimensional limit of the web (and the relative integration measure) since the integral is now finite. The next step is to handle the remaining hard-collinear divergences present in the integral over γ in Eq. (2.27). Unlike the case of soft radiation, the exponentiation of the unresolved hardcollinear emissions is more delicate in that every hard-collinear emission could potentially change the colour charge felt by subsequent radiation. However, the treatment of hard-collinear radiation is much simplified by observing that, owing to rIRC safety, only a fixed number of hard-collinear emissions is to be considered at a given logarithmic order. Therefore, instead of proceeding as in the soft case, it is convenient to start from the integral over the anomalous dimension γ in the virtual corrections (2.27) and split it into two pieces at the collinear scale of the resolution variable, that is found by setting the rapidity η ( ) to its maximum (i.e. ln(Q/k t )) in Eq. (2.7), which yields V sc (k) ∼ k a+b t . Inspired by this, the integral over γ then can be split at k = v 1/(a+b ) Q and becomes (2.29) Next, we expand the exponential of the second integral in the r.h.s. of the above equation considering only a fixed number of terms in its expansion as where, in our case of emitting quarks, the leading-order anomalous dimension is given by The first non-trivial order in the expansion must be included at NNLL, the second at N 3 LL (together with the squared of the first), and so forth. The divergences of these terms will cancel order-by-order in perturbation theory against those of the hard-collinear emissions in the real radiation. The last step to obtain a NNLL expression for Σ(v) is to handle the squared matrix element for real emissions M 2 in Eq. (2.8). At NLL accuracy, rIRC safety ensures that resolved radiation contains no hard-collinear emissions, and the real matrix element squared is approximated by its soft approximation M 2 s . Moreover, the squared amplitude at this order reduces to the product of n independent, soft-collinear emission probabilities. In fact, and C the colour factor of leg , C F for a quark and C A for a gluon. In order to achieve NNLL accuracy, it is sufficient to correct the products of independently emitted single-particle clusters with the insertion of a single tree-level correlated cluster of two soft and collinear emissionsM 2 s,0 (k a , k b ), and of the one-loop correction to the single-emission clusterM 2 s,1 (k). Moreover, beyond NLL, a finite number of hard-collinear emissions must be considered. In particular, at NNLL, it is sufficient to allow one single emission to be hard and collinear. When combined with Eq. (2.30), this leads to a finite, logarithmically enhanced, left over, as it will be shown shortly. In such configurations, at NNLL, the remaining soft radiation consists of an arbitrary number of single-emission clusters.
With the above decomposition, the NNLL resummed cross section Σ(v) of Eq. (2.8) takes the form where the squared amplitude M 2 s (k 1 , ..., k n ) is approximated by where V sc (k a +k b ) is defined as in Eq. (2.7), and k ( ) t , η ( ) , φ ( ) are the transverse momentum, rapidity and azimuth of the four-vector k a +k b with respect to leg . The cancellation of infrared and collinear divergences in the first term of Eq. (2.33) can be easily handled with a simple subtraction scheme as outlined in Ref. [54], which allows for a numerical evaluation in d = 4 dimensions. The cancellation of collinear singularities in the second term still requires the use of dimensional regularisation. In order to make the second term suitable for a numerical evaluation, we add and subtract the following and recast Eq. (2.33) as follows (2.36) The integral in round brackets of the last line of the above equation can be evaluated analytically as follows. For each leg = 1, 2, we expand the Θ(v − V sc (k hc )) function in the last line of Eq. (2.36) as where we neglect N 3 LL corrections in the expansion. With the above expansion it is sufficient to use the azimuthally averaged splitting function in d = 4 − 2 dimensions to construct the hard-collinear squared matrix element, since in Eq. (2.37) the only term that involves a non-trivial φ ( ) dependence is finite in d = 4 dimensions. In this approximation, the hard-collinear emission probability relative to each leg is given by where we introduced the hard-collinear radiator defined as ( 2.40) The hard-collinear constant C hc, is given by Finally, the constant part of virtual corrections at NNLL is given by where Each of the terms in the resolved contribution to Eq. (2.39) can be further decomposed into a finite set of corrections so that the NNLL cross section Σ NNLL (v) can be parametrised with the following master formula (we define λ = α s (Q)β 0 ln(1/v)) where the functions F NLL and δF NNLL have a general expression for any rIRC safe observable [18,54,55] and can be efficiently evaluated numerically in d = 4 dimensions. The NNLL function is decomposed as follows where each term has a well-defined physical origin.
• The corrections δF sc , δF wa , δF correl and δF clust have soft origin and they all originate from the first term in the curly brackets of Eq. (2.39).
The function δF sc accounts for running coupling corrections to the real emissions in the CMW scheme [61] as well as for the correct rapidity boundary for a single soft-collinear emission.
For event shapes variables this correction is particularly simple [54] owing to the fact that the rapidity dependence of the observable can be always handled analytically. For observables with a more complicated rapidity dependence, such as jet rates [55], the running coupling correction and the correction due to the rapidity boundary must be treated separately.
The function δF wa accounts for the difference between the observable and its soft-collinear parametrisation for a single soft-non-collinear (wide-angle) emission accompanied by many soft-collinear gluons.
At NLL all resolved emissions are strongly ordered in angle, and thus emitted independently. The matrix element used to compute the function F NLL is simply given by a product of an arbitrary number of single-gluon emission squared amplitudesM 2 sc (k i ), in the decomposition of Eq. (2.17). However, starting from NNLL two or more resolved emissions can become close in angle. In this type of configurations, the squared amplitude is given by an abelian term (defined by the product of n single emission probabilities) and by non-abelian, correlated clusters of two or more particles (see Eq. (2.17)). At NNLL it is sufficient to account for the effect of only two emissions getting close in angle, while the others can be considered far apart. This induces two types of corrections: δF correl and F clust .
The correlated correction δF correl accounts for the insertion in the resolved ensemble of softcollinear, independently emitted gluons of a single double-soft clusterM 2 s,0 (k 1 , k 2 ) that is defined as the non-abelian part of the square of the double-soft current [54]. More details on this correction will be given in Sec. 3.4.
The clustering correction δF clust (defined in Ref. [55]), on the other hand, accounts for the contribution of two independently emitted gluons that become close in angle. Due to the nature of most event shapes, this correction normally vanishes and it becomes different from zero only when the observable has a non-trivial dependence on the rapidity of the emissions, as for instance in the case of jet rates [55].
• The corrections δF hc and δF rec originate from the second term in the curly brackets of Eq. (2.39), and have a hard-collinear nature. The emission of a hard collinear parton induces two types of corrections: at the level of the squared amplitude (encoded in δF hc ), and at the level of the kinematics, due to the recoil of the whole event against the hard-collinear emission (encoded in δF rec ).
• Finally, the term arises from the first and third term in the curly brackets of Eq. (2.39), where N 3 LL corrections were neglected. The function F NLL is purely NLL [18], while the multiplying constants H (1) and C hc, induce NNLL corrections.
Since the detailed formulation of the functions F NLL and δF NNLL is given in refs. [18,54,55], we do not report their expressions here, and refer to the original publications.
In the next section we perform the calculation of the Sudakov radiator at NNLL accuracy. Before we proceed, it is important to stress that in the definition of the unresolved soft radiation given in Eq. (2.24) one clearly has some freedom in deciding precisely how the resolution variable is defined. In particular, instead of V sc , one could use any observable V res that shares the same leading logarithms as the full observable V that is being resummed. This is in fact the only requisite for this method to be applied. The choice of V sc is mainly due to computational convenience, as all ingredients in the Sudakov radiator can be computed analytically for any rIRC safe observable. Choosing another resolution variable would change the expression in the Sudakov radiator, and consequently the expression of the functions F NLL and δF NNLL .
A particularly important aspect of the definition of the resolution variable concerns the way V sc in Eq. (2.24) is evaluated on the total momentum k clust. of a cluster of more than one particle. Although the cluster has a non-zero invariant mass, V sc (2.7) does not depend on the mass, and hence, in the definition resolution scale, the cluster is treated as if it were massless. This will lead to great simplifications in evaluating the Sudakov radiator, where the integral over the invariant mass can be evaluated analytically as it will be shown in Sec. 3. The prescription of treating the cluster as massless also impacts the definition of the correlated correction δF correl , that will be the subject of Sec. 3.4.

The Sudakov radiator at NNLL accuracy
In this section we explicitly compute the Sudakov radiator defined in Eqs. (2.28) and (2.40) to NNLL.

The soft radiator
With the choice of resolution variable as in Eq. (2.7), the soft radiator reads where the phase space is given in Eq. (2.10). Notice here that the phase-space measure contains the massive rapidity y of the web as defined in Eq. (2.10), while the observable is expressed in terms of the rapidity of a massless parton η ( ) . Now we wish to write Eq. (3.1) in a way in which LL, NLL and NNLL contributions are separated. The first step to achieve this is to isolate the dependence on d g (φ) by expanding the observable constraint in Eq. (3.1) as follows 2) The first term in the above equation starts at LL accuracy, the second at NLL accuracy, and so on. This gives

(3.4)
We now concentrate on R (v). The kinematic boundary for the rapidity integral is ln(Q/ k 2 t + m 2 ). Instead of computing directly the integral in Eq. (3.4), we split it into the sum of two terms as where R 0 (v) is defined as in Eq. (3.4) but with a massless rapidity boundary, i.e. ln(Q/k t ). This defines a massless radiator, in which η ( ) coincides with y, i.e.
that starts at LL accuracy. The function δR (v) defines a mass correction, which accounts for the correct rapidity bound. By inspecting the phase space constraints due to the physical rapidity bound and the observable, one finds that the rapidity integral is bounded by ln(Q/ k 2 t + m 2 ) only when k t > v 1 a+b Q. This leads to the following expression for the mass correction to the soft radiator The separation of the radiator as in Eq. (3.5) has a physical justification. If one ignores the running of the coupling constant, then the massless radiator R 0 (v) only contains double logarithmic terms, while the mass correction δR (v) is purely single logarithmic.
Let us now focus of the massless radiator at NNLL. Eq. (3.6) can be easily evaluated by observing that the resolution variable does not depend on the mass of the web. Therefore, in Eq. (3.6) we can freely integrate over this variable. The integral of the web w(m 2 , k 2 t + m 2 ) over its invariant mass defines a generalisation of the physical CMW coupling [61]: The set of constants K (n) is perturbatively calculable and, once identified, the massless radiator R 0 (v) is fully determined for any observable, and given by At NNLL accuracy, in the expression of α phys s one needs to include only K (1) and K (2) , whose expressions are obtained by integrating the web up to order α 3 s . This requires contributions up to the triple-soft current at tree level [62], the single-soft two loop current [63], and the double-soft current at one loop. One obtains Note that K (1) is proportional to the two-loop cusp anomalous dimension, but this is not true any more starting from K (2) . We stress that this is not the case for rIRC unsafe cases, such as threshold resummation. In this case, one needs to perform the integration over the web transverse momentum in d dimensions and subtract the residual collinear singularity in a given factorisation scheme. The divergent integral over k t then gives an extra contribution that enters at same order as K (2) so that the coefficient of the α 3 s L 2 term in the Sudakov coincides with the cusp anomalous dimension. For the computation of the mass correction at NNLL, we only need to consider the web up to α 2 s . In fact, the only non-vanishing contribution arises from the double-emission soft block, M 2 s,0 (k a , k b ), that can be found in Appendix A. Using the rescaled variable µ 2 = m 2 /k 2 t , one finally gets (3.12) From the above derivation we observe that the O(α 3 s ) correction to the physical coupling (3.9) is the only place where the α 3 s web enters in a NNLL resummation of any rIRC safe observable. Indeed the resolved real corrections only involve correlated blocks with up to two soft partons. Therefore, the physical coupling defined in Eq. (3.9) constitutes a universal ingredient to account for the triplecorrelated soft contribution (and relative virtual corrections due to the double-correlated soft at one loop as well as single-correlated soft at two loops) at NNLL accuracy. In particular, it defines a building block of a parton shower algorithm at this order, that will be relevant in the context of the current efforts that aim at improving the accuracy of these algorithms [64][65][66][67][68][69][70][71].

The hard-collinear radiator
The hard-collinear part of the radiator, defined in Eq. (2.40), starts at NLL accuracy. Up to NNLL accuracy, its expression is In our case, the coefficients γ (0) and γ (1) are the coefficients of the δ(1 − x) piece of the P qq (x) splitting functions with an overall minus sign. In particular, γ (0) is given in Eq. (2.31), and (3.14)

The radiator up to NNLL accuracy
The computation of the radiator proceeds by integrating the equations of the above sections with a running coupling. In particular, at NNLL accuracy, this is given by the renormalisation group equation where β 0 is given in Eq. (2.22), and the other coefficients of the beta function, in the MS scheme, are given by For resummation purposes, it suffices to solve Eq. (3.15) using the following ansatz Plugging the above in Eq. (3.15), one finds and we can neglect the contributions of f 4 (t) and beyond, since they start to matter from N 3 LL accuracy.
It is customary to express the radiator in terms of λ = α s β 0 ln(1/v), with α s = α s (Q). We further parametrize the radiator in terms of functions of λ, in such a way as to separate LL, NLL and NNLL contributions: , . (3.30) The corresponding results for b = 0 can be obtained by taking the limit of the above expressions for b → 0. We also include here various derivatives of the massless radiator that appear in the evaluation of the correction functions: . (3.33)

The correlated correction δF correl
The definition of the correlated correction δF correl is given by [ where k a and k b are the two soft emissions close in angle and collinear to the same leg, while the remaining soft-collinear emissions k i have very disparate angles, and hence are emitted independently. The invariant mass of the web is denoted by m 2 = (k a + k b ) 2 . The configuration in which k a and k b are collinear to different Born legs requires the parent gluon to be emitted with a large angle, and hence gives at most a N 3 LL contribution. The observable simply denotes the soft-collinear approximation of the full observable, as this is the only limit relevant here. In order to make sense of the difference between the two Θ functions, we should first discuss how the cluster of two soft partons is treated in the ARES algorithm.
The two-particle correlated soft blockM 2 s,0 (k a , k b ) diverges when k a and k b are collinear, i.e. when the invariant mass of the cluster tends to zero. Such a divergence is entirely cancelled by the one-loop correction to the single-emission clusterM 2 s,1 (k), that should be included at NNLL. In the evaluation of the Sudakov radiator, as shown in Sec. 3, the cancellation of the above collinear singularity is performed analytically. This is because we chose to define a resolution variable that does not depend on the invariant mass of the cluster, and hence the corresponding integral becomes straightforward. This inclusive integration ofM 2 s,0 (k a , k b ) andM 2 s,1 (k) is responsible for the presence of the physical coupling in the radiator, see Eq. (3.9).
In the resolved radiation, however, the situation is more complicated, as in the full observable we are not allowed to integrate over the invariant mass of the cluster inclusively. We can, however, define the following subtraction scheme to cancel the collinear singularity arising fromM 2 s,0 (k a , k b ) in four dimensions. We first treat the {k a , k b } cluster inclusively, as done in the definition of the radiator. This can be done by considering only the total momentum k a + k b when evaluating the contribution of the cluster to the observable, and treat it as if it were a massless (lightlike) momentum in the computation of the observable. This once again allows us to combine it with the one loop correction to the single-emission clusterM 2 s,1 (k) analytically. This contribution is encoded in δF sc which features the coefficient K (1) of Eq. (3.11a), as explained in refs. [54,55]. As a second step, we consider the difference between the full observable, where the {k a , k b } cluster is treated exclusively, and its inclusive approximation that we considered above to cancel the singularity against the virtual correction. This is represented by the difference in the two Θ functions in Eq. (3.34). It is important to bear in mind that in the second Θ function the observable V sc treats the momentum k a + k b as if it were massless, in order to exactly match our convention for the cancellation of real and virtual corrections in δF sc . This is implemented by the limit in the second Θ function in Eq. (3.34).
There are many ways to parametrise the phase space for k a and k b to keep, in δF correl , only NNLL contributions that correspond to configurations in which k a and k b are collinear to the same Born leg. One possible parametrisation was presented in Ref. [54]. Instead, here we adopt the notation of appendix A, that is the same we used to compute the NNLL radiator. We define the rescaled invariant mass of the {k a , k b } cluster µ 2 = m 2 /k 2 t , and we introduce the pseudo-parent parton k of k a and k b with transverse momentum k t and observable fraction ζ defined as Using Eqs. (A.6) and (A.8), the squared amplitude for a double-soft correlated emission reads with C ab (µ, z, φ) = C(k a , k b ) given by Eq. (A.9), and µ 2 ∈ [0, ∞), z ∈ [0, 1], φ ∈ [0, 2π). The matrix element squared and phase space for the pseudo-parent k is given in Eq. (2.32). Actually, due to the fact that the pseudo-parent has a non-zero invariant-mass, the integral over its rapidity η ( ) should have the boundary |η ( ) | < ln(Q/ k 2 t + m 2 ). However, following what is done in the computation of the NLL function F NLL [18], we observe that the exact position of the rapidity integration bound in the resolved radiation enters at one logarithmic order higher. Therefore, in order to neglect all N 3 LL corrections and obtain a result that is purely NNLL, we replace the actual rapidity integration limit with the massless one, as done in Eq. (2.32). The integral over η ( ) can be evaluated analytically [54] and the correlated correction takes the following simple form (3.38) Naturally, for analytic calculations the parametrisation for the phase space and matrix element should be chosen in order to simplify the integrals for any given observable. The present choice will make the integrations in the next section simpler, while an alternative parametrisation was reported in Ref. [54].

Additive observables
A particularly interesting case is that of additive observable, for which all NNLL corrections admit a simple analytic form, reported in the appendix of Ref. [54]. The correlated correction for such observables can be simplified considerably. The additivity of the observable implies that (3.39) We now introduce f ( ) correl (µ, z, φ) as: This gives (3.41) Following the derivation in the appendix of Ref. [54], we can now rescale the momenta k 1 , . . . , k n in two ways. In Eq. (3.38), in the term involving the step function of V sc ({p}, k a , k b , k 1 , . . . , k n ), we construct rescaled momentak 1 , . . . ,k n such that In the term with the step function of V sc ({p}, k a +k b , k 1 , . . . , k n ) we construct another set of rescaled momentak 1 , . . . ,k n such that Performing similar formal manipulations as in Ref. [54], we obtain After performing the ζ and φ ( ) integrations analytically, we obtain the factorised form in which the correlated correction reduces to a number that multiplies the NLL function F NLL . This result will be used in Sec. 4 to obtain the NNLL resummation for angularities and moments of energy-energy correlation.

Fully worked out examples: angularities and moments of EEC
In this section we apply the resummation procedure described in the previous sections to interesting observables in e + e − annihilation, namely angularities and moments of energy-energy correlation. Angularities are defined with respect to some reference axis, usually the thrust [72] or the winner-take-all (WTA) [17,73] axis. They depend on a parameter x, as follows where the sum runs over all hadrons in the event, (E i , q i ) is the four-momentum of hadron i, and θ i is the angle between hadron i and the reference axis. In Ref. [18], another class of observables was introduced, the fractional moments of energyenergy correlation (EEC), defined by where, as before, the sums run over all hadrons in the event, θ ij denotes the angle between hadrons i and j, and n T is the thrust axis. Note that similar variables have attracted interest due to their discriminating power between quark-and gluon-initiated jets [73]. For instance, in jet studies for e + e − collisions, one considers where one considers the particles i, j within a given jet. At hadron colliders the definition of C (β) 1 [73] involves the transverse momentum and the angular distance R 2 ij = ∆y 2 ij + ∆φ 2 ij between final state particles. The global component of the resummed cross section for these observables has analogous resummation properties as the observables F C 2−β studied here. However, in this case, the cross section receives a non-global logarithmic correction starting at NLL. Both angularities with respect to the WTA axis and moments of EEC have the property that, in the presence of multiple soft and collinear emissions k 1 , . . . , k n , they are always additive, i.e.
Angularities with respect to the thrust axis are additive as long as x < 1 [18]. For these observables, x < 1 is the range of values of x that we will implicitly consider in the following. The NNLL resummed distribution is given by Eq. (2.45). Our task is to compute each ingredient of that formula. First, we observe that both τ x and F C x are infrared safe, and collinear safe for x < 2. With a single soft and collinear emission, and in the range of values of x appropriate for each observable, we have Therefore, the soft-collinear radiator R s (v), the hard-collinear radiator R hc (v) and the hard-collinear constant C (1) hc, are obtained by computing Eqs. (3.3), (3.13) and (2.41) respectively, with a = 1, b = 1−x, and d = g (φ) = 1. In the following subsections, we compute the corrections due to real radiation.
The results we will present below for the WTA-axis angularities have been found to be in complete agreement with the findings of Ref. [33], that have been obtained in a SCET framework.

Soft-collinear corrections
Since the observables we consider are additive, they fall into the category studied in appendix C of Ref. [54]. This gives Here, for notational convenience, we introduced R NLL = R NLL,1 + R NLL,2 , and similarly for R NNLL and R . Note that both expressions in Eqs. (4.6) and (4.7) do not depend explicitly on the parameter x, but this dependence is implicit in the functions R NLL , R NNLL and R .

Hard-collinear and recoil corrections
If we add to an ensemble of soft and collinear emissions k 1 , . . . , k n a single hard emission k, collinear to either p 1 or p 2 , our observables behave as follows For the angularities with respect to the thrust axis, we obtain whereas if we compute angularities with respect to the WTA axis, we obtain (4.10) Finally, for fractional moments of EEC, we get If we extrapolate f ( ) hc (z, φ ( ) ) for z → 0 we obtain the same result for all observables: The function f ( ) sc is the only one needed to compute the correction δF hc according to the procedure described in appendix C of Ref. [54], which leads to sc . Specialising the formulae of appendix C of Ref. [54] to the present case, for the angularities with respect to the thrust axis, we obtain (4.14) If we consider τ x with respect to the WTA axis, we obtain Finally, for the moments of EEC we obtain (4.16)

Soft wide-angle corrections
If we add to an ensemble of soft and collinear emissions, k 1 , . . . , k n , a single soft emission k, at an angle θ with respect to the thrust axis that is much larger than that of all other emissions, we have | sin θ k1 | | sin θ k2 | | sin θ| , (4.17) where θ k is the angle between k andp , with = 1, 2. Also, for any appropriate value of x, since k is the emission at the largest angle, its transverse momentum with respect to its emitter is the same as that with respect to the thrust axis. Therefore, for all considered observables, we have (4.18) Therefore, V wa ({p}, k, k 1 , . . . , k n ) = V sc ({p}, k, k 1 , . . . , k n ), and δF wa = 0.

Correlated corrections
Since all the observables we consider are the same in the soft and collinear limit, whenever we have any two soft emissions k a , k b , collinear to the same leg , together with an ensemble of soft-collinear emissions k 1 , . . . , k n , we obtain In terms of the variables defined in appendix A, we have Using the rescaled variables µ 2 ≡ m 2 /k 2 t and u i ≡ q i /k t , and using the notation of section 3.4.1, we obtain (4.22) Note that, for x = 0, which is the same as 1 − T , we have f correl (z, µ, φ) = 1 This implies that, for x = 0, one has Therefore, only for x = 0 are the considered observables fully inclusive with respect to multiple collinear splittings. This result generalises to an arbitrary number of soft and collinear emissions. Now we can compute δF correl for any value of x using the general formula in Eq. (3.44). We obtain where (4.26) In the above equation, 2S and H g are defined in Eqs. (A.11a) and (A.11b) respectively, and H g is defined in Eq. (A.11c). We have computed lnf correl C A and lnf correl n f numerically as a function of x, and the result can be found in Fig. 1. For x = 0, lnf correl C A and lnf correl n f can be computed analytically, which gives (4.27)

Matching and issues with Sudakov shoulders for F C x
It is interesting to study the matching to fixed order for the moments of energy energy correlation (4.2). Such observables feature a Sudakov shoulder [74], whose position can get dangerously close to the Sudakov peak for certain values of x. To examine this feature we now match the resummed NNLL distributions to NLO fixed-order differential cross sections obtained with EVENT2 [75]. Although we only analyse F C x below, the procedure discussed in the following applies to all observables considered in this article. The matching is performed according to the log-R scheme (see for instance [25,76]). As it is customary in resummed calculations, to probe the size of subleading logarithmic terms we introduce a rescaling constant x V as and expand the cross section around ln x V /v neglecting subleading terms. 6 Eventually we modify the resummed logarithm ln x V /v in order to impose that the total cross section is reproduced at the kinematical endpoint v max (4.29) Here, p denotes a positive number which controls how quickly the logarithms are switched off close to the endpoint. Since in the following we do not perform a phenomenological study, we simply set p = 1 and v max = 1 for the sake of simplicity.
To obtain our central predictions we set µ R = Q, with Q being the centre-of-mass energy of the hard scattering, corresponding to α s (M Z ) = 0.118, and x V = 1. We then construct the uncertainty bands by varying µ R and x V individually by a factor of two in either direction. The relevant formulae for the scale dependence are reported in the appendix of Ref. [54].  Figure 2 shows the NLO differential distribution matched to both NLL and NNLL for the cases x = 1/2 and x = 1. From the plots one can appreciate the following two interesting features.
The first is that the size of NNLL corrections increases with x. This can be explained by inspecting the parametrisation of the observable in the soft and collinear limit (4.5). The characteristic transverse momentum of the soft radiation is k t ∼ Qv, while the hard-collinear radiation occurs at scales k t ∼ Qv 1 2−x , with x < 2. The soft scale is therefore lower than the collinear scales for x < 1, the two coincide for x = 1, and the situation is inverted for 2 > x > 1. For a given value v of the observable, the typical size of the soft logarithms (and hence of the soft corrections) does not depend on the moment parameter x. Conversely, the size of the hard-collinear logarithms increases with x, hence leading to larger subleading corrections. One also expects that corrections beyond NNLL become more sizeable as x increases, as it is reflected by the scale uncertainty band in Figure 2. For x > 1 the subleading corrections grow very large as the collinear scale becomes smaller than the soft one, which corresponds to a badly convergent logarithmic series. A consequence of this fact is that for x > 1 the abscissa of the Landau pole moves towards larger values of v, and hence the differential distribution becomes non-perturbative at moderate values of v.
A second interesting observation is that the relative distance between the Sudakov peak and the shoulder decreases for increasing x. This implies that there is a value of x for which the two overlap. Such a situation can be observed in the left plot of Figure 3 for x = 3/2, where the curves are obtained at the Z resonance and for central values of the scales.
In this case the solution provided by the resummation in the two-jet limit is obviously unphysical. The position of the Sudakov peak represents the bulk of the soft and collinear radiation probability in the two-jet configuration, and it coincides with the kinematic endpoint for the threejet configuration, above which the distribution is again dominated by soft and collinear emissions. This phenomenon is due to the violation of momentum conservation in the formulation of the resummed calculation, in which the exact kinematics in the presence of an extra hard parton is ignored. In such a situation one should perform a simultaneous resummation of the Sudakov  logarithms treated here together with the logarithms that originate at the shoulder. This is currently out of reach at the logarithmic order analysed in this article.
While in this case a matching to fixed order results in an unphysical prediction, the right plot of Figure 3 shows that the situation improves at higher collider energies. As can be seen from this plot, for higher collider energies the position of the Sudakov peak moves towards smaller values of the observable, driven by the smaller coupling constant, while the position of the shoulder does note depend on Q. At these scales the resummed result is physical and can be matched to the fixed order. An example is reported in Figure 4, where the matched distribution for Q = 1 TeV is shown.

Conclusions
In this paper we have completed the study of jet observables at NNLL accuracy in e + e − annihilation, that started in refs. [54,55]. These results constitute the core of the ARES method for the seminumerical resummation of jet observables at NNLL accuracy, which generalises the NLL procedure of refs. [18][19][20]. This involves the calculation of an observable dependent Sudakov form factor, the radiator, which encodes the all-order cancellation of infrared singularities between real and virtual contributions, and that we have computed at NNLL accuracy for a generic rIRC safe observable.
As a byproduct, we have defined a generalisation of the well known CMW physical coupling in the soft limit, and given a closed expression for its relation to the MS coupling up to O(α 3 s ).
This quantity is a universal ingredient for all resummations of rIRC safe observables, and as such it constitutes one of the main ingredients for a NNLL accurate parton shower algorithm.
As an application, we have computed NNLL resummed distributions for angularities and fractional moments of EEC, for all allowed values of the parameter x they depend on. We have also presented a very basic phenomenology of moments of EEC, highlighting their main features. A particularly severe issue is that, for x > 1, the Sudakov peak of the differential distribution, where the observable should be dominated by multiple soft-collinear emissions, becomes dangerously close to the edge of the phase space for real emissions, so that the core approximation underlying softcollinear resummations breaks down. This situation is severe at LEP energies and prevents us to make any sense of resummed prediction matched to fixed order for certain values of x > 1. This is not the case at a future e + e − collider with centre-of-mass energy of 1 TeV, where the Sudakov peak moves towards lower values of the observable, while the position of the kinematical boundary stays unchanged. This is a feature that should be considered when using these observables for phenomenology.
We stress that the procedure outlined in this paper is only an example on how to construct a Sudakov radiator. The same calculation could for instance be performed at all orders using the methods of Soft-Collinear Effective Theory (SCET), along the lines of what has been shown for thrust in Ref. [77]. To achieve a full resummation, the radiator should be supplemented by appropriate NLL and NNLL corrections due to the resolved real radiation. In this paper, we have defined the Sudakov radiator in such a way that all NLL, and most NNLL corrections are the same as in refs. [54,55]. The only exception is the NNLL correction δF correl , which we had to redefine to ensure that the radiator could be computed analytically for a generic rIRC safe observable.
We also note that this very same radiator could in principle be used also for processes with incoming hadrons, for reactions with two hard emitting legs at the Born level. In that case, all corrections due to soft radiation, F NLL , δF sc , δF wa , δF correl stay unchanged, and one has to evaluate parton distribution functions at the factorisation scales of order v 1/(a+b ) Q, and recompute only the hard-collinear contributions C (1) hc, , δF hc , δF rec , as done for instance in Ref. [51,78] for certain classes of observables. For processes with more than two legs, all terms in the master formula (2.45) must be redefined in order to account for the structure of the wide-angle soft radiation. This will be left for future work.
In summary, this work completes the formulation of a general method for to the calculation of any jet observable in processes with two legs, that can be systematically generalised to more complicated cases. We hope that the results presented here will define a solid starting point for future systematic studies of jet observables at all perturbative orders.