# Anti-$k_{T}$ jet function at next-to-next-to-leading order

Hao-Yu Liu<EMAIL_ADDRESS>Center of Advanced Quantum Studies, Department
of Physics, Beijing Normal University, Beijing 100875, China  Xiaohui Liu
<EMAIL_ADDRESS>Center of Advanced Quantum Studies, Department of Physics,
Beijing Normal University, Beijing 100875, China  Center for High Energy
Physics, Peking University, Beijing 100871, China  Sven-Olaf Moch  sven-
<EMAIL_ADDRESS>II. Institut für Theoretische Physik Universität Hamburg,
Luruper Chaussee 149, D-22761 Hamburg, Germany

###### Abstract

Jets constructed via clustering algorithms (e.g., anti-$k_{T}$, soft-drop)
have been proposed for many precision measurements, such as the strong
coupling $\alpha_{s}$ and the nucleon intrinsic dynamics. However, the
theoretical accuracy is affected by missing QCD corrections at higher orders
for the jet functions in the associated factorization theorems. Their
calculation is complicated by the jet clustering procedure. In this work, we
propose a method to evaluate jet functions at higher orders in QCD. The
calculation involves the phase space sector decomposition with suitable soft
subtractions. As a concrete example, we present the quark-jet function using
the anti-$k_{T}$ algorithm with E-scheme recombination at next-to-next-to-
leading order.

††preprint: DESY 21-032

##  I Introduction

Jets have been the center of attention at many frontiers of high-energy
particle and nuclear physics. At the Large Hadron Collider (LHC), for many
years, jets have served as indispensable tools for precision test of the
Standard Model (SM) and new physics searches at the TeV scale. Measurements at
the LHC and the Relativistic Heavy Ion Collider (RHIC) demonstrate that jets
can also be unique probes of non-perturbative dynamics, such as collinear
parton distribution functions (PDFs) Accardi:2016ndt  as well as transverse
momentum dependent ones (TMDPDFs) Angeles-Martinez:2015sea , the intrinsic
spin of the nucleon Boer:2014lka ; Aschenauer:2016our ; Adamczyk:2017wld  or
the hot medium effects of the quark-gluon plasma Chien:2016led ;
Connors:2017ptx ; Qiu:2019sfj . These efforts are supported by theory
developments of novel jet substructures Larkoski:2014wba ; Larkoski:2017jix ,
and enable accurate extractions of the SM parameters Larkoski:2017jix ;
Kardos:2020ppl  and the three-dimensional tomographic images of nucleons
Gutierrez-Reyes:2018qez ; Liu:2018trl ; Liu:2020dct ; Arratia:2020nxw ;
Kang:2020fka ; Kang:2020xez ; Chien:2020hzh  out of the jets. All these
studies will receive a further boost at the future Electron-Ion Collider (EIC)
Boer:2011fh ; Accardi:2012qut , including the possibility to consider
polarization degrees of freedom.

The theoretical foundation of such precision measurements roots in the
factorization of the cross section $d\sigma$ with $N$ exclusive jets with
small jet radius $R$, whose schematic form is Becher:2015hka ; Becher:2016mmh
; Liu:2017pbb

$\displaystyle d\sigma={\cal F}_{a}{\cal
F}_{b}\operatorname{Tr}[HS_{G}]\,\prod_{c}^{N}\sum_{m}\operatorname{Tr}[J^{c}_{m}\otimes_{\Omega}S^{c}_{cs,m}]\,,$
(1)

where ${\cal F}$ encodes the parton distributions of the colliding nucleons,
such as the collinear PDFs or TMDPDFs, in case transverse momentum dependence
is kept. The traces are over the color matrices. $H$ is the hard function that
describes the hard interaction which initiates the process while $S_{G}$
includes the wide angle radiations with long wavelength. The jet function
containing $m$ partons is given by Becher:2015hka ; Becher:2016mmh

$\displaystyle J_{m}^{c}$ $\displaystyle=$ $\displaystyle\int_{\sum
p^{+,\perp}_{i}=p^{+,\perp}_{J}}\prod_{i=1}^{m}\frac{dE_{i}E_{i}^{d-3}}{2(2\pi)^{d-1}}$
(2) $\displaystyle\times$ $\displaystyle\Theta_{\text{jet alg.}}\sum_{\rm
spins}|{\cal M}_{m}\rangle\langle{\cal M}_{m}|\,,$

where $|{\cal M}_{m}\rangle$ is the amplitude for the collinear field with
momentum $p^{\mu}_{J}\approx\frac{p^{+}_{J}}{2}n^{\mu}_{J}$, with $n_{J}$ a
light-like vector along the jet axis, to split into $m$ particles with momenta
$\\{p_{1}\,,\dots\,,p_{m}\\}$ and $\Theta_{\text{jet alg.}}$ enforces that all
the $m$ collinear particles are clustered within one jet.

The measurements of jets and jet substructures are formulated in the angular
convolution $\otimes_{\Omega}$ of the energetic collinear jet function $J_{m}$
and the collinear-soft function $S_{cs,m}$, where $m$ is the energetic
collinear parton multiplicity. The angular convolution guarantees the proper
inclusion of the non-global logarithms (NGLs) Dasgupta:2001sh  if present. We
remark that the multi-emission feature of the jet cross section in Eq. (1) is
also shared by the forward scattering in the small-$x$ framework Liu:2020mpy .
The necessity for increasing theoretical precision for jet cross sections,
both at fixed order perturbation theory and with resummation, requires the
evaluation of each component in the factorization Eq. (1) at least to next-to-
next-to-leading order (NNLO).

On the other hand, many theoretical predictions of the jet measurements are
currently based on the simplified formula Liu:2012sz ; Liu:2013hba ;
Sun:2014gfa ; Liu:2018trl ; Liu:2018ktv ; Liu:2020dct ; Liu:2020jjv ;
Arratia:2020nxw

$\displaystyle d\sigma={\cal F}_{a}{\cal
F}_{b}\operatorname{Tr}[HS_{G}]\,\prod_{c}^{N}\,J^{c}\,S^{c}_{cs}\,e^{L_{\rm
ngl}}\,,$ (3)

where one decouples the angular correlation between $J_{m}$ and $S_{cs,m}$ and
averages each function over its solid angle in Eq. (2) to have
$J=\sum_{m}\langle J_{m}\rangle_{\Omega}$ and $S_{cs}=\langle
S_{cs,1}\rangle_{\Omega}$. As indicated by the factor $e^{L_{\rm ngl}}$, the
leading NGLs are resummed through the exponential. In order to reliably
justify the logarithmic accuracy of the simplified form, a direct NNLO
calculation of the factorization is also desired.

For many processes, the hard and the soft functions are already known to two
loops or even beyond, see, e.g. heinrich2020collider . However the two-loop
information for the jet function is still missing and therefore limits the
theoretical accuracy of the resummation to the next-to-leading logarithm
(NLL). Compared with the other components in the factorization, the
computation of the jet function beyond next-to-leading order (NLO) is
dramatically complicated by the presence of the recursive clustering
procedure, which is the essence of the widely used algorithms for jets and jet
substructure, such as anti-$k_{T}$ Cacciari_2008  with either the E-scheme or
the winner-take-all (WTA) recombination scheme Bertolini:2013iqa ;
Neill:2016vbi  and soft-drop grooming Larkoski:2014wba .

In this work, we propose a method to evaluate $J_{m}$ or $J$ at NNLO
efficiently using sector decomposition for the phase space together with
suitable designed soft subtractions. We demonstrate the feasibility of the
method by presenting the first anti-$k_{T}$ quark-jet function at ${\cal
O}(\alpha_{s}^{2})$ using E-scheme recombination as a concrete example.

##  II NLO

We take the well-known NLO calculation Liu:2012sz ; Arratia:2020nxw  of the
quark-jet function as a warm up to illustrate our approach. The NLO
calculation involves evaluating the integration of the matrix element for
$q_{a}\to q_{i}g_{j}$ over the collinear two-body phase space, which reads in
$d=4-2\epsilon$ dimensions

$\displaystyle J_{bare}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\frac{1}{(2\pi)^{d-1}}\frac{2\pi^{1-\epsilon}}{\Gamma(1-\epsilon)}\int\mathrm{d}z\mathrm{d}s_{ij}s_{ij}^{-\epsilon}(z{\bar{z}})^{1-\epsilon}$
(4)
$\displaystyle\times\frac{8\pi\alpha_{s}Z_{\alpha}\,\mu^{2\epsilon}e^{\gamma_{E}\epsilon}}{(4\pi)^{\epsilon}s_{ij}}C_{F}\left[\frac{1+{\bar{z}}^{2}}{z}-\epsilon
z\right]\theta(R^{2}-\Delta R_{ij}^{2})\,,$

where $Z_{\alpha}$ is the $\alpha_{s}$ renormalization factor, which at NLO is

$\displaystyle
Z_{\alpha}=1-\frac{\alpha_{s}}{4\pi\epsilon}\beta_{0}=1-\frac{\alpha_{s}}{4\pi\epsilon}\left(\frac{11}{3}C_{A}-\frac{4}{3}T_{F}N_{F}\right)\,,$
(5)

$s_{ij}=2p_{i}\cdot p_{j}$ and $z$ is the momentum fraction carried by the
gluon $g_{j}$ and ${\bar{z}}=1-z$. The first line of Eq. (4) is the phase
space measure in the collinear limit and we note that the matrix element in
second line of Eq. (4) is nothing but the leading order (LO) splitting
function $\frac{8\pi\alpha_{s}\mu^{2\epsilon}}{s_{ij}}P_{q\to qg}$. The
$\theta$-function originates from the requirements of the jet algorithm to
cluster two partons into one single jet with radius $R$. Here

$\displaystyle\Delta
R_{ij}^{2}=\Delta\eta_{ij}^{2}+\Delta\phi_{ij}^{2}\approx\frac{2p_{i}\cdot
p_{j}}{p_{i,T}p_{j,T}}=\frac{s_{ij}}{z{\bar{z}}p_{T}^{2}}\,,$ (6)

where $\Delta\eta_{ij}$ and $\Delta\phi_{ij}$ are the rapidity and azimuthal
angle differences between parton $i$ and $j$, respectively. $p_{T}$ is the jet
transverse momentum and $p_{i,T}$ and $p_{j,T}$ are the transverse momenta for
partons $i$ and $j$. The narrow jet approximation $R\ll 1$ has been used.

If we introduce the variables

$\displaystyle
x_{1}\equiv\tilde{s}_{ij}=\frac{s_{ij}}{z{\bar{z}}(p_{T}R)^{2}}\leq 1\,,\quad
x_{2}\equiv z\leq 1\,,$ (7)

we can write

$\displaystyle J_{bare}^{(1)}$ $\displaystyle=$ $\displaystyle e^{2\epsilon
L}\frac{\alpha_{s}C_{F}}{2\pi}\frac{Z_{\alpha}\,e^{\gamma_{E}\epsilon}}{\Gamma(1-\epsilon)}\int_{0}^{1}\mathrm{d}x_{1}\mathrm{d}x_{2}x_{1}^{-1-\epsilon}x_{2}^{-1-2\epsilon}$
(8) $\displaystyle\times$
$\displaystyle(1-x_{2})^{-2\epsilon}\,\left[1+(1-x_{2})^{2}-\epsilon\,x_{2}^{2}\right]\,,$

