aainstitutetext: Rudolf Peierls Centre for Theoretical Physics, Clarendon
Laboratory, Parks Road, University of Oxford, Oxford OX1 3PU,
UKbbinstitutetext: Lancaster-Manchester-Sheffield Consortium for Fundamental
Physics, Department of Physics & Astronomy, University of Manchester,
Manchester M13 9PL, United Kingdomccinstitutetext: CERN, Theoretical Physics
Department, CH-1211 Geneva 23, Switzerland

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

Melissa van Beekveld b,c Mrinal Dasgupta b Basem Kamal El-Menoufi a Jack
Helliwell c and Pier Francesco Monni<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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
${\mathcal{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.

††preprint: CERN-TH-2023-143, OUTP-23-08P

##  1 Introduction and motivation

In this paper and in a forthcoming article vanBeekveld:2023? , we initiate a
formulation of the methods of jet calculus Konishi:1979cb ; Bassetto:1983mvz ;
Dokshitzer:1991wu ; Ellis:1996mzs  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. Campbell:2022qmc  for an overview).

More specifically this work addresses the following points:

  1. 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 Dasgupta:2014yra ; Kang:2016mcy ,
fragmentation Jain:2011xz ; Chang:2013rca ; Ritzmann:2014mka ; Elder:2017bkd ;
Chen:2020adz ; Neill:2020bwv ; Chen:2021gdk ; Neill:2021std ; Chen:2022jhb ;
Chen:2022muj ; Chen:2023zlx  or non-global observables Dasgupta:2001sh ;
Banfi:2002hw ; Hatta:2013iba ; Caron-Huot:2015bja ; Larkoski:2015zka ;
Neill:2016stq ; AngelesMartinez:2018cfz ; Hatta:2020wre ; Banfi:2021xzn ;
Banfi:2021owj ; Becher:2023vrh ) 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 Banfi:2001bz ;
Banfi:2014sua ; Banfi:2016zlc ). 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 Konishi:1979cb ; Bassetto:1983mvz ;
Dokshitzer:1991wu ; Ellis:1996mzs , 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. 2.

The recent development of new techniques to construct parton- shower
algorithms with demonstrable higher logarithmic accuracy (see e.g. Refs.
Li:2016yez ; Hoche:2017hno ; Dasgupta:2018nvj ; Dulat:2018vuy ; Bewick:2019rbu
; Forshaw:2019ver ; Dasgupta:2020fwr ; Forshaw:2020wrq ; Hamilton:2020rcu ;
Nagy:2020rmk ; Nagy:2020dvz ; Hamilton:2021dyz ; Karlberg:2021kwr ;
Herren:2022jej ; vanBeekveld:2022zhl ; vanBeekveld:2022ukn ; Hamilton:2023dwb
; vanBeekveld:2023lfu ; Assi:2023rbu ; Ravasio:2023anp ) 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 Dasgupta:2020fwr
; Hamilton:2020rcu ; Hamilton:2021dyz ; Karlberg:2021kwr ; vanBeekveld:2022zhl
; vanBeekveld:2022ukn ; Hamilton:2023dwb ; vanBeekveld:2023lfu ;
Ravasio:2023anp  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.111A related description
of the evolution of soft radiation away from the collinear limit using the
same mathematical language was discussed in Refs. Banfi:2021xzn ;
Banfi:2021owj  in the context of non-global resummations in the large$-N_{c}$
limit.

