Collinear fragmentation at NNLL: generating functionals, groomed correlators and angularities

Jet calculus offers a unique mathematical technique to bridge the area of QCD resummation with Monte Carlo parton showers. With the ultimate goal of constructing next-to-next-to-leading logarithmic (NNLL) parton showers we study, using the language of generating functionals, the collinear fragmentation of final-state partons. In particular, we focus on the definition and calculation of the Sudakov form factor, which physically describes the no-emission probability in an ordered branching process. We review recent results for quark jets and compute the Sudakov form factor for the collinear fragmentation of gluon jets at NNLL. The NNLL corrections are encoded in a $z$ dependent two-loop anomalous dimension $B_2(z)$, with $z$ being a suitably defined longitudinal momentum fraction. This is obtained from the integration of the relevant $1\to 3$ collinear splitting kernels combined with the one-loop corrections to the $1\to 2$ counterpart. This work provides the missing ingredients to extend the methods of jet calculus in the collinear limit to NNLL and gives an important element of the next generation of NNLL parton shower algorithms. As an application we derive new NNLL results for both the fractional moments of energy-energy correlation $FC_x$ and the angularities $\lambda_x$ measured on mMDT/Soft-Drop ($\beta=0$) groomed jets.


Introduction and motivation
In this paper and in a forthcoming article [1], we initiate a formulation of the methods of jet calculus [2][3][4][5] beyond the NLL order to describe the dynamics of collinear fragmentation.Jet calculus techniques, and specifically the generating functional method, stand out as a mathematical language to formulate self-similar branching processes.For this reason, they are particularly suitable to bridge the field of QCD resummation to that of parton shower Monte Carlo algorithms.This is an important direction to pursue in the context of improving the logarithmic accuracy of such algorithms for QCD phenomenology and is receiving widespread attention in the collider physics context (see e.g.Ref. [6] for an overview).
More specifically this work addresses the following points: 1.A large class of QCD resummations cannot be handled using purely analytic techniques since they do not admit a closed analytic form.This can be either due to the non-linear nature of the evolution equations (like for instance the resummation of microjet observables [7,8], fragmentation [9][10][11][12][13][14][15][16][17][18][19] or non-global observables [20][21][22][23][24][25][26][27][28][29][30]) or to the high complexity of the observable which may not easily factorise in a suitable conjugate space (like for instance some complex event shapes or jet resolution parameters [31][32][33]).In such cases it is necessary to formulate the resummation in a way that can be solved accurately using numerical methods.The methods of jet calculus [2][3][4][5], notably Generating Functionals (GFs), offer a powerful mathematical tool to describe the resummation of logarithmic radiative corrections in collider observables, and constitute one avenue to achieve a numerical solution via Monte Carlo techniques.
2. The recent development of new techniques to construct parton-shower algorithms with demonstrable higher logarithmic accuracy (see e.g.Refs.[34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53]) offers an opportunity to further increase systematically their precision.The connection between these new developments in the area of parton showers and resummations is based on a well-defined set of accuracy criteria that are based directly on QCD dynamics (i.e.approximating correctly the real and virtual matrix elements in specific kinematical limits).Recent parton-shower algorithms [40,42,45,46,[48][49][50][51]53] are constructed in such a way that the above criteria are satisfied, hence achieving NLL accuracy for broad classes of observable at once.GFs help us push this important correspondence to a deeper level, by deriving evolution equations which provide an analytic connection between parton showers and all-order resummations.As a result, one can unify parton showers and resummations within a single framework, which allows for systematic developments in both areas simultaneously.
In this work we focus on the class of observables which feature only collinear sensitivity, and on the description of final state fragmentation, leaving the extension to initial state radiation to future work. 1 Such observables are by definition only sensitive to non-soft small-angle emissions, where the longitudinal momentum fraction z of the splitting satisfies z ∼ (1−z) ∼ O (1).In this regime, the relative angle and normalised transverse momentum of each splitting are commensurate, and one can expand one about the other systematically.Examples of observables of this kind are moments of energy-energy correlations [54] or angularities [55] measured on mMDT/soft-drop (β = 0) groomed jets [56,57], which we will study here, or fragmentation functions [58][59][60][61][62].
Here we calculate the analytical ingredients needed in the formulation of NNLL resummation within the GFs method.As we will demonstrate in forthcoming work [1], these ingredients can then also be used in Monte Carlo algorithms to correctly include virtual corrections, 2 which are not captured by the shower's unitary approximation.For this reason, ingredients of this type (such as the soft physical coupling scheme [63][64][65]) have to be computed in dimensional regularisation.
The aforementioned analytical ingredients amount to the anomalous dimensions defining the NNLL Sudakov form factor, which encodes the no-emission probability in the branching process, and is the back bone of both the GFs method and parton-shower algorithms.After reviewing the recent work of Ref. [66], which provides the NNLL Sudakov for quark jets, here we frame these results in the context of the GFs method and furthermore define and compute the NNLL Sudakov form factor for the case of gluon fragmentation.We also discuss applications of our results in the context of QCD resummation, by deriving new NNLL accurate results for specific groomed jet observables related closely to those measured at the LHC.For these observables, the evolution equations can be solved in closed form and allows us to derive analytic results.
The layout of the paper is as follows.In Sec. 2 we introduce the Generating Functionals method as a way to understand the role of such analytic ingredients in a Monte Carlo calculation.We then discuss the extension to NNLL order in the collinear limit and define the NNLL Sudakov form factor and the anomalous dimension B 2 (z).In Secs.4-7 we calculate B g 2 (z) for gluon jets, which complements previous studies done for quark jets and discuss the important differences between these two cases.In Sec. 8 we provide concrete applications to collider phenomenology by deriving new results for the NNLL resummation of moments of energy-energy correlation [54] and angularities [55] measured on groomed jets with the mMDT/Soft-Drop procedure [56,57].We conclude in Sec. 9 and provide a number of appendices with additional technical details about the method and the calculations performed in this paper.The results for the B 2 (z) anomalous dimensions for quarks and gluons are given in electronic format with the arXiv submission of this article.

Generating functionals for collinear fragmentation
In this section we introduce the GF method to resum logarithmic corrections of collinear nature.We start by reviewing the NLL case and then we discuss the extension to NNLL.mathematical language was discussed in Refs.[28,29] in the context of non-global resummations in the large−Nc limit.
2 Specifically, by unitary approximation we mean that in a parton shower virtual corrections simply amount to (minus) the integral of the real (soft and/or collinear) matrix element.This implies a unitary effect on the total cross section.