where $L=\log\frac{\mu}{p_{T}R}$. All the singularities are given by $x_{i}\to
0$ and have been isolated in the factor
$x_{1}^{-1-\epsilon}x_{2}^{-1-2\epsilon}$ so that the remaining terms in the
second line of the above equation are finite as $x_{i}\to 0$.

With this manipulation, to evaluate the integral, we can first expand
$x_{i}^{-1-a_{i}\epsilon}$ using the Laurent expansion

$\displaystyle
x_{i}^{-1-a_{i}\epsilon}=-\frac{1}{a_{i}\epsilon}\delta(x_{i})+\sum_{n=0}\frac{(-a_{i}\epsilon)^{n}}{n!}\left[\frac{\log^{n}x_{i}}{x_{i}}\right]_{+}\,,\quad$
(9)

and then perform the integration either numerically or analytically to find
the coefficients of the $\epsilon$-poles and the finite terms. This reproduces
the well-known NLO bare quark-jet function

$\displaystyle J_{bare}^{(1)}$ $\displaystyle=$ $\displaystyle e^{2\epsilon
L}Z_{\alpha}\,\frac{\alpha_{s}}{2\pi}C_{F}\left(\frac{1}{\epsilon^{2}}+\frac{3}{2\epsilon}+\frac{13}{2}-\frac{3\pi^{2}}{4}\right.$
(10)
$\displaystyle+\left[26-\frac{9\pi^{2}}{8}-\frac{49}{3}\zeta_{3}\right]\epsilon$
$\displaystyle\left.+\left[104-\frac{39}{8}\pi^{2}-\frac{49}{2}\zeta_{3}-\frac{11}{32}\pi^{4}\right]\epsilon^{2}\right)\,.\quad$

We note that at NLO, all $k_{T}$-type jet algorithms give the same results.

Although we present here only the $x_{i}$-integrated jet function, our
approach is actually fully exclusive in $x_{i}$. There are no difficulties to
generate distributions differential in $x_{i}$, as can be seen from Eq. (9),
especially in the angular variable $x_{1}=\tilde{s}_{ij}$. The implications of
the possibility to provide exclusive results, especially their role in
resumming non-global logarithms, will be discussed later.

A similar strategy will be used to compute the NNLO jet functions but the
calculation is more involved, as we will explained in detail in the rest of
the work.

##  III NNLO

###  III.1 Real-virtual contribution

Now we proceed to calculate the jet function at NNLO. We start with the real-
virtual corrections. The phase space of the real-virtual corrections is
identical to the NLO calculation. The matrix element for the jet function is
given by the one-loop correction to the splitting kernel, which is well
documented in Kosower:1999xi ; Kosower:1999rx ; Bern:1999ry ; Sborlini:2013jba
; Ritzmann:2014mka . The calculation for the real-virtual correction is
straightforward and gives

$\displaystyle J^{(2)}_{rv}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{s}^{2}e^{4\epsilon
L}}{(2\pi)^{2}}C_{F}\Big{(}C_{F}{\cal K}^{rv}_{C_{F}}+C_{A}{\cal
K}^{rv}_{C_{A}}\Big{)}\,,$ (11)

with

$\displaystyle{\cal K}^{rv}_{C_{F}}$ $\displaystyle=$
$\displaystyle\left(-\frac{5}{4}+\frac{\pi^{2}}{3}\right)\frac{1}{\epsilon^{2}}+\left(-\frac{31}{2}+\frac{\pi^{2}}{2}+22\zeta_{3}\right)\frac{1}{\epsilon}$
(12)
$\displaystyle-\frac{575}{4}+\frac{137}{24}\pi^{2}+33\zeta_{3}+\frac{10}{9}\pi^{4}\,,$

and

$\displaystyle{\cal K}^{rv}_{C_{A}}$ $\displaystyle=$
$\displaystyle-\frac{1}{4\epsilon^{4}}-\frac{3}{4\epsilon^{3}}+\left(-5+\frac{11\pi^{2}}{24}\right)\frac{1}{\epsilon^{2}}$
(13)
$\displaystyle+\left(-\frac{63}{2}+\frac{13\pi^{2}}{8}+\frac{26}{3}\zeta_{3}\right)\frac{1}{\epsilon}$
$\displaystyle-\frac{781}{4}+11\pi^{2}+\frac{85}{2}\zeta_{3}-\frac{67}{1440}\pi^{4}\,.\qquad$

###  III.2 Real-real contribution

####  III.2.1 Matrix elements

We turn to the calculation of the real-real corrections. The matrix element
involved is nothing but the tree level $a\to ijk$ splitting kernel,

$\displaystyle|{\cal
M}|^{2}=\left(\frac{\mu^{2}e^{\gamma_{E}}}{4\pi}\right)^{2\epsilon}\,\frac{64\pi^{2}\alpha_{s}^{2}}{s^{2}_{ijk}}P_{a\to
ijk}(z_{i},z_{j},z_{k})\,,$ (14)

where $s_{ijk}=s_{ij}+s_{ik}+s_{jk}$ and $z_{i}$ is the momentum fraction
carried by parton $i$, with $z_{i}+z_{j}+z_{k}=1$. The explicit form of the
splitting function $P_{a\to ijk}(z_{i},z_{j},z_{k})$ can be found in
Catani:1999ss  (see also Campbell:1997hg ). For the quark-jet we need

$\displaystyle P_{{\bar{q}}_{1}^{\prime}q_{2}^{\prime}q_{3}}\,,\quad
P_{{\bar{q}}_{1}q_{2}q_{3}}^{({\rm id})}\,,\quad P_{g_{1}g_{2}q_{3}}^{({\rm
ab})}\,,\quad\text{and}\quad P_{g_{1}g_{2}q_{3}}^{(\rm nab)}\,,$ (15)

and explicit expressions are listed in the appendix A.

The matrix element $|{\cal M}|^{2}$ contains structures

$\displaystyle|{\cal
M}|^{2}\supset\frac{1}{s_{ijk}^{\alpha}s_{ij}^{\beta}}\,,\frac{1}{s_{ij}s_{jk}}\,,\text{permutations
in $i$, $j$, $k$}\,,\quad$ (16)

which will develop singularities in $\epsilon$ if one or more of the
$s_{ij}$’s vanish upon integration over the phase space. To extract the real-
real corrections, our strategy is to write the phase space integration as

$\displaystyle\int d\Phi_{3}|{\cal
M}|^{2}=\int_{0}^{1}\prod_{i}\mathrm{d}x_{i}\,x_{i}^{-1-a_{i}\epsilon}\,F(\\{x_{i}\\},\epsilon)\,,$
(17)

where $F$ is regular as $x_{i}\to 0$. Hence we can first perform the Laurent
expansion in $x_{i}^{-1-a_{i}\epsilon}$ to isolate the $\epsilon$-poles and
compute their coefficients and the remaining constant terms either numerically
or analytically. We emphasize again that Eq. (17) allows not only for the
integrated jet function but also the differential distributions in $x_{i}$’s.
To arrive at the expression in Eq. (17), we need a proper parameterization of
the phase space.

####  III.2.2 Phase Space Parameterization

In the collinear limit, the three-body phase space can be shown to take the
form GehrmannDeRidder:1997gf ; Ritzmann:2014mka

$\displaystyle d\Phi_{3}$ $\displaystyle=$ $\displaystyle
4\frac{ds_{ij}ds_{ik}ds_{jk}dz_{i}dz_{j}}{(4\pi)^{5-2\epsilon}\Gamma(1-2\epsilon)}\,\Delta^{-\frac{1}{2}-\epsilon}\,,$
(18)

where the Gram determinant is given by

$\displaystyle\Delta=4z_{i}z_{j}s_{ik}s_{jk}-(z_{k}s_{ij}-z_{i}s_{jk}-z_{j}s_{ik})^{2}>0\,,$
(19)

with $z_{k}=1-z_{i}-z_{j}$.

Now we introduce the angular variable between two arbitrary partons $a$ and
$b$

$\displaystyle\tilde{s}_{ab}=\frac{1}{(p_{T}R)^{2}}\frac{s_{ab}}{z_{a}z_{b}}$
(20)

and further make a variable transformation using one of the $\tilde{s}_{ab}$

$\displaystyle\tilde{s}_{ik}=(\sqrt{\tilde{s}_{jk}}-\sqrt{\tilde{s}_{ij}})^{2}+4\sqrt{\tilde{s}_{ij}\tilde{s}_{jk}}t\,,$
(21)

with $t\in[0,1]$ to find

$\displaystyle d\Phi_{3}$ $\displaystyle=$
$\displaystyle(2p_{T}R)^{4-4\epsilon}\frac{d\tilde{s}_{ij}d\tilde{s}_{jk}dtdz_{i}dz_{j}}{2(4\pi)^{5-2\epsilon}\Gamma(1-2\epsilon)}(z_{i}z_{j}z_{k})^{1-2\epsilon}$
(22)
$\displaystyle\times\left(\tilde{s}_{ij}\tilde{s}_{jk}\right)^{-\epsilon}t^{-\frac{1}{2}-\epsilon}\left(1-t\right)^{-\frac{1}{2}-\epsilon}\,.$

If $\tilde{s}_{ik}$ happens to be in the denominator of the matrix element,
the variable change in Eq. (21) will develop a line singularity, which can be
mapped to $\tilde{s}_{ij}=\tilde{s}_{jk}$ by a non-linear transformation
Anastasiou:2003gr ; Boughezal:2015eha

$\displaystyle
x_{5}^{\prime}=\frac{(\sqrt{\tilde{s}_{ij}}-\sqrt{\tilde{s}_{jk}})^{2}(1-t)}{\tilde{s}_{ik}}\,,$
(23)

and thus

$\displaystyle\tilde{s}_{ik}=(\tilde{s}_{ij}-\tilde{s}_{jk})^{2}\left((\sqrt{\tilde{s}_{ij}}-\sqrt{\tilde{s}_{jk}})^{2}+4\sqrt{\tilde{s}_{ij}\tilde{s}_{jk}}\,x_{5}^{\prime}\right)^{-1}\,,$

which gives

$\displaystyle d\Phi_{3}$ $\displaystyle=$
$\displaystyle(2p_{T}R)^{4-4\epsilon}\,\frac{\pi\,d\tilde{s}_{ij}d\tilde{s}_{jk}dz_{i}dz_{j}dx_{5}}{2(4\pi)^{5-2\epsilon}\Gamma(1-2\epsilon)}(z_{i}z_{j}z_{k})^{1-2\epsilon}$
(25) $\displaystyle\times$
$\displaystyle\left(\tilde{s}_{ij}\tilde{s}_{jk}\right)^{-\epsilon}\,{x_{5}^{\prime}}^{-\epsilon}(1-x_{5}^{\prime})^{-\epsilon}\,|\tilde{s}_{ij}-\tilde{s}_{jk}|^{1-2\epsilon}$
$\displaystyle\times$
$\displaystyle\left((\sqrt{\tilde{s}_{ij}}-\sqrt{\tilde{s}_{jk}})^{2}+4\sqrt{\tilde{s}_{ij}\tilde{s}_{jk}}\,x_{5}^{\prime}\right)^{-1+2\epsilon}\,.$