Here we calculate the analytical ingredients needed in the formulation of NNLL
resummation within the GFs method. As we will demonstrate in forthcoming work
vanBeekveld:2023? , these ingredients can then also be used in Monte Carlo
algorithms to correctly include virtual corrections,222Specifically, 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.
which are not captured by the shower’s unitary approximation. For this reason,
ingredients of this type (such as the soft physical coupling scheme
Catani:1990rr ; Banfi:2018mcq ; Catani:2019rvy ) 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.
Dasgupta:2021hbh , 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 $\mathcal{B}_{2}(z)$. In Secs. 4-7 we calculate
$\mathcal{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 Banfi:2004yd  and angularities Berger:2003iw  measured on groomed
jets with the mMDT/Soft-Drop procedure Dasgupta:2013ihk ; Larkoski:2014wba .
Finally, 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 $\mathcal{B}_{2}(z)$ anomalous dimensions
for quarks and gluons are given in electronic format with the arXiv submission
of this article.

##  2 Generating functionals for collinear fragmentation

In the present 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.

###  2.1 Review of NLL resummation

We consider the collinear fragmentation of a parton of initial energy $E$
resolved at an angular resolution $\theta\ll 1$. The all-order resummation at
NLL resums single logarithms $L=-\ln\theta$ of the form $\alpha_{s}^{n}L^{n}$.
At this order, the fragmentation can be formulated as a branching process in
which emissions are strongly ordered in angle.333Angular 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. At NNLL
we aim at resumming corrections of order $\alpha_{s}^{n}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

$t_{i}=\int_{\theta_{i}^{2}}^{1}\frac{d\theta^{2}}{\theta^{2}}\,\frac{\alpha_{s}(E^{2}g^{2}(z)\,\theta^{2})}{2\pi}\,,$
(1)

where $\alpha_{s}$ is the $\overline{\rm 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$ (see for example Ref. Dasgupta:2014yra ). In what follows, the use of
the $\beta(\alpha_{s})$ function that drives the running of $\alpha_{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\in\\{q,g\\}$ and momentum fraction $x$ 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\in\\{q,g\\}$ and momentum fraction $x$
starting at an evolution time $t$ is given by Konishi:1979cb ;
Dokshitzer:1991wu

$\int dP^{(f)}_{m}=\left.\frac{1}{m!}\frac{\delta^{m}}{\delta
u^{m}}\,G_{f}(x,t)\right|_{\\{u\\}=0}.$ (2)

The quantity $u\equiv 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) 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 Dasgupta:2014yra 444An equivalent differential
form can be easily obtained by dividing $G_{f}$ by $\Delta_{f}$ and
subsequently taking the $t$ derivative. Note that Ref. Dasgupta:2014yra
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.

$\displaystyle G_{q}(x,t)$ $\displaystyle=u\,\Delta_{q}(t)+\int_{t}^{t_{0}}\
dt^{\prime}\int_{z_{0}}^{1-z_{0}}dz\,P_{qq}(z)\,G_{q}(x\,z,t^{\prime})\,G_{g}(x\,(1-z),t^{\prime})\frac{\Delta_{q}(t)}{\Delta_{q}(t^{\prime})}\,,$
$\displaystyle G_{g}(x,t)$
$\displaystyle=u\,\Delta_{g}(t)+\int_{t}^{t_{0}}dt^{\prime}\int_{z_{0}}^{1-z_{0}}dz\,\bigg{[}P_{gg}(z)\,G_{g}(x\,z,t^{\prime})\,G_{g}(x\,(1-z),t^{\prime})$
$\displaystyle+P_{qg}(z)\,G_{q}(x\,z,t^{\prime})\,G_{q}(x\,(1-z),t^{\prime})\bigg{]}\frac{\Delta_{g}(t)}{\Delta_{g}(t^{\prime})}\,,$
(3)

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 $\Delta_{f}$ describe the no-emission probability and
can be derived by imposing the unitarity condition

$\left.G_{f}(x,t)\right|_{u=1}=1\,,$ (4)

from which one obtains

$\displaystyle\ln\Delta_{q}(t)$
$\displaystyle=-\int_{t}^{t_{0}}dt^{\prime}\int_{z_{0}}^{1-z_{0}}dz\,P_{qq}(z)\,,$
(5) $\displaystyle\ln\Delta_{g}(t)$
$\displaystyle=-\int_{t}^{t_{0}}dt^{\prime}\int_{z_{0}}^{1-z_{0}}dz\left(P_{gg}(z)+P_{qg}(z)\right)\,.$
(6)

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\sigma^{(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

$\displaystyle d\sigma^{(f)}=\sigma_{0}\,C(\alpha_{s})\otimes
J^{(f)}(\alpha_{s},v)\,,$ (7)

where $\sigma_{0}$ is the Born cross section for the hard process under
consideration. The jet distribution $J^{(f)}$ is simply obtained by
integrating Eq. (2) with the measurement function of the observable for any
final-state multiplicity $m$, that is

$J^{(f)}(\alpha_{s},v)=\sum_{m=1}^{\infty}\int\,dP^{(f)}_{m}\,\delta(v-V(\\{k\\}_{m}))\,.$
(8)

The $\otimes$ 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(\alpha_{s})$
admits a fixed-order perturbative expansion in powers of the strong coupling
constant $C(\alpha_{s})=1+{\cal O}(\alpha_{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\to 0$.
Specifically, at NNLL $C(\alpha_{s})$ is required at the one-loop order, which
entails the difference between the full ${\cal O}(\alpha_{s})$ QCD calculation
in the logarithmic limit (i.e. $v\to 0$) for the observable $v$ and the
expansion of the jet function at the same order.

###  2.2 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.1) to ${\cal
O}(\alpha_{s}^{2})$. To understand what the calculation entails, we observe
that Eqs. (2.1) describe the fragmentation by a sequence of angular ordered
branchings. The resulting emission probability is the iteration of $1\to 2$
splitting kernels, an independent emission pattern, which correctly describes
the NLL strongly ordered regime, i.e. $t_{i}\ll t_{i+1}$. At NNLL one needs to
account for unordered corrections, i.e. $t_{i}\sim t_{i+1}$, which are
described by the full set of $1\to 3$ splitting kernels Campbell:1997hg ;
Catani:1998nv ; Braun-White:2022rtg . Such a correction will generate an extra
term in Eqs. (2.1) which contains a product of three GFs (e.g. the splitting
$q\to 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\to 2$ splitting
kernels Sborlini:2013jba , which will have the same GFs structure as the NLL
kernel (2.1). Finally, to avoid double counting with the ${\cal
O}(\alpha_{s}^{2})$ iteration of the NLL evolution, we must subtract the
latter from the r.h.s. of Eqs. (2.1).

The aforementioned calculation can be consistently performed in the
dimensional regularisation scheme in $d=4-2\epsilon$ 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 $\epsilon\to 0$ at the integrand level. In order to make the cancellation
of $\epsilon$ divergences manifest, we include a local subtraction term with
the goal of regularising both the $1\to 3$ real and $1\to 2$ virtual
corrections. The subtraction can be built directly from the $1\to 3$ real
corrections, by integrating them with the same $1\to 2$ GFs structure present
in the NLL kernel (2.1). 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

$\displaystyle G_{q}(x,t)=u\,\Delta_{q}(t)\,+$
$\displaystyle\int_{t}^{t_{0}}dt^{\prime}\int_{z_{0}}^{1-z_{0}}dz\,\,G_{q}(x\,z,t^{\prime})\,G_{g}(x\,(1-z),t^{\prime})\frac{\Delta_{q}(t)}{\Delta_{q}(t^{\prime})}\,{\mathcal{P}}_{q}(z,\theta)$
(9) $\displaystyle+{\mathbb{K}}_{q}^{\rm finite}[G_{q},G_{g}]\,.$

In the above equation we defined the inclusive emission probability
${\mathcal{P}}_{q}(z,\theta)$ as 555The precise form of $g(z)$ in the
coefficient of ${\mathcal{B}}^{q}_{2}(z)$ is not relevant at NNLL, but we keep
it in for simplicity.

$\displaystyle{\mathcal{P}}_{q}(z,\theta)\equiv$
$\displaystyle\frac{2\,C_{F}}{1-z}\left(1+\frac{\alpha_{s}(E^{2}g^{2}(z)\theta^{2})}{2\pi}K^{(1)}\right)$
$\displaystyle+{\mathcal{B}}_{1}^{q}(z)+\frac{\alpha_{s}(E^{2}g^{2}(z)\theta^{2})}{2\pi}\left({\mathcal{B}}^{q}_{2}(z)+{\mathcal{B}}_{1}^{q}(z)b_{0}\ln
g^{2}(z)\right)\,,$ (10)

where $\theta$ is the angle between the final state quark and gluon (set by
the evolution time $t^{\prime}$). The inclusive emission probability
${\mathcal{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 vanBeekveld:2023? . The quantity $K^{(1)}$ is
the ratio of the two-loop to the one-loop cusp anomalous dimension

$K^{(1)}=\left(\frac{67}{18}-\frac{\pi^{2}}{6}\right)\,C_{A}-\frac{10}{9}\,T_{R}\,n_{f}\equiv
K^{(1),C_{A}}\,C_{A}+K^{(1),n_{f}}\,T_{R}\,n_{f}\,,$ (11)

and $b_{0}$ is the first coefficient of the QCD beta function

$\displaystyle b_{0}=\frac{11}{6}\,C_{A}-\frac{2}{3}\,T_{R}\,n_{f}\equiv
b_{0}^{(C_{A})}\,C_{A}+b_{0}^{(n_{f})}\,T_{R}\,n_{f}\,.$ (12)

We have further decomposed the splitting function $P_{qq}(z)$ into its soft
part ($z\simeq 1$) and the hard-collinear left over, where the latter is given
by

${\mathcal{B}}_{1}^{q}(z)\equiv-C_{F}\,(1+z)\,.$ (13)

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^{\prime}$ in Eq. (14),
such that the quantity ${\mathcal{B}}_{2}^{q}(z)$ is independent of the choice
of $g(z)$ and matches the definition and calculation given in Ref.
Dasgupta:2021hbh . The definition of ${\mathcal{P}}_{q}(z,\theta)$ implicitly
relies upon a scheme to define the variables $z$ and $\theta$ while
integrating over a second emission (and adding the corresponding virtual
corrections). This scheme must be IRC safe and effectively
${\mathcal{P}}_{q}(z,\theta)$ can be thought of as a next-to-leading-order
correction to a $1\to 2$ collinear splitting.

The functional ${\mathbb{K}}_{q}^{\rm finite}$ contains 4-dimensional terms
which can be handled efficiently with Monte Carlo methods.666The application
of this formalism to collinear unsafe quantities, such as fragmentation
functions, entails considerable conceptual subtleties vanBeekveld:2023? .
These terms satisfy a unitarity condition $\left.{\mathbb{K}}_{q}^{\rm
finite}\right|_{u=1}=0$, they are at most NNLL ($\alpha_{s}^{n}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. Banfi:2014sua ; Banfi:2016zlc
; Banfi:2018mcq ; Anderle:2020mxj ; Dasgupta:2022fim . On the other hand, for
purely collinear problems, such as the dynamics of small-$R$ jets
Dasgupta:2014yra , ${\mathbb{K}}_{q}^{\rm finite}$ contributes as a whole. The
precise definition of ${\mathbb{K}}_{q}^{\rm finite}[G_{q},G_{g}]$ goes hand-
in-hand with the definition of the function ${\mathcal{B}}_{2}^{q}$, 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. (9) and its gluonic counterpart.

The function ${\mathcal{B}}^{q}_{2}(z)$ was calculated recently in Ref.
Dasgupta:2021hbh  and will be reviewed in the next section. One can now define
the Sudakov $\Delta_{q}(t)$ at NNLL from Eq. (9). We can set $u=1$ and obtain

$\displaystyle\ln\Delta_{q}(t)$
$\displaystyle=-\int_{t}^{t_{0}}dt^{\prime}\int_{z_{0}}^{1-z_{0}}dz\,{\mathcal{P}}_{q}(z,\theta)\,.$
(14)

The goal of this article is to define and calculate the analogous quantity for
gluon jets ${\mathcal{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\simeq 1$) it is fixed to Catani:1990rr

$\lim_{z\to 1}\frac{g(z)}{1-z}=1\,,$ (15)

so that the coupling is evaluated at the relative transverse momentum of the
branching. Beyond the soft limit, the form of $g(z)$ and
${\mathcal{B}}_{2}(z)$ are linked and emerge from an ${\cal
O}(\alpha_{s}^{2})$ calculation. For the sake of concreteness, in the rest of
this article we take

$g(z)=1-z\,.$ (16)

This defines unambiguously the terms in ${\mathcal{B}}_{2}(z)$ that are
proportional to the QCD $\beta(\alpha_{s})$ function at one loop (i.e.
$b_{0}$).

##  3 $\mathcal{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.
Dasgupta:2021hbh . 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.2), the intensity of
radiation in the soft limit ($z\simeq 1$) can be encoded in a physical scheme
of the strong coupling constant Catani:1990rr , in such a way that the
emission probability is defined as

$\displaystyle\left.dt\,dz\,{\mathcal{P}}_{q}(z,\theta)\right|_{\rm
soft}\simeq
C_{F}\frac{d\theta^{2}}{\theta^{2}}dz\,\frac{2}{1-z}\,\frac{\alpha_{s}^{\text{CMW}}(E^{2}(1-z)^{2}\theta^{2})}{2\pi}\
,$ (17)

where the r.h.s. amounts to using the first line of Eq. (2.2) and the CMW
coupling is related to the $\overline{\rm MS}$ one by the following relation

$\displaystyle\alpha_{s}^{\text{CMW}}(\mu^{2})=\alpha_{s}(\mu^{2})\left(1+\frac{\alpha_{s}(\mu^{2})}{2\pi}K^{(1)}\right)\,.$
(18)

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 $\theta$.
While generalisations of Eq. (18) to three loops (encoded in the coefficient
$K^{(2)}$) in the soft limit (relevant for double-logarithmic resummations)
were obtained in Refs. Banfi:2018mcq ; Catani:2019rvy , 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 Davies:1984hs ;
deFlorian:2001zd ; deFlorian:2004mp ; Banfi:2018mcq . 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. Davies:1984hs ; deFlorian:2001zd ; deFlorian:2004mp ; Banfi:2018mcq
)

$\displaystyle B_{2}^{q}=-\gamma_{q}^{(2)}+b_{0}X\ ,$ (19)

where $\gamma_{q}^{(2)}$ is the endpoint of the singlet DGLAP splitting kernel
at two loops (e.g. Ellis:1996mzs )

$\displaystyle\gamma_{q}^{(2)}=C_{F}^{2}\left(\frac{3}{8}-\frac{\pi^{2}}{2}+6\zeta_{3}\right)+C_{F}C_{A}\left(\frac{17}{24}+\frac{11\pi^{2}}{18}-3\zeta_{3}\right)-T_{R}n_{f}\left(\frac{1}{6}+\frac{2\pi^{2}}{9}\right)\,,$
(20)

and the quantity $X$ is observable dependent.

Eq. (14) however suggests that we can encode this observable dependence into
the observable-independent inclusive emission probability
${\mathcal{P}}_{q}(z,\theta)$, which encapsulates a differential anomalous
dimension ${\cal B}^{q}_{2}(z)$. The observable dependence of $B_{2}$ thus
emerges entirely from the integration over $z$, and specifically from
integrating the Sudakov (14) with the single-emission measurement function
parameterised in terms of $z$ and $\theta$, i.e.777We note that some
definitions of the $B^{q}_{2}$ coefficient (and thus of $z$ and $\theta$) in
the literature contain an extra contribution arising from single-logarithmic
soft physics (see e.g. Ref. Banfi:2018mcq ). These terms do not contribute to
$B^{q}_{2}$ in our scheme for quark jets, and they emerge from the integration
of the ${\mathbb{K}}_{q}^{\rm finite}[G_{q},G_{g}]$ contribution. We will
return to this point when discussing the case of gluon jets.

$\Theta(V(z,\theta)-v)\,.$ (21)

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. Dasgupta:2021hbh  presented a two loop
calculation of the quantity ${\mathcal{B}}^{q}_{2}(z)$ in Eq. (14) for the
quark case, where the corresponding $1\to 3$ splitting kernels are integrated
by fixing $z$ and $\theta$ to the momentum fraction and angle of either the
first emission, i.e. the one at larger angles for the $C_{F}^{2}$ channel, or
that of the radiated pair (for the $C_{F}C_{A}$, $T_{R}n_{f}$ and
$C_{F}\,(C_{F}-C_{A}/2)$ channels).888The definition of the radiated pair is
ambiguous in the $C_{F}\,(C_{F}-C_{A}/2)$ colour channel due to the symmetry
of the splitting kernel, however this ambiguity does not affect the form of
${\mathcal{B}}^{q}_{2}(z)$. Its expression is given in Appendix B.

The differential anomalous dimension ${\mathcal{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. Dasgupta:2021hbh ).
Specifically, for the definition of $z$ and $\theta$ adopted in Ref.
Dasgupta:2021hbh  (and reported in Appendix B) one obtains

$\displaystyle\int_{0}^{1}dz\,{\cal B}^{q}_{2}(z)\equiv
B_{2,\theta^{2}}^{q}=-\gamma_{q}^{(2)}+b_{0}X^{q}_{\theta^{2}}\ ,$ (22)

where

$\displaystyle
X^{q}_{\theta^{2}}=C_{F}\left(\frac{2\pi^{2}}{3}-\frac{13}{2}\right)\,.$ (23)

An analogous integration taking into account the observable constraint (21)
would produce the observable-specific constant $X$. Explicit examples of this
will be considered in Sec. 8.

##  4 ${\mathcal{B}}_{2}^{g}$ in the gluon case: definitions and computational
strategy

To calculate ${\mathcal{B}}^{g}_{2}(z)$, one needs an IRC safe definition of
$z$ and $\theta$ such that it projects the $1\to 3$ phase space $\Phi_{3}$
onto the underlying $1\to 2$ kinematics $\Phi_{2}$. The definition of
${\mathcal{B}}^{g}_{2}(z)$ is then uniquely specified by a kinematical map

${\cal M}:\Phi_{3}\to\Phi_{2}\,,$ (24)

that provides a definition of the longitudinal momentum fraction $z$ and the
angle $\theta$ of the first branching in the $\Phi_{2}$ phase space. This also
defines in a unique manner the inclusive emission probability
${\mathcal{P}}(z,\theta)$, where one integrates over a second emission while
keeping the variables $z$ and $\theta$ fixed. A property of the definition of
$z$ and $\theta$ is that this projection reproduces $K^{(1)}$ in the soft
limit. The generalisation of the calculation of the quark case
Dasgupta:2021hbh  to the gluonic case is non-trivial and it entails two
conceptual subtleties:

  1. 1.

In the case of quark jets it was possible to define $z$ and $\theta$ based on
the colour structure of the triple collinear matrix elements. Specifically,
independent (i.e. $C_{F}^{2}$) 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 $\theta$ such that in the
abelian channel ($C_{F}^{2}$) 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_{A}^{2}$ 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 $\theta$ in all gluonic channels.

  2. 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.1). 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\to 2$ splitting $g\to
g_{a}(z)g_{b}(1-z)$ followed by, e.g., $g_{b}\to q\bar{q}$ the gluon $g_{a}$
which does not branch can still bring a soft divergence (this time at $z\simeq
0$). This behaviour is absent in the quark case (as $z\to 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\to 2$
branching that is safe in this limit.

We give below a definition of $z$ and $\theta$ that addresses the above
features and can be used to derive ${\mathcal{B}}^{g}_{2}(z)$. The procedure
follows a strategy inspired by the mMDT/Soft-Drop $(\beta=0)$ algorithm
Dasgupta:2013ihk ; Larkoski:2014wba  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 $\theta_{ij}\to 0$.

  * •

In each sector the angle $\theta$ 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\equiv z_{k}$. We then symmetrise the result by adding the permutation $z\to
1-z$ and multiply by $1/2!$ if the Born $1\to 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}\to 0$, the definition of the angle
$\theta$ becomes unsafe. Therefore, when $z_{k}$ is below a certain threshold,
say $z_{\text{cut}}\ll 1$, we discard the soft branch containing $z_{k}$ and
fix the angle $\theta\equiv\theta_{ij}$. In addition the longitudinal momentum
is chosen as $z\equiv z_{i}/(z_{i}+z_{j})$ which must necessarily pass the
$z_{\text{cut}}$ condition in order to define a hard-collinear splitting.

  * •

Finally, we explicitly take the limit $z_{\rm cut}\to 0$.

In practice we consider the fragmentation of a gluon jet defined as one
hemisphere in the two-jet limit of the $H\to gg$ decay in the heavy-top-mass
limit. The corresponding NNLL resummed doubly-differential distribution in $z$
and $\theta$ is given in Eq. (138) of Appendix D. The extraction of
${\mathcal{B}}^{g}_{2}(z)$ then can proceed by simply relating the second
order calculation for
$\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma}{d\theta^{2}dz}$ to the second
order expansion of Eq. (138). In general our procedure is not exactly the same
as mMDT 999An exception is given by the $C_{A}^{2}$ channel, where our
procedure is equivalent to the mMDT algorithm. 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 ${\cal B}^{q}_{2}(z)$
yields the simple form given in Eq. (22), the richer singularity structure
present in the gluon case (specifically in the $C_{A}^{2}$ channel) and
reflected in the emergence of multiple singular sectors as outlined above will
lead to the general form

$\int_{0}^{1}{\cal B}^{g}_{2}(z)\equiv
B_{2,\theta^{2}}^{g}=-\gamma_{g}^{(2)}+b_{0}X^{g}_{\theta^{2}}+\mathcal{F}_{\text{clust.}}^{C_{A}^{2}}\,,$
(25)

where Ellis:1996mzs

$-\gamma_{g}^{(2)}=\frac{4}{3}C_{A}T_{R}n_{f}+C_{F}T_{R}n_{f}-C_{A}^{2}\left(\frac{8}{3}+3\zeta_{3}\right),$
(26)

and $X^{g}_{\theta^{2}}$ is a constant that we shall determine in Sec. 7.
Finally, the extra term $\mathcal{F}_{\text{clust.}}^{C_{A}^{2}}$ emerges from
the non-trivial sectorisation that becomes necessary in the $C_{A}^{2}$
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 $\mathcal{P}_{g}(z,\theta)$. In
particular, in order to calculate $\mathcal{P}_{g}(z,\theta)$ we start by
computing the doubly-differential distribution
$\theta^{2}d^{2}\sigma/d\theta^{2}dz$, and then in Sec. 7 we will extract
$\mathcal{P}_{g}(z,\theta)$. We will discuss separately the contribution of
the one-loop correction to the collinear $1\to 2$ splitting kernel and of the
tree-level $1\to 3$ splitting kernels.

##  5 Virtual corrections to $1\to 2$ collinear splitting

We start by discussing the ${\cal O}(\alpha_{s}^{2})$ real-virtual corrections
to the $\theta^{2}d^{2}\sigma/d\theta^{2}dz$ distribution. We parameterise
them as follows:

$\displaystyle\left(\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma_{\mathcal{V}}^{(2)}}{d\theta^{2}\,dz}\right)=\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\mathcal{V}(\theta^{2},z,\epsilon)\
,$ (27)

where, in this case, the kinematical variables $(\theta,z)$ are uniquely
defined as those of a $1\to 2$ collinear splitting. Our results will be
expressed in terms of a renormalised $\overline{\text{MS}}$ coupling related
to the bare strong coupling by

$\mu^{2\epsilon}\alpha_{s}=S^{-1}_{\epsilon}\mu_{R}^{2\epsilon}\alpha_{s}(\mu_{R}^{2})\left(1-\frac{b_{0}}{\epsilon}\frac{\alpha_{s}(\mu_{R})}{2\pi}+\mathcal{O}(\alpha_{s}^{2})\right)\,,$
(28)

where

$S_{\epsilon}=(4\pi)^{\epsilon}e^{-\epsilon\gamma_{E}}\,.$ (29)

Finally, we choose $\mu_{R}=E$ (the energy of the hard parton initiating the
fragmentation) to present the results, and define
$\alpha_{s}\equiv\alpha_{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 (28), and reads:

$\mathcal{V}_{\text{renorm.}}(\theta^{2},z,\epsilon)=-\frac{b_{0}}{\epsilon}\left(z(1-z)\theta\right)^{-2\epsilon}\left(C_{A}\,p_{gg}(z)+T_{R}n_{f}\,p_{qg}(z,\epsilon)\right)\,.\\\
$ (30)

The second is the one-loop correction to the $1\to 2$ splitting function
itself. These corrections are taken from Ref. Sborlini:2013jba , in the CDR
scheme.101010In the $T_{R}^{2}n_{f}^{2}$ channel, we noticed a typographical
error in Eq. (5.19) of the journal version of Ref. Sborlini:2013jba , which we
fix to recover the divergent structure predicted by Catani’s one-loop formula
Catani:1998bh . 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. Sborlini:2013jba . Explicitly, the one-loop correction to the
collinear $g\to q\bar{q}$ splitting is given by

$\mathcal{V}^{(1)}_{g\to
q\bar{q}}(\theta^{2},z,\epsilon)=(z(1-z))^{-3\epsilon}\theta^{-4\epsilon}\,T_{R}n_{f}\,p_{qg}(z,\epsilon)\left(\frac{1}{\epsilon^{2}}-\frac{2\pi^{2}}{3}\right)\times\\\
\times\bigg{[}T_{R}n_{f}\,f_{T_{R}n_{f}}^{g\to
q\bar{q}}+C_{F}\,f_{C_{F}}^{g\to q\bar{q}}+C_{A}\,f_{C_{A}}^{g\to
q\bar{q}}\bigg{]}\ ,$ (31)

where

$\displaystyle f_{T_{R}n_{f}}^{g\to q\bar{q}}$
$\displaystyle=-\frac{4\epsilon}{3}-\frac{20\epsilon^{2}}{9}+\mathcal{O}\left(\epsilon^{3}\right)\
,$ (32) $\displaystyle f_{C_{F}}^{g\to q\bar{q}}$
$\displaystyle=-2-3\epsilon-8\epsilon^{2}+\mathcal{O}\left(\epsilon^{3}\right)\
,$ (33) $\displaystyle f_{C_{A}}^{g\to q\bar{q}}$
$\displaystyle=1+\frac{11\epsilon}{3}+\frac{76\epsilon^{2}}{9}+\epsilon\ln(z(1-z))+\epsilon^{2}\,\text{Li}_{2}\left(\frac{z-1}{z}\right)+\epsilon^{2}\,\text{Li}_{2}\left(\frac{z}{z-1}\right)+\mathcal{O}(\epsilon^{3})\,.$
(34)

For the one-loop correction to the collinear $g\to gg$ splitting we have

$\displaystyle\mathcal{V}^{(1)}_{g\to
gg}(\theta^{2},z,\epsilon)=(z(1-z))^{-3\epsilon}\theta^{-4\epsilon}C^{2}_{A}\,p_{gg}(z,\epsilon)\left(\frac{1}{\epsilon^{2}}-\frac{2\pi^{2}}{3}\right)\,f_{C_{A}}^{g\to
gg}+\frac{1}{6}C_{A}^{2}-\frac{1}{3}C_{A}T_{R}n_{f}\ ,$ (35)

where

$\displaystyle f_{C_{A}}^{g\to
gg}=-1+\epsilon\ln(z(1-z))+\epsilon^{2}\,\text{Li}_{2}\left(\frac{z-1}{z}\right)+\epsilon^{2}\,\text{Li}_{2}\left(\frac{z}{z-1}\right)+\mathcal{O}\left(\epsilon^{3}\right)\,.$
(36)

We finally note that if we were to compute the full
$\mathcal{O}(\alpha_{s}^{2})$ 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\to gg$ process in
the heavy-top effective theory considered in this article it reads (with
$\mu_{R}=E$)

$\displaystyle\mathcal{V}^{(1)\,{\rm Born}}_{H\to gg}(\theta^{2},z,\epsilon)$
$\displaystyle=\left(C_{A}\,p_{gg}(z)+T_{R}\,n_{f}\,p_{qg}(z,\epsilon)\right)\,\left(1-\frac{\pi^{2}}{12}\epsilon^{2}\right)(z^{2}(1-z)^{2}\theta^{2})^{-\epsilon}$
$\displaystyle\times\left(\left(-\frac{4^{-\epsilon}}{\epsilon^{2}}-\frac{11}{6\epsilon}+\frac{7}{12}\pi^{2}\right)\,C_{A}+\frac{2}{3\epsilon}T_{R}\,n_{f}\right)\,.$
(37)

In the above expression, the process dependence is embodied in the ${\cal
O}(\epsilon^{0})$ terms of the above equation, and it cancels in the
extraction of ${\cal B}^{g}_{2}(z)$ when taking the difference between our
second-order calculation and the expansion of the resummation formula given in
Eq. (138). The function $\mathcal{V}(\theta^{2},z,\epsilon)$ in Eq. (27) then
reads

$\mathcal{V}(\theta^{2},z,\epsilon)=\mathcal{V}_{\text{renorm.}}(\theta^{2},z,\epsilon)+\mathcal{V}^{(1)}_{g\to
q\bar{q}}(\theta^{2},z,\epsilon)+\mathcal{V}^{(1)}_{g\to
gg}(\theta^{2},z,\epsilon)+\mathcal{V}^{(1)\,{\rm Born}}_{H\to
gg}(\theta^{2},z,\epsilon)\,.$ (38)

##  6 Real $1\to 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
$\theta$ 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. For an initial gluon there are three
distinct splitting functions to consider organised by colour structure. In the
notation of Ref. Catani:1998nv  these are:

  * •

$\langle\hat{P}_{g_{1}g_{2}g_{3}}\rangle$ corresponding to a $1\to 3$ gluon
splitting into three gluons.

  * •

$\langle\hat{P}^{\textrm{(ab)}}_{g_{1}q_{2}\bar{q}_{3}}\rangle$ corresponding
to the abelian channel, $C_{F}T_{R}\,n_{f}$, of a gluon splitting into a
$q\bar{q}$ pair, with subsequent emission of a gluon.

  * •

$\langle\hat{P}^{\textrm{(nab)}}_{g_{1}q_{2}\bar{q}_{3}}\rangle$ 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\bar{q}$ pair.

The splitting functions in Ref. Catani:1998nv  are expressed in terms of
invariants, $s_{ij}$, and energy fractions $z_{i}$. In the relevant collinear
limit, we have that $s_{ij}=(p_{i}+p_{j})^{2}\simeq
E^{2}z_{i}z_{j}\theta_{ij}^{2}$ where $E$ denotes the energy of the initial
gluon. The triple-collinear phase space in $d=4-2\epsilon$ dimensions may be
expressed in the form Gehrmann-DeRidder:1997fom

$\displaystyle
d\Phi_{3}=\frac{1}{\pi}\frac{E^{4-4\epsilon}}{(4\pi)^{4-2\epsilon}\Gamma(1-2\epsilon)}dz_{2}dz_{3}d\theta_{13}^{2}d\theta_{23}^{2}d\theta_{12}^{2}(z_{1}z_{2}z_{3})^{1-2\epsilon}\Delta^{-1/2-\epsilon}\,\Theta(\Delta)\
\ ,$ (39)

where the Gram determinant is given by

$\displaystyle\Delta=4\theta_{ik}^{2}\theta_{jk}^{2}-\left(\theta_{ij}^{2}-\theta_{ik}^{2}-\theta_{jk}^{2}\right)^{2},\quad
i\neq j\neq k\,,$ (40)

and

$\sum_{i=1}^{3}z_{i}=1\,.$ (41)

As explained before, we are interested in the double-differential distribution
in $(z,\theta)$, which is defined according to the procedure of Sec. 4 and
will be illustrated in the following:

  * •

$C_{A}T_{R}\,n_{f}$: This colour channel exhibits a simple collinear structure
because the sole collinear singularity appears as $\theta_{23}\to 0$, see Fig.
1. Therefore, we fix the energy fraction of the gluon to be $z_{1}=z$ and the
angle $\theta=\theta_{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_{\text{cut}}$, the angle will be
defined as $\theta=\theta_{23}$ and the longitudinal momentum defined as
$z=z_{2}$ as shown in Fig. 2.

$\theta_{g}=\theta$$z_{2}=(1-z)z_{p}$$\theta_{23}$$z_{3}=(1-z)(1-z_{p})$$z_{1}=z>z_{\text{cut}}$
Figure 1: The Feynman diagram representing gluon emission followed by its
subsequent decay to a $q\bar{q}$ pair.
$z_{1}<z_{\text{cut}}$$z=z_{2}$$\theta_{23}=\theta$$1-z=z_{3}$ Figure 2: The
Feynman diagram representing gluon emission followed by its subsequent decay
to a $q\bar{q}$ pair, where $g_{1}$ goes below the energy threshold and is
allowed to fly off at wide angles.

  * •

$C_{F}T_{R}\,n_{f}$: This colour channel exhibits two collinear singularities
as $\theta_{13}\to 0$ and $\theta_{12}\to 0$, see Fig. 3. Therefore, we
partition the phase space into two sectors $\theta_{13}<\theta_{12}$ and
$\theta_{12}<\theta_{13}$. In the former, we fix the quark energy $z_{2}=1-z$
and the angle $\theta=\theta_{2,13}$. In the latter we fix the anti-quark
energy $z_{3}=z$ and the angle $\theta=\theta_{3,12}$. The configurations in
which $z<z_{\text{cut}}$ ($z>1-z_{\text{cut}}$) in sector 1 (2) correspond to
power-suppressed contributions in $z_{\text{cut}}$ (due to a soft quark) and
are therefore neglected.

$z_{2}=1-z$$\theta_{2,13}$$z_{3}=z\,z_{p}$$\theta_{13}$$z_{1}=z(1-z_{p})$
Figure 3: The Feynman diagram representing the $q\bar{q}$ emission, and the
energy fraction parameterisation is suitable for the region
$\theta_{13}<\theta_{12}$.

  * •

$C_{A}^{2}$: This channel is the most complicated since there is a collinear
pole as any angle $\theta_{ij}\to 0$. Therefore, we partition the phase space
into three distinct sectors defined by the smallest angle in each. For
example, in the sector $\text{min}\\{\theta_{ij}\\}=\theta_{12}$ we fix the
angle between the parent of the clustered pair and the remaining gluon, i.e.
$\theta=\theta_{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\to ggg$ splitting function we require that $z_{p}\geq 1/2$ in the
parameterisation of Fig. 4.111111We observe that any choice of the angle that
is kept fixed to $\theta$ which is equivalent in the singular limits will
produce the same result for ${\cal B}^{g}_{2}(z)$, e.g. one could also choose
to fix $\theta=\theta_{13}$, and would only differ by further subleading
(N3LL) corrections. The energy fraction being fixed is $z_{3}=z$. When $z$
falls below $z_{\text{cut}}$, then $\theta=\theta_{12}$ and $z=z_{2}$.

$z_{3}=z$$\theta_{12,3}$$z_{1}=(1-z)z_{p}$$\theta_{12}$$z_{2}=(1-z)(1-z_{p})$
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 $\text{min}\\{\theta_{ij}\\}=\theta_{12}$ is
shown.

We integrate the splitting functions over the three-body phase space in
$d=4-2\epsilon$ dimensions to obtain real emission contributions. The
integrals we carry out are generically of the form

$\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma}{dz\,d\theta^{2}}=\int\text{d}\Phi_{3}(z_{i},\theta_{ij})\frac{\left(8\pi\alpha_{s}\mu^{2\epsilon}\right)^{2}}{s^{2}_{123}}\,\langle\hat{P}\rangle\,\theta^{2}\,\delta\left(\theta^{2}-\theta^{2}\left(z_{i},\theta_{ij}\right)\right)\delta(z-z(z_{i}))\,\Theta_{\text{cut}}(\theta_{ij})\,,$
(42)

where $\Theta_{\text{cut}}(\theta_{ij})$ denotes the angular cuts due to
sectoring that we described above. We evaluate the above integrals by
expanding out the singular structure as a Laurent series in distributions, and
then handle the remaining integrations with Mathematica.

###  6.1 The $C_{F}T_{R}\,n_{f}$ channel

We start with the $C_{F}T_{R}\,n_{f}$ channel. The splitting function to be
integrated is the unpolarised quantity
$\langle\hat{P}_{g_{1}q_{2}\bar{q}_{3}}\rangle$ from Ref. Catani:1998nv . 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. Anderle:2020mxj ;
Dasgupta:2021hbh .

Following our definition of $z$ and $\theta$, we consider two sectors:
$\theta_{13}<\theta_{12}$, which contains the collinear singularity along the
antiquark, and $\theta_{12}<\theta_{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 $\theta\equiv\theta_{2,13}$ in the first sector, and
$\theta\equiv\theta_{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\to 1-z$. We report below the real emission contribution for the
sector $\theta_{13}<\theta_{12}$ which can be expressed as:

$\left(\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma_{\mathcal{R}}^{(2)}}{d\theta^{2}\,dz}\right)^{C_{F}T_{R}n_{f},\mathrm{sec.1}}=C_{F}T_{R}n_{f}\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\bigg{(}\frac{H_{\text{soft-
coll.}}^{\mathrm{sec.1}}(\theta^{2},z,\epsilon)}{\epsilon^{2}}+\frac{H_{\mathrm{coll.}}^{\mathrm{sec.1}}(\theta^{2},z,\epsilon)}{\epsilon}\\\
+\frac{H_{\mathrm{soft}}^{\mathrm{sec.1}}(\theta^{2},z,\epsilon)}{\epsilon}+H^{\mathrm{sec.1}}_{\text{fin.}}(z)\bigg{)}\,,$
(43)

where the functions are labeled according to the origin of the singular
behaviour and $H^{\mathrm{sec.1}}_{\text{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_{F}^{2}$
quark jet case, and reported in Refs. Anderle:2020mxj ; Dasgupta:2021hbh ,
viz.

$\begin{split}H^{\mathrm{sec.1}}_{\text{soft-
coll.}}(\theta^{2},z,\epsilon)&=z^{-4\epsilon}(1-z)^{-2\epsilon}\theta^{-4\epsilon}\,p_{qg}(z,\epsilon)\left(1-\frac{\pi^{2}}{6}\epsilon^{2}+\mathcal{O}\left(\epsilon^{3}\right)\right)\
,\\\
H_{\text{coll.}}^{\mathrm{sec.1}}(\theta^{2},z,\epsilon)&=z^{-4\epsilon}(1-z)^{-2\epsilon}\theta^{-4\epsilon}\,p_{qg}(z,\epsilon)\left(\frac{3}{2}+\frac{13}{2}\epsilon-\frac{2\pi^{2}}{3}\epsilon+\mathcal{O}\left(\epsilon^{2}\right)\right)\
,\\\ H_{\text{soft}}^{\mathrm{sec.1}}(\theta^{2},z,\epsilon)&=0\ .\end{split}$
(44)

The result in the second sector can be obtained quite easily from the first by
sending $z\to 1-z$. In particular,

$\displaystyle
H^{\mathrm{sec.2}}_{\text{fin.}}(z)=H^{\mathrm{sec.1}}_{\text{fin.}}(1-z)\ .$
(45)

We plot the total finite function, $H^{C_{F}T_{R}}_{\text{fin.}}\equiv
H^{\mathrm{sec.1}}_{\text{fin.}}+H^{\mathrm{sec.2}}_{\text{fin.}}$ in Fig. 5.
The divergent behavior of $H^{C_{F}T_{R}}_{\text{fin.}}(z)$ as $z\to 0$ or
$z\to 1$ is only logarithmic, and thus it is integrable over $z\in[0,1]$. The
integral read

$\displaystyle\int_{0}^{1}dz\,H^{C_{F}T_{R}}_{\text{fin.}}(z)=-3\ .$ (46)

Figure 5: The function $H^{C_{F}T_{R}}_{\text{fin.}}(z)$.

###  6.2 The $C_{A}T_{R}\,n_{f}$ channel

We start with the configuration depicted in Fig. 1. For this channel channel,
a convenient parametrisation of the three body phase space can be obtained in
terms of web variables as discussed in Ref. Dasgupta:2021hbh . The only
collinear singularity in this channel emerges from the collinear $g\to
q\bar{q}$ splitting, i.e. when $\theta_{23}\to 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_{\text{cut}}$
threshold. In the first case, the longitudinal momentum kept fixed is
$z_{1}=z$ (with $1-z_{\text{cut}}>z>z_{\text{cut}}$). We fix the angle
$\theta$ to be that of the parent of the $(23)$ pair,
$\theta=\theta_{23,1}=\theta_{g}$, which reads

$\displaystyle\theta_{g}^{2}=\frac{z_{2}}{z_{2}+z_{3}}\theta^{2}_{12}+\frac{z_{3}}{z_{2}+z_{3}}\theta^{2}_{13}-\frac{z_{2}z_{3}}{(z_{2}+z_{3})^{2}}\theta^{2}_{23}\
.$ (47)

Finally, we symmetrise in $z\leftrightarrow 1-z$ and multiply by $1/2!$ to
obtain:

$\displaystyle\left(\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma_{\mathcal{R}}^{(2)}}{d\theta^{2}\,dz}\right)_{z_{1}>z_{\text{cut}}}^{C_{A}T_{R}}$
$\displaystyle=\frac{1}{2!}C_{A}T_{R}n_{f}\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\,z^{-3\epsilon}(1-z)^{-3\epsilon}\,\theta^{-4\epsilon}\bigg{[}p_{gg}(z)\left(-\frac{8}{3\epsilon}-\frac{40}{9}\right)$
$\displaystyle-\frac{4}{3}(1+z)\,\ln
z-\frac{4}{3}(2-z)\,\ln(1-z)+\frac{26}{9}\left(z^{2}+(1-z)^{2}-\frac{1}{z(1-z)}\right)+\frac{10}{3}\bigg{]}$
$\displaystyle\times\Theta(z-z_{\text{cut}})\,\Theta(1-z_{\text{cut}}-z)\ .$
(48)

In the second case (see Fig. 2), the first gluon fails the $z_{\text{cut}}$
condition. Accordingly, $z$ and $\theta$ are now fixed by the kinematics of
the hard $g\to q\bar{q}$ splitting. We note that, in contrast to the
$C_{F}T_{R}\,n_{f}$ channel, the failure of the $z_{\text{cut}}$ condition is
now associated with a soft gluon divergence. In the $z_{1}\to 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\to 2$ splitting
function. In this configuration we obtain:

$\left(\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma_{\mathcal{R}}^{(2)}}{d\theta^{2}\,dz}\right)_{z_{1}<z_{\text{cut}}}^{C_{A}T_{R}}=\left(\frac{\alpha_{s}}{2\pi}\right)^{2}C_{A}T_{R}n_{f}\left(2\,z(1-z)\right)^{-2\epsilon}\theta^{-2\epsilon}p_{qg}(z,\epsilon)\times\\\
\times\left(-\frac{z_{\mathrm{cut}}^{-2\epsilon}}{\epsilon}\ln\frac{4}{\theta^{2}}-\frac{\pi^{2}}{6}-\frac{1}{2}\ln^{2}\frac{4}{\theta^{2}}\right)\Theta(z-z_{\text{cut}})\,\Theta(1-z_{\text{cut}}-z)\,.$
(49)

This result is similar to what was obtained for quark jets in the $C_{F}^{2}$
channel reported in Eq. (27) of Ref. Anderle:2020mxj , modulo the replacement
$p_{qq}\to p_{qg}$. The final result is obtained as the sum of the two
contributions Eqs. (6.2) and (49).

###  6.3 The $C_{A}^{2}$ channel

As explained in Sec. 4, we work in the region
$\textrm{min}\\{\theta_{ij}\\}=\theta_{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 $\theta_{12}\to 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_{\text{cut}}$ condition. In the first scenario, the angular variable to be
fixed is $\theta=\theta_{12,3}=\theta_{g}$ now defined as

$\displaystyle\theta_{g}^{2}=\frac{z_{1}}{z_{1}+z_{2}}\theta^{2}_{13}+\frac{z_{2}}{z_{1}+z_{2}}\theta^{2}_{23}-\frac{z_{1}z_{2}}{(z_{1}+z_{2})^{2}}\theta^{2}_{12}\
,$ (50)

and $z_{3}=z$ (with $1-z_{\text{cut}}>z>z_{\text{cut}}$). We obtain:

$\left(\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma_{\mathcal{R}}^{(2)}}{d\theta^{2}\,dz}\right)_{z_{3}>z_{\text{cut}}}^{C_{A}^{2}}=C_{A}^{2}\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\bigg{(}\frac{H_{\text{soft-
coll.}}(\theta^{2},z,\epsilon)}{\epsilon^{2}}+\frac{H_{\mathrm{coll.}}(\theta^{2},z,\epsilon)}{\epsilon}\\\
+\frac{H_{\mathrm{soft}}(\theta^{2},z,\epsilon)}{\epsilon}+H^{C_{A}^{2}}_{\text{fin.}}(z)\bigg{)}\Theta(z-z_{\text{cut}})\,\Theta(1-z_{\text{cut}}-z)\,,$
(51)

where we have symmetrised over $z\leftrightarrow 1-z$. The functions appearing
in the above equation are given by

$\displaystyle H_{\text{soft-coll.}}(\theta^{2},z,\epsilon)$
$\displaystyle=\,z^{-2\epsilon}(1-z)^{-2\epsilon}\,\left(z^{-2\epsilon}+(1-z)^{-2\epsilon}\right)\theta^{-4\epsilon}4^{2\epsilon}\,p_{gg}(z)\left(1-\frac{\pi^{2}}{6}\epsilon^{2}\right)\,,$
$\displaystyle H_{\text{coll.}}(\theta^{2},z,\epsilon)$
$\displaystyle=-z^{-2\epsilon}(1-z)^{-2\epsilon}\left(z^{-2\epsilon}+(1-z)^{-2\epsilon}\right)\theta^{-4\epsilon}4^{\epsilon}\,p_{gg}(z)$
$\displaystyle\times\left(-\frac{11}{6}+2\ln
2+\epsilon\left(-\frac{67}{9}+\frac{2\pi^{2}}{3}+2\ln^{2}2\right)\right)\,,$
(52) $\displaystyle H_{\text{soft}}(\theta^{2},z,\epsilon)$
$\displaystyle=-z^{-2\epsilon}(1-z)^{-2\epsilon}\left(z^{-2\epsilon}+(1-z)^{-2\epsilon}\right)\theta^{-4\epsilon}4^{\epsilon}p_{gg}(z)\left(2\ln
2+h^{\rm poles}_{\rm pass}\,\epsilon\right)\,,$

where $h^{\rm poles}_{\rm pass}$ is a numerical constant that reads

$h^{\rm poles}_{\rm pass}\simeq 1.32644693(2)\,.$ (53)

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^{C_{A}^{2}}_{\text{fin.}}(z)$ has soft divergences
as $z\to 0$ and $z\to 1$ which can be singled out, leading to

$\displaystyle H^{C_{A}^{2}}_{\text{fin.}}(z)=-\frac{h_{\rm pass}^{\rm
fin}}{z(1-z)}+G_{z_{3}>z_{\text{cut}}}(z)\,.$ (54)

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_{\rm pass}^{\rm fin}\simeq 3.31674336(8)\,.$ (55)

The quantity $G_{z_{3}>z_{\text{cut}}}(z)$ is fully regular and can be
integrated over $z\in[0,1]$

$\displaystyle\int_{0}^{1}dz\,G_{z_{3}>z_{\text{cut}}}(z)\simeq 16.947\pm
0.004\,.$ (56)

The function $G_{z_{3}>z_{\text{cut}}}(z)$ for this colour channel is
displayed in Fig. 6.

Figure 6: The function $G_{z_{3}>z_{\text{cut}}}(z)$.

The decomposition in Eq. (54) will be important later in the extraction of the
NNLL hard collinear coefficient ${\cal B}_{2}^{g}(z)$, which requires a
consistent subtraction of the soft contributions.

We now consider the second scenario in which $z_{3}<z_{\text{cut}}$. In the
$C_{A}^{2}$ channel, this contribution contains a physical difference from the
corresponding result quoted in Eq. (49). The difference comes from the fact
that in the $z_{3}\to 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,
$\theta_{13}\sim\theta_{23}\gg\theta_{12}$. We can then organise the result as
the sum of a term analogous to Eq. (49), given below

$\left(\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma_{\mathcal{R}}^{(2)}}{d\theta^{2}\,dz}\right)_{z_{3}<z_{\text{cut}}}^{C_{A}^{2},\,(A)}=\left(\frac{\alpha_{s}}{2\pi}\right)^{2}C_{A}^{2}\left(2\,z(1-z)\right)^{-2\epsilon}\theta^{-2\epsilon}\,p_{gg}(z)\times\\\
\times\left(-\frac{z_{\mathrm{cut}}^{-2\epsilon}}{\epsilon}\ln\frac{4}{\theta^{2}}-\frac{\pi^{2}}{6}-\frac{1}{2}\ln^{2}\frac{4}{\theta^{2}}\right)\Theta(z-z_{\text{cut}})\,\Theta(1-z_{\text{cut}}-z)\,,$
(57)

and a new contribution stemming from the non-independent (i.e. correlated)
emission contribution. This is evaluated numerically and obtain

$\left(\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma_{\mathcal{R}}^{(2)}}{d\theta^{2}\,dz}\right)_{z_{3}<z_{\text{cut}}}^{C_{A}^{2},\,(B)}=h^{p_{gg}}_{\rm
fail}\,p_{gg}(z)\,\Theta(z-z_{\text{cut}})\,\Theta(1-z_{\text{cut}}-z)+h^{\delta}_{\rm
fail}\,\delta(z)\,,$ (58)

where

$h^{p_{gg}}_{\rm fail}\simeq 0.731081807(5)\,,\quad h^{\delta}_{\rm
fail}\simeq-8.42858916(9)\,.$ (59)

##  7 Extraction of $\mathcal{B}_{2}^{g}(z)$

Now we use the results obtained in the previous two sections to extract the
function $\mathcal{B}_{2}^{g}(z)$. As we stressed before
$\mathcal{B}_{2}^{g}(z)$ reflects the NNLL dynamics that does not arise from
strongly-ordered physics. We follow the procedure outlined in Sec. 4 to
extract $\mathcal{B}_{2}^{g}(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.2) for quark jets we write down our defining
equation for $\mathcal{B}_{2}^{g}(z)$ starting from the NNLL evolution
equation for $G_{g}(x,t)$

$\displaystyle G_{g}(x,t)$
$\displaystyle=u\,\Delta_{g}(t)+\int_{t}^{t_{0}}dt^{\prime}\int_{z_{0}}^{1-z_{0}}dz\,\bigg{[}{\mathcal{P}}_{gg}(z,\theta)\,G_{g}(x\,z,t^{\prime})\,G_{g}(x\,(1-z),t^{\prime})$
$\displaystyle+{\mathcal{P}}_{qg}(z,\theta)\,G_{q}(x\,z,t^{\prime})\,G_{q}(x\,(1-z),t^{\prime})\bigg{]}\frac{\Delta_{g}(t)}{\Delta_{g}(t^{\prime})}+{\mathbb{K}}_{g}^{\rm
finite}[G_{q},G_{g}]\,,$ (60)

where ${\mathbb{K}}_{g}^{\rm finite}[G_{q},G_{g}]$ is given in Appendix C and
the inclusive emission probabilities for gluon jets are given by

$\displaystyle{\mathcal{P}}_{gg}(z,\theta)$
$\displaystyle\equiv\,\frac{C_{A}}{1-z}\left(1+\frac{\alpha_{s}(E^{2}(1-z)^{2}\theta^{2})}{2\pi}K^{(1)}\right)$
$\displaystyle+\frac{\alpha_{s}(E^{2}z^{2}\theta^{2})}{\alpha_{s}(E^{2}(1-z)^{2}\theta^{2})}\,\frac{C_{A}}{z}\left(1+\frac{\alpha_{s}(E^{2}z^{2}\theta^{2})}{2\pi}K^{(1)}\right)$
$\displaystyle+{\mathcal{B}}_{1}^{gg}(z)+\frac{\alpha_{s}(E^{2}(1-z)^{2}\theta^{2})}{2\pi}\left({\mathcal{B}}^{gg}_{2}(z)+{\mathcal{B}}_{1}^{gg}(z)b_{0}\ln(1-z)^{2}\right)\,,$
$\displaystyle{\mathcal{P}}_{qg}(z,\theta)$
$\displaystyle\equiv\,{\mathcal{B}}_{1}^{qg}(z)+\frac{\alpha_{s}(E^{2}(1-z)^{2}\theta^{2})}{2\pi}\left({\mathcal{B}}^{qg}_{2}(z)+{\mathcal{B}}_{1}^{qg}(z)b_{0}\ln(1-z)^{2}\right)\,,$
(61)

where the LO anomalous dimensions read

$\displaystyle{\mathcal{B}}_{1}^{gg}(z)=C_{A}\left(z(1-z)-2\right)\,,\quad{\mathcal{B}}_{1}^{qg}(z)=T_{R}n_{f}\left(z^{2}+(1-z)^{2}\right)\,.$
(62)

The ratio of strong couplings in the second line of
${\mathcal{P}}_{gg}(z,\theta)$ has the role of restoring the correct scale of
the coupling corresponding to the soft singularity as $z\to 0$ as opposed to
the one at $z\to 1$ that is encoded in the evolution time (1). This feature is
of course present exclusively for gluon jets, and it is absent in the quark
case (2.2). The quantity $\mathcal{B}_{2}^{g}(z)$ is then simply the sum of
$\mathcal{B}_{2}^{gg}(z)$ and $\mathcal{B}_{2}^{qg}(z)$. A subtle aspect of
the above decomposition between terms that are interpreted as corrections to
either the $g\to gg$ or the $g\to q\bar{q}$ channel is that they both receive
a contribution from the $C_{A}T_{R}$ colour factor. The separation between
such contributions to $\mathcal{B}_{2}^{g}(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\to gg$. We thus write:

$\displaystyle{\mathcal{B}}^{gg}_{2}(z)$ $\displaystyle\equiv
C_{A}\,T_{R}\,n_{f}\,\mathcal{B}_{2}^{g,\,C_{A}T_{R}}(z)+C_{A}^{2}\,\mathcal{B}_{2}^{g,\,C_{A}^{2}}(z)\,,$
$\displaystyle{\mathcal{B}}^{qg}_{2}(z)$ $\displaystyle\equiv
T_{R}^{2}n_{f}^{2}\,\mathcal{B}_{2}^{g,\,T_{R}^{2}}(z)+C_{F}\,T_{R}\,n_{f}\,\mathcal{B}_{2}^{g,\,C_{F}\,T_{R}}(z)\,.$
(63)

The Sudakov form factor for gluon jets, given at NLL in Eq. (5), at NNLL then
reads

$\displaystyle\ln\Delta_{g}(t)$
$\displaystyle=-\int_{t}^{t_{0}}dt^{\prime}\int_{z_{0}}^{1-z_{0}}dz\,\left({\mathcal{P}}_{gg}(z,\theta)+{\mathcal{P}}_{qg}(z,\theta)\right)\,.$
(64)

In the following we carry out the calculation of $\mathcal{B}_{2}^{g}(z)$ in
each of the above colour channels.

###  7.1 The $T_{R}^{2}\,n_{f}^{2}$ channel

This channel is distinct in that it has no double real contribution. From Eqs.
(30), (31), (5) and after subtracting the NLL contribution that emerges from
Eq. (138) we obtain

$\displaystyle\mathcal{B}^{g,T_{R}^{2}}_{2}(z)=\,p_{qg}(z)\left(\frac{1}{3}+\frac{4}{3}\ln(z(1-z))\right)\,,$
(65)

which integrates to

$\displaystyle\int_{0}^{1}dz\,\mathcal{B}^{g,T_{R}^{2}}_{2}(z)=-\frac{46}{27}\,.$
(66)

###  7.2 The $C_{F}T_{R}n_{f}$ channel

The extraction of $\mathcal{B}_{2}^{g}(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. (43) and its mirror
symmetric obtained by the swap $z\leftrightarrow 1-z$, and virtual, Eq. (31),
corrections to find

$\displaystyle\mathcal{B}_{2}^{g,C_{F}T_{R}}(z)=p_{qg}(z)\left(\ln^{2}\left(\frac{z}{1-z}\right)-\frac{\pi^{2}}{3}+5\right)+H^{C_{F}T_{R}}_{\textrm{fin.}}(z)\,,$
(67)

whose integral reads

$\displaystyle\int_{0}^{1}dz\,\mathcal{B}^{g,C_{F}T_{R}}_{2}(z)=1\,.$ (68)

###  7.3 The $C_{A}T_{R}n_{f}$ channel

We combine the real and virtual terms Eqs. (6.2), (49), (30), (31), (35) and
(5) and then subtract the NLL contribution that emerges from Eq. (138), to
obtain

$\mathcal{B}_{2}^{g,C_{A}T_{R}}(z)=-p_{qg}(z)\left(\ln^{2}z+\ln^{2}(1-z)\right)+\frac{1}{9}(28-41z+41z^{2})\\\
+\ln
z\left(\frac{4}{3(1-z)}-\frac{26}{3}z^{2}+8z-7\right)+\ln(1-z)\left(\frac{4}{3z}-\frac{26}{3}(1-z)^{2}+8(1-z)-7\right),$
(69)

which integrates to

$\int_{0}^{1}dz\,\mathcal{B}_{2}^{g,C_{A}T_{R}}(z)=\frac{593}{54}-\frac{4\pi^{2}}{9}.$
(70)

Before we discuss the more involved $C_{A}^{2}$ channel, it is instructive to
compare the integral of the contributions to $\mathcal{B}_{2}^{g}(z)$ computed
so far to the expectation given in Eq. (25). In particular, given that the
term $\mathcal{F}_{\text{clust.}}^{C_{A}^{2}}$ is a pure $C_{A}^{2}$
contribution, the combination of Eqs. (66), (70) allows us to extract the
constant $X^{g}_{\theta^{2}}$ corresponding to the observable
$\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma}{d\theta^{2}dz}$ used in the
calculation. This yields

$X^{g}_{\theta^{2}}=-C_{A}\left(\frac{67}{9}-\frac{2\pi^{2}}{3}\right)+\frac{23}{9}T_{R}n_{f}\,,$
(71)

which can be used as a cross check in the $C_{A}^{2}$ channel below.

###  7.4 The $C_{A}^{2}$ channel

The procedure for extracting $\mathcal{B}_{2}(z)$ is the same as in the
previous channel. The relevant real emission terms are given in Eqs. (51),
(57), (58) which can then be combined with the virtual contributions emerging
from Eqs. (30), (35), (5). After removal of NLL terms the result can be
expressed as

$\mathcal{B}_{2}^{g,C_{A}^{2}}(z)=\mathcal{B}_{2}^{g,C_{A}^{2},\text{analytic}}(z)+G_{z_{3}>z_{\rm
cut}}(z)+\mathcal{B}_{2}^{\text{endpoint}}(z).$ (72)

The term $\mathcal{B}_{2}^{g,C_{A}^{2},\text{analytic}}(z)$ is computed
analytically while the contribution $G_{z_{3}>z_{\rm cut}}(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_{\mathrm{cut}}$. The term
$\mathcal{B}_{2}^{g,C_{A}^{2},\text{analytic}}(z)$ reads:

$\mathcal{B}_{2}^{g,C_{A}^{2},\text{analytic}}(z)=p_{gg}(z)\left(-2\ln
z\ln(1-z)-\frac{11}{3}\ln z-\frac{11}{3}\ln(1-z)\right)\\\
+\frac{11}{3}\left(\frac{\ln
z}{z}+\frac{\ln(1-z)}{1-z}\right)-\frac{265}{18}+\frac{\pi^{2}}{3}(2-z(1-z))+4\,h_{\rm
pass}^{\rm poles}-\frac{44}{3}\ln 2+8\ln^{2}2\\\
-2z(1-z)\left(-\frac{67}{18}+h_{\rm pass}^{\rm
poles}+2\ln^{2}2-\frac{11}{3}\ln 2\right)+h_{\rm
fail}^{p_{gg}}\,(-2+z(1-z))\,.$ (73)

Its integral is straightforward and gives:

$\int
dz\,\mathcal{B}_{2}^{g,C_{A}^{2},\text{analytic}}(z)=-\frac{535}{27}+\frac{11\pi^{2}}{9}+\frac{11}{3}h_{\rm
pass}^{\rm poles}-\frac{121}{9}\ln
2+\frac{22}{3}\ln^{2}2-4\zeta(3)-\frac{11}{6}\,h_{\rm fail}^{p_{gg}}.$ (74)

The end-point contribution is simply given by

$\mathcal{B}_{2}^{\text{endpoint}}(z)=h_{\rm fail}^{\delta}\,\delta(z)\,.$
(75)

The constants $h_{\rm pass}^{\rm poles}$, $h_{\rm fail}^{p_{gg}}$, and $h_{\rm
fail}^{\delta}$ are defined in Eqs. (53), (59) and the function
$G_{z_{3}>z_{\rm cut}}(z)$ is defined in Eq. (54). The total integral of
$\mathcal{B}_{2}^{g,C_{A}^{2}}(z)$ then yields

$\int_{0}^{1}dz\,\mathcal{B}_{2}^{g,C_{A}^{2}}(z)\simeq-6.314\pm 0.004\,.$
(76)

As a cross check of the above result we can compare to the expected integral
in Eq. (25), where $X_{\theta^{2}}^{g}$ is given in Eq. (71). For this check
we still have to determine the function
$\mathcal{F}_{\text{clust.}}^{C_{A}^{2}}$. For the $C_{A}^{2}$ channel our
procedure for the extraction of $\mathcal{B}_{2}^{g,C_{A}^{2}}(z)$ agrees
exactly with the mMDT/Soft drop grooming with the Cambridge/Aachen algorithm,
for which $\mathcal{F}_{\text{clust.}}^{C_{A}^{2}}$ 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
$\mathcal{F}_{\text{clust.}}^{C_{A}^{2}}$ to the known result for quark jets
(see e.g. Anderle:2020mxj ) (cf. also the discussion in Appendix E). Replacing
$C_{F}\to C_{A}$ in the quark jet result we get the following result the
$C_{A}^{2}$ component reads:

$\mathcal{F}_{\text{clust.}}^{C_{A}^{2}}=\frac{1}{2}C^{2}_{A}\left(\frac{4\pi}{3}\text{Cl}_{2}\left(\frac{\pi}{3}\right)+\,h^{C_{A}}_{\text{clust.}}\right)\,,$
(77)

where

$h^{C_{A}}_{\text{clust.}}\simeq-1.16363257(4)\,,$ (78)

and $\text{Cl}_{2}$ denotes the Clausen function. The overall factor of $1/2$
is multiplied to obtain the result for a single leg. This gives

$(-\gamma_{g}^{(2)}+b_{0}X_{\theta^{2}}^{g}+\mathcal{F}_{\text{clust.}})_{C_{A}^{2}}\simeq-6.31426325(8)\,C_{A}^{2},$
(79)

in good agreement with our result.

We conclude this section with a remark on the endpoint contribution
$\mathcal{B}_{2}^{\text{endpoint}}(z)$ given in Eq. (75). As discussed above,
this originates from double-soft configurations and is a consequence of the
scheme used here to define $z$ and $\theta$. 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 $\theta$ 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_{A}^{2}$ channel. The same double-soft origin is shared by the clustering
correction (77). Terms of double-soft origin are discussed further in Appendix
E.

##  8 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
($\beta=0$) procedure Dasgupta:2013ihk ; Larkoski:2014wba . These classes of
jet substructure observables have received widespread attention in the
literature, with several applications both at hadron and lepton colliders
Frye:2016aiz ; Frye:2016okc ; Larkoski:2017bvj ; Marzani:2017kqd ; CMS:2018vzn
; ATLAS:2019mgf ; Marzani:2019evv ; Kardos:2020gty ; Kardos:2020ppl ;
CMS:2021vsp ; Cal:2021fla ; Caletti:2021oor ; Reichelt:2021svh ;
Hannesdottir:2022rsl . As a concrete example we consider the processes $Z\to
q\bar{q}$ and $H\to 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 Banfi:2004yd

$\displaystyle FC^{\mathcal{H}}_{x}=\frac{2^{-x}}{E^{2}}\sum_{i\neq
j}E_{i}E_{j}|\sin\theta_{ij}|^{x}(1-|\cos\theta_{ij}|)^{1-x}\ ,$ (80)

where the sum runs over all particles within a given hemisphere,
$\mathcal{H}\equiv\mathcal{H}_{R}$ or $\mathcal{H}\equiv\mathcal{H}_{L}$, and
$E$ denotes the total energy in the hemisphere. We also consider the
angularities Berger:2003iw  (the corresponding NNLL resummation in the
ungroomed case is given in Refs. Larkoski:2014uqa ; Bell:2018gce ;
Banfi:2018mcq ; Procura:2018zpn ; Bauer:2020npd ) defined w.r.t. the Winner-
Take-All (WTA) axis as

$\displaystyle\lambda^{\mathcal{H}}_{x}=\frac{2^{1-x}}{E}\sum_{i}E_{i}\,|\sin\theta_{i}|^{x}(1-|\cos\theta_{i}|)^{1-x}\,,$
(81)

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:

$\displaystyle
FC_{x}\equiv{\max}\left\\{FC^{\mathcal{H}_{R}}_{x},FC^{\mathcal{H}_{L}}_{x}\right\\},\,\lambda_{x}\equiv{\max}\left\\{\lambda^{\mathcal{H}_{R}}_{x},\lambda^{\mathcal{H}_{L}}_{x}\right\\}\,.$
(82)

A NNLL calculation for WTA angularities for quark jets has been recently
presented in Ref. Dasgupta:2022fim ,121212Notice that here we adopt a
different normalisation of the observable, compared to Ref. Dasgupta:2022fim ,
in order to match the corresponding hadron collider jet definitions. 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 $\Sigma(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. (7). In the following we work in the limit
$v\ll z_{\rm cut}\ll 1$, that is commonly considered for these types of
groomed observables. This regime has the advantage that one can neglect power
corrections in $z_{\rm 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 ${\mathbb{K}}_{q}^{\rm finite}[G_{q},G_{g}]$
and ${\mathbb{K}}_{g}^{\rm finite}[G_{q},G_{g}]$ functions (cf. Eqs. (9), (7)
and Appendix C). With a little abuse of notation we have used the same
parameter $z_{\rm cut}$ for the definition of groomed observables and in the
calculation of ${\cal B}^{g}_{2}(z)$. However, we stress that (cf. Sec. 4)
${\cal B}^{g}_{2}(z)$ is defined strictly by taking the limit $z_{\rm cut}\to
0$ while in the observables considered here one can opt to retain finite
$z_{\rm cut}$ effects to improve the accuracy of the calculation. The quantity
${\cal B}^{g}_{2}(z)$ we have derived applies to both cases with and without
finite $z_{\rm 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. (9), (7), and following similar steps to
those outlined in Appendix D, we obtain the following NNLL results for quark
and gluon jets, respectively

$\displaystyle\Sigma^{q}(v)$ $\displaystyle=\sigma^{Z\to
q\bar{q}}_{0}\left(1+\frac{\alpha_{s}(E^{2})}{2\pi}C_{v}^{q(1)}(z_{\rm
cut})\right)e^{-2R^{q}_{v}(v,z_{\rm
cut})}\left(1+\frac{\alpha^{2}_{s}(E^{2})}{(2\pi)^{2}}\,2{\cal
F}^{q}_{\text{clust}}(v)\right)\,,$ $\displaystyle\Sigma^{g}(v)$
$\displaystyle=\sigma^{H\to
gg}_{0}\left(1+\frac{\alpha_{s}(E^{2})}{2\pi}C_{v}^{g(1)}(z_{\rm
cut})\right)e^{-2R^{g}_{v}(v,z_{\rm
cut})}\left(1+\frac{\alpha^{2}_{s}(E^{2})}{(2\pi)^{2}}\,2{\cal
F}^{g}_{\text{clust}}(v)\right)\,.$ (83)

The quantity $R_{v}(v,z_{\rm cut})$ denotes the Sudakov radiator, which can be
obtained by integrating the inclusive emission probability in Eqs. (2.2), (7)
with the measurement function of the observables (80), (81) for a single
collinear $a\to bc$ splitting, namely

$\Theta(V_{x}(z,\theta)-v)\,,$ (84)

with $V_{x}=\\{FC_{x},\lambda_{x}\\}$ and

$\displaystyle FC_{x}(z,\theta)\simeq$ $\displaystyle\,z(1-z)\theta^{2-x}\ \
,$ (85) $\displaystyle\lambda_{x}(z,\theta)\simeq$
$\displaystyle\,\,\text{min}(z,1-z)\theta^{2-x}\,,$ (86)

where $z$ is the energy fraction of $b$. This gives (see also Ref.
Dasgupta:2022fim ) (with $f\in\\{q,g\\}$)

$\displaystyle R^{f}_{v}(v,z_{\rm
cut})\equiv-g^{f}_{1}(\lambda_{v},\lambda_{z_{\mathrm{cut}}})-g_{2}^{f}(\lambda_{v},\lambda_{z_{\mathrm{cut}}})-h^{f}_{2}(\lambda_{v})-\frac{\alpha_{s}}{\pi}h^{f}_{3}(\lambda_{v};B_{2,v}^{f})\,,$
(87)

where $\lambda_{v}\equiv\alpha_{s}\beta_{0}\ln 1/v$, $\lambda_{z_{\rm
cut}}\equiv\alpha_{s}\beta_{0}\ln 1/z_{\rm cut}$,
$\alpha_{s}\equiv\alpha_{s}(E^{2})$, and

$\beta_{0}=\frac{b_{0}}{2\pi}\,.$ (88)

The functions in Eq. (87) read (we use the notation $C_{R}^{q}\equiv C_{F}$
and $C_{R}^{g}\equiv C_{A}$)

$\displaystyle
g^{f}_{1}(\lambda_{v},\lambda_{z_{\mathrm{cut}}})=\frac{C_{R}^{f}}{\pi\alpha_{s}\beta_{0}^{2}}\lambda_{z_{\mathrm{cut}}}\ln\left(1-\frac{2\lambda_{v}}{2-x}\right)\,,$
(89) $\displaystyle g^{f}_{2}(\lambda_{v},\lambda_{z_{\mathrm{cut}}})$
$\displaystyle=\frac{C_{R}^{f}}{\pi\alpha_{s}\beta_{0}^{2}}\lambda_{z_{\mathrm{cut}}}^{2}\left(\frac{1-2\lambda_{v}}{2-x-2\lambda_{v}}\right)+\frac{C_{R}^{f}\beta_{1}}{\pi\beta_{0}^{3}}\lambda_{z_{\mathrm{cut}}}\frac{2\lambda_{v}+(2-x)\ln\left(1-\frac{2\lambda_{v}}{2-x}\right)}{2-x-2\lambda_{v}}$
$\displaystyle-\frac{C_{R}^{f}}{\pi^{2}\beta_{0}^{2}}\lambda_{z_{\mathrm{cut}}}\frac{\lambda_{v}\,K^{(1)}}{2-x-2\lambda_{v}}\,,$
(90) $\displaystyle h^{f}_{2}(\lambda_{v})$
$\displaystyle=\frac{B_{1}^{f}}{2\pi\beta_{0}}\ln\left(1-\frac{2\lambda_{v}}{2-x}\right)\
\ ,$ (91) $\displaystyle h^{f}_{3}(\lambda_{v};B_{2,v}^{f})$
$\displaystyle=\frac{B_{1}^{f}\,\beta_{1}}{2\beta_{0}^{2}(2-x-2\lambda_{v})}\left((2-x)\ln\left(1-\frac{2\lambda_{v}}{2-x}\right)+2\lambda_{v}\right)$
$\displaystyle-\frac{B_{2,v}^{f}}{2\pi\beta_{0}(2-x-2\lambda_{v})}\lambda_{v}\,,$
(92)

where

$\displaystyle\beta_{1}=\frac{17\,C_{A}^{2}-10\,C_{A}T_{R}n_{f}-6\,C_{F}T_{R}n_{f}}{24\pi^{2}}\,.$
(93)

For the observables of interest $B_{2,v}^{f}$ emerges from the integration of
Eqs. (2.2), (7) with the phase space constraint in Eq. (84) and it is given by

$\displaystyle B_{2,FC_{x}}^{q}$
$\displaystyle=B^{q}_{2,\theta^{2}}+C_{F}\,b_{0}\frac{2}{2-x}\left(3-\frac{\pi^{2}}{3}\right)=-\gamma^{(2)}_{q}+b_{0}X_{FC_{x}}^{q}\,,$
(94) $\displaystyle B_{2,FC_{x}}^{g}$
$\displaystyle=B_{2,\theta^{2}}^{g}+\frac{b_{0}}{2-x}\left(C_{A}\left(\frac{67}{9}-\frac{2\pi^{2}}{3}\right)-\frac{26}{9}T_{R}n_{f}\right)$
$\displaystyle=-\gamma^{(2)}_{g}+b_{0}X_{FC_{x}}^{g}+\frac{1}{2}C_{A}^{2}\left(\frac{4\pi}{3}\text{Cl}_{2}\left(\frac{\pi}{3}\right)+\,h^{C_{A}}_{\text{clust.}}\right)\,,$
(95)

for the fractional moments of EEC, and

$\displaystyle B_{2,\lambda_{x}}^{q}$
$\displaystyle=B^{q}_{2,\theta^{2}}+C_{F}\,b_{0}\frac{(9-\pi^{2}+9\ln
2)}{3(2-x)}=-\gamma^{(2)}_{q}+b_{0}X_{\lambda_{x}}^{q}\,,$ (96) $\displaystyle
B_{2,\lambda_{x}}^{g}$
$\displaystyle=B_{2,\theta^{2}}^{g}+\frac{b_{0}}{2-x}\left(C_{A}\left(\frac{137}{36}-\frac{\pi^{2}}{3}+\frac{44\ln
2}{12}\right)-T_{R}n_{f}\left(\frac{29}{18}+\frac{4\ln 2}{3}\right)\right)$
$\displaystyle=-\gamma^{(2)}_{g}+b_{0}X_{\lambda_{x}}^{g}+\frac{1}{2}C_{A}^{2}\left(\frac{4}{3}\pi\text{Cl}_{2}\left(\frac{\pi}{3}\right)+\,h^{C_{A}}_{\text{clust.}}\right)\,,$
(97)

for the WTA angularities. The quantities $B_{2,\theta^{2}}^{q}$ and
$B_{2,\theta^{2}}^{g}$ are defined in Eqs. (22) and (25) while $X^{f}_{v}$
denotes the observable dependent constants for which we obtain

$\displaystyle X_{FC_{x}}^{q}$
$\displaystyle=C_{F}\,\frac{2}{2-x}\left(3-\frac{\pi^{2}}{3}\right)+X^{q}_{\theta^{2}}\,,$
$\displaystyle X_{FC_{x}}^{g}$
$\displaystyle=\frac{1}{2-x}\left(C_{A}\left(\frac{67}{9}-\frac{2\pi^{2}}{3}\right)-\frac{26}{9}T_{R}n_{f}\right)+X^{g}_{\theta^{2}}\,,$
(98)

and

$\displaystyle X_{\lambda_{x}}^{q}$
$\displaystyle=C_{F}\,\frac{(9-\pi^{2}+9\ln 2)}{3(2-x)}+X^{q}_{\theta^{2}}\,,$
$\displaystyle X_{\lambda_{x}}^{g}$
$\displaystyle=\frac{1}{2-x}\left(C_{A}\left(\frac{137}{36}-\frac{\pi^{2}}{3}+\frac{44\ln
2}{12}\right)-T_{R}n_{f}\left(\frac{29}{18}+\frac{4\ln
2}{3}\right)\right)+X^{g}_{\theta^{2}}\,,$ (99)

where the constants $X^{q}_{\theta^{2}}$ and $X^{g}_{\theta^{2}}$ are given in
Eqs. (23), (71). Finally, $h^{C_{A}}_{\text{clust.}}$ is given in Eq. (78). We
note that the result for $B_{2,\lambda_{x}}^{q}$ has been previously obtained
in Ref. Dasgupta:2022fim .

The clustering corrections ${\cal F}_{\text{clust.}}(v)$ in Eq. (8) arise from
the finite terms ${\mathbb{K}}_{q}^{\rm finite}[G_{q},G_{g}]$ and
${\mathbb{K}}_{g}^{\rm finite}[G_{q},G_{g}]$ in the GF Eqs. (9), (7) (cf.
Appendix C), which give rise to a correction to the Sudakov in Eq. (8) and are
given by (see also Ref. Dasgupta:2022fim )

$\displaystyle{\cal F}^{q}_{\text{clust.}}(v)$
$\displaystyle=C_{F}\left(C_{F}\frac{4\pi}{3}\text{Cl}_{2}\left(\frac{\pi}{3}\right)+C_{A}\,h^{C_{A}}_{\text{clust.}}+T_{R}n_{f}\,h^{T_{R}n_{f}}_{\text{clust.}}\right)\frac{\ln
v}{2-x-2\lambda_{v}}\,,$ $\displaystyle{\cal F}^{g}_{\text{clust.}}(v)$
$\displaystyle=C_{A}T_{R}n_{f}\,h^{T_{R}n_{f}}_{\text{clust.}}\frac{\ln
v}{2-x-2\lambda_{v}}\,,$ (100)

where

$h^{T_{R}n_{f}}_{\text{clust.}}\simeq-1.75559363(5)\,.$ (101)

A comment on the difference between quark and gluon jets is in order. This
difference can be traced back to the projection (24) used in the definition of
the inclusive emission probability and hence of ${\cal B}_{2}^{f}(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_{2,v}^{f}$ given in Eqs. (94), (96). Here we can see that the gluon jet
result $B_{2,v}^{g}$, unlike $B_{2,v}^{q}$, includes a part of the clustering
correction, specifically in the $C_{A}^{2}$ channel (see discussion in Sec.
7.4).

Finally, we comment on the coefficient functions $C_{v}^{f(1)}(z_{\rm cut})$
($f=\\{q,g\\}$) in Eq. (8). These are matching coefficients between the
collinear approximation of the generating functionals and the full ${\cal
O}(\alpha_{s})$ calculation in the limit $v\ll z_{\rm cut}\ll 1$. For the
processes and observables considered here they have the general structure
Banfi:2018mcq

$\displaystyle C_{v}^{q(1)}(z_{\rm cut})$
$\displaystyle=H^{q(1)}-2X_{v}^{q}+C_{F}\left(8\ln 2\ln
z_{\mathrm{cut}}+\,6\ln 2-\frac{\pi^{2}}{3}\right)\,,$ $\displaystyle
C_{v}^{g(1)}(z_{\rm cut})$ $\displaystyle=H^{g(1)}-2X_{v}^{g}+C_{A}\left(8\ln
2\ln z_{\mathrm{cut}}-\frac{\pi^{2}}{3}\right)\,,$ (102)

where $X_{v}^{f}$ are the observable dependent constants of hard-collinear
origin given in Eqs. (8), (8), and the remaining constants are the the (only)
process-dependent ingredients of the calculation. Specifically, $H^{f(1)}$
accounts for the hard-virtual corrections at one loop order

$\displaystyle H^{q(1)}$ $\displaystyle=\pi^{2}-8\,,$ $\displaystyle H^{g(1)}$
$\displaystyle=\pi^{2}\,.$ (103)

The presence of $X_{v}^{f}$ both in Eqs. (8) and in Eqs. (94), (96) indicates
that the GF solution is defined in a resummation scheme in which the running
of $\alpha_{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.
Banfi:2014sua ; Banfi:2018mcq ; Dasgupta:2022fim ) the coupling multiplying
these constants is evaluated at the observable-dependent collinear (low)
scale, which would remove $X_{v}^{f}$ from Eqs. (94), (96). The physical
result at a given logarithmic order is of course resummation-scheme invariant
(cf. also footnote 17 in Appendix D).

As a check of our results, we compare the resummed predictions for $FC_{x}$
and $\lambda_{x}$ against the formalism of Refs. Banfi:2014sua ; Banfi:2018mcq
adapted to groomed observables in Refs. Anderle:2020mxj ; Dasgupta:2022fim .
Specifically, as observed in Ref. Dasgupta:2022fim , in the case of groomed
quark jets the integrated quantities $B_{2,\lambda_{x}}^{q}$ and
$B_{2,FC_{x}}^{q}$ can be extracted from Sec. 3 of Ref. Banfi:2018mcq , by
combining the endpoint of the two loop DGLAP anomalous dimension
$\gamma^{(2)}_{f}$ with the running of the one-loop hard-collinear constant
$C^{(1)}_{\rm hc}$ and the recoil correction $\delta{\cal F}_{\rm
rec}$,131313The hard-collinear correction $\delta{\cal F}_{\rm hc}$ would also
contribute for ungroomed angularities and fractional moments of EEC. which
agrees with our findings in Eqs. (94), (96). Similarly, an independent
calculation of the clustering corrections for quark angularities can be found
in Ref. Dasgupta:2022fim . We also carried out analogous checks for the
gluonic results, based on a straightforward extension of the above formalism
to $H\to gg$. 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).

##  9 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 triple-collinear 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 ${\cal 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. 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. Dasgupta:2014yra . We will address some of these aspects in a
future publication vanBeekveld:2023?  together with an algorithmic solution to
the problem of collinear fragmentation.

  2. 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 Frye:2016aiz ; Frye:2016okc ;
Larkoski:2017bvj ; Marzani:2017kqd ; Marzani:2019evv ; Kardos:2020gty ;
Kardos:2020ppl ; Caletti:2021oor ; Reichelt:2021svh ; Hannesdottir:2022rsl .
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. 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. Banfi:2021xzn ; Banfi:2021owj ).
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 wide classes of observable. 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.

## Acknowledgements

We thank Ludovic Scyboz, Grégory Soyez and Gavin Salam for useful discussions
through the course of this work and comments on the manuscript. Furthermore,
we are grateful to Gavin Salam for collaboration on the related forthcoming
Ref. vanBeekveld:2023? . We would also like to thank other colleagues on the
PanScales collaboration whose thoughts and efforts towards pushing forward the
logarithmic accuracy of parton showers have directly motivated this work. This
work has been partly funded by the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program (grant agreement
No 788223) (MvB, MD, BKE and JH) and by the U.K.’s Science and Technologies
Facilities Council under grant ST/T001038 (MD). The work of PM is funded by
the European Union (ERC, grant agreement No. 101044599, JANUS). Views and
opinions expressed are however those of the authors only and do not
necessarily reflect those of the European Union or the European Research
Council Executive Agency. Neither the European Union nor the granting
authority can be held responsible for them.

##  Appendix A Leading order splitting functions

In this appendix we report the well known expression of the unregularised
tree-level splitting kernels used in the main text. We factor out the colour
factors using the notation

$\displaystyle P_{qq}(z,\epsilon)=C_{F}\,p_{qq}(z,\epsilon)\,,\quad
P_{gg}(z)=C_{A}\,p_{gg}(z)\,,\quad
P_{qg}(z,\epsilon)=T_{R}\,n_{f}\,p_{qg}(z,\epsilon)\,,$ (104)

where the splitting functions in $4-2\epsilon$ dimensions read

$\displaystyle p_{qq}(z,\epsilon)$
$\displaystyle=\frac{1+z^{2}}{1-z}-\epsilon(1-z)\,,$ $\displaystyle p_{gg}(z)$
$\displaystyle=\frac{1}{1-z}+\frac{1}{z}-2+z(1-z)\,,$ $\displaystyle
p_{qg}(z,\epsilon)$
$\displaystyle=\frac{z^{2}+(1-z)^{2}-\epsilon}{1-\epsilon}\,.$ (105)

For the sake of simplicity, in the text we also use the following notation for
the $\epsilon=0$ case

$\displaystyle p_{qq}(z)$ $\displaystyle\equiv p_{qq}(z,\epsilon=0)\,,$
$\displaystyle p_{qg}(z)$ $\displaystyle\equiv p_{qg}(z,\epsilon=0)\,.$ (106)

##  Appendix B Expression of ${\mathcal{B}}_{2}(z)$ for quark fragmentation

The expression of ${\mathcal{B}}^{q}_{2}(z)$ can be organised as
follows:141414The function ${\mathcal{B}}^{q}_{2}(z)$ computed in Ref.
Dasgupta:2021hbh  is defined as the ${\mathcal{B}}^{q}_{2}(z)$ used here
multiplied by $\left(\alpha_{s}/(2\pi)\right)^{2}$.

$\displaystyle{\mathcal{B}}^{q}_{2}(z)$
$\displaystyle=C^{2}_{F}\,\mathcal{B}_{2}^{q,C_{F}^{2}}(z)+C_{F}C_{A}\,\mathcal{B}_{2}^{q,\,C_{F}C_{A}}(z)+C_{F}T_{R}n_{f}\,\mathcal{B}_{2}^{q,\,C_{F}T_{R}}(z)+C_{F}\left(C_{F}-\frac{C_{A}}{2}\right)\mathcal{B}_{2}^{q,\,\textrm{id.}}(z)\
.$ (107)

The above functions read

$\displaystyle\mathcal{B}_{2}^{q,\,\textrm{id.}}(z)$
$\displaystyle=4z-\frac{7}{2}+\frac{5z^{2}-2}{2(1-z)}\ln
z+\frac{1+z^{2}}{1-z}\left(\frac{\pi^{2}}{6}-\ln
z\ln(1-z)-\text{Li}_{2}(z)\right)\,,$ (108)
$\displaystyle\mathcal{B}_{2}^{q,\,C_{F}T_{R}}(z)$
$\displaystyle=-b_{0}^{(n_{f})}p_{qq}(z)\ln
z+b_{0}^{(n_{f})}(1-z)-K^{(1),n_{f}}(1+z)+2\,b_{0}^{(n_{f})}(1+z)\ln(1-z)\,,$
(109) $\displaystyle\mathcal{B}_{2}^{q,\,C_{F}C_{A}}(z)$
$\displaystyle=-b_{0}^{(C_{A})}p_{qq}(z)\ln
z+b_{0}^{(C_{A})}(1-z)+\frac{3}{2}\frac{z^{2}\ln z}{1-z}+\frac{1}{2}(2z-1)$
(110)
$\displaystyle+2\,b_{0}^{(C_{A})}(1+z)\ln(1-z)+p_{qq}(z)\left(\ln^{2}z+\mathrm{Li}_{2}\left(\frac{z-1}{z}\right)+2\,\mathrm{Li}_{2}(1-z)\right)-K^{(1),C_{A}}(1+z)\,,$
$\displaystyle\mathcal{B}_{2}^{q,C_{F}^{2}}(z)$
$\displaystyle=p_{qq}(z)\left(-3\ln z-2\ln
z\ln(1-z)+2\,\text{Li}_{2}\left(\frac{z-1}{z}\right)\right)-1+H^{\rm
fin.}(z)\,,$ (111)

where $H^{\rm fin.}(z)$ is given by a 1-fold integral (cf. Figure 4 of ref.
Dasgupta:2021hbh ), that is provided in Mathematica format as an ancillary
file with the arXiv preprint of this article. The function
${\mathcal{B}}^{q}_{2}(z)$ is regular in the soft limit $z\to 1$ and is thus
fully integrable over $z\in[0,1]$.

##  Appendix C The NNLL ${\mathbb{K}}^{\rm finite}$ kernel

In this appendix we report the functions ${\mathbb{K}}_{f}^{\rm finite}$ (with
$f\in\\{q,g\\}$) entering the NNLL evolution equation for the quark generating
functionals given in Eqs. (9), (7). 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 ${\mathbb{K}}_{q}^{\rm finite}$ as a difference between two
terms which encode the (subtracted) real corrections to the first of Eqs.
(2.1) and its double counting with the iteration of the NLL kernel,
respectively. That is:

${\mathbb{K}}_{q}^{\rm finite}[G_{q},G_{g}]\equiv{\mathbb{K}}_{q}^{\rm
R}[G_{q},G_{g}]-{\mathbb{K}}_{q}^{\rm DC}[G_{q},G_{g}]\,.$ (112)

The difference of ${\mathbb{K}}_{q}^{\rm R}$ and ${\mathbb{K}}_{q}^{\rm DC}$
ensures that the quantity ${\mathbb{K}}_{q}^{\rm finite}$ is infrared finite
and purely NNLL.

$\theta_{12,3}$$z_{1}=(1-z)z_{p}$$\theta_{12}$$z_{2}=(1-z)(1-z_{p})$$z_{3}=z$
Figure 7: The diagram representing gluon decay to a $q\bar{q}$ pair, where the
quark from the gluon decay is either identical or non-identical to the
initiating quark.
$z_{1}=1-z$$\theta_{1,23}$$z_{2}=z(1-z_{p})$$\theta_{23}$$z_{3}=zz_{p}$ Figure
8: The diagram representing the gluon emission $C_{F}^{2}$ channel.

Following the definition of the inclusive emission probability (2.2) we
obtain:

$\displaystyle{\mathbb{K}}_{q}^{\rm R}$
$\displaystyle[G_{q},G_{g}]=\sum_{(A)}\frac{1}{S_{2}}\int
d\Phi^{(A)}_{3}\,P^{(A)}_{1\to
3}\,\bigg{\\{}G_{f_{1}}(x\,z_{p}\,(1-z),t_{1,2})\,G_{f_{2}}(x\,(1-z_{p})\,(1-z),t_{1,2})$
$\displaystyle\times
G_{q}(x\,z,t_{12,3})-G_{f_{12}}(x\,(1-z),t_{12,3})\,G_{q}(x\,z,t_{12,3})\bigg{\\}}\frac{\Delta_{q}(t)}{\Delta_{q}(t_{1,2})}$
$\displaystyle+\int d\Phi^{(B)}_{3}\,P^{(B)}_{1\to
3}\,\bigg{\\{}G_{g}(x\,(1-z),t_{1,23})\,G_{g}(x\,z\,(1-z_{p}),t_{2,3})$
$\displaystyle\times
G_{q}(x\,z\,z_{p},t_{2,3})-G_{g}(x\,(1-z),t_{1,23})\,G_{q}(x\,z,t_{1,23})\bigg{\\}}\frac{\Delta_{q}(t)}{\Delta_{q}(t_{2,3})}\,\Theta(t_{2,3}-t_{1,3})\,,$
(113)

where we have used the notation $t_{i,j}$ to indicate the value of the
evolution time (1) corresponding to the angle $\theta_{i,j}$. Moreover, we
have parameterised the integrands as 151515We notice that the production of a
$q\bar{q}q$ final state with identical flavours contributing to the
$C_{F}(C_{F}-C_{A}/2)$ is finite, and therefore does not factorise into the
product of two splitting functions in the strong angular ordered limit.

$\displaystyle P_{1\to 3}^{(A)}$
$\displaystyle\equiv\frac{(8\pi)^{2}}{s_{123}^{2}}\alpha^{2}_{s}(E^{2}g^{2}(z)\,\theta_{12,3}^{2})\langle\hat{P}\rangle_{\tiny{C_{F}C_{A}\,,C_{F}T_{R}n_{f}\,,C_{F}(C_{F}-C_{A}/2)}}\,,$
(114) $\displaystyle P_{1\to 3}^{(B)}$
$\displaystyle\equiv\frac{(8\pi)^{2}}{s_{123}^{2}}\alpha^{2}_{s}(E^{2}g^{2}(z)\,\theta_{1,3}^{2})\langle\hat{P}\rangle_{\tiny{C_{F}^{2}}}\,.$
(115)

Here, $P_{1\to 3}^{(A)}$ is parameterised according to the phase space
depicted in Fig. 7, while $P_{1\to 3}^{(B)}$ is parameterised according to
Fig. 8. The corresponding phase space measures, denoted by $d\Phi^{(A,B)}_{3}$
in (C.1) are obtained from Eq. (39) by setting $\epsilon=0$ and performing the
change of variable (with the corresponding Jacobian) of Figs. 7, 8. The sum in
Eq. (C.1) runs over all the $(A)$ colour channels of the type $q\to
f_{1}f_{2}q(\bar{q})$ defined in Eq. (114), 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}\geq\\{z,z_{p}\\}\geq z_{0}$. All the calculations are then meant to
be performed by taking the limit of $t_{0}\to\infty$ and $z_{0}\to 0$ after
combining ${\mathbb{K}}_{q}^{\rm R}$ and ${\mathbb{K}}_{q}^{\rm DC}$. We
observe that each term in curly brackets in ${\mathbb{K}}_{q}^{\rm R}$ 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-collinear-
safe definition of the branching variables $\theta$ and $z$, that uniquely
specifies the definition of ${\mathcal{B}}^{q}_{2}(z)$. Therefore, the scheme
that defines ${\mathcal{B}}^{q}_{2}(z)$ is tightly connected to the specific
form of ${\mathbb{K}}_{q}^{\rm finite}[G_{q},G_{g}]$ in Eq. (9), which
specifically adopts the definition of $z$ and $\theta$ of Ref.
Dasgupta:2021hbh  outlined in Sec. 3. An alternative scheme will modify the
expression of both ${\mathcal{B}}^{q}_{2}(z)$ and ${\mathbb{K}}_{q}^{\rm
finite}[G_{q},G_{g}]$ by modifying the precise definition of $\theta$ 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

$\displaystyle{\mathbb{K}}_{q}^{\rm DC}$
$\displaystyle[G_{q},G_{g}]=\sum_{g\to
f\bar{f}}\int_{t}^{t_{0}}dt_{12,3}\,dt_{1,2}\int_{z_{0}}^{1-z_{0}}dz\,dz_{p}P_{qq}(z)P_{fg}(z_{p})\,\bigg{\\{}G_{f}(x\,z_{p}\,(1-z),t_{1,2})$
$\displaystyle\times
G_{f}(x\,(1-z_{p})\,(1-z),t_{1,2})G_{q}(x\,z,t_{12,3})-G_{g}(x\,(1-z),t_{12,3})\,G_{q}(x\,z,t_{12,3})\bigg{\\}}$
$\displaystyle\times\frac{\Delta_{q}(t)}{\Delta_{q}(t_{1,2})}\,\Theta(t_{1,2}-t_{12,3})$
$\displaystyle+\int_{t}^{t_{0}}dt_{1,23}\,dt_{2,3}\int_{z_{0}}^{1-z_{0}}dz\,dz_{p}P_{qq}(z)P_{qq}(z_{p})\,\bigg{\\{}G_{g}(x\,(1-z),t_{1,23})\,G_{g}(x\,z\,(1-z_{p}),t_{2,3})$
$\displaystyle\times
G_{q}(x\,z\,z_{p},t_{2,3})-G_{g}(x\,(1-z),t_{1,23})\,G_{q}(x\,z,t_{1,23})\bigg{\\}}\frac{\Delta_{q}(t)}{\Delta_{q}(t_{2,3})}\,\Theta(t_{2,3}-t_{1,23})\,,$
(116)

where the sum runs over the $g\to gg$ and $g\to q\bar{q}$ splitting channels.
Eq. (C.1) is obtained by calculating the contribution to Eq. (C.1) due to the
iteration of the NLL kernel in Eq. (2.1). We notice that the difference
between Eqs. (C.1) and (C.1) (i.e. ${\mathbb{K}}_{q}^{\rm
finite}[G_{q},G_{g}]$) only contributes in the regime where all angles are
commensurate $t_{2,3}\simeq t_{1,23}\simeq t_{1,3}$. Conversely,
${\mathbb{K}}_{q}^{\rm finite}[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}\simeq t_{1,23}\simeq t_{12,3}\simeq
t_{1,3}$ neglecting higher-order, N3LL terms.

###  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 ${\mathcal{B}}^{g}_{2}(z)$ used in the main text for its
computation. As we discussed at length in Sec. 4, the definition of
${\mathcal{B}}^{g}_{2}(z)$ relies on a specific projection (24) 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 $\theta$ used in the calculation of ${\mathbb{K}}^{\rm
R}[G_{q},G_{g}]$. The latter now reads:

${\mathbb{K}}_{g}^{\rm R}[G_{q},G_{g}]={\mathbb{K}}_{g}^{\rm
R,\,C_{A}\,T_{R}}[G_{q},G_{g}]+{\mathbb{K}}_{g}^{\rm
R,\,C_{F}\,T_{R}}[G_{q},G_{g}]+{\mathbb{K}}_{g}^{\rm
R,\,C_{A}^{2}}[G_{q},G_{g}]\,,$ (117)

where (using the parameterisation of Figs. 1 and 3)

$\displaystyle{\mathbb{K}}_{g}^{\rm R,\,C_{A}\,T_{R}}$
$\displaystyle[G_{q},G_{g}]=\int d\Phi^{(A)}_{3}\,P^{C_{A}\,T_{R}}_{1\to
3}\,\bigg{\\{}G_{q}(x\,z_{p}\,(1-z),t_{2,3})\,G_{q}(x\,(1-z_{p})\,(1-z),t_{2,3})$
$\displaystyle\times
G_{g}(x\,z,t_{1,23})-G_{g}(x\,(1-z),t_{1,23})\,G_{g}(x\,z,t_{1,23})\bigg{\\}}\frac{\Delta_{g}(t)}{\Delta_{g}(t_{2,3})}\,,$
(118) $\displaystyle{\mathbb{K}}_{g}^{\rm R,\,C_{F}\,T_{R}}$
$\displaystyle[G_{q},G_{g}]=\int d\Phi^{(B)}_{3}\,P^{C_{F}\,T_{R}}_{1\to
3}\,\Bigg{[}\bigg{\\{}G_{q}(x\,(1-z),t_{2,13})\,G_{g}(x\,z\,(1-z_{p}),t_{1,3})$
$\displaystyle\times
G_{q}(x\,z\,z_{p},t_{1,3})-G_{q}(x\,(1-z),t_{2,13})\,G_{q}(x\,z,t_{2,13})\bigg{\\}}\frac{\Delta_{g}(t)}{\Delta_{g}(t_{1,3})}\,\Theta(t_{1,3}-t_{1,2})$
$\displaystyle+\\{2\leftrightarrow 3,\,z\leftrightarrow 1-z\\}\bigg{]}\,.$
(119)

Finally, for the $C_{A}^{2}$ channel we must distinguish between the case in
which $z>z_{\rm cut}$ and $z<z_{\rm cut}$, following the scheme used to define
${\cal B}_{2}^{g}(z)$. However, since the latter is defined in the $z_{\rm
cut}\to 0$ limit, we can entirely neglect the region $z<z_{\rm cut}$ in the
calculation of ${\mathbb{K}}_{g}^{\rm R,\,C_{A}^{2}}[G_{q},G_{g}]$. Using the
phase space parameterisation of Fig. 4 we then find

$\displaystyle{\mathbb{K}}_{g}^{\rm
R,\,C_{A}^{2}}[G_{q},G_{g}]={\mathbb{K}}_{g,\\{12,3\\}}^{\rm
R,\,C_{A}^{2}}[G_{q},G_{g}]+\\{{\rm 2\,cyclic\,permutations}\\}\,,$ (120)

where we defined

$\displaystyle{\mathbb{K}}_{g,\\{12,3\\}}^{\rm R,\,C_{A}^{2}}$
$\displaystyle[G_{q},G_{g}]=\frac{1}{2!}\int
d\Phi^{(A)}_{3}\,P^{C_{A}^{2}}_{1\to
3}\,\bigg{\\{}G_{g}(x\,z_{p}\,(1-z),t_{1,2})\,G_{g}(x\,(1-z_{p})\,(1-z),t_{1,2})$
$\displaystyle\times
G_{g}(x\,z,t_{12,3})-G_{g}(x\,(1-z),t_{12,3})\,G_{g}(x\,z,t_{12,3})\bigg{\\}}\frac{\Delta_{g}(t)}{\Delta_{g}(t_{1,2})}\,\Theta(t_{1,2}-\max\\{t_{1,3},t_{2,3}\\})\,.$
(121)

The additional combinatorial factor $1/2!$ accounts for the remaining two-fold
symmetry in the $g\to ggg$ splitting function.

The second important difference between the gluon and quark cases concerns the
double counting term ${\mathbb{K}}^{\rm 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.1)).
We write it as

${\mathbb{K}}_{g}^{\rm DC}[G_{q},G_{g}]={\mathbb{K}}_{g}^{\rm
DC,\,C_{A}\,T_{R}}[G_{q},G_{g}]+{\mathbb{K}}_{g}^{\rm
DC,\,C_{F}\,T_{R}}[G_{q},G_{g}]+{\mathbb{K}}_{g}^{\rm
DC,\,C_{A}^{2}}[G_{q},G_{g}]\,,$ (122)

where (using the parameterisation of Figs. 1, 3 and 4, respectively)161616We
note that the $C_{A}\,T_{R}$ contribution to the double counting term
${\mathbb{K}}_{g}^{\rm DC,\,C_{A}\,T_{R}}[G_{q},G_{g}]$ is symmetric in
$z\leftrightarrow 1-z$, while ${\mathbb{K}}_{g}^{\rm
R,\,C_{A}\,T_{R}}[G_{q},G_{g}]$ given above is not. This symmetrisation in
${\mathbb{K}}_{g}^{\rm R,\,C_{A}\,T_{R}}[G_{q},G_{g}]$ is not necessary but it
could be performed (i.e. adding the same integrand with $z\leftrightarrow 1-z$
and dividing by 2) to mirror the iteration of the NLL kernel.

$\displaystyle{\mathbb{K}}_{g}^{\rm DC,\,C_{A}\,T_{R}}$
$\displaystyle[G_{q},G_{g}]=\int_{t}^{t_{0}}dt_{1,23}\,dt_{2,3}\int_{z_{0}}^{1-z_{0}}dz\,dz_{p}P_{gg}(z)P_{qg}(z_{p})\,\bigg{[}\bigg{\\{}G_{q}(x\,z_{p}\,(1-z),t_{2,3})$
$\displaystyle\times
G_{q}(x\,(1-z_{p})\,(1-z),t_{2,3})G_{g}(x\,z,t_{1,23})-G_{g}(x\,(1-z),t_{1,23})\,G_{g}(x\,z,t_{1,23})\bigg{\\}}$
$\displaystyle\times\frac{\Delta_{g}(t)}{\Delta_{g}(t_{2,3})}\,\Theta(t_{2,3}-t_{1,23})+\\{z\leftrightarrow
1-z\\}\bigg{]}\,,$ (123) $\displaystyle{\mathbb{K}}_{g}^{\rm
DC,\,C_{F}\,T_{R}}$
$\displaystyle[G_{q},G_{g}]=\int_{t}^{t_{0}}dt_{2,13}\,dt_{1,3}\int_{z_{0}}^{1-z_{0}}dz\,dz_{p}P_{qg}(z)P_{qq}(z_{p})\,\bigg{[}\bigg{\\{}G_{q}(x\,(1-z),t_{2,13})$
$\displaystyle\times
G_{g}(x\,z\,(1-z_{p}),t_{1,3})G_{q}(x\,z\,z_{p},t_{1,3})-G_{q}(x\,(1-z),t_{2,13})\,G_{q}(x\,z,t_{2,13})\bigg{\\}}$
$\displaystyle\times\frac{\Delta_{g}(t)}{\Delta_{g}(t_{1,3})}\,\Theta(t_{1,3}-t_{2,13})+\\{2\leftrightarrow
3,\,z\leftrightarrow 1-z\\}\bigg{]}\,,$ (124)
$\displaystyle{\mathbb{K}}_{g}^{\rm DC,\,C_{A}^{2}}$
$\displaystyle[G_{q},G_{g}]=\int_{t}^{t_{0}}dt_{12,3}\,dt_{1,2}\int_{z_{0}}^{1-z_{0}}dz\,dz_{p}P_{gg}(z)P_{gg}(z_{p})\,\bigg{[}\bigg{\\{}G_{g}(x\,z_{p}\,(1-z),t_{1,2})$
$\displaystyle\times
G_{g}(x\,(1-z_{p})\,(1-z),t_{1,2})G_{g}(x\,z,t_{12,3})-G_{g}(x\,(1-z),t_{12,3})\,G_{g}(x\,z,t_{12,3})\bigg{\\}}$
$\displaystyle\times\frac{\Delta_{g}(t)}{\Delta_{g}(t_{1,2})}\,\Theta(t_{1,2}-t_{12,3})+\\{z\leftrightarrow
1-z\\}\bigg{]}\,.$ (125)

We note that in Eqs. (C.2) and (C.2) the symmetrisation in $z\leftrightarrow
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)$.

##  Appendix D Derivation of
$\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma}{d\theta^{2}dz}$

In this appendix we derive the double-differential distribution in Eq. (138).
It is instructive to start with the NLL case, for which the distribution
$\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma}{d\theta^{2}dz}$ is defined as
(cf. Eq. (7))

$\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma^{(f)}}{d\theta^{2}dz}=\sum_{m=1}^{\infty}\int\,dP^{(f)}_{m}\,\theta^{2}\delta(\theta^{2}-\theta^{2}_{\rm
pass})\delta(z-z_{\rm pass})\,,$ (126)

where $\theta_{\rm pass}$ and $z_{\rm pass}$ denote the angle and momentum
fraction of the largest-angle branching with $1-z_{\rm cut}>z_{\rm
pass}>z_{\rm cut}$. The probabilities $dP^{(f)}_{m}$ are computed with Eq. (2)
using the NLL evolution equation (2.1). 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
$\theta$ and $z$ to the first one that passes the $z_{\rm cut}$ (with $z_{\rm
cut}\ll 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

$G_{g}(x,t)\simeq u\,.$ (127)

With this simplification, the first few terms of the series read

$\displaystyle\int dP^{(q)}_{1}$ $\displaystyle=\Delta_{q}(t)\,,$ (128)
$\displaystyle\int dP^{(q)}_{2}$
$\displaystyle=\Delta_{q}(t)\,\int_{t}^{t_{0}}dt_{1}\int_{z_{0}}^{1-z_{0}}dz_{1}\,P_{qq}(z_{1})\,,$
(129) $\displaystyle\int dP^{(q)}_{3}$
$\displaystyle=\Delta_{q}(t)\,\int_{t}^{t_{0}}dt_{1}\int_{t_{1}}^{t_{0}}dt_{2}\int_{z_{0}}^{1-z_{0}}dz_{1}dz_{2}\,P_{qq}(z_{1})\,P_{qq}(z_{2})\,,$
(130) $\displaystyle...$ $\displaystyle=...$ (131)

The final state with a single parton does not contribute to the distribution
(more precisely it contributes a $\delta(\theta^{2})$ term), while the final
state with two partons (described by $dP_{2}^{(q)}$) contributes according to
the measurement function

$\theta^{2}\delta(\theta^{2}-\theta^{2}_{1})\delta(z-z_{1})\Theta(z_{1}-z_{\rm
cut})\,.$ (132)

Similarly, the state with three partons will involve the measurement function

$\theta^{2}\left(\delta(\theta^{2}-\theta^{2}_{1})\delta(z-z_{1})\Theta(z_{1}-z_{\rm
cut})+\delta(\theta^{2}-\theta^{2}_{2})\delta(z-z_{2})\Theta(z_{\rm
cut}-z_{1})\Theta(z_{2}-z_{\rm cut})\right)\,,$ (133)

where the first term corresponds to a configuration in which the first,
largest-angle emission (declustering) passes the $z_{\rm cut}$ condition,
while the second term corresponds to a configuration in which the first
emission fails the condition and the observable is then defined by the second
branching provided it passes the cut. The configuration in which neither of
the declusterings passes the cut only gives a $\delta(\theta^{2})$
contribution to the observable, and can be discarded for the differential
distribution.

One can then evaluate explicitly the infinite sum in Eq. (126) by factoring
out the contribution of the branching the passes the cut and summing
inclusively over the remaining branchings. If we trade the angular ordering
for a $1/n!$ combinatorial factor, we obtain the formula

$\displaystyle\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma^{(f)}}{d\theta^{2}dz}$
$\displaystyle=\Delta_{q}(t)\,\frac{\alpha_{s}(E^{2}g(z)^{2}\theta^{2})}{2\pi}\,P_{qq}(z)$
$\displaystyle\times\exp\left\\{\int_{t}^{t_{0}}dt^{\prime}\int_{z_{0}}^{1-z_{0}}dz\,P_{qq}(z)\,\left(\Theta(z_{\rm
cut}-z)+\Theta(\theta-\theta^{\prime})\Theta(z-z_{\rm
cut})\right)\right\\}\,\,.$ (134)

We can now evaluate the integral in the exponent and combine it with the
Sudakov $\Delta_{q}(t)$ given in Eq. (5). After taking the limits
$t_{0}\to\infty$ and $z_{0}\to 0$ we obtain

$\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma^{(f)}}{d\theta^{2}dz}=\frac{\alpha_{s}(E^{2}g(z)^{2}\theta^{2})}{2\pi}\,P_{qq}(z)\,e^{-R^{\rm
NLL}_{q}(\theta,z_{\rm cut})}\,,$ (135)

where $t_{\theta}$ denotes the evolution time corresponding to the angle
$\theta$ and

$R^{\rm NLL}_{q}(\theta,z_{\rm
cut})\equiv\int_{0}^{t_{\theta}}dt^{\prime}\int_{z_{\rm
cut}}^{1}dz\,P_{qq}(z)$ (136)

The quantity in the exponent defines the Sudakov radiator featuring in Eq.
(138). The NNLL version of Eq. (135) can be obtained following the same
procedure. We start from Eq. (7), where now the matching coefficient needs to
be evaluated at one loop, that is

$C(\alpha_{s})=1+\frac{\alpha_{s}(E^{2})}{2\pi}C_{q}^{(1)}+{\cal
O}(\alpha_{s}^{2})\,.$ (137)

The constant $C_{q}^{(1)}$ contains hard-virtual corrections to the Born
process, as well as non-logarithmic terms of soft and/or collinear origin that
contribute in the limit $\theta\to 0$.171717Importantly, the coupling in Eq.
(137) is evaluated at the hard scale of the process (i.e. $\mu^{2}\sim
E^{2}$), which formally corresponds to working in a resummation scheme
Catani:2000vq  in which such non-logarithmic terms are evaluated at the hard
scale. In an alternative resummation scheme one may want to calculate some of
these non-logarithmic terms (for instance those of soft and/or collinear
origin) at the low scales (soft or collinear) of the problem, so that their
running generates logarithmic terms at higher perturbative orders.
Consistently with a resummation scheme transformation, this implies that the
higher-order logarithmic terms generated by the running of the coupling must
be subtracted from the anomalous dimension ${\mathcal{B}}^{f}_{2}(z)$ computed
here.

As a second step, we consider the functional derivatives of the NNLL quark
generating functional, starting from the term in the first line of Eq. (9).
The effect of this term is simply to replace $P_{qq}$ in Eq. (135) with the
full inclusive emission probability ${\mathcal{P}}_{q}(z,\theta)$ given in
Eqs. (9), (2.2). We then consider the contribution of the
${\mathbb{K}}_{q}^{\rm finite}$ term, given in Appendix C. By inspecting the
structure of ${\mathbb{K}}_{q}^{\rm finite}$ we observe that the local
counterterms in Eqs. (C.1), (C.1) are defined precisely by means of the same
projection between the $1\to 3$ and the $1\to 2$ phase spaces that was used in
the definition and calculation of ${\mathcal{B}}^{f}_{2}(z)$. Therefore the
observable is defined in exactly the same way for the real contributions as
well as for the corresponding counterterms, giving a zero correction.
Therefore, the contribution of ${\mathbb{K}}_{q}^{\rm finite}$ to the
differential distribution
$\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma^{(f)}}{d\theta^{2}dz}$
vanishes trivially, and one is left with the general form in Eq. (138). The
same derivation holds for gluon jets, leading to the expression

$\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma}{d\theta^{2}dz}=\left(1+\frac{\alpha_{s}(E^{2})}{2\pi}C_{g}^{(1)}\right)e^{-R^{\rm
NNLL}_{g}(\theta,z_{\rm
cut})}\,\frac{\alpha_{s}(E^{2}g(z)^{2}\theta^{2})}{2\pi}\,\left({\cal
P}_{gg}(z,\theta)+{\cal P}_{qg}(z,\theta)\right)\,,$ (138)

with

$R^{\rm NNLL}_{g}(\theta,z_{\rm
cut})\equiv\int_{0}^{t_{\theta}}dt^{\prime}\int_{z_{\rm
cut}}^{1}dz\,\left({\cal P}_{gg}(z,\theta)+{\cal P}_{qg}(z,\theta)\right)\,,$
(139)

where the inclusive emission probabilities are given in Eq. (7). In the
extraction of ${\mathcal{B}}^{g}_{2}(z)$ we have chosen to work with the
fragmentation of a gluon jet defined as one hemisphere in the two-jet limit of
the $H\to gg$ decay in the heavy-top-mass limit. In this case, the one-loop
matching coefficient is given by

$C_{g}^{(1)}=C_{A}\left(4\ln 2\ln
z_{\mathrm{cut}}+\frac{67}{9}-\frac{\pi^{2}}{3}\right)-T_{R}n_{f}\frac{23}{9}\,.$
(140)

##  Appendix E Double-soft end-point contributions to ${\cal B}_{2}^{g}(z)$

In this appendix we discuss how to devise an alternative scheme for
${\mathcal{B}}^{g,C_{A}^{2}}_{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. (25), and
likewise the end-point contribution in Eq. (72), is due to the specific
definition of $z$ and $\theta$ in the kinematic projection that we used to
define ${\mathcal{B}}^{g,C_{A}^{2}}_{2}(z)$. To make the physics clear it is
useful to recapitulate the case of quark jets Dasgupta:2021hbh , in which the
definition of ${\mathcal{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_{F}^{2}$ 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 $\theta$, 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 Dasgupta:2021hbh . An equally valid, albeit tedious, procedure for
defining $z$ and $\theta$ 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
${\mathcal{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_{A}^{2}$ 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. (25). 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 ${\mathcal{B}}^{g}_{2}(z)$ will directly integrate to

$\int_{0}^{1}dz\,{\mathcal{B}}^{g}_{2}(z)=-\gamma^{(2)}_{g}+b_{0}X_{\theta}^{2}\,.$
(141)

A scheme of this type has the advantage of eliminating the end-point
contribution present in Eq. (72), but presents a slight disadvantage in that
it makes the calculation of ${\mathbb{K}}_{q}^{\rm finite}$ more cumbersome
due to the different projections adopted in the double-soft and in the hard-
collinear limits.

Below we present the quantity ${\mathcal{B}}^{g}_{2}(z)$ in this alternate
scheme, which differs from the one given in the article exclusively in the
$C_{A}^{2}$ colour channel. In order to carry out the computation of
$\mathcal{B}_{2}^{g,C_{A}^{2}}(z)$ in the new scheme, we start by separating
the double-soft limit of the splitting function, which satisfies Casimir
scaling and is therefore identical (up to an overall quadratic Casimir
operator) to the quark case. Let gluons $(g_{i},g_{j})$ be soft, i.e.
$z_{i},z_{j}\ll 1$, while parton $k$ is hard $z_{k}\to 1$. In the triple-
collinear limit we have

$\displaystyle\langle\hat{P}_{g_{i}g_{j};p_{k}}\rangle=C_{k}^{2}\langle\hat{P}^{\rm
ind.}_{g_{i}g_{j};p_{k}}\rangle+C_{k}C_{A}\langle\hat{P}^{\rm
corr.}_{g_{i}g_{j};p_{k}}\rangle\ ,$ (142)

where $C_{k}$ is the Casimir of the hard parton $k$. The various functions
read 181818Notice here in particular that we do not send $z_{k}\to 1$ in the
double-soft splitting functions, because the triple-collinear phase space will
be exactly retained when the computation is performed. This choice is made to
ensure the subtraction in Eq. (145) is achieved locally in phase space.

$\displaystyle\langle\hat{P}^{\rm
ind.}_{g_{i}g_{j};p_{k}}\rangle=\frac{4}{s_{ik}s_{jk}}\frac{z_{k}}{z_{i}z_{j}}\
,$ (143)

and

$\displaystyle\langle\hat{P}^{\rm corr.}_{g_{i}g_{j};p_{k}}\rangle$
$\displaystyle=\frac{(1-\epsilon)}{4z_{k}(s_{ik}+s_{jk})^{2}s_{ij}^{2}}\left(2\frac{z_{i}s_{jk}-z_{j}s_{ik}}{z_{i}+z_{j}}\right)^{2}+\frac{1}{s_{ij}s_{ik}}\left(\frac{1}{z_{j}}+\frac{1}{z_{i}+z_{j}}\right)-\frac{1}{2s_{ik}s_{jk}}\frac{z_{k}}{z_{i}z_{j}}$
$\displaystyle+\frac{1}{2(s_{ik}+s_{jk})s_{ij}}\left(\frac{1}{z_{i}}+\frac{1}{z_{j}}-\frac{8}{z_{i}+z_{j}}\right)-\frac{1}{(s_{ik}+s_{jk})s_{ik}}\frac{z_{k}}{z_{i}(z_{i}+z_{j})}+(i\leftrightarrow
j)\ .$ (144)

In the case of three identical gluons, we have to properly symmetrise the
double-soft function in order to account for any pair of gluons becoming soft.
Therefore we define the following subtracted splitting function

$\displaystyle\langle\hat{P}^{\rm
sub.}_{g_{1}g_{2}g_{3}}\rangle\equiv\frac{1}{s_{123}^{2}}\langle\hat{P}_{g_{1}g_{2}g_{3}}\rangle-\langle\hat{P}_{g_{1}g_{2};g_{3}}\rangle-\langle\hat{P}_{g_{1}g_{3};g_{2}}\rangle-\langle\hat{P}_{g_{2}g_{3};g_{1}}\rangle\
\ .$ (145)

Having removed the double-soft limit, we are now ready to repeat our standard
computations using Eq. (145) as our integrand, for which we follow the
sectorisation introduced in Sec. 4.

The results presented in Sec. 6.3 are then modified as follows:

$\left(\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma_{\mathcal{R}}^{(2)}}{d\theta^{2}\,dz}\right)_{z_{3}>z_{\text{cut}}}^{C_{A}^{2},\,{\rm
sub.}}=\left(\frac{C_{A}\alpha_{s}}{2\pi}\right)^{2}\left(z(1-z)\theta^{2}\right)^{-2\epsilon}\bigg{(}H^{\text{sub.}}_{\text{soft-
coll.}}(z,\epsilon)+H^{\text{sub.}}_{\text{coll.}}(z,\epsilon)\\\
+H^{\text{sub.}}_{\text{soft}}(z,\epsilon)+H^{C_{A}^{2},\,\text{sub.}}_{\text{fin.}}(z)\bigg{)}\
,$ (146)

where

$\displaystyle\begin{split}H^{\text{sub.}}_{\text{soft-
coll.}}(z,\epsilon)=&4^{\epsilon}\,\left(z^{-2\epsilon}+(1-z)^{-2\epsilon}\right)\,\left(-2+z(1-z)\right)\left(\frac{1}{\epsilon^{2}}+\frac{2\ln
2}{\epsilon}+4\ln^{2}2-\frac{\pi^{2}}{2}\right)\ ,\\\
H^{\text{sub.}}_{\text{coll.}}(z,\epsilon)=\,&4^{\epsilon}\,p_{gg}(z)\,\left(z^{-2\epsilon}+(1-z)^{-2\epsilon}\right)\bigg{[}\frac{(11-12\ln
2)}{6\epsilon}+\,\left(\frac{67}{9}-\frac{\pi^{2}}{3}-4\ln^{2}2\right)\bigg{]}\\\
&-\frac{11}{6\epsilon}4^{\epsilon}\bigg{(}\frac{(1-z)^{-2\epsilon}}{1-z}+\frac{z^{-2\epsilon}}{z}\bigg{)}+\,4^{\epsilon}\left(z^{-2\epsilon}+(1-z)^{-2\epsilon}\right)\frac{2\ln
2}{z(1-z)\epsilon}\\\
&-\frac{1}{z(1-z)}\bigg{(}\frac{67}{9}-\frac{2\pi^{2}}{3}-8\ln^{2}2\bigg{)}\,,\\\
H^{\text{sub.}}_{\text{soft}}(z,\epsilon)=&-\left(z^{-2\epsilon}+(1-z)^{-2\epsilon}\right)(-2+z(1-z))\bigg{(}\frac{2\ln(2)}{\epsilon}-\frac{h^{\text{sub.}}_{\rm
pass}}{2}\bigg{)}\,,\end{split}$ (147)

where the constant $h^{\text{sub.}}_{\rm soft}$ is the result of a 2-fold
integration and reads

$\displaystyle h^{\text{sub.}}_{\rm soft}\simeq-6.49651797(5)\,.$ (148)

The function $H^{C_{A}^{2},\,\text{sub.}}_{\text{fin.}}$ has the usual
decomposition (54)

$\displaystyle
H^{C_{A}^{2},\,\text{sub.}}_{\text{fin.}}(z)=-\frac{11\ln(2)}{3}\frac{1}{z(1-z)}+G^{\text{sub.}}_{z_{3}>z_{\text{cut}}}(z)\,,$
(149)

and the remainder function $G^{\text{sub.}}_{z_{3}>z_{\text{cut}}}$, see Fig.
9, is integrable in $z\in[0,1]$ and its integral reads

$\displaystyle\int_{0}^{1}dz\,G^{\text{sub.}}_{z_{3}>z_{\text{cut}}}(z)\simeq
6.974\pm 0.001\,.$ (150)

Figure 9: Plot of the function $G^{\text{sub.}}_{z_{3}>z_{\text{cut}}}(z)$.

Similarly, below the soft threshold $z_{3}<z_{\rm cut}$ we have

$\left(\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma_{\mathcal{R}}^{(2)}}{d\theta^{2}\,dz}\right)_{z_{3}<z_{\text{cut}}}^{C_{A}^{2},\,\text{sub.}}=\left(\frac{C_{A}\alpha_{s}}{2\pi}\right)^{2}\left(-2+z(1-z)\right)\Theta\left(z-z_{\rm
cut}\right)\Theta\left(1-z_{\rm cut}-z\right)\times\\\ \times
4^{-\epsilon}(z(1-z))^{-2\epsilon}\,(\theta^{2})^{-\epsilon}\,\
\bigg{(}\frac{z^{-2\epsilon}_{\rm
cut}}{-\epsilon}\,\ln\frac{4}{\theta^{2}}-\frac{\pi^{2}}{6}-\frac{1}{2}\ln^{2}\frac{4}{\theta^{2}}+h^{p_{gg}}_{\rm
fail}\bigg{)}\,,$ (151)

where the constant $h^{p_{gg}}_{\rm fail}$ is given in Eq. (59). The most
important feature of eq. (151) is that its $z$-dependence is fully regular
over $z\in[0,1]$ unlike Eq. (58) 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. Dasgupta:2021hbh  by enforcing the
appropriate limit. Starting with the correlated portion of the double-soft
function, eq. (E), we find

$\displaystyle\left(\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma_{\mathcal{R}}^{(2)}}{d\theta^{2}\,dz}\right)^{C_{A}^{2},\,\text{corr.}}=$
$\displaystyle\left(\frac{C_{A}\alpha_{s}}{2\pi}\right)^{2}\left(z(1-z)\theta^{2}\right)^{-2\epsilon}\left(\frac{1}{\epsilon^{2}}+\frac{11}{6\epsilon}-\frac{2\pi^{2}}{3}+\frac{67}{18}\right)\times$
$\displaystyle\times\left(\frac{(1-z)^{-2\epsilon}}{1-z}+\frac{z^{-2\epsilon}}{z}\right)\Theta\left(z-z_{\rm
cut}\right)\Theta\left(1-z_{\rm cut}-z\right)\,,$ (152)

while for the independent emission contribution we get191919The third line of
Eq. (E) emerges when the gluon at larger angle fails the $z_{\rm cut}$
condition. Although this contribution can not be extracted directly from Ref.
Dasgupta:2021hbh , it entails a straightforward computation.

$\displaystyle\left(\frac{\theta^{2}}{\sigma_{0}}\frac{d^{2}\sigma_{\mathcal{R}}^{(2)}}{d\theta^{2}\,dz}\right)^{C_{A}^{2},\,\text{ind.}}=$
$\displaystyle\left(\frac{C_{A}\alpha_{s}}{2\pi}\right)^{2}\left(z(1-z)\theta^{2}\right)^{-2\epsilon}\left(\frac{1}{\epsilon^{2}}-\frac{5\pi^{2}}{6}\right)\left(\frac{z^{-2\epsilon}}{1-z}+\frac{(1-z)^{-2\epsilon}}{z}\right)\times$
$\displaystyle\times\Theta\left(z-z_{\rm cut}\right)\Theta\left(1-z_{\rm
cut}-z\right)$ $\displaystyle+$
$\displaystyle\left(\frac{C_{A}\alpha_{s}}{2\pi}\right)^{2}\frac{1}{z(1-z)}\Theta\left(z-z_{\rm
cut}\right)\Theta\left(1-z_{\rm cut}-z\right)\times$ $\displaystyle\times
4^{-\epsilon}(z(1-z))^{-2\epsilon}\,(\theta^{2})^{-\epsilon}\,\
\bigg{(}\frac{z^{-2\epsilon}_{\rm
cut}}{-\epsilon}\,\ln\frac{4}{\theta^{2}}-\frac{\pi^{2}}{6}-\frac{1}{2}\ln^{2}\frac{4}{\theta^{2}}\bigg{)}.$
(153)

Adding all the results in Eqs. (146), (151), (E), (E), and the relevant
virtual corrections, the alternative scheme yields

$\displaystyle\mathcal{B}_{2}^{g,C_{A}^{2},\,\text{alt.}}\left(z\right)=G^{\text{sub.}}_{z_{3}>z_{\text{cut}}}(z)+\mathcal{B}^{g,\,C_{A}^{2},\,\rm{analytic,\,alt.}}_{2}\left(z\right)\,,$
(154)

where

$\displaystyle\mathcal{B}^{g,\,C_{A}^{2},\,\rm{analytic,\,alt.}}_{2}\left(z\right)$
$\displaystyle=\frac{1}{18}\left(2-z(1-z)\right)\left(-18\,h^{\text{sub.}}_{\rm
soft}-18\,h^{p_{gg}}_{\rm fail}+6\pi^{2}-132\ln 2-72\ln^{2}2\right)$
$\displaystyle+\frac{-265+134z-134z^{2}}{18}+\frac{11\left(1-3z+2z^{2}-z^{3}\right)}{6(1-z)}\ln
z$
$\displaystyle+\frac{11\left(-1+2z-z^{2}+z^{3}\right)}{6z}\ln(1-z)-\frac{11}{6}\frac{\ln
z}{1-z}-\frac{11}{6}\frac{\ln(1-z)}{z}$ $\displaystyle-2\,p_{gg}(z)\ln(1-z)\ln
z+\frac{11}{6}(2-z(1-z))\ln(z(1-z))\,.$ (155)

The integral of $\mathcal{B}_{2}^{g,C_{A}^{2},\,\text{alt.}}(z)$ then reads

$\displaystyle\int_{0}^{1}dz\,\mathcal{B}_{2}^{g,C_{A}^{2},\,\text{alt.}}(z)\simeq-7.858\pm
0.001\,.$ (156)

It is easy to verify that eq. (156) agrees with the sum rule given in Eq.
(141).

## References

  * (1) M. van Beekveld, M. Dasgupta, B. El-Menoufi, J. Helliwell, P. Monni and G. P. Salam, _In preparation_ , .
  * (2) K. Konishi, A. Ukawa and G. Veneziano, _Jet Calculus: A Simple Algorithm for Resolving QCD Jets_ , _Nucl. Phys. B_ 157 (1979) 45–107.
  * (3) A. Bassetto, M. Ciafaloni and G. Marchesini, _Jet Structure and Infrared Sensitive Quantities in Perturbative QCD_ , _Phys. Rept._ 100 (1983) 201–272.
  * (4) Y. L. Dokshitzer, V. A. Khoze, A. H. Mueller and S. I. Troian, _Basics of perturbative QCD_.  1991\.
  * (5) R. K. Ellis, W. J. Stirling and B. R. Webber, _QCD and collider physics_ , vol. 8.  Cambridge University Press, 2, 2011.  10.1017/CBO9780511628788.
  * (6) J. M. Campbell et al., _Event Generators for High-Energy Physics Experiments_ , in _Snowmass 2021_ , 3, 2022.  2203.11110.
  * (7) M. Dasgupta, F. Dreyer, G. P. Salam and G. Soyez, _Small-radius jets to all orders in QCD_ , _JHEP_ 04 (2015) 039, [1411.5182].
  * (8) Z.-B. Kang, F. Ringer and I. Vitev, _The semi-inclusive jet function in SCET and small radius resummation for inclusive jet production_ , _JHEP_ 10 (2016) 125, [1606.06732].
  * (9) A. Jain, M. Procura and W. J. Waalewijn, _Parton Fragmentation within an Identified Jet at NNLL_ , _JHEP_ 05 (2011) 035, [1101.4953].
  * (10) H.-M. Chang, M. Procura, J. Thaler and W. J. Waalewijn, _Calculating Track-Based Observables for the LHC_ , _Phys. Rev. Lett._ 111 (2013) 102002, [1303.6637].
  * (11) M. Ritzmann and W. J. Waalewijn, _Fragmentation in Jets at NNLO_ , _Phys. Rev. D_ 90 (2014) 054029, [1407.3272].
  * (12) B. T. Elder, M. Procura, J. Thaler, W. J. Waalewijn and K. Zhou, _Generalized Fragmentation Functions for Fractal Jet Observables_ , _JHEP_ 06 (2017) 085, [1704.05456].
  * (13) H. Chen, I. Moult and H. X. Zhu, _Quantum Interference in Jet Substructure from Spinning Gluons_ , _Phys. Rev. Lett._ 126 (2021) 112003, [2011.02492].
  * (14) D. Neill and F. Ringer, _Soft Fragmentation on the Celestial Sphere_ , _JHEP_ 06 (2020) 086, [2003.02275].
  * (15) H. Chen, I. Moult and H. X. Zhu, _Spinning gluons from the QCD light-ray OPE_ , _JHEP_ 08 (2022) 233, [2104.00009].
  * (16) D. Neill, F. Ringer and N. Sato, _Leading jets and energy loss_ , _JHEP_ 07 (2021) 041, [2103.16573].
  * (17) H. Chen, I. Moult, J. Sandor and H. X. Zhu, _Celestial blocks and transverse spin in the three-point energy correlator_ , _JHEP_ 09 (2022) 199, [2202.04085].
  * (18) H. Chen, M. Jaarsma, Y. Li, I. Moult, W. J. Waalewijn and H. X. Zhu, _Collinear Parton Dynamics Beyond DGLAP_ , 2210.10061.
  * (19) W. Chen, J. Gao, Y. Li, Z. Xu, X. Zhang and H. X. Zhu, _NNLL Resummation for Projected Three-Point Energy Correlator_ , 2307.07510.