Review of NLL resummation
We consider the collinear fragmentation of a parton of initial energy E resolved at an angular resolution θ ≪ 1.The all-order resummation at NLL resums single logarithms L = − ln θ of the form α n s L n .At this order, the fragmentation can be formulated as a branching process in which emissions are strongly ordered in angle. 3At NNLL we aim at resumming corrections of order α n s L n−1 , where now emissions can be unordered and have commensurate angles.To start, we therefore define an angular resolution scale t (evolution time) as where E p is the energy of the parton branching at angular resolution t i .In Eq. (2.1) α s is the MS coupling, and g(z) is a function of the longitudinal momentum fraction 1 − z carried by the i-th emission.At NLL the precise form of g(z) is irrelevant and one can set g(z) = 1 as well as E p = E (see for example Ref. [7]).In what follows, the use of the β(α s ) function that drives the running of α s at either one-loop or two-loop order, depending on whether the target accuracy is NLL or NNLL is understood.
We introduce the GF G f (x, t), which encodes the probability of resolving a fixed number of emissions in the time-like fragmentation of an initial parton of flavour f ∈ {q, g} and momentum fraction x (i.e.E p ≡ x E) below a resolution angle set by an initial evolution time t.The GFs are defined in such a way that the probability of exclusively resolving m partons in the collinear fragmentation of a parton of flavour f ∈ {q, g} and momentum fraction x starting at an evolution time t is given by [2,4] 2) The quantity u ≡ u(x, t; f ) is the probing function and has the role of tagging a real emission in the functional derivative of G f .Eq. (2.2) can be taken as the definition of the generating functional G f .The evolution of the GFs with the resolution scale t is described, at NLL, by the coupled system of integral equations [7] 4 3 Angular ordering ensures the full coverage of the relevant phase space at NLL.Beyond NLL, emissions can have commensurate angles and therefore the precise definition of angular ordering (i.e. with respect to a specific reference direction) has to be specified in order to guarantee the coverage of the full phase space. 4An equivalent differential form can be easily obtained by dividing G f by ∆ f and subsequently taking the t derivative.Note that Ref. [7] adopts a slightly different definition of the GFs which in this reference describe the fragmentation of a parton from a starting time t = 0 down to a resolution angle set by the time t.This convention is complementary to the one adopted in this paper, but their difference is irrelevant at the level of physical results.
where t 0 is a collinear cutoff at which the evolution stops, while z 0 is an infrared cutoff on the energy fraction of each emission, which must be taken to zero in the calculation of an IRC safe observable (unlike for z 0 , the value of t 0 is bounded by the presence of the Landau pole).The standard leading-order splitting functions are given in Appendix A.
The above set of equations can be used to derive collinear resummations with NLL (single logarithmic) accuracy for final state radiation.The Sudakov form factors ∆ f describe the no-emission probability and can be derived by imposing the unitarity condition from which one obtains An obvious boundary condition is then G f (x, t 0 ) = u, indicating a 100% probability of finding a parton f with momentum x.With the above definition of the GFs, the resummed distribution (or equivalently cumulative distribution) dσ (f ) for a given observable v = V ({k i }) (with {k i } denoting the final state momenta produced in the fragmentation of a jet of initial flavour f ) has the general form where σ 0 is the Born cross section for the hard process under consideration.The jet distribution J (f ) is simply obtained by integrating Eq. (2.2) with the measurement function of the observable for any final-state multiplicity m, that is The ⊗ operation is observable dependent.It is usually a regular product, but for some specific observables (e.g.fragmentation functions) it can take the form of a convolution over the longitudinal momentum fraction.The process-and observable-dependent matching coefficient C(α s ) admits a fixed-order perturbative expansion in powers of the strong coupling constant C(α s ) = 1 + O(α s ), and it accounts for constant terms stemming from the matching of the jet distribution (defined by the GFs evolution equation) to the fixed order QCD calculation in the limit v → 0. Specifically, at NNLL C(α s ) is required at the one-loop order, which entails the difference between the full O(α s ) QCD calculation in the logarithmic limit (i.e.v → 0) for the observable v and the expansion of the jet function at the same order.