Here we have introduced the variable $x_{5}$ through

$\displaystyle x_{5}^{\prime}=\sin^{2}\left(\frac{\pi}{2}x_{5}\right)\,.$ (26)

Without loss of generality, we can assume $z_{i}\leq z_{j}$ and
$\tilde{s}_{ij}\leq\tilde{s}_{jk}$, and in this case we can then introduce
111It is understood that if $x_{i}$ is not in $[0,1]$ but $[0,a]$ instead, we
can do $x_{i}\to ax_{i}$ to map it onto $[0,1]$

$\displaystyle\tilde{s}_{jk}=x_{1}\,,\quad\tilde{s}_{ij}=x_{1}x_{2}\,,\quad
z_{j}=x_{3}\,,\quad z_{i}=x_{3}x_{4}\,,\quad$ (27)

to reach the phase space parameterization

$\displaystyle d\Phi_{3}$ $\displaystyle=$
$\displaystyle(2p_{T}R)^{-4\epsilon}\frac{\pi\,dx_{1}dx_{2}dx_{3}dx_{4}dx_{5}}{2(4\pi)^{5-2\epsilon}\Gamma(1-2\epsilon)}\times
z_{k}$ $\displaystyle\times$
$\displaystyle\left[z_{k}^{2}\,\,x^{\prime}_{5}(1-x^{\prime}_{5})\,\left((\sqrt{x_{2}}-1)^{2}+4\sqrt{x_{2}}\,x_{5}^{\prime}\right)^{-2}\right]^{-\epsilon}$
$\displaystyle\times$ $\displaystyle
x_{1}^{-1-2\epsilon}\,x_{2}^{-1-\epsilon}\,(1-x_{2})^{-1-2\epsilon}\,x_{3}^{-1-4\epsilon}\,\,x_{4}^{-1-2\epsilon}\,$
$\displaystyle\times$
$\displaystyle\left[x_{1}^{2}x_{2}(1-x_{2})^{2}x_{3}^{4}x_{4}^{2}((1-\sqrt{x_{2}})^{2}+4\sqrt{x_{2}}x_{5}^{\prime})^{-1}\right]\,.$

Other orders can be obtained in the same manner. In the third line of the
phase space parameterization, the Laurent expansion in Eq. (9) can be
performed to extract the $\epsilon$-poles and the finite contributions. On the
other hand, the last line of Eq. (III.2.2), when multiplied with the matrix
element, will be free of any singularities as $x_{i}\to 0$ or $x_{2}\to 1$. It
is suitable for a numerical evaluation, as long as appropriate subtraction
terms are constructed as will be explained in the following section.

####  III.2.3 Jet algorithm and subtraction terms

In this work, we are always interested in $k_{T}$-type jet algorithms, which
involves the comparison of the metrics

$\displaystyle\rho_{ij}=\min\left[p_{T,i}^{-2\alpha},p_{T,j}^{-2\alpha}\right]\frac{\Delta
R_{ij}^{2}}{R^{2}}\,,\quad\rho_{i}=p_{T,i}^{-2\alpha}\,,$ (29)

for some parameter $\alpha$ and $\alpha=1$ defines the anti-$k_{T}$ jet
algorithm.

In the small jet radius limit, or in the collinear limit, we can approximate
$p_{T,i}=z_{i}\,p_{T}$ up to power corrections. Furthermore, if we demand all
three partons to be grouped into the same jet, we will require in the last
stage of the clustering,

$\displaystyle\Delta R_{\tilde{ij},k}$ $\displaystyle=$
$\displaystyle\Delta\eta_{\tilde{ij},k}^{2}+\Delta\phi_{\tilde{ij},k}^{2}$
(30) $\displaystyle\approx$ $\displaystyle
2(\cosh\Delta\eta_{\tilde{ij},k}-\cos\Delta\phi_{\tilde{ij},k})\leq R^{2}\,,$

where we have assumed that the partons $i$ and $j$ are clustered first into
the pseudo-particle $\tilde{ij}$, which then clusters with parton $k$ to form
the final jet.

Within the E-scheme recombination, the four-momentum of the pseudo particle
$\tilde{ij}$ is given by $p^{\mu}_{\tilde{ij}}=p^{\mu}_{i}+p^{\mu}_{j}$ and
thus

$\displaystyle 2p_{\tilde{ij}}\cdot
p_{k}=2p_{k,T}(m_{\tilde{ij},T}\cosh\Delta\eta_{\tilde{ij},k}-p_{\tilde{ij},T}\cos\Delta\phi_{\tilde{ij},k})\,$
(31) $\displaystyle=$ $\displaystyle
2p_{k,T}\left(\sqrt{p^{2}_{\tilde{ij},T}+s_{ij}}\cosh\Delta\eta_{\tilde{ij},k}-p_{\tilde{ij},T}\cos\Delta\phi_{\tilde{ij},k}\right)$
$\displaystyle=$ $\displaystyle
2p_{k,T}p_{\tilde{ij},T}(\cosh\Delta\eta_{\tilde{ij},k}-\cos\Delta\phi_{\tilde{ij},k})+p_{k,T}\frac{s_{ij}}{p_{\tilde{ij},T}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\times\left(1+{\cal
O}\left(s_{ij}/p^{2}_{{\tilde{ij}},T}\right)\right)$ $\displaystyle\approx$
$\displaystyle
2p_{T}^{2}(z_{i}+z_{j})z_{k}(\cosh\Delta\eta_{\tilde{ij},k}-\cos\Delta\phi_{\tilde{ij},k})$
$\displaystyle+\frac{z_{k}z_{i}z_{j}}{(z_{i}+z_{j})}\frac{s_{ij}}{z_{i}z_{j}p_{T}^{2}R^{2}}p_{T}^{2}R^{2}\,.$

From this, it is easy to see that the requirement in Eq. (30) of all three
partons being clustered into one jet can be written as

$\displaystyle c_{ijk}\equiv\left(z_{i}\tilde{s}_{ik}+z_{j}\tilde{s}_{jk}\leq
z_{i}+z_{j}+\frac{z_{i}z_{j}}{z_{i}+z_{j}}\tilde{s}_{ij}\right)\,,$ (32)

valid in the small-$R$ limit. Here again,
$\tilde{s}_{ij}\equiv\frac{s_{ij}}{z_{i}z_{j}}\frac{1}{p_{T}^{2}R^{2}}$.

In the earlier history of the clustering, in order to have $i$ and $j$ to be
clustered first, we need

$\displaystyle\rho_{ij}<\min\left[\rho_{ik}\,,\rho_{jk}\,,\rho_{i}\,,\rho_{j}\,,\rho_{k}\right]\,,$
(33)

as required by the jet algorithm. This maps to the conditions for the
following cases:

  * •

case I:
$\tilde{s}_{ij}<z_{j}^{2\alpha}z_{k}^{-2\alpha}\min[\tilde{s}_{ik},\tilde{s}_{jk},1]$,
for $z_{i}\leq z_{j}\leq z_{k}$;

  * •

case II:
$\tilde{s}_{ij}<\min[z_{k}^{-2\alpha}z_{j}^{2\alpha}\tilde{s}_{ik},\tilde{s}_{jk},1]$,
for $z_{i}\leq z_{k}\leq z_{j}$;

  * •

case III:
$\tilde{s}_{ij}<\min[z_{j}^{2\alpha}z_{i}^{-2\alpha}\tilde{s}_{ik},\tilde{s}_{jk},1]$,
for $z_{k}\leq z_{i}\leq z_{j}$.

Other cases can be obtained by switching the order between $i$ and $j$.

The jet clustering conditions can be expressed as constraints in the phase
space integrals. For case I, we have

$\displaystyle
d\Phi^{\text{I}}_{3}=d\Phi_{3}\,\theta\left(\tilde{s}_{ij}<\left(\frac{z_{j}}{z_{k}}\right)^{2\alpha}\tilde{s}_{jk}\right)\theta\left(\tilde{s}_{ij}<\left(\frac{z_{j}}{z_{k}}\right)^{2\alpha}\tilde{s}_{ik}\right)\,$
$\displaystyle\times{\theta\left(\tilde{s}_{ij}<\left(\frac{z_{j}}{z_{k}}\right)^{2\alpha}\right)}\theta\left(c_{ijk}\right)\theta(z_{i}\leq
z_{j}\leq z_{k})\,,$ (34)

where $d\Phi_{3}$ is the three-body phase space given in the previous section
and $c_{ijk}$ is given in Eq. (32). We note that the $\theta$-functions in the
first line due to the clustering condition are dangerous, in the sense that
they will change the pre-determined fractional power of $\epsilon$ in the
exponent in Eq. (III.2.2), associated with the $z_{j}\to 0$ poles and
therefore invalidate the Laurent expansion in Eq. (9). To understand the
issue, we consider a simple example by comparing the integrals

$\displaystyle
I_{1}=\int_{0}^{1}\mathrm{d}x_{3}\,\mathrm{d}x_{2}\,x_{2}^{-1-a_{2}\epsilon}x_{3}^{-1-a_{3}\epsilon}\,,$
(35)

and

$\displaystyle
I_{2}=\int_{0}^{1}\mathrm{d}x_{3}\mathrm{d}x_{2}\,x_{2}^{-1-a_{2}\epsilon}x_{3}^{-1-a_{3}\epsilon}\theta(x_{3}^{\alpha}-x_{2})\,.$
(36)

Apparently we can safely Laurent-expand the integrand
$x_{2}^{-1-a_{2}\epsilon}x_{3}^{-1-a_{3}\epsilon}$ in $I_{1}$ to perform the
integration. This naive expansion is not allowed in $I_{2}$, since the correct
power of $\epsilon$ for $x_{3}$ is not $-a_{3}\epsilon$ but
$-(a_{3}+a_{2}\alpha)\epsilon$, which can be easily seen by integrating out
$x_{2}$. To avoid this complication, we introduce a subtraction term, to
remove all possible the $z_{j}\to 0$ singularities by taking the $z_{j}\to 0$
limit in the clustering condition to find

$\displaystyle d\Phi^{\text{I}}_{3,sub.}$ $\displaystyle=$ $\displaystyle
d\Phi_{3}\,\theta\left(\tilde{s}_{ij}<z_{j}^{2\alpha}\tilde{s}_{jk}\right)\,\theta\left(\tilde{s}_{jk}\leq
1\right)\theta\left(\tilde{s}_{ik}\leq 1\right)$ (37)
$\displaystyle\times\theta(z_{i}\leq z_{j}\leq z_{k})\,,$

where we have used the fact that as $z_{j}\to 0$, $\tilde{s}_{ij}\to 0$ as
forced by the $\theta$-functions, and therefore
$\tilde{s}_{jk}\to\tilde{s}_{ik}$. With the subtraction, the integral
$(d\Phi_{3}^{I}-d\Phi^{I}_{3,sub.})|{\cal M}|^{2}$ will be free of the
$z_{j}\to 0$ poles. The term $d\Phi^{I}_{3,sub.}\,|{\cal M}|^{2}$ alone still
contains the $z_{j}\to 0$ singularity and therefore the potential $\epsilon$
fractional power modification due to
$\theta\left(\tilde{s}_{ij}<z_{j}^{2\alpha}\tilde{s}_{jk}\right)$. But we will
see later that when combined with other subtraction terms introduced soon, all
such dangerous $\theta$-functions go away.

We have, for case II and case III, the phase space integrals

$\displaystyle d\Phi^{\text{II}}_{3}$ $\displaystyle=$ $\displaystyle
d\Phi_{3}\left[1-\theta\left(\tilde{s}_{ij}<\left(\frac{z_{j}}{z_{k}}\right)^{2\alpha}\tilde{s}_{ik}\right)\right]\theta(\tilde{s}_{ik}<1)$
(38)
$\displaystyle\times\theta(\tilde{s}_{ik}\leq\tilde{s}_{jk})\theta\left(c_{ikj}\right)\theta(z_{i}\leq
z_{j}\leq z_{k})\,,$

and

$\displaystyle d\Phi_{3}^{\text{III}}$ $\displaystyle=$ $\displaystyle
d\Phi_{3}\times\left[1-\theta\left(\tilde{s}_{ij}<\left(\frac{z_{j}}{z_{k}}\right)^{2\alpha}\tilde{s}_{jk}\right)\right]\theta(\tilde{s}_{jk}<1)$
(39)
$\displaystyle\times\theta(\tilde{s}_{jk}<\tilde{s}_{ik})\theta\left(c_{jki}\right)\theta(z_{i}\leq
z_{j}\leq z_{k})\,,$

respectively. Here for later use, we have relabeled subscripts $i$, $j$ and
$k$ to stick to the case that $z_{i}\leq z_{j}\leq z_{k}$, but with $i$, $k$
clustered first in case II and $j$ and $k$ first in case III. This is always
allowed since eventually we sum over all possible orders of $z_{i}$’s and all
possible clustering histories. Following the same logic as in case I, we
construct the subtraction terms for case II

$\displaystyle
d\Phi^{\text{II}}_{3,sub.}\,=\,d\Phi_{3}\left[1-\theta\left(\tilde{s}_{ij}<z_{j}^{2\alpha}\tilde{s}_{jk}\right)\right]\theta(\tilde{s}_{ik}\leq\tilde{s}_{jk})$
$\displaystyle\times\,\theta(\tilde{s}_{ik}<1)\theta\left(\tilde{s}_{jk}<1\right)\theta(z_{i}<z_{j}<z_{k})\,,\qquad$
(40)

and case III

$\displaystyle
d\Phi_{3,sub.}^{\text{III}}\,=\,d\Phi_{3}\left[1-\theta\left(\tilde{s}_{ij}<z_{j}^{2\alpha}\tilde{s}_{jk}\right)\right]\theta(\tilde{s}_{jk}<\tilde{s}_{ik})$
$\displaystyle\times\,\theta(\tilde{s}_{jk}<1)\theta\left(\tilde{s}_{ik}<1\right)\times\theta(z_{i}<z_{j}<z_{k})\,,\qquad$
(41)

respectively, to remove the $z_{j}\to 0$ singularities.

The phase space integral with the jet algorithm implemented reads

$\displaystyle d\Phi_{\rm alg.}|{\cal
M}|^{2}=\sum_{i=\text{I}}^{\text{III}}\left(d\Phi_{3}^{i}-d\Phi_{3,sub.}^{i}\right)|{\cal
M}|^{2}+\,d\Phi_{3,sub.}^{i}|{\cal M}|^{2}$ $\displaystyle\hskip
21.52771pt+\left(\text{all $z_{i}$, $z_{j}$, $z_{k}$ orders}\right)\,,$ (42)

and we note that by construction, it is easy to check that

$\displaystyle\sum_{i={\text{I}}}^{\text{III}}d\Phi_{3,sub.}^{i}=d\Phi_{3}\theta(\tilde{s}_{jk}<1)\theta\left(\tilde{s}_{ik}<1\right)\,\theta(z_{i}<z_{j}<z_{k})\,,$

where all $\theta$-function constraints involving powers of $z_{j}^{\alpha}$
cancel in the sum. Therefore, the phase space integration does not induce any
modifications of the fractional powers in $\epsilon$ as $z_{j}\to 0$.

However at this stage, the phase space parameterization in Eq. (III.2.2) is
still not sufficient to isolate all singularities within $|{\cal M}|^{2}$.
Particularly in $\sum_{i={\text{I}}}^{\text{III}}d\Phi_{3,sub.}^{i}$, a
further ordering between the variables $x_{2}$ and $x_{4}$ is required to
fully disentangle the poles generated by
$s_{ijk}=z_{i}z_{k}\tilde{s}_{ik}+z_{j}z_{k}\tilde{s}_{jk}+z_{i}z_{j}\tilde{s}_{ij}\to
0$. To avoid such an ordering, we introduce instead a subtraction term in the
matrix element to regulate the single soft limit, as $x_{4}\to 0$, of $|{\cal
M}|^{2}$. Thus, a natural choice for the subtraction term $|{\cal
M}|_{sub.}^{2}$ is

$\displaystyle|{\cal M}|_{sub.}^{2}=\lim_{x_{4}\to 0}|{\cal M}|^{2}=T_{j}\cdot
T_{k}\,\frac{s_{jk}}{s_{ij}s_{ik}}P^{(0)}_{a\to jk}\,,$ (44)

which is nothing but the product of the LO eikonal factor and the LO splitting
kernel.

The final phase space integration suitable for a numerical evaluation is thus
given by

$\displaystyle d\Phi_{\rm alg.}|{\cal M}|^{2}$ $\displaystyle=$
$\displaystyle\sum_{i=\text{I}}^{\text{III}}d\Phi_{3,sub.}^{i}|{\cal
M}|_{sub.}^{2}$ (45)
$\displaystyle+\left(d\Phi_{3}^{i}-d\Phi_{3,sub.}^{i}\right)|{\cal M}|^{2}$
$\displaystyle+\,\sum_{i=\text{I}}^{\text{III}}d\Phi_{3,sub.}^{i}(|{\cal
M}|^{2}-|{\cal M}|_{sub.}^{2})$ $\displaystyle+\left(\text{all $z_{i}$,
$z_{j}$, $z_{k}$ orders}\right)\,.$

Similar subtraction counter terms can be constructed for the WTA scheme
recombination as well as other jet algorithms such as the Cambridge-Aachen
algorithm and thus used to compute the soft-drop groomed jet function. Further
applications to obtain the semi-inclusive jet function Kang:2016mcy ;
Dai:2016hzf  at ${\cal O}(\alpha_{s}^{2})$ are also straightforward.

####  III.2.4 Results and validations

With the setups developed in the previous sections, we find the double-real
correction to the anti-$k_{T}$ jet function to be

$\displaystyle J_{rr}^{(2)}=\,\frac{\alpha_{s}^{2}e^{4\epsilon
L}}{(2\pi)^{2}}C_{F}\Big{(}C_{F}{\cal K}^{rr}_{C_{F}}+C_{A}{\cal
K}^{rr}_{C_{A}}+N_{F}T_{F}{\cal K}^{rr}_{N_{F}T_{F}}\Big{)}\,,$

where we find

$\displaystyle{\cal K}^{rr}_{C_{F}}$ $\displaystyle=$
$\displaystyle\frac{1}{2\epsilon^{4}}+\frac{3}{2\epsilon^{2}}-\frac{1.8171(3)}{\epsilon^{2}}-\frac{{21.272(3)}}{\epsilon}-{76.42(2)}\,,$
$\displaystyle{\cal K}^{rr}_{C_{A}}$ $\displaystyle=$
$\displaystyle\frac{1}{4\epsilon^{4}}+\frac{1.20833}{\epsilon^{3}}+\frac{1.5484(2)}{\epsilon^{2}}$
(47) $\displaystyle-\frac{{7.941(4)}}{\epsilon}-{75.425(2)}\,,$

and

$\displaystyle{\cal K}^{rr}_{N_{F}T_{F}}$ $\displaystyle=$
$\displaystyle-\frac{1}{6\epsilon^{3}}-\frac{7}{9\epsilon^{2}}+\frac{0.1067(3)}{\epsilon}+{17.230(2)}\,.\qquad$
(48)

Here and below we always include in brackets the numerical error on the last
digit.

To validate our results, we have performed independent calculations and found
full agreement. Also we note that parts of our calculations associated with
$d\Phi_{3,sub.}$ can be done (semi-)analytically and for these parts, we have
checked our numerical computation against the (semi-)analytic results to find
full agreement. Some of these checks are presented in the appendix C.

####  III.2.5 Renormalization group equation

The renormalization group equation (RGE) for the quark-jet function serves as
strong independent check on the calculation. The leading poles up to
$\epsilon^{-2}$ of the double-real correction in Eq. (III.2.4) combined with
the real-virtual one in Eq. (11) can be predicted by solving the RGE up to the
logarithmic term $\alpha_{s}^{2}L^{2}$ and taking the NLO quark-jet function
from Eq. (10) as input. Details are given in the appendix B.

As a result, the RGE predicts the leading three poles normalized to
$\alpha_{s}^{2}/(2\pi)^{2}$ for the $C_{F}^{2}$ term as

$\displaystyle{\cal K}^{rr}_{C_{F}}\biggr{|}_{\epsilon^{-4}}$ $\displaystyle=$
$\displaystyle\frac{1}{2\epsilon^{4}}\,,$ (49) $\displaystyle{\cal
K}^{rr}_{C_{F}}\biggr{|}_{\epsilon^{-3}}$ $\displaystyle=$
$\displaystyle\frac{3}{2\epsilon^{3}}\,,$ (50) $\displaystyle\big{(}{\cal
K}^{rv}_{C_{F}}+{\cal K}^{rr}_{C_{F}}\big{)}\biggr{|}_{\epsilon^{-2}}$
$\displaystyle=$
$\displaystyle\left(\frac{61}{8}-\frac{3\pi^{2}}{4}\right)\frac{1}{\epsilon^{2}}\approx\frac{0.222797}{\epsilon^{2}}$
(51) $\displaystyle=$ $\displaystyle\frac{0.2228(3)}{\epsilon^{2}}\,,\qquad$

which all agree exactly with Eqs. (12) and (III.2.4), in the case of the
$\epsilon^{-2}$ pole within numerical errors as shown in the second line in
Eq. (51).

For the $C_{F}N_{F}T_{F}$ contribution, solving the RGE we find the leading
poles

$\displaystyle{\cal K}^{rr}_{N_{F}T_{F}}\biggr{|}_{\epsilon^{-3}}$
$\displaystyle=$ $\displaystyle-\frac{1}{6\epsilon^{3}}\,,$ (52)
$\displaystyle{\cal K}^{rr}_{N_{F}T_{F}}\biggr{|}_{\epsilon^{-2}}$
$\displaystyle=$ $\displaystyle-\frac{7}{9\epsilon^{2}}\,,$ (53)

which fully agrees with what we obtained in Eq. (48).

As for the $C_{F}C_{A}$ term, the naive RGE alone can predict the
$\epsilon^{-4}$ and $\epsilon^{-3}$ poles correctly, however the
$\epsilon^{-2}$ term is affected by the non-global contribution, for which we
find that an additional $-\frac{\pi^{2}}{12\epsilon^{2}}$ correction is needed
to produce the correct pole and the logarithm at $\alpha_{s}^{2}L^{2}$. The
same situation has been encountered in the Sterman-Weinberg cone jet in
Becher:2015hka  and for the hemisphere mass distribution in Schwartz:2014wha .

Once every piece is taken into account, the poles predicted for the
$C_{F}C_{A}$ contribution are

$\displaystyle\big{(}{\cal K}^{rv}_{C_{A}}+{\cal
K}^{rr}_{C_{A}}\big{)}\biggr{|}_{\epsilon^{-4}}$ $\displaystyle=$
$\displaystyle 0\,,$ (54) $\displaystyle\big{(}{\cal K}^{rv}_{C_{A}}+{\cal
K}^{rr}_{C_{A}}\big{)}\biggr{|}_{\epsilon^{-3}}$ $\displaystyle=$
$\displaystyle\frac{11}{24\epsilon^{3}}\approx\frac{0.458333}{\epsilon^{3}}$
(55) $\displaystyle=$ $\displaystyle\frac{0.45833(5)}{\epsilon^{3}}\,,\quad$
$\displaystyle\big{(}{\cal K}^{rv}_{C_{A}}+{\cal
K}^{rr}_{C_{A}}\big{)}\biggr{|}_{\epsilon^{-2}}$ $\displaystyle=$
$\displaystyle\left(\frac{83}{36}-\frac{\pi^{2}}{8}\right)\frac{1}{\epsilon^{2}}\approx\frac{1.07186}{\epsilon^{2}}$
(56) $\displaystyle=$ $\displaystyle\frac{1.0720(2)}{\epsilon^{2}}\,,$

where the agreement is within numerical errors. For the $\epsilon^{-2}$ term
it is at ${\cal O}(1\mbox{\textperthousand})$ as shown in the second line in
Eq. (56).

##  IV Two-loop jet function and discussions

Now we present the anti-$k_{T}$ bare quark-jet function at two loops

$\displaystyle J_{bare}$ $\displaystyle=$ $\displaystyle
1+J_{bare}^{(1)}+\frac{\alpha_{s}^{2}e^{4\epsilon
L}}{4\pi^{2}}C_{F}\Big{(}C_{F}{\cal J}_{C_{F}}$ (57) $\displaystyle\hskip
38.74988pt+C_{A}{\cal J}_{C_{A}}+N_{F}T_{F}{\cal J}_{N_{F}T_{F}}\Big{)}\,,$

where $J^{(1)}_{bare}$ has been given in Eq. (10). The new two-loop result
reads

$\displaystyle{\cal J}_{C_{F}}$ $\displaystyle=$
$\displaystyle\frac{1}{2\epsilon^{4}}+\frac{3}{2\epsilon^{3}}+\frac{1}{\epsilon^{2}}\left(\frac{61}{8}-\frac{3\pi^{2}}{4}\right)$
(58) $\displaystyle-\frac{{5.392(3)}}{\epsilon}-{15.93(2)}\,,$
$\displaystyle{\cal J}_{C_{A}}$ $\displaystyle=$
$\displaystyle\frac{11}{24\epsilon^{3}}+\frac{1}{\epsilon^{2}}\left(\frac{83}{36}-\frac{\pi^{2}}{8}\right)-\frac{{12.985(4)}}{\epsilon}-{115.55(2)}\,,$

and

$\displaystyle{\cal
J}_{N_{F}T_{F}}=-\frac{1}{6\epsilon^{3}}-\frac{7}{9\epsilon^{2}}+\frac{0.1067(3)}{\epsilon}+{17.230(2)}\,,$
(60)

with additional ${\cal O}(\alpha_{s}^{2})$ contributions coming from the
$\alpha_{s}$ renormalization of $J_{bare}^{(1)}$ in Eq. (10). The renormalized
jet function $J$ is then found to be

$\displaystyle J$ $\displaystyle=$ $\displaystyle
1+\frac{\alpha_{s}}{2\pi}C_{F}\left(\frac{13}{2}-\frac{3\pi^{2}}{4}\right)+\frac{\alpha_{s}^{2}}{4\pi^{2}}\left({-1.78(2)}C_{F}^{2}\right.$
(61)
$\displaystyle\left.{-106.87(2)}C_{A}C_{F}+{14.072(2)}C_{F}N_{F}T_{F}\right)\,,\qquad$

where we have set the scale $\mu=p_{T}R$. The associated numerical errors
originate entirely from the NNLO real-real corrections.

Before concluding we make two comments on our calculation and the result.
First, although we present here the jet function $J$ in the simplified
factorization in Eq. (3), there is no difficulty to provide also the angular
dependent jet function $J_{m}$ in Eq. (1) with explicit dependence on the
collinear parton multiplicity $m$ up to $m=3$. In fact, the results of
$J_{m=2}^{(0)}$, $J_{m=2}^{(1)}$ and $J_{m=3}^{(0)}$ are already encoded in
our calculation, which can be easily seen from the fully differential nature
of the phase space sector decomposition and the relation that
$J=\sum_{m}\langle J_{m}\rangle_{\Omega}$. To demonstrate such feasibility, we
plot in fig. 1, the distribution in $\tilde{s}_{ij}$ for the $\epsilon^{0}$
part of the jet function $J_{m=3}^{(0)}$ with the $P^{({\rm
id})}_{{\bar{q}}_{1}q_{2}q_{3}}$ contribution only. A detailed study of the
angular-dependent jet function will be left for future work.

Figure 1: The distribution $dJ^{({\rm id})}_{{\bar{q}}qq}/d\tilde{s}_{ij}$ for
the $\epsilon^{0}$ part of the jet function $J_{m=3}^{(0)}$ with $P^{({\rm
id})}_{{\bar{q}}_{1}q_{2}q_{3}}$ term only.

The results for the individual $J_{m}$ will be more suitable for resumming
NGLs and for improving the logarithmic accuracy of the predictions.

Second, our two-loop result confirms that the factorization theorem in Eq.
(3), as used in many phenomenological studies of jets, TMDs and in spin
physics, is only valid up to the single logarithmic level $({\cal
O}\left(e^{\alpha_{s}^{n}L^{n}}\right))$ due to the existence of the NGLs,
which is generic for many exclusive jet processes. Beyond this order, the
factorization in Eq. (1) or alternative approaches Nagy:2020rmk  should be
used for predictions. Our calculational framework can be useful for realizing
the resummation in Eq. (1). We leave this for future studies.

##  V Conclusions

In this work, we have developed a method to calculate the jet functions, which
arise in the factorization of cross sections with exclusive jets of small
radius $R$, with a jet clustering procedure up to NNLO. Explicit results at
NNLO have been presented, for the first time, for the quark-jet function with
the anti-$k_{T}$ jet algorithm. The new two-loop contribution to the
anti-$k_{T}$ quark-jet function provides the missing input for the evaluation
of the factorization theorem in Eq. (1) to NNLO and for resummation beyond NLL
accuracy. This will substantially improve the precision of predictions for
observables with jets as hard probes. The computational framework is not
limited to anti-$k_{T}$ E-scheme jets but readily applicable to many other jet
observables up to NNLO, such as the semi-inclusive jet function, the WTA jet
and the soft-drop groomed jet radius $R_{g}$ and the jet momentum fraction
$z_{g}$. In addition, the obtained two-loop quark-jet function will serve as
an important building block when using the factorization in Eq. (1) for the
development of a colorful $q_{T}$-subtraction scheme, which is applicable to
both, calculations at fixed order in perturbation theory up to NNLO and to
parton shower simulations.

In summary, given the broad interest in jets from both the high-energy
particle and the nuclear communities, we expect the results to have numerous
phenomenological applications at colliders, such as the LHC and a future EIC
for the study of jet cross sections, TMDs and in addressing the spin structure
of the nucleon.

## Acknowledgments

We thank Ding-Yu Shao for clarifications on the winner-take-all jet. X.L. is
supported by National Natural Science Foundation of China under Grant
No.11775023. S.M. is supported by Deutsche Forschungsgemeinschaft (DFG)
through the Research Unit FOR 2926, “Next Generation pQCD for Hadron
Structure: Preparing for the EIC”, project MO 1801/5-1. We thank Kyle Lee,
Jürg Haag, Luca Buonocore, Luca Rottoli and Massimiliano Grazzini for pointing
out the inappropriate approximation in the early version of the draft. The
method proposed in this work are not affected, while the numerics in the
concrete example are corrected.

##  Appendix A The $a\to ijk$ splitting functions

The splitting functions used in our calculation are taken from Catani:1999ss .
For the case in which the quark-anti-quark pair
${\bar{q}}_{1}^{\prime}q_{2}^{\prime}$ has a different flavor from the quark
$q_{3}$, we use

$\displaystyle
P_{{\bar{q}}_{1}^{\prime}q_{2}^{\prime}q_{3}}=\frac{1}{2}C_{F}T_{F}\frac{s_{123}}{s_{12}}\left[-\frac{t_{12,3}^{2}}{s_{12}s_{123}}+\frac{4z_{3}+(z_{1}-z_{2})^{2}}{z_{1}+z_{2}}+(1-2\epsilon)\left(z_{1}+z_{2}-\frac{s_{12}}{s_{123}}\right)\right]\,,$
(62)

with the abbreviation

$\displaystyle t_{ij,k}\equiv
2\frac{z_{i}s_{jk}-z_{j}s_{ik}}{z_{i}+z_{j}}+\frac{z_{i}-z_{j}}{z_{i}+z_{j}}s_{ij}\,.$
(63)

For the final state with identical quark flavors ${\bar{q}}_{1}q_{2}q_{3}$, we
need

$\displaystyle
P_{{\bar{q}}_{1}q_{2}q_{3}}=\Big{[}P_{{\bar{q}}_{1}^{\prime}q_{2}^{\prime}q_{3}}+(2\leftrightarrow
3)\Big{]}+P_{{\bar{q}}_{1}q_{2}q_{3}}^{(\rm id)}\,,$ (64)

where $P_{{\bar{q}}_{1}q_{2}q_{3}}^{(\rm id)}$ accounts for the interference
term

$\displaystyle P_{{\bar{q}}_{1}q_{2}q_{3}}^{(\rm id)}$ $\displaystyle=$
$\displaystyle
C_{F}\left(C_{F}-\frac{C_{A}}{2}\right)\,\Bigg{\\{}(1-\epsilon)\left(\frac{2s_{23}}{s_{12}}-\epsilon\right)+\frac{s_{123}}{s_{12}}\left[\frac{1+z_{1}^{2}}{1-z_{2}}-\frac{2z_{2}}{1-z_{3}}-\epsilon\left(\frac{(1-z_{3})^{2}}{1-z_{2}}+1+z_{1}-\frac{2z_{2}}{1-z_{3}}\right)\right.$
(65)
$\displaystyle\left.-\epsilon^{2}(1-z_{3})\right]-\frac{s_{123}^{2}}{s_{12}s_{13}}\frac{z_{1}}{2}\left[\frac{1+z_{1}^{2}}{(1-z_{2})(1-z_{3})}-\epsilon\left(1+2\frac{1-z_{2}}{1-z_{3}}\right)-\epsilon^{2}\right]\Bigg{\\}}+(2\leftrightarrow
3)\,.$

The splitting function for the quark decay into $g_{1}g_{2}q_{3}$ is given by

$\displaystyle P_{g_{1}g_{2}q_{3}}=C_{F}^{2}P_{g_{1}g_{2}q_{3}}^{({\rm
ab})}+C_{F}C_{A}P_{g_{1}g_{2}q_{3}}^{({\rm nab})}\,,$ (66)

where the Abelian contribution is

$\displaystyle P_{g_{1}g_{2}q_{3}}^{({\rm ab})}$ $\displaystyle=$
$\displaystyle\Bigg{\\{}\frac{s_{123}^{2}}{2s_{13}s_{23}}z_{3}\left[\frac{1+z_{3}^{2}}{z_{1}z_{2}}-\epsilon\frac{z_{1}^{2}+z_{2}^{2}}{z_{1}z_{2}}-\epsilon(1+\epsilon)\right]+(1-\epsilon)\left[\epsilon-(1-\epsilon)\frac{s_{23}}{s_{13}}\right]$
(67)
$\displaystyle+\frac{s_{123}}{s_{13}}\left[\frac{z_{3}(1-z_{1})+(1-z_{2})^{3}}{z_{1}z_{2}}+\epsilon^{2}(1+z_{3})-\epsilon(z_{1}^{2}+z_{1}z_{2}+z_{2}^{2})\frac{1-z_{2}}{z_{1}z_{2}}\right]\Bigg{\\}}+(1\leftrightarrow
2)\,,$

and the non-Abelian one reads

$\displaystyle P_{g_{1}g_{2}q_{3}}^{({\rm nab})}$ $\displaystyle=$
$\displaystyle\Bigg{\\{}(1-\epsilon)\left(\frac{t_{12,3}^{2}}{4s_{12}^{2}}+\frac{1}{4}-\frac{\epsilon}{2}\right)+\frac{s_{123}^{2}}{2s_{12}s_{13}}\left[\frac{(1-z_{3})^{2}(1-\epsilon)+2z_{3}}{z_{2}}+\frac{z_{2}^{2}(1-\epsilon)+2(1-z_{2})}{1-z_{3}}\right]$
(68)
$\displaystyle-\frac{s_{123}^{2}}{4s_{13}s_{23}}\,z_{3}\left[\frac{(1-z_{3})^{2}(1-\epsilon)+2z_{3}}{z_{1}z_{2}}+\epsilon(1-\epsilon)\right]$
$\displaystyle+\frac{s_{123}}{2s_{12}}\left[(1-\epsilon)\frac{z_{1}(2-2z_{1}+z_{1}^{2})-z_{2}(6-6z_{2}+z_{2}^{2})}{z_{2}(1-z_{3})}+2\epsilon\frac{z_{3}(z_{1}-2z_{2})-z_{2}}{z_{2}(1-z_{3})}\right]$
$\displaystyle+\frac{s_{123}}{2s_{13}}\left[(1-\epsilon)\frac{(1-z_{2})^{3}+z_{3}^{2}-z_{2}}{z_{2}(1-z_{3})}-\epsilon\left(\frac{2(1-z_{2})(z_{2}-z_{3})}{z_{2}(1-z_{3})}-z_{1}+z_{2}\right)\right.$
$\displaystyle\hskip
34.44434pt\left.-\frac{z_{3}(1-z_{1})+(1-z_{2})^{3}}{z_{1}z_{2}}+\epsilon(1-z_{2})\left(\frac{z_{1}^{2}+z_{2}^{2}}{z_{1}z_{2}}-\epsilon\right)\right]\Bigg{\\}}+(1\leftrightarrow
2)\,.$

##  Appendix B The leading $\epsilon$-poles and large logarithms from the RGE

Up to non-global contributions, the jet function $J$ satisfies the RGE,

$\displaystyle\mu\frac{\mathrm{d}J}{\mathrm{d}\mu}=\left(2\Gamma[\alpha_{s}]L+\gamma_{J_{q}}\right)J\,,$
(69)

where $\Gamma$ is the cusp anomalous dimension, cf. e.g. Moch:2004pa ,

$\displaystyle\Gamma=\sum_{n=0}\left(\frac{\alpha_{s}}{4\pi}\right)^{n+1}\Gamma_{n}\,,\quad\quad\mbox{with}\quad\quad\Gamma_{0}=4C_{F}\,,\quad\quad\Gamma_{1}=4C_{F}\left[\left(\frac{67}{9}-\frac{\pi^{2}}{3}\right)C_{A}-\frac{20}{9}T_{F}N_{F}\right]\,,$
(70)

and $\gamma_{J_{q}}$ is the anomalous dimension of the quark-jet,

$\displaystyle\gamma_{J_{q}}=\sum_{n=0}\left(\frac{\alpha_{s}}{4\pi}\right)^{n+1}\gamma_{n}\,,\quad\quad\mbox{with}\quad\quad\gamma_{0}=6C_{F}\,.$
(71)

Using the fact that
$\mathrm{d}\log\mu=\beta^{-1}[\alpha_{s}]\,\mathrm{d}\alpha_{s}$, with
$\beta[\alpha_{s}]=-2\alpha_{s}\sum_{n=0}\beta_{n}\left(\frac{\alpha_{s}}{4\pi}\right)^{n+1}$,
the RGE can be solved to find

$\displaystyle
J(\mu)=J(\mu=p_{T}R)\exp\left[2\int_{\alpha_{s}(p_{T}R)}^{\alpha_{s}(\mu)}\mathrm{d}\alpha_{s}\frac{\Gamma[\alpha_{s}]}{\beta[\alpha_{s}]}\int_{\alpha_{s}(p_{T}R)}^{\alpha_{s}}\frac{\mathrm{d}\alpha_{s}^{\prime}}{\beta[\alpha_{s}^{\prime}]}\,+\int_{\alpha_{s}(p_{T}R)}^{\alpha_{s}(\mu)}\mathrm{d}\alpha_{s}\frac{\gamma_{J_{q}}[\alpha_{s}]}{\beta[\alpha_{s}]}\right]\,.$
(72)

Expanding the solution to ${\cal O}(\alpha_{s}^{2})$, we find the logarithmic
terms in the jet function at $\alpha_{s}^{2}$ contain

$\displaystyle J^{(2)}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{s}^{2}}{(2\pi)^{2}}\left[\frac{\Gamma_{0}^{2}}{8}L^{4}+\Gamma_{0}\left(\frac{\beta_{0}}{6}+\frac{\gamma_{0}}{4}\right)L^{3}+\frac{1}{8}\left(2\beta_{0}\gamma_{0}+\gamma_{0}^{2}+4\Gamma_{0}J^{(1)}_{0,ren}+2\Gamma_{1}\right)L^{2}\right.$
(73)
$\displaystyle\left.+\frac{1}{4}\left(2\gamma_{0}J^{(1)}_{0,ren}+4\beta_{0}J^{(1)}_{0,ren}+\gamma_{1}\right)L+J^{(2)}_{0,ren}\right]\,,$