Extension to the NNLL case
The extension of the above formulation to NNLL order requires the generalisation of the r.h.s. of the evolution equations (2.3) to O(α 2 s ).To understand what the calculation entails, we follow an analogous derivation of that presented in Refs.[28,29] in the context of softgluon evolution in the planar limit.We observe that Eqs.(2.3) describe the fragmentation by a sequence of angular ordered branchings.The resulting emission probability is the iteration of 1 → 2 splitting kernels, an independent emission pattern, which correctly describes the NLL strongly ordered regime, i.e. t i ≪ t i+1 .At NNLL one needs to account for unordered corrections, i.e. t i ∼ t i+1 , which are described by the full set of 1 → 3 splitting kernels [67][68][69].Such a correction will generate an extra term in Eqs.(2.3) which contains a product of three GFs (e.g. the splitting q → qgg will be proportional to G q G g G g ).At the same time one has to consider the virtual one-loop corrections to the 1 → 2 splitting kernels [70], which will have the same GFs structure as the NLL kernel (2.3).Finally, to avoid double counting with the O(α 2 s ) iteration of the NLL evolution, we must subtract the latter from the r.h.s. of Eqs.(2.3).
The aforementioned calculation can be consistently performed in the dimensional regularisation scheme in d = 4 − 2ϵ dimensions.However, in order to exploit the flexibility of Monte Carlo integration, we need to bring the integral equations into a form that is manifestly finite so that we can take ϵ → 0 at the integrand level.In order to make the cancellation of ϵ divergences manifest, we include a local subtraction term with the goal of regularising both the 1 → 3 real and 1 → 2 virtual corrections.The subtraction can be built directly from the 1 → 3 real corrections, by integrating them with the same 1 → 2 GFs structure present in the NLL kernel (2.3).In the reals, this effectively plays the role of a virtual correction obtained by unitarity.This procedure results in evolution equations that are manifestly finite in four space time dimensions.In the case of quark jets, the NNLL evolution equation takes the form In the above equation we defined the inclusive emission probability P q (z, θ) as5 where θ is the angle between the final state quark and gluon (set by the evolution time t ′ ).The inclusive emission probability P q (z) (and its gluonic counterpart) will be also central to the construction of an algorithmic solution to the NNLL problem, as we shall demonstrate in forthcoming work [1].The quantity K (1) is the ratio of the two-loop to the one-loop cusp anomalous dimension and b 0 is the first coefficient of the QCD beta function (2.12) We have further decomposed the splitting function P qq (z) into its soft part (z ≃ 1) and the hard-collinear left over, where the latter is given by The term proportional to b 0 ln 2 g(z) balances the g(z) dependence of the argument of the overall coupling encoded in dt ′ in Eq. (2.14), such that the quantity B q 2 (z) is independent of the choice of g(z) and matches the definition and calculation given in Ref. [66].The definition of P q (z, θ) implicitly relies upon a scheme to define the variables z and θ while integrating over a second emission (and adding the corresponding virtual corrections).This scheme must be IRC safe and effectively P q (z, θ) can be thought of as a next-to-leadingorder correction to a 1 → 2 collinear splitting.
The functional K finite q contains 4-dimensional terms which can be handled efficiently with Monte Carlo methods. 6These terms satisfy a unitarity condition K finite q u=1 = 0, they are at most NNLL (α n s L n−1 ) and contain structures of the G i G j G k type.Although in the most general case they enter in NNLL resummations, for event shapes and jet rates they only enter in the soft limit and are responsible for correlated and clustering corrections in the language of Refs.[32,33,64,71,72].On the other hand, for purely collinear problems, such as the dynamics of small-R jets [7], K finite q contributes as a whole.The precise definition of K finite q [G q , G g ] goes hand-in-hand with the definition of the function B q 2 , of which it specifies the scheme.Nevertheless, any physical prediction is independent of such a scheme, and for any IRC safe observables the scheme change can be performed directly in d = 4 dimensions.We report the expressions of such terms for quark and gluon jets in Appendix C, while in the body of this article we focus on the first line of Eq. (2.9) and its gluonic counterpart.
The function B q 2 (z) was calculated recently in Ref. [66] and will be reviewed in the next section.One can now define the Sudakov ∆ q (t) at NNLL from Eq. (2.9).We can set u = 1 and obtain ln ∆ q (t) = − (2.14) The goal of this article is to define and calculate the analogous quantity for gluon jets B g 2 (z), which paves the way to obtain the Sudakov for the gluon fragmentation and the complete set of integral equations which generalises jet calculus in the collinear limit to NNLL.
Finally, the function g(z) fixes the precise scale of the coupling that becomes relevant at NNLL.An important constraint on g(z) is that in the soft limit (z ≃ 1) it is fixed to [63] so that the coupling is evaluated at the relative transverse momentum of the branching.Beyond the soft limit, the form of g(z) and B 2 (z) are linked and emerge from an O(α 2 s ) calculation.For the sake of concreteness, in the rest of this article we take (2.16) This defines unambiguously the terms in B 2 (z) that are proportional to the QCD β(α s ) function at one loop (i.e.b 0 ).
3 B q 2 (z) in the quark case and the physical coupling scheme In this section we give a concise review of the results for quark jets in Ref. [66].This helps elucidate the main messages of the current work, in addition to showcasing some differences to the case of gluon jets.In the inclusive emission probability for quark jets (2.10), the intensity of radiation in the soft limit (z ≃ 1) can be encoded in a physical scheme of the strong coupling constant [63], in such a way that the emission probability is defined as where the r.h.s.amounts to using the first line of Eq. (2.10) and the CMW coupling is related to the MS one by the following relation The above scheme is customarily used in QCD resummations and is a crucial analytical ingredient for NLL parton shower algorithms.The extension of this picture to higher logarithmic accuracy is not unique, and depends on the definition of the kinematic variables z and θ.While generalisations of Eq. (3.2) to three loops (encoded in the coefficient K (2) ) in the soft limit (relevant for double-logarithmic resummations) were obtained in Refs.[64,65], much less is known about its generalisation to the hard-collinear limit.In resummation literature this information is encoded, at the integrated level, in the observable-dependent collinear anomalous dimension commonly known as B 2 , which nowadays has been calculated for several observables [64,[73][74][75].When used in resummations of double logarithmic observables such as event shapes or transverse momentum of a colour singlet, B 2 always takes the form (see e.g.Refs.[64,[73][74][75]) where γ (2) q is the endpoint of the singlet DGLAP splitting kernel at two loops (e.g.[5]) and the quantity X q v is observable dependent.
Eq. (2.14) however suggests that we can encode this observable dependence into the observable-independent inclusive emission probability P q (z, θ), which encapsulates a differential anomalous dimension B q 2 (z).The observable dependence of B 2 thus emerges entirely from the integration over z, and specifically from integrating the Sudakov (2.14) with the single-emission measurement function parameterised in terms of z and θ, i.e. 7   Θ(V (z, θ) − v) . (3.5) This offers a natural path to incorporate this information in Monte Carlo parton showers, and it constitutes a crucial step towards the development of NNLL algorithms.Recently, Ref. [66] presented a two loop calculation of the quantity B q 2 (z) in Eq. (2.14) for the quark case, where the corresponding 1 → 3 splitting kernels are integrated by fixing z and θ to the momentum fraction and angle of either the first emission, i.e. the one at larger angles for the C 2 F channel, or that of the radiated pair (for the The differential anomalous dimension B q 2 (z) is sufficient to derive B 2 for any rIRC safe global observable defined on quark jets (see also the related discussions in Section 3 of Ref. [66]).Specifically, for the definition of z and θ adopted in Ref. [66] (and reported in Appendix B) one obtains where An analogous integration taking into account the observable constraint (3.5) would produce the observable-specific constant X q v .Explicit examples of this will be considered in Sec. 8.

B g 2 in the gluon case: definitions and computational strategy
To calculate B g 2 (z), one needs an IRC safe definition of z and θ such that it projects the 1 → 3 phase space Φ 3 onto the underlying 1 → 2 kinematics Φ 2 .The definition of B g 2 (z) is then uniquely specified by a kinematical map that provides a definition of the longitudinal momentum fraction z and the angle θ of the first branching in the Φ 2 phase space.This also defines in a unique manner the inclusive 7 We note that some definitions of the B q 2 coefficient (and thus of z and θ) in the literature contain an extra contribution arising from single-logarithmic soft physics (see e.g.Ref. [64]).These terms do not contribute to B q 2 in our scheme for quark jets, and they emerge from the integration of the K finite q [Gq, Gg] contribution.We will return to this point when discussing the case of gluon jets. 8The definition of the radiated pair is ambiguous in the CF (CF − CA/2) colour channel due to the symmetry of the splitting kernel, however this ambiguity does not affect the form of B q 2 (z).
emission probability P(z, θ), where one integrates over a second emission while keeping the variables z and θ fixed.A property of the definition of z and θ is that this projection reproduces K (1) in the soft limit.The generalisation of the calculation of the quark case [66] to the gluonic case is non-trivial and it entails two conceptual subtleties: 1.In the case of quark jets it was possible to define z and θ based on the colour structure of the triple collinear matrix elements.Specifically, independent (i.e.C 2 F ) and correlated (i.e.C F C A , C F T R , and C F (C F − C A /2)) channels are separated by their colour structure.This allows one to choose z and θ such that in the abelian channel (C 2 F ) these correspond to the kinematics of the first emission in an angular ordered picture, while in the remaining channels these variables correspond to the kinematics of the parent emitter of the final state pair of real partons (see e.g.Fig. 7).In the gluonic case, and specifically for the C 2 A piece, this correspondence between independent and correlated emissions and different colour channels is no longer valid, in that independent and correlated radiation are mixed under the same colour structure.As a result, we need to formulate an IRC safe procedure to define z and θ in all gluonic channels.

2.
A second subtlety concerns the flavour structure and the divergences in the gluonic case, which is more involved than in the quark case as reflected in the evolution equation for the gluonic GF (2.3).Unlike in the quark case, the gluon GF encodes two different splitting kernels already at NLL.Of these, p gg (z) is such that in an iterated 1 → 2 splitting g → g a (z)g b (1 − z) followed by, e.g., g b → q q the gluon g a which does not branch can still bring a soft divergence (this time at z ≃ 0).This behaviour is absent in the quark case (as z → 0 would correspond to a soft quark configuration), and ensuring the cancellation of all IRC divergences in the gluon case requires a definition of the underlying 1 → 2 branching that is safe in this limit.
We give below a definition of z and θ that addresses the above features and can be used to derive B g 2 (z).The procedure follows a strategy inspired by the mMDT/Soft-Drop (β = 0) algorithm [56,57] which was designed to identify hard-collinear splittings in a jet.It can be summarised in the following steps: • For each colour channel, decompose the phase space into angular sectors, such that each sector contains exactly one collinear singularity.Evidently, the number of required sectors matches the number of singular collinear configurations in each colour channel.
• For each sector, we cluster the two partons, say i and j, which give rise to the collinear singularity as θ ij → 0.
• In each sector the angle θ is identified with the angle between the clustered pair (ij) and the remaining parton k.Similarly, the longitudinal momentum fraction z coincides with that of parton k, i.e. z ≡ z k .We then symmetrise the result by adding the permutation z → 1 − z and multiply by 1/2! if the Born 1 → 2 splitting to which we are projecting involves identical particles.The procedure follows the angular-ordered pattern of the generating functional evolution equation.
• In the region of phase space where z k → 0, the definition of the angle θ becomes unsafe.Therefore, when z k is below a certain threshold, say z cut ≪ 1, we discard the soft branch containing z k and fix the angle θ ≡ θ ij .In addition the longitudinal momentum is chosen as z ≡ z i /(z i +z j ) which must necessarily pass the z cut condition in order to define a hard-collinear splitting.
• Finally, we explicitly take the limit z cut → 0.
In practice we consider the fragmentation of a gluon jet defined as one hemisphere in the two-jet limit of the H → gg decay in the heavy-top-mass limit.The corresponding NNLL resummed doubly-differential distribution in z and θ is given in Eq. (D.13) of Appendix D. The extraction of B g 2 (z) then can proceed by simply relating the second order calculation for θ 2 σ 0 d 2 σ dθ 2 dz to the second order expansion of Eq. (D.13).In general our procedure is not exactly the same as mMDT9 since it does not correspond to the exact C/A declustering sequence.Yet, our clustering sequence is sufficient to capture all the divergent structure and obtain an IRC safe effective emission probability.
In contrast to the quark case, in which the integral of B q 2 (z) yields the simple form given in Eq. (3.6), the richer singularity structure present in the gluon case (specifically in the C 2 A channel) and reflected in the emergence of multiple singular sectors as outlined above will lead to the general form where [5] −γ (2) and X g θ 2 is a constant that we shall determine in Sec. 7. Finally, the extra term emerges from the non-trivial sectorisation that becomes necessary in the C 2 A channel (cf.Sec.7.4).
In the next sections we discuss how to use the above algorithm to define and calculate the inclusive emission probability P g (z, θ).In particular, in order to calculate P g (z, θ) we start by computing the doubly-differential distribution θ 2 d 2 σ/dθ 2 dz, and then in Sec.7 we will extract P g (z, θ).We will discuss separately the contribution of the one-loop correction to the collinear 1 → 2 splitting kernel and of the tree-level 1 → 3 splitting kernels.

Virtual corrections to 1 → 2 collinear splitting
We start by discussing the O(α 2 s ) real-virtual corrections to the θ 2 d 2 σ/dθ 2 dz distribution.We parameterise them as follows: where, in this case, the kinematical variables (θ, z) are uniquely defined as those of a 1 → 2 collinear splitting.Our results will be expressed in terms of a renormalised MS coupling related to the bare strong coupling by where Finally, we choose µ R = E (the energy of the hard parton initiating the fragmentation) to present the results, and define α s ≡ α s (E 2 ).There are three kinds of one-loop corrections that are involved in the calculation, and they can be separated according to their origin.
The first is the one-loop UV counter-term that is related to the renormalisation of the bare coupling (5.2), and reads: The second is the one-loop correction to the 1 → 2 splitting function itself.These corrections are taken from Ref. [70], in the CDR scheme. 10Explicitly, the one-loop correction to the collinear g → q q splitting is given by where (5.7) (5.8) 10 In the T 2 R n 2 f channel, we noticed a typographical error in Eq. (5.19) of the journal version of Ref. [70], which we fix to recover the divergent structure predicted by Catani's one-loop formula [76].Specifically, one should replace −3 with +3 in the denominator of the first term on the r.h.s. of Eq. (5.19) in the journal version of Ref. [70].
For the one-loop correction to the collinear g → gg splitting we have where (5.10) We finally note that if we were to compute the full O(α 2 s ) result for a given observable then the one-loop virtual corrections to the Born process would also have to be included.In practice, for the double-differential calculation at second order that we carry out here only the product of such one-loop virtual corrections with the single collinear splitting function is needed.For the H → gg process in the heavy-top effective theory considered in this article it reads (with In the above expression, the process dependence is embodied in the O(ϵ 0 ) terms of the above equation, and it cancels in the extraction of B g 2 (z) when taking the difference between our second-order calculation and the expansion of the resummation formula given in Eq. (D.13).The function V(θ 2 , z, ϵ) in Eq. (5.1) then reads g→q q(θ 2 , z, ϵ) + V (1)  g→gg (θ 2 , z, ϵ) + V (1) Born H→gg (θ 2 , z, ϵ) .(5.12) 6 Real 1 → 3 collinear splitting In this section we discuss the integration of the triple-collinear splitting functions over the three-body phase space, while keeping fixed z and θ as discussed in Sec. 4. We start by introducing the spin-averaged triple-collinear splitting functions relevant for gluon fragmentation.By considering the spin-averaged case, we implicitly focus on observables that are insensitive to spin correlations. 12For an initial gluon there are three distinct splitting functions to consider organised by colour structure.In the notation of Ref. [68] these are: • ⟨ Pg 1 g 2 g 3 ⟩ corresponding to a 1 → 3 gluon splitting into three gluons.
• ⟨ P (ab) g 1 q 2 q3 ⟩ corresponding to the abelian channel, C F T R n f , of a gluon splitting into a q q pair, with subsequent emission of a gluon.
• ⟨ P (nab) g 1 q 2 q3 ⟩ corresponding to the non-abelian channel, C A T R n f , of a gluon splitting into two gluons, one of which branches further into a q q pair.The splitting functions in Ref. [68] are expressed in terms of invariants, s ij , and energy fractions z i .In the relevant collinear limit, we have that ij where E denotes the energy of the initial gluon.The triple-collinear phase space in d = 4 − 2ϵ dimensions may be expressed in the form [77] where the Gram determinant is given by and As explained before, we are interested in the double-differential distribution in (z, θ), which is defined according to the procedure of Sec. 4 and will be illustrated in the following: This colour channel exhibits a simple collinear structure because the sole collinear singularity appears as θ 23 → 0, see Fig. 1.Therefore, we fix the energy fraction of the gluon to be z 1 = z and the angle θ = θ g shown in Fig. 1.We map the set (z 1 , z 2 , z 3 ) into an independent set (z, z p ) defined as shown in the figure.Instead when z falls below z cut , the angle will be defined as θ = θ 23 and the longitudinal momentum defined as z = z 2 as shown in Fig. 2.
The Feynman diagram representing gluon emission followed by its subsequent decay to a q q pair.
The Feynman diagram representing gluon emission followed by its subsequent decay to a q q pair, where g 1 goes below the energy threshold and is allowed to fly off at wide angles.
Figure 3: The Feynman diagram representing the q q emission, and the energy fraction parameterisation is suitable for the region θ 13 < θ 12 .
The configurations in which z < z cut (z > 1−z cut ) in sector 1 (2) correspond to powersuppressed contributions in z cut (due to a soft quark) and are therefore neglected.
• C 2 A : This channel is the most complicated since there is a collinear pole as any angle θ ij → 0. Therefore, we partition the phase space into three distinct sectors defined by the smallest angle in each.For example, in the sector min{θ ij } = θ 12 we fix the angle between the parent of the clustered pair and the remaining gluon, i.e. θ = θ 12,3 as shown in Fig. 4. The above sectoring removes a 1/3 combinatorial factor, and in order to remove the remaining twofold symmetry in the g → ggg splitting function we require that z p ≥ 1/2 in the parameterisation of Fig. 4. 13 The energy fraction being fixed is z 3 = z.When z falls below z cut , then θ = θ 12 and z = z 2 .
We integrate the splitting functions over the three-body phase space in d = 4 − 2ϵ dimensions to obtain real emission contributions.The integrals we carry out are generically of the form ) where, s 123 = s 12 + s 13 + s 23 and Θ cut (θ ij ) denotes the angular cuts due to sectoring that we described above.We evaluate the above integrals by expanding out the singular Figure 4: The Feynman diagram representing gluon emission followed by its subsequent decay to a gg pair, where the energy fraction parameterisation suitable for the angular region min{θ ij } = θ 12 is shown.
structure as a Laurent series in distributions, and then handle the remaining integrations with Mathematica.
We start with the C F T R n f channel.The splitting function to be integrated is the unpolarised quantity ⟨ Pg 1 q 2 q3 ⟩ from Ref. [68].The approach to the calculation of the corresponding phase-space integrals, as well as other computational details, are discussed in the corresponding calculation for the quark fragmentation given in Refs.[66,71].Following our definition of z and θ, we consider two sectors: θ 13 < θ 12 , which contains the collinear singularity along the antiquark, and θ 12 < θ 13 , which contains the collinear singularity along the quark.In the first sector we parameterise the momenta as if the gluon is emitted from the antiquark so that z 2 = (1 − z), z 3 = zz p and z 1 = z(1 − z p ) as shown in Fig. 3.In the second sector we instead parameterise the energy fractions as if the gluon is emitted from the quark, i.e. z 3 = z, z 1 = (1 − z p )(1 − z) and z 2 = z p (1 − z).We also fix the angle θ ≡ θ 2,13 in the first sector, and θ ≡ θ 3,12 in the second sector.We can just perform the calculation in the first sector and obtain the answer in the second sector by sending z → 1 − z.We report below the real emission contribution for the sector θ 13 < θ 12 which can be expressed as: where the functions are labeled according to the origin of the singular behaviour and H sec.1 fin.(z) is a finite function that we cast as a 1-fold integral over z p .As the gluon is emitted from the (anti)-quark, the singular structure is identical to that for the C 2 (z).jet case, and reported in Refs.[66,71], viz.
The result in the second sector can be obtained quite easily from the first by sending z → 1 − z.In particular, We plot the total finite function, H C F T R fin.
≡ H sec.1 fin.+H sec.2 fin. in Fig. 5.The divergent behavior of H C F T R fin.
(z) as z → 0 or z → 1 is only logarithmic, and thus it is integrable over z ∈ [0, 1].The integral read (z) = −3 . (6.8) We start with the configuration depicted in Fig. 1.For this channel, a convenient parametrisation of the three body phase space can be obtained in terms of web variables as discussed in Ref. [66].The only collinear singularity in this channel emerges from the collinear g → q q splitting, i.e. when θ 23 → 0. Therefore, following our procedure outlined in Sec. 4, there is a single sector and thus we can integrate inclusively over the angular phase space.In this channel we need to consider the two situations depicted in Figs. 1 and 2, in which the branching with a larger angle either passes or does not pass the z cut threshold.In the first case, the longitudinal momentum kept fixed is z 1 = z (with 1 − z cut > z > z cut ).We fix the angle θ to be that of the parent of the (23) pair, θ = θ 23,1 = θ g , which reads Finally, we symmetrise in z ↔ 1 − z and multiply by 1/2! to obtain: In the second case (see Fig. 2), the first gluon fails the z cut condition.Accordingly, z and θ are now fixed by the kinematics of the hard g → q q splitting.We note that, in contrast to the C F T R n f channel, the failure of the z cut condition is now associated with a soft gluon divergence.In the z 1 → 0 soft limit, the spin averaged splitting kernel factorises into a product of an eikonal factor (associated to the emission of g 1 ) and a 1 → 2 splitting function.In this configuration we obtain: This result is similar to what was obtained for quark jets in the C 2 F channel reported in Eq. ( 27) of Ref. [71], modulo the replacement p qq → p qg .The final result is obtained as the sum of the two contributions Eqs.(6.10) and (6.11).

The C 2
A channel As explained in Sec. 4, we work in the region min{θ ij } = θ 12 , see Fig. 4.This is completely general as the other sectors have an identical structure due to the three-fold symmetry present in the splitting kernel.To account for the full 1/3!symmetry factor, one can simply choose z p < 1 − z p .Clearly in this sector the only collinear singularity arises when θ 12 → 0. As in the C A T R n f channel, also in this case we have to consider two scenarios, i.e. when the gluon at larger angle either passes or fails the z cut condition.In the first scenario, the angular variable to be fixed is θ = θ 12,3 = θ g now defined as and z 3 = z (with 1 − z cut > z > z cut ).We obtain: where we have symmetrised over z ↔ 1 − z.The functions appearing in the above equation are given by where h poles pass is a numerical constant that reads The uncertainty in the above number is on the last quoted digit, which has be rounded as indicated by the bracket notation.We will use the same notation for all numerical constants quoted in the following, unless an error is explicitly quoted.As expected from the structure of the gluon splitting functions, the function H The coefficient of the soft divergences at z = 1 and z = 0 in the above expression is the result of a two-fold integration which can be performed with very high numerical precision.We obtain h fin pass ≃ 3.31674336 (8) .(6.17) The quantity G z 3 >zcut (z) is fully regular and can be integrated over z ∈ [0, 1] The function G z 3 >zcut (z) for this colour channel is displayed in Fig. 6.The decomposition in Eq. (6.16) will be important later in the extraction of the NNLL hard collinear coefficient B g 2 (z), which requires a consistent subtraction of the soft contributions.We now consider the second scenario in which z 3 < z cut .In the C 2 A channel, this contribution contains a physical difference from the corresponding result quoted in Eq. (6.11).The difference comes from the fact that in the z 3 → 0 limit, the matrix element does not  factorise into the product of independent emissions matrix elements, except for the limit in which the gluon pair that passes the cut is well separated in angle from the gluon that fails the cut, i.e. with the notation of Fig. 4, θ 13 ∼ θ 23 ≫ θ 12 .We can then organise the result as the sum of a term analogous to Eq. (6.11), given below A , (A) and a new contribution stemming from the non-independent (i.e.correlated) emission contribution.This is evaluated numerically and obtain where h pgg fail ≃ 0.731081807(5) , h δ fail ≃ −8.42858916( 9) .(6.21) 7 Extraction of B g 2 (z) Now we use the results obtained in the previous two sections to extract the function B g 2 (z).As we stressed before B g 2 (z) reflects the NNLL dynamics that does not arise from stronglyordered physics.We follow the procedure outlined in Sec. 4 to extract B g 2 (z), and after adding real and virtual corrections for each colour channel we need to remove pieces pertaining to NLL physics.In analogy with Eq. (2.10) for quark jets we write down our defining equation for B g 2 (z) starting from the NNLL evolution equation for G g (x, t) where K finite g [G q , G g ] is given in Appendix C and the inclusive emission probabilities for gluon jets are given by where the LO anomalous dimensions read The ratio of strong couplings in the second line of P gg (z, θ) has the role of restoring the correct scale of the coupling corresponding to the soft singularity as z → 0 as opposed to the one at z → 1 that is encoded in the evolution time (2.1).This feature is of course present exclusively for gluon jets, and it is absent in the quark case (2.10).The quantity B g 2 (z) is then simply the sum of B gg 2 (z) and B qg 2 (z).A subtle aspect of the above decomposition between terms that are interpreted as corrections to either the g → gg or the g → q q channel is that they both receive a contribution from the C A T R colour factor.The separation between such contributions to B g 2 (z) is not unique and the ambiguity is immaterial as one can decide to assign the whole correction to either of the two flavour channels.Conventionally we include it as a correction to g → gg.We thus write: The Sudakov form factor for gluon jets, given at NLL in Eq. (2.5), at NNLL then reads In the following we carry out the calculation of B g 2 (z) in each of the above colour channels.
7.1 The T 2 R n 2 f channel This channel is distinct in that it has no double real contribution.From Eqs. (5.4), (5.5), (5.11) and after subtracting the NLL contribution that emerges from Eq. (D.13) we obtain which integrates to The extraction of B g 2 (z) in this channel is very simple since it starts at NNLL and hence there is no need to subtract NLL contributions from the total result.Thus we add real, Eq. (6.5) and its mirror symmetric obtained by the swap z ↔ 1 − z, and virtual, Eq. (5.5), corrections to find (z) , whose integral reads We combine the real and virtual terms Eqs.(6.10), (6.11), (5.4), (5.5), (5.9) and (5.11) and then subtract the NLL contribution that emerges from Eq. (D.13), to obtain which integrates to Before we discuss the more involved C 2 A channel, it is instructive to compare the integral of the contributions to B g 2 (z) computed so far to the expectation given in Eq. (4.2).In particular, given that the term A contribution, the combination of Eqs.(7.7), (7.11) allows us to extract the constant X g θ 2 corresponding to the observable dθ 2 dz used in the calculation.This yields which can be used as a cross check in the C 2 A channel below.

The C 2
A channel The procedure for extracting B 2 (z) is the same as in the previous channel.The relevant real emission terms are given in Eqs.(6.13), (6.19), (6.20) which can then be combined with the virtual contributions emerging from Eqs. (5.4), (5.9), (5.11).After removal of NLL terms the result can be expressed as The term B g,C 2 A ,analytic 2 (z) is computed analytically while the contribution G z 3 >zcut (z) is determined numerically.Finally there is an endpoint contribution which originates purely from the double-soft limit.This includes the clustering correction as well as the numerically computed contributions from the region where z 3 < z cut .The term B g,C 2  A ,analytic 2 (z) reads: As a cross check of the above result we can compare to the expected integral in Eq. (4.2), where X g θ 2 is given in Eq. (7.12).For this check we still have to determine the function F clust. is known to stem from the double-soft limit (within the triple-collinear approximation).Therefore, due to Casimir scaling it is simple to relate the clustering correction to the known result for quark jets (see e.g.[71]) (cf.also the discussion in Appendix E).Replacing C F → C A in the quark jet result we get the following result the C 2 A component reads: where h C A clust.≃ −1.16363257( 4) , (7.19) and Cl 2 denotes the Clausen function.The overall factor of 1/2 is multiplied to obtain the result for a single leg.This gives in good agreement with our result.We conclude this section with a remark on the endpoint contribution B endpoint 2 (z) given in Eq. (7.16).As discussed above, this originates from double-soft configurations and is a consequence of the scheme used here to define z and θ.We observe that such endpoint terms mark an important difference to the case of quark jets (cf.Appendix B), where separate definitions of z and θ can be adopted for correlated and independent emission contributions.In the quark case these are separated by colour factors, while in the gluon case they are mixed together in the C 2 A channel.The same double-soft origin is shared by the clustering correction (7.18).Terms of double-soft origin are discussed further in Appendix E.

Moments of EEC and angularities in groomed jets at NNLL
In this section we use the calculations presented in this article to derive NNLL results for the moments of energy-energy correlation (EEC) and angularities measured on jets groomed according to the mMDT/Soft drop (β = 0) procedure [56,57].These classes of jet substructure observables have received widespread attention in the literature, with several applications both at hadron and lepton colliders [78][79][80][81][82][83][84][85][86][87][88][89][90][91].As a concrete example we consider the processes Z → q q and H → gg to analyse both quark and gluon jets.Specifically, we calculate, for the first time, the groomed fractional moments of EEC, which are defined as [54] where the sum runs over all particles within a given hemisphere, H ≡ H R or H ≡ H L , and E denotes the total energy in the hemisphere.We also consider the angularities [55] (the corresponding NNLL resummation in the ungroomed case is given in Refs.[64,[92][93][94][95]) defined w.r.t. the Winner-Take-All (WTA) axis as where once again the sum involves all particles in a given hemisphere.For both observables the parameter x is constrained by IR safety to be x < 2. Finally, we define the observables as follows: A NNLL calculation for WTA angularities for quark jets has been recently presented in Ref. [72], 14 and below we present new results for fractional moments of EEC for quark and gluon jets as well as for angularities measured on gluon jets.These results allow for a complete phenomenological analysis of moments of EEC and angularities on groomed jets at hadron and lepton colliders.The master formula for the NNLL cumulative cross section Σ(v) for a groomed observable v (that here denotes either a moment of EEC or an angularity) in the two processes considered here can be easily derived by applying the GFs method to Eq. (2.7).In the following we work in the limit v ≪ z cut ≪ 1, that is commonly considered for these types of groomed observables.This regime has the advantage that one can neglect power corrections in z cut and hence the evolution equations for the GFs can be solved analytically.Specifically, in this limit we can restrict ourselves to taking the soft limit of the K finite q [G q , G g ] and K finite g [G q , G g ] functions (cf.Eqs.(2.9), (7.1) and Appendix C).With a little abuse of notation we have used the same parameter z cut for the definition of groomed observables and in the calculation of B g 2 (z).However, we stress that (cf.Sec. 4) B g 2 (z) is defined strictly by taking the limit z cut → 0 while in the observables considered here one can opt to retain finite z cut effects to improve the accuracy of the calculation.The quantity B g 2 (z) we have derived applies to both cases with and without finite z cut effects, since these would be captured by the full (numerical) solution of the GFs equations (or equivalently by a NNLL accurate parton shower algorithm).
Using the evolution equation in Eqs.(2.9), (7.1), and following similar steps to those outlined in Appendix D, 15 we obtain the following NNLL results for quark and gluon jets, respectively with σ Z→q q,H→gg 0 denoting the leading-order cross sections.The quantity R v (v, z cut ) denotes the Sudakov radiator, which can be obtained by integrating the inclusive emission probability in Eqs.(2.10), (7.2) with the measurement function of the observables (8.1), (8.2) for a single collinear a → bc splitting, namely with where z is the energy fraction of b.This gives (see also Ref. [72]) (with f ∈ {q, g}) where λ v ≡ α s β 0 ln 1/v, λ zcut ≡ α s β 0 ln 1/z cut , α s ≡ α s (E 2 ), and The functions in Eq. (8.8) read (we use the notation ) where For the observables of interest B f 2,v emerges from the integration of Eqs.(2.10), (7.2) with the phase space constraint in Eq. (8.5) and it is given by for the fractional moments of EEC, and for the WTA angularities.The quantities B q 2,θ 2 and B g 2,θ 2 are defined in Eqs.(3.6) and (4.2) while X f v denotes the observable dependent constants for which we obtain and where the constants X q θ 2 and X g θ 2 are given in Eqs.(3.7), (7.12).Finally, h C A clust. is given in Eq. (7.19).We note that the result for B q 2,λx has been previously obtained in Ref. [72].The clustering corrections F clust.(v) in Eq. (8.4) arise from the finite terms K finite q [G q , G g ] and K finite g [G q , G g ] in the GF Eqs.(2.9), (7.1) (cf.Appendix C), which give rise to a correction to the Sudakov in Eq. (8.4) and are given by (see also Ref. [72]) where h A comment on the difference between quark and gluon jets is in order.This difference can be traced back to the projection (4.1) used in the definition of the inclusive emission probability and hence of B f 2 (z).As we discussed extensively in the paper, we adopt a different projection for quark and gluon jets due to the more involved structure of collinear singularities in the latter case.This translates into the final anomalous dimensions B f 2,v given in Eqs.(8.15), (8.17).Here we can see that the gluon jet result B g 2,v , unlike B q 2,v , includes a part of the clustering correction, specifically in the C 2  A channel (see discussion in Sec.7.4).
Finally, we comment on the coefficient functions C (8.4).These are matching coefficients between the collinear approximation of the generating functionals and the full O(α s ) calculation in the limit v ≪ z cut ≪ 1.For the processes and observables considered here they have the general structure [64] where X f v are the observable dependent constants of hard-collinear origin given in Eqs.(8.19), (8.20), and the remaining constants are the (only) process-dependent ingredients of the calculation.Specifically, H f (1) accounts for the hard-virtual corrections at one loop order The presence of X f v both in Eqs.(8.23) and in Eqs.(8.15), (8.17) indicates that the GF solution is defined in a resummation scheme in which the running of α s multiplying the constants of collinear origin is absorbed into the anomalous dimensions which define the inclusive emission probability.In alternative approaches to this class of resummations (e.g. that of Refs.[32,64,72]) the coupling multiplying these constants is evaluated at the observable-dependent collinear (low) scale, which would remove X f v from Eqs. (8.15), (8.17).The physical result at a given logarithmic order is of course resummation-scheme invariant (cf.also footnote 20 in Appendix D).
For quark jets, as a check of our results, we compare the resummed predictions for F C x and λ x against Event2 in Appendix F, and against the formalism of Refs.[32,64] adapted to groomed observables in Refs.[71,72].Specifically, as observed in Ref. [72], in the case of groomed quark jets the integrated quantities B q 2,λx and B q 2,F Cx can be extracted from Sec. 3 of Ref. [64], by combining the endpoint of the two loop DGLAP anomalous dimension γ f with the running of the one-loop hard-collinear constant C (1) hc and the recoil correction δF rec , 16 which agrees with our findings in Eqs.(8.15), (8.17).Similarly, an independent calculation of the clustering corrections for quark angularities can be found in Ref. [72].We also carried out analogous checks for the gluonic results, based on a straightforward extension of the above formalism to H → gg.In particular, we checked that the result for F C 0 , which coincides with the heavy hemisphere mass, agrees with the result for the latter of Ref. [78].We reiterate that, due to the process-independence of the collinear limit, our results can be used for hadron-collider jets after supplying the appropriate process-dependent analogues of the constants in Eqs.(8.23).

Conclusions and Outlook
In this paper we have presented a generating functional formulation of the NNLL resummation of collinear logarithms produced by multiple timelike collinear parton splittings.In particular, we have addressed one of the key elements of the generating functionals method, namely the Sudakov form factor, which appears in the formalism as the no-branching probability for a given fragmenting parton.We have pointed out that for general NNLL accuracy in the collinear limit, alongside correcting the real emission matrix-element using the triplecollinear splitting functions, the Sudakov form factor also needs to be augmented to the two-loop order.Here we focussed on gluon fragmentation, and have provided the necessary extension by computing the anomalous dimension B g 2 (z) which governs the intensity of collinear radiation off a gluon as a function of a suitably defined longitudinal momentum fraction z.We envisage that the results presented here will be instrumental to address the following classes of resummation problems: 1. Problems that do not admit a closed-form solution such as, for instance, the differential fragmentation of a jet into a system of microjets.This class of problems is sensitive to the whole non-linear structure of the generating functionals evolution equations provided in this article, which reflects the recursive nature of collinear fragmentation.As such, in the general case only a numerical solution of these equations can be envisioned, which can be achieved either using discretisation techniques or Monte Carlo methods.In particular, our results are necessary for the extension to NNLL of the NLL results of Ref. [7].We will address some of these aspects in a future publication [1] together with an algorithmic solution to the problem of collinear fragmentation.
2. The second class of collinear problems that can be addressed with the formalism outlined here is that of semi-inclusive observables that are not sensitive to the full non-linear structure of the evolution equations, which can now be solved with analytic methods [78-81, 84-86, 89-91].Examples of these problems, considered in this article, are groomed event shapes and moments of EEC.We derived new NNLL results for a family of groomed angularities defined w.r.t. the Winner-Take-All axis and fractional moments of energy-energy correlators, which open up a range of interesting phenomenological studies at hadron colliders such as the LHC.
3. Another interesting direction is the extension to more general observables that also exhibit sensitivity to soft physics, including the important cases of ungroomed global event shapes and non-global observables (in the latter case a GFs formulation is given in Refs.[28,29]).Although many of the above NNLL resummations (e.g. for global event shapes) are known from the literature, their formulation in the GFs language would pave the way for a single resummation tool capable of achieving NNLL accuracy across a wide class of observables.Similar considerations apply to the space-like collinear branching of initial-state partons, which will be a crucial future step.Another key consideration will be the inclusion of spin correlations at NNLL via the use of polarised splitting kernels in the GF evolution equations.
An important additional aspect of the formalism presented here is that it allows one to formulate collinear resummation in a language that resembles that of parton showers, hence offering insight on the design of future NNLL algorithms.For instance, the Sudakov form factor calculated in this article constitutes a crucial ingredient to reach this perturbative accuracy and it is therefore a key element of future NNLL parton showers.Finally, the results derived in this article are distributed in Mathematica format with the arXiv submission of this paper.
The above functions read where H fin. (z) is given by a 1-fold integral (cf. Figure 4 of ref. [66]), that is provided in Mathematica format as an ancillary file with the arXiv preprint of this article.The function B q 2 (z) is regular in the soft limit z → 1 and is thus fully integrable over z ∈ [0, 1].

C The NNLL K finite kernel
In this appendix we report the functions K finite f (with f ∈ {q, g}) entering the NNLL evolution equation for the quark generating functionals given in Eqs.(2.9), (7.1).We will start by presenting the result for quark jets due to its simpler structure, and later give the gluon counterpart.

C.1 Quark fragmentation
We can express K finite q as a difference between two terms which encode the (subtracted) real corrections to the first of Eqs.(2.3) and its double counting with the iteration of the NLL kernel, respectively.That is: The difference of K R q and K DC q ensures that the quantity K finite q is infrared finite and purely NNLL.
The diagram representing gluon decay to a q q pair, where the quark from the gluon decay is either identical or non-identical to the initiating quark.Following the definition of the inclusive emission probability (2.10) we obtain: where we have used the notation t i,j to indicate the value of the evolution time (2.1) corresponding to the angle θ i,j , and E p is set to the energy of the parent parton in the 1 → 3 branching.Moreover, we have parameterised the integrands as18 Here, P 1→3 is parameterised according to the phase space depicted in Fig. 7, while P (B) 1→3 is parameterised according to Fig. 8.The corresponding phase space measures, denoted by dΦ (A,B) 3 in (C.2) are obtained from Eq. (6.1) by setting ϵ = 0 and performing the change of variable (with the corresponding Jacobian) of Figs. 7, 8.The sum in Eq. (C.2) runs over all the (A) colour channels of the type q → f 1 f 2 q(q) defined in Eq. (C.3), and f i denotes the flavour of parton i.The factor 1/S 2 is the symmetry factor to be applied in the case of identical particles (notably S 2 = 2! in the C F C A and in the C F (C F − C A /2) channels).
The phase space integrals are now intended to be regulated by the usual cutoff procedure used so far, that implies an upper cut t 0 on evolution time (i.e. a lower cut on angles) and a cut on the energy fractions 1 − z 0 ≥ {z, z p } ≥ z 0 .All the calculations are then meant to be performed by taking the limit of t 0 → ∞ and z 0 → 0 after combining K R q and K DC q .We observe that each term in curly brackets in K R q is obtained by subtracting from each real configuration a differential counterterm defined by projecting the real kinematics into a given underlying Born phase space point.Crucially, this requires an infrared-and-collinearsafe definition of the branching variables θ and z, that uniquely specifies the definition of B q 2 (z).Therefore, the scheme that defines B q 2 (z) is tightly connected to the specific form of K finite q [G q , G g ] in Eq. (2.9), which specifically adopts the definition of z and θ of Ref. [66] outlined in Sec. 3.An alternative scheme will modify the expression of both B q 2 (z) and K finite q [G q , G g ] by modifying the precise definition of θ and z, but the prediction of the evolution equation is clearly invariant under any scheme change of this type.Following similar considerations for the double-counting term we find where the sum runs over the g → gg and g → q q splitting channels.Eq. (C.5) is obtained by calculating the contribution to Eq. (C.2) due to the iteration of the NLL kernel in Eq. (2.3).We notice that the difference between Eqs. (C.2) and (C.5) (i.e.K finite q [G q , G g ]) only contributes in the regime where all angles are commensurate t 2,3 ≃ t 1,23 ≃ t 1,3 .Conversely, K finite q [G q , G g ] vanishes in the strongly ordered limit.This implies that one can always approximate the angles in the Sudakov and in the GFs according to t 2,3 ≃ t 1,23 ≃ t 12,3 ≃ t 1,3 neglecting higher-order, N 3 LL terms.This property proves useful when performing a (semi-)analytic calculation using this formalism.

C.2 Gluon fragmentation
We now extend the result of the previous section to the case of gluon jets.The derivation proceeds through analogous arguments to the quark case, with two main differences.The first difference is related to the scheme for the anomalous dimension B g 2 (z) used in the main text for its computation.As we discussed at length in Sec. 4, the definition of B g 2 (z) relies on a specific projection (4.1) from the three-particle to the two-particle phase space, which is reflected in definition of the inclusive emission probability as well as in the definition of z and θ used in the calculation of K R [G q , G g ].The latter now reads: where (using the parameterisation of Figs. 1 and 3) Finally, for the C 2 A channel we must distinguish between the case in which z > z cut and z < z cut , following the scheme used to define B g 2 (z).However, since the latter is defined in the z cut → 0 limit, we can entirely neglect the region z < z cut in the calculation of Using the phase space parameterisation of Fig. 4 we then find where we defined The additional combinatorial factor 1/2! accounts for the remaining two-fold symmetry in the g → ggg splitting function, while the factor 1/3 in Eq. (C.9) can be removed by choosing one of the permutations.
The second important difference between the gluon and quark cases concerns the double counting term K DC [G q , G g ], which in the gluonic case is more involved due to the fact that two different splitting channels contribute to the NLL evolution equation for G g (cf.Eq. (2.3)).We write it as where (using the parameterisation of Figs. 1, 3 and 4, respectively) 19 C.12) 19 We note that the CA TR contribution to the double counting term is not necessary but it could be performed (i.e.adding the same integrand with z ↔ 1 − z and dividing by 2) to mirror the iteration of the NLL kernel.
We note that in Eqs.(C.12) and (C.14) the symmetrisation in z ↔ 1−z is not accompanied by a factor of 1/2! since this is already included in the definition of the leading-order splitting function P gg (z).

D Derivation of
In this appendix we derive the double-differential distribution in Eq. (D.13).It is instructive to start with the NLL case, for which the distribution θ 2 σ 0 d 2 σ dθ 2 dz is defined as (cf.Eq. (2.7)) where θ pass and z pass denote the angle and momentum fraction of the largest-angle branching with 1 − z cut > z pass > z cut .The probabilities dP (f ) m are computed with Eq. (2.2) using the NLL evolution equation (2.3).For the sake of establishing our argument, we consider the simpler case of the fragmentation of a quark.For a given number m of final state particles, the observable is defined by considering all sequential primary declusterings with respect to the hard leg that initiates the fragmentation in an angular ordered picture, and setting θ and z to the first one that passes the z cut (with z cut ≪ 1) condition.At NLL, each declustering amounts to a single primary emission off the initial hard leg, which then fragments inclusively.Therefore, we can ignore the secondary branchings of the gluons and approximate their generating functional with the first order expansion, that is E Double-soft end-point contributions to B g 2 (z) In this appendix we discuss how to devise an alternative scheme for B g,C 2 A 2 (z) in which the clustering and endpoint contributions, whose origin lies in the double-soft limit, completely disappear.The emergence of the clustering correction in Eq. (4.2), and likewise the endpoint contribution in Eq. (7.13), is due to the specific definition of z and θ in the kinematic projection that we used to define B g,C 2 A 2 (z).To make the physics clear it is useful to recapitulate the case of quark jets [66], in which the definition of B q 2 (z) does not involve such terms.The procedure given in Sec. 4 is based upon decomposing, for each colour channel, the three-particle phase space into angular sectors each containing exactly one collinear singularity.In the case of quark jets, this task is transparent.For correlated channels, there is a unique collinear singularity between the offspring of the gluon decay, and thus there is a single sector defined according to our procedure.For the C 2 F channel the only collinear singularities arise when emissions g 1 or g 2 are collinear to the quark q 3 , and thus there are two sectors in our procedure.These features lead to a natural definition of z and θ, in all colour channels, that is based on a naive clustering procedure in which pairs of partons which do not develop a collinear singularity never get clustered together [66].An equally valid, albeit tedious, procedure for defining z and θ would be to follow a strict C/A algorithm and define z as the momentum fraction between the two branches produced by the C/A declustering.This procedure would differ from the previous one in finite angular regions resulting in an alternative scheme for B q 2 (z), which would contain a clustering correction.
Let us now move to the case of gluon jets.As discussed in the main text, the C 2 A channel has a richer structure of collinear singularities which involves all three partons (gluons) on an equal footing.In particular, it is not possible to separate out correlated and independent emissions.In this case one cannot use the naive clustering argument adopted for quark jets and for this reason we had to resort to strict C/A procedure outlined in the article.This ultimately leads to the appearance of the clustering correction in Eq. (4.2).However, in a situation in which two of the three final state gluons are soft the parent gluon can be uniquely identified and one can separate the double-soft squared amplitude into the sum of correlated and independent emission terms which are identical to those in the quark case up to Casimir scaling.One could then introduce an alternative scheme in which the double-soft limit is treated as in the quark case, while the hard-collinear leftover is treated as explained in Sec. 4, i.e. using strict C/A declustering.Given that the clustering correction originates from the double-soft region of phase space, this alternative B g 2 (z) will directly integrate to A scheme of this type has the advantage of eliminating the end-point contribution present in Eq. (7.13), but presents a slight disadvantage in that it makes the calculation of K finite q more cumbersome due to the different projections adopted in the double-soft and in the hard-collinear limits.where the constant h sub.soft is the result of a 2-fold integration and reads h sub.soft ≃ −6.49651797 (5) . (E.8) The function H + G sub. z 3 >zcut (z) , (E.9) and the remainder function G sub. z 3 >zcut , see Fig. 9, is integrable in z ∈ [0, 1] and its integral reads  where the constant h pgg fail is given in Eq. (6.21).The most important feature of eq.(E.11) is that its z-dependence is fully regular over z ∈ [0, 1] unlike Eq. (6.20) which contains an end-point contribution.
Finally we compute the pure double-soft contributions, which can be read off from the analytic results of Ref. [66] by enforcing the appropriate limit.Starting with the correlated portion of the double-soft function, eq.(E.4), we find F Comparison to Event2 for quark-jet observables In this appendix, we present the comparison between the O(α 2 s ) expansions of our resummed results with the numerical fixed-order code Event2.In order to focus on the NNLL terms, and improve the precision on the pure O(α 2 s ) contribution, we define the difference between observable distributions (as done e.g. in Ref. [32] for ungroomed event shapes) If our resummation correctly captures all NNLL terms, then we would expect the quantity ∆ to tend to zero as ln 1/v grows large, in all colour channels.For groomed observables, however, the resummation we have carried out is valid in the regime v ≪ z cut ≪ 1 and hence neglects log-enhanced terms which are power-suppressed in z cut .Therefore, we expect to find agreement with Event2 in the limit of z cut → 0. Since taking this limit is numerically challenging, we show plots for three different values of z cut , showing convergence between our results and Event2 as this limit is approached.The plots shown in Fig. 10 are obtained with 3 × 10 12 events, although residual statistical instabilities remain in the Event2 results.We note that in all colour channels, the agreement with Event2 substantially improves with decreasing z cut , strongly validating our resummed predictions.

Figure 5 :
Figure 5: The function H C F T R fin.

C 2 Ag,C 2 A 2 ( 2 A
clust. .For the C 2A channel our procedure for the extraction of B z) agrees exactly with the mMDT/Soft drop grooming with the Cambridge/Aachen algorithm, for which F C

Figure 10 :
Figure 10: Validation of resummation against Event2 for quark jets, using three different values of z cut .See text for details.