where $J^{(l)}_{0,ren}$ are the renormalized jet functions at NLO ($l=1$) and
at NNLO ($l=2$) with $\mu=p_{T}R$ and $\epsilon\to 0$.

On the other hand, if we assume the bare jet function have the form

$\displaystyle J_{bare}=1+Z_{\alpha}\frac{\alpha_{s}}{2\pi}e^{2\epsilon
L}\left(\frac{\Gamma_{0}}{4\epsilon^{2}}+\frac{\gamma_{0}}{4\epsilon}+J_{0}^{(1)}+J_{1}^{(1)}\epsilon\right)+\frac{\alpha_{s}^{2}}{(2\pi)^{2}}e^{4\epsilon
L}\left(\frac{J^{(2)}_{-4}}{\epsilon^{4}}+\frac{J^{(2)}_{-3}}{\epsilon^{3}}+\frac{J^{(2)}_{-2}}{\epsilon^{2}}+\frac{J^{(2)}_{-1}}{\epsilon^{1}}+J^{(2)}_{0}\right)\,,$
(74)

and hence the NLO renormalization factor of the jet function is

$\displaystyle
Z_{1}=-\frac{\alpha_{s}}{2\pi}\left(\frac{\Gamma_{0}}{4\epsilon^{2}}+\frac{\gamma_{0}}{4\epsilon}+\frac{\Gamma_{0}\,L}{2\epsilon}\right)\,.$
(75)

It is straightforward to check that $J_{0}^{(1)}=J_{0,ren}^{(1)}$. From these,
the logarithmic contributions can be obtained via

$\displaystyle J_{bare}+Z_{1}+Z_{1}J_{bare}^{(1)}$ $\displaystyle\supset$
$\displaystyle\frac{\alpha_{s}^{2}}{(2\pi)^{2}}\left[\left(\frac{32J^{(2)}_{-4}}{3}-\frac{5\Gamma_{0}^{2}}{24}\right)L^{4}+\left(\frac{32J^{(2)}_{-3}}{3}-\frac{\beta_{0}\Gamma_{0}}{6}-\frac{5\gamma_{0}\Gamma_{0}}{12}\right)L^{3}\right.$
(76)
$\displaystyle\left.+\left(8J^{(2)}_{-2}-\frac{1}{8}(\gamma_{0}+2\beta_{0})\gamma_{0}-\frac{3}{2}\Gamma_{0}J^{(1)}_{0}\right)L^{2}+\left(4J^{(2)}_{-1}-\frac{1}{2}(\gamma_{0}+2\beta_{0})J^{(1)}_{0}-\Gamma_{0}J^{(1)}_{1}\right)L\right]\,.\qquad$

The coefficients of the $\epsilon$-poles $J_{-n}^{(2)}$ can be obtained by
relating Eqs. (73) and (76). The terms $J_{-n}^{(2)}$ with $n=4,3,2$ serve as
a check, although due to the existence of the non-global logarithms, the
$\epsilon^{-2}$-pole of $C_{A}C_{F}$ contribution will be affected and will
deviate from the RGE prediction. The term $J_{-1}^{(2)}$ allows for the
extraction of the two-loop expression for the anomalous dimension of the
quark-jet $\gamma_{1}$, which can be determined as

$\displaystyle\gamma_{1}={11.17(5)}\,C_{F}^{2}{-181.30(6)}\,C_{F}C_{A}{-7.916(5)}\,C_{F}N_{F}T_{F}\,.$
(77)

With quark-jet function computed to NNLO in Eq. (61) and all anomalous
dimensions known to the relevant order we can also predict the logarithmic
structure up to the next-to-next-to-leading logarithm at ${\cal
O}(\alpha_{s}^{3})$ based on the naive RGE, which reads

$\displaystyle J^{(3)}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{s}^{3}}{(2\pi)^{3}}\Bigg{\\{}\frac{\Gamma_{0}^{3}}{48}L^{6}+\frac{1}{48}\left(4\beta_{0}+3\gamma_{0}\right)\Gamma_{0}^{2}\,L^{5}+\frac{1}{48}\Gamma_{0}\left(4\beta_{0}^{2}+10\beta_{0}\gamma_{0}+3\gamma_{0}^{2}+6J^{(1)}_{0,ren}\Gamma_{0}+6\Gamma_{1}\right)L^{4}$
(78)
$\displaystyle+\frac{1}{48}\left(8\beta_{0}^{2}\gamma_{0}+6\beta_{0}\gamma_{0}^{2}+\gamma_{0}^{3}+4\beta_{1}\Gamma_{0}+6\Gamma_{0}\gamma_{1}+6\gamma_{0}(2J^{(1)}_{0,ren}\Gamma_{0}+\Gamma_{1})+8\beta_{0}(4J^{(1)}_{0,ren}\Gamma_{0}+\Gamma_{1})\right)L^{3}$
$\displaystyle+\left(J^{(1)}_{0,ren}\beta_{0}^{2}+\frac{3}{4}J^{(1)}_{0,ren}\beta_{0}\gamma_{0}+\frac{\gamma_{0}\beta_{1}}{8}+\frac{\gamma_{0}^{2}}{8}J^{(1)}_{0,ren}\,+\frac{\Gamma_{0}}{2}J^{(2)}_{0,ren}\,+\frac{\gamma_{1}\beta_{0}}{4}+\frac{\gamma_{0}\gamma_{1}}{8}+\frac{\Gamma_{1}}{4}\,J^{(1)}_{0,ren}+\frac{\Gamma_{2}}{8}\right)L^{2}$
$\displaystyle+{\cal O}(L)\Bigg{\\}}\,,$

where sub-leading logarithmic terms require additional inputs from higher
order calculations.

Finally, we note that in dimensional regularization the $\epsilon$-poles of
the quark-jet function display the a pattern of exponentiation analogous to
the quark form-factor, cf. e.g. Moch:2005id , with the cusp anomalous
dimensions determining the leading poles and the anomalous dimension of the
quark-jet $\gamma_{J_{q}}$ governing the sub-leading ones.

##  Appendix C Analytic results

The analytic computation of some of the phase space integrals associated with
Eq. (III.2.3)

$\displaystyle
d\Phi_{incl.}\equiv\sum_{i={\text{I}}}^{\text{III}}d\Phi_{3,sub.}^{i}\times(p_{T}R)^{-4}=d\Phi_{3}\times(p_{T}R)^{-4}\times\theta(\tilde{s}_{jk}<1)\theta\left(\tilde{s}_{ik}<1\right)\,\theta(z_{i}<z_{j}<z_{k})\,,$
(79)

is straightforward, specifically for integrals like

$\displaystyle(a)=d\Phi_{incl.}\,\frac{1}{\tilde{s}_{ik}\tilde{s}_{jk}}\,,\quad\quad\text{and}\quad\quad(b)=d\Phi_{incl.}\,\frac{1}{\tilde{s}_{ik}\tilde{s}_{ij}}\,.$
(80)

For $(a)$ we directly obtain

$\displaystyle(a)$ $\displaystyle=$ $\displaystyle
2^{2-4\epsilon}(p_{T}R)^{-4\epsilon}\frac{\int_{0}^{1}dz_{i}dz_{j}}{(4\pi)^{5-2\epsilon}\Gamma(1-2\epsilon)}(z_{i}z_{j}z_{k})^{1-2\epsilon}\,\int
d\tilde{s}_{ik}d\tilde{s}_{jk}\left(\tilde{s}_{ik}\tilde{s}_{jk}\right)^{-1-\epsilon}\int_{0}^{1}dt\,t^{-\frac{1}{2}-\epsilon}\left(1-t\right)^{-\frac{1}{2}-\epsilon}$
(81) $\displaystyle=$
$\displaystyle\frac{(p_{T}R)^{-4\epsilon}}{\epsilon^{2}\Gamma(1-\epsilon)^{2}}\frac{\int_{0}^{1}dz_{i}dz_{j}}{(4\pi)^{4-2\epsilon}}(z_{i}z_{j}z_{k})^{1-2\epsilon}\theta(z_{i}<z_{j}<z_{k})\,,$

while for $(b)$, we find

$\displaystyle(b)$ $\displaystyle=$ $\displaystyle
2^{2-4\epsilon}(p_{T}R)^{-4\epsilon}\frac{\int_{0}^{1}dz_{i}dz_{j}}{(4\pi)^{5-2\epsilon}\Gamma(1-2\epsilon)}(z_{i}z_{j}z_{k})^{1-2\epsilon}$
(82) $\displaystyle\times\int
d\tilde{s}_{ik}d\tilde{s}_{jk}\tilde{s}_{ik}^{-1-\epsilon}\tilde{s}_{jk}^{-\epsilon}\int_{0}^{1}dt\left((\sqrt{\tilde{s}_{jk}}-\sqrt{\tilde{s}_{ik}})^{2}+4\sqrt{\tilde{s}_{jk}\tilde{s}_{ik}}t\right)^{-1}\,t^{-\frac{1}{2}-\epsilon}\left(1-t\right)^{-\frac{1}{2}-\epsilon}$
$\displaystyle=$ $\displaystyle
2^{2-4\epsilon}(p_{T}R)^{-4\epsilon}\,\frac{\Gamma\left(\frac{1}{2}-\epsilon\right)^{2}}{\Gamma(1-2\epsilon)}\,\frac{\int_{0}^{1}dz_{i}dz_{j}}{(4\pi)^{5-2\epsilon}\Gamma(1-2\epsilon)}(z_{i}z_{j}z_{k})^{1-2\epsilon}$
$\displaystyle\times\int
d\tilde{s}_{ik}d\tilde{s}_{jk}\tilde{s}_{ik}^{-1-\epsilon}\tilde{s}_{jk}^{-\epsilon}\,\left(\tilde{s}_{ik}+\tilde{s}_{jk}-2\sqrt{\tilde{s}_{ik}\tilde{s}_{jk}}\right)^{-1}{}_{2}F_{1}\left[1,\frac{1}{2}-\epsilon,1-2\epsilon,-\frac{4\sqrt{\tilde{s}_{ik}\tilde{s}_{jk}}}{\tilde{s}_{ik}+\tilde{s}_{jk}-2\sqrt{\tilde{s}_{ik}\tilde{s}_{jk}}}\right]\,.$

By splitting up the integral into two regions,
$\tilde{s}_{ik}\leq\tilde{s}_{jk}$ and $\tilde{s}_{ik}>\tilde{s}_{jk}$, we
have for $\tilde{s}_{ik}\leq\tilde{s}_{jk}$, we let
$\tilde{s}_{ik}=\tilde{s}_{jk}u^{2}$,

$\displaystyle(b_{1})$ $\displaystyle=$ $\displaystyle
2^{2-4\epsilon}(p_{T}R)^{-4\epsilon}\frac{\Gamma\left(\frac{1}{2}-\epsilon\right)^{2}}{\Gamma(1-2\epsilon)}\
\frac{\int_{0}^{1}dz_{i}dz_{j}}{(4\pi)^{5-2\epsilon}\Gamma(1-2\epsilon)}(z_{i}z_{j}z_{k})^{1-2\epsilon}$
(83) $\displaystyle\times
2\int_{0}^{1}d\tilde{s}_{jk}\tilde{s}_{jk}^{-1-2\epsilon}\,\int_{0}^{1}duu^{-1-2\epsilon}\,\left(1-u\right)^{-2}{}_{2}F_{1}\left[1,\frac{1}{2}-\epsilon,1-2\epsilon,\frac{-4u}{(1-u)^{2}}\right]$
$\displaystyle=$ $\displaystyle
2^{2-4\epsilon}(p_{T}R)^{-4\epsilon}\,\frac{\Gamma\left(\frac{1}{2}-\epsilon\right)^{2}}{\Gamma(1-2\epsilon)}\,\frac{\int_{0}^{1}dz_{i}dz_{j}}{(4\pi)^{5-2\epsilon}\Gamma(1-2\epsilon)}(z_{i}z_{j}z_{k})^{1-2\epsilon}$
$\displaystyle\times
2\int_{0}^{1}d\tilde{s}_{jk}\,\tilde{s}_{jk}^{-1-2\epsilon}\int_{0}^{1}du\,u^{-1-2\epsilon}\,{}_{2}F_{1}\left[1,1+\epsilon,1-\epsilon,u^{2}\right]$
$\displaystyle=$
$\displaystyle(p_{T}R)^{-4\epsilon}\frac{\int_{0}^{1}dz_{i}dz_{j}}{(4\pi)^{4-2\epsilon}\Gamma(1-2\epsilon)}(z_{i}z_{j}z_{k})^{1-2\epsilon}\left(\frac{3}{4\epsilon^{2}}-\frac{\pi^{2}}{24}-\frac{5\zeta_{3}}{2}\epsilon-\frac{3\pi^{4}}{32}\epsilon^{2}\right)\,,$

where we have used

$\displaystyle{}_{2}F_{1}\left[1,\frac{1}{2}-\epsilon,1-2\epsilon,\frac{-4u}{(1-u)^{2}}\right]=(1-u)^{2}{}_{2}F_{1}\left[1,1+\epsilon,1-\epsilon,u^{2}\right]\,.$
(84)

For $\tilde{s}_{jk}<\tilde{s}_{ik}$, we let
$\tilde{s}_{jk}=u^{2}\tilde{s}_{ik}$ and obtain

$\displaystyle(b_{2})$ $\displaystyle=$
$\displaystyle(p_{T}R)^{-4\epsilon}\frac{\int_{0}^{1}dz_{i}dz_{j}}{(4\pi)^{4-2\epsilon}\Gamma(1-2\epsilon)}(z_{i}z_{j}z_{k})^{1-2\epsilon}\left(\frac{1}{4\epsilon^{2}}+\frac{\pi^{2}}{24}+\frac{\zeta_{3}}{2}\epsilon+\frac{19\pi^{4}}{1440}\epsilon^{2}\right)\,.$
(85)

In summary, this gives us

$\displaystyle(b)=(b_{1})+(b_{2})=\frac{(p_{T}R)^{-4\epsilon}}{\Gamma(1-2\epsilon)}\frac{\int_{0}^{1}dz_{i}dz_{j}}{(4\pi)^{4-2\epsilon}}(z_{i}z_{j}z_{k})^{1-2\epsilon}\,\theta(z_{i}<z_{j}<z_{k})\left(\frac{1}{\epsilon^{2}}-2\zeta_{3}\epsilon-\frac{29\pi^{4}}{360}\epsilon^{2}\right)\,.$
(86)

The results of $(a)$ and $(b)$ in Eqs. (81) and (86) can be used to validate
the numerical calculations. For instance, for the term proportional to
$\frac{1}{s_{12}s_{13}}$ in the splitting function
$P_{{\bar{q}}_{1}q_{2}q_{3}}^{({\rm id})}$,

$\displaystyle P^{({\rm
id})}_{{\bar{q}}_{1}q_{2}q_{3}}\left[\frac{1}{s_{12}s_{13}}\right]=4(4\pi)^{2}\alpha_{s}^{2}\left(\frac{\mu^{2}e^{\gamma_{E}}}{4\pi}\right)^{2\epsilon}\frac{-1}{2s_{12}s_{13}}z_{1}\left[\frac{1+z_{1}^{2}}{(1-z_{2})(1-z_{3})}-\epsilon\left(1+2\frac{1-z_{2}}{1-z_{3}}\right)-\epsilon^{2}\right]\,,$
(87)

we have the following result for $z_{1}<z_{2}<z_{3}$, by using $(b)$ in Eq.
(86) and after integrating over $z_{1}$ and $z_{2}$,

$\displaystyle\frac{\alpha_{s}^{2}}{(2\pi)^{2}}\,e^{4\epsilon
L}\,\left(\frac{-0.189353}{\epsilon^{2}}-\frac{1.70334}{\epsilon}-8.02839\right)\,,$
(88)

which agrees with the numerical calculation,

$\displaystyle\frac{\alpha_{s}^{2}}{4\pi^{2}}\,e^{4\epsilon
L}\left(\frac{-0.18935(1)}{\epsilon^{2}}-\frac{1.70335(7)}{\epsilon}-8.0283(5)\right)\,.$
(89)

For $z_{2}<z_{3}<z_{1}$, using $(a)$ in Eq. (81) and integrating over $z_{2}$
and $z_{3}$, we find

$\displaystyle\frac{\alpha_{s}^{2}}{4\pi^{2}}\,e^{4\epsilon
L}\left(\frac{-0.0911635}{\epsilon^{2}}-\frac{0.668211}{\epsilon}-2.3708\right)\,,$
(90)

in agreement with the numerical calculation, which returns

$\displaystyle\frac{\alpha_{s}^{2}}{4\pi^{2}}e^{4\epsilon
L}\left(\frac{-0.091162(2)}{\epsilon^{2}}-\frac{0.66821(2)}{\epsilon}-2.3709(1)\right)\,.$
(91)

The same check can be done for the splitting function $P_{g_{1}g_{2}q_{3}}$.
For instance for the contribution proportional to
$\frac{1}{s_{12}s_{13}}+(1\leftrightarrow 2)$ in the non-Abelian term
$P_{g_{1}g_{2}q_{3}}^{({\rm nab})}$

$\displaystyle P_{g_{1}g_{2}q_{3}}^{({\rm
nab})}\left[\frac{1}{s_{12}s_{13}}+(1\leftrightarrow 2)\right]$
$\displaystyle=$ $\displaystyle
4(4\pi)^{2}\alpha_{s}^{2}\left(\frac{\mu^{2}e^{\gamma_{E}}}{4\pi}\right)^{2\epsilon}\,\Bigg{\\{}\,\frac{1}{2s_{12}s_{13}}\left[\frac{2z_{3}+(1-\epsilon)(1-z_{3})^{2}}{z_{2}}+\frac{2(1-z_{2})+(1-\epsilon)z_{2}^{2}}{1-z_{3}}\right]\Bigg{\\}}$
(92) $\displaystyle+(1\leftrightarrow 2)\,,$

we find for $z_{1}<z_{2}<z_{3}$,

$\displaystyle\frac{\alpha_{s}^{2}}{4\pi^{2}}e^{4\epsilon
L}\left(\frac{3}{8\epsilon^{4}}+\frac{21+48\log
2}{32\epsilon^{3}}+\frac{2.23062}{\epsilon^{2}}-\frac{1.7911}{\epsilon}+0.703757\right)\,,$
(93)

again in complete agreement with the numerical calculations

$\displaystyle\frac{\alpha_{s}^{2}}{4\pi^{2}}e^{4\epsilon
L}\left(\frac{0.375}{\epsilon^{4}}+\frac{1.69596(5)}{\epsilon^{3}}+\frac{2.2305(2)}{\epsilon^{2}}-\frac{1.791(8)}{\epsilon}+0.703(3)\right)\,.$
(94)

Similar other checks have also been performed but are not listed here.

## References

  * (1) A. Accardi et al.,  Eur. Phys. J. C 76, 471 (2016), 1603.08906.
  * (2) R. Angeles-Martinez et al.,  Acta Phys. Polon. B 46, 2501 (2015), 1507.05267.
  * (3) D. Boer and C. Pisano,  Phys. Rev. D 91, 074024 (2015), 1412.5556.
  * (4) E.-C. Aschenauer et al.,  (2016), 1602.03922.
  * (5) STAR, L. Adamczyk et al.,  Phys. Rev. D 97, 032004 (2018), 1708.07080.
  * (6) Y.-T. Chien and I. Vitev,  Phys. Rev. Lett. 119, 112301 (2017), 1608.07283.
  * (7) M. Connors, C. Nattrass, R. Reed, and S. Salur,  Rev. Mod. Phys. 90, 025005 (2018), 1705.01974.
  * (8) J.-W. Qiu, F. Ringer, N. Sato, and P. Zurita,  Phys. Rev. Lett. 122, 252301 (2019), 1903.01993.
  * (9) A. J. Larkoski, S. Marzani, G. Soyez, and J. Thaler,  JHEP 05, 146 (2014), 1402.2657.
  * (10) A. J. Larkoski, I. Moult, and B. Nachman,  Phys. Rept. 841, 1 (2020), 1709.04464.
  * (11) A. Kardos, A. J. Larkoski, and Z. Trócsányi,  Phys. Rev. D 101, 114034 (2020), 2002.05730.
  * (12) D. Gutierrez-Reyes, I. Scimemi, W. J. Waalewijn, and L. Zoppi,  Phys. Rev. Lett. 121, 162001 (2018), 1807.07573.
  * (13) X. Liu, F. Ringer, W. Vogelsang, and F. Yuan,  Phys. Rev. Lett. 122, 192003 (2019), 1812.08077.
  * (14) X. Liu, F. Ringer, W. Vogelsang, and F. Yuan,  Phys. Rev. D 102, 094022 (2020), 2007.12866.
  * (15) M. Arratia, Z.-B. Kang, A. Prokudin, and F. Ringer,  Phys. Rev. D 102, 074015 (2020), 2007.07281.
  * (16) Z.-B. Kang, X. Liu, S. Mantry, and D. Y. Shao,  Phys. Rev. Lett. 125, 242003 (2020), 2008.00655.
  * (17) Z.-B. Kang, K. Lee, D. Y. Shao, and J. Terry,  (2020), 2008.05470.
  * (18) Y.-T. Chien et al.,  Phys. Lett. B 815, 136124 (2021), 2005.12279.
  * (19) D. Boer et al.,  (2011), 1108.1713.
  * (20) A. Accardi et al.,  Eur. Phys. J. A 52, 268 (2016), 1212.1701.
  * (21) T. Becher, M. Neubert, L. Rothen, and D. Y. Shao,  Phys. Rev. Lett. 116, 192001 (2016), 1508.06645.
  * (22) T. Becher, M. Neubert, L. Rothen, and D. Y. Shao,  JHEP 11, 019 (2016), 1605.02737,  [Erratum: JHEP 05, 154 (2017)].
  * (23) X. Liu, S.-O. Moch, and F. Ringer,  Phys. Rev. Lett. 119, 212001 (2017), 1708.04641.
  * (24) M. Dasgupta and G. P. Salam,  Phys. Lett. B 512, 323 (2001), hep-ph/0104277.
  * (25) H.-Y. Liu, Z.-B. Kang, and X. Liu,  Phys. Rev. D 102, 051502 (2020), 2004.11990.
  * (26) X. Liu and F. Petriello,  Phys. Rev. D 87, 014018 (2013), 1210.1906.
  * (27) X. Liu and F. Petriello,  Phys. Rev. D 87, 094027 (2013), 1303.4405.
  * (28) P. Sun, C. P. Yuan, and F. Yuan,  Phys. Rev. Lett. 113, 232001 (2014), 1405.1105.
  * (29) X. Liu, S.-O. Moch, and F. Ringer,  Phys. Rev. D 97, 056026 (2018), 1801.07284.
  * (30) X. Liu, F. Ringer, W. Vogelsang, and F. Yuan,  Phys. Rev. D 102, 114012 (2020), 2008.03666.
  * (31) G. Heinrich,  (2020), 2009.00516.
  * (32) M. Cacciari, G. P. Salam, and G. Soyez,  Journal of High Energy Physics 2008, 063 (2008).
  * (33) D. Bertolini, T. Chan, and J. Thaler,  JHEP 04, 013 (2014), 1310.7584.
  * (34) D. Neill, I. Scimemi, and W. J. Waalewijn,  JHEP 04, 020 (2017), 1612.04817.
  * (35) D. A. Kosower,  Nucl. Phys. B 552, 319 (1999), hep-ph/9901201.
  * (36) D. A. Kosower and P. Uwer,  Nucl. Phys. B 563, 477 (1999), hep-ph/9903515.
  * (37) Z. Bern, V. Del Duca, W. B. Kilgore, and C. R. Schmidt,  Phys. Rev. D 60, 116001 (1999), hep-ph/9903516.
  * (38) G. F. R. Sborlini, D. de Florian, and G. Rodrigo,  JHEP 01, 018 (2014), 1310.6841.
  * (39) M. Ritzmann and W. J. Waalewijn,  Phys. Rev. D 90, 054029 (2014), 1407.3272.
  * (40) S. Catani and M. Grazzini,  Nucl. Phys. B 570, 287 (2000), hep-ph/9908523.
  * (41) J. M. Campbell and E. W. N. Glover,  Nucl. Phys. B 527, 264 (1998), hep-ph/9710255.
  * (42) A. Gehrmann-De Ridder and E. W. N. Glover,  Nucl. Phys. B 517, 269 (1998), hep-ph/9707224.
  * (43) C. Anastasiou, K. Melnikov, and F. Petriello,  Phys. Rev. D 69, 076010 (2004), hep-ph/0311311.
  * (44) R. Boughezal, X. Liu, and F. Petriello,  Phys. Rev. D 91, 094035 (2015), 1504.02540.
  * (45) Z.-B. Kang, F. Ringer, and I. Vitev,  JHEP 10, 125 (2016), 1606.06732.
  * (46) L. Dai, C. Kim, and A. K. Leibovich,  Phys. Rev. D 94, 114023 (2016), 1606.07411.
  * (47) M. D. Schwartz and H. X. Zhu,  Phys. Rev. D 90, 065004 (2014), 1403.4949.
  * (48) Z. Nagy and D. E. Soper,  (2020), 2011.04773.
  * (49) S. Moch, J. A. M. Vermaseren, and A. Vogt,  Nucl. Phys. B 688, 101 (2004), hep-ph/0403192.
  * (50) S. Moch, J. A. M. Vermaseren, and A. Vogt,  JHEP 08, 049 (2005), hep-ph/0507039.