aainstitutetext: Institut für Theoretische Teilchenphysik und Kosmologie, RWTH
Aachen University,
D-52056 Aachen, Germanybbinstitutetext: Institute of Nuclear Physics, Polish
Academy of Sciences, Radzikowskiego 152,
31-342 Kraków, Poland

# Complete collection of one-loop triple-collinear splitting operators for
dimensionally-regulated QCD

Michał Czakon<EMAIL_ADDRESS>b and Sebastian Sapeta
<EMAIL_ADDRESS>
###### Abstract

We provide results for the one-loop triple-collinear color/spin-space
splitting operators for the five possible processes, $q\to
qq^{\prime}\bar{q}^{\prime}$, $q\to qq\bar{q}$, $q\to qgg$, $g\to gq\bar{q}$
and $g\to ggg$. The expressions are exact in dimensionally-regulated massless
QCD up to a single integral, which we expand to second order in the
dimensional-regularisation parameter. We also evaluate the related splitting
functions. Our results are both sufficient and indispensable for the
construction of subtraction and integrated-subtraction terms for triple-
collinear singularities of one-loop double-real-emission cross-section
contributions as part of a next-to-next-to-next-to leading order subtraction
scheme.

###### Keywords:

QCD, Scattering Amplitudes, Higher-Order Perturbative Calculations

††preprint: P3H-22-042, TTK-22-15, IFJPAN-IV-2022-7

##  1 Introduction

The study of the singular behaviour of massless-gauge-theory scattering
amplitudes in multiparticle soft/collinear limits has a long history. One of
the several important applications of the results of the analyses is the
construction of efficient methods for the evaluation of cross sections in
perturbative QCD defined within dimensional regularisation. Of particular
interest for the present work are subtraction schemes, see e.g. Refs.
Catani:1996vz ; Frixione:1995ms  for classic work at the next-to leading order
(NLO). The quest for predictions at ever higher orders of perturbation theory
has led to the construction of a plethora of schemes. The current frontier
lies between the next-to-next-to (NNLO) and next-to-next-to-next-to (N3LO)
order, see the review Ref. Caola:2022ayt  and references therein.

The known singular limits of scattering amplitudes can be classified by the
order of perturbation theory of the involved matrix elements, and by the
number of unresolved partons whose presence cannot be detected either because
they have vanishingly small (soft) energies or their momenta are parallel
(collinear) to the momentum of a single parton. A complete set of formulae at
tree-level for up to three unresolved partons has been amassed in Refs.
campbell ; Catani:1998nv ; Catani:1999ss ; DelDuca:1999iql ; Kosower:2002su ;
Catani:2019nqv ; DelDuca:2019ggv ; DelDuca:2020vst ; Braun-White:2022rtg . At
the one-loop level, complete information is available for a single unresolved
parton Bern:1994zx ; Bern:1998sc ; Kosower:1999rx ; Bern:1999ry ;
Catani:2000pi , while the case of two unresolved partons has attracted much
attention Catani:2003vu ; Sborlini:2014mpa ; Sborlini:2014kla ; Badger:2015cxa
; Zhu:2020ftr ; Catani:2021kcy . Finally, at the two-loop level only single
unresolved limits have been analysed extensively Bern:2004cz ; Badger:2004uk ;
Duhr:2014nda ; Li:2013lsa ; Duhr:2013msa .

In the present publication, we are concerned with the triple-collinear limit
of one-loop amplitudes. It is surprising that the known results for this case
do not even cover all possible splitting processes. While Ref. Badger:2015cxa
provides results for a gluon splitting into three partons, the much earlier
Ref. Catani:2003vu  contains the case of a quark splitting into a quark and a
quark-anti-quark pair of different flavor, but the result is incomplete111The
authors provide the antisymmetric part of the results with respect to the
exchange of the quark and the anti-quark momenta. The tree-level splitting is
symmetric under this transformation.. The previously cited Refs.
Sborlini:2014mpa ; Sborlini:2014kla  are rather concerned with processes
involving a photon. Thus, results for quark splittings, $q\to
qq^{\prime}\bar{q}^{\prime}$, $q\to qq\bar{q}$, $q\to qgg$, are either
incomplete or not available. In view of this situation, the following quote
from the almost twenty years old Ref. Catani:2003vu  is truly ironic:
”For the sake of brevity, we have limited ourselves, in this letter, to
presenting a few explicit results for the one-loop triple collinear splitting.
These results have mainly an illustrative purpose. The method and the tools
(in particular, the one-loop integrals) used to obtain them are sufficient and
can be applied straightforwardly to evaluate the one-loop splitting matrix of
any splitting process $a\to a_{1}+a_{2}+a_{3}$.”
There is yet another problem with the available results for the one-loop
triple-collinear splittings. They are restricted to an expansion in the
dimensional-regularisation parameter up to finite terms. Below, we explain
that this is insufficient for an application to the construction of a
subtraction scheme. This fact was, unsurprisingly, known to at least the
authors of Ref. Badger:2015cxa  as can be read in the conclusions to that
work.

In the present publication, we resolve the aforementioned issues and provide
the complete set of results for triple-collinear splittings, $q\to
qq^{\prime}\bar{q}^{\prime}$, $q\to qq\bar{q}$, $q\to qgg$, $g\to gq\bar{q}$
and $g\to ggg$, in the form necessary for the construction of an N3LO
subtraction scheme.

We aim at a self-contained but concise publication. In the first Section, we
provide the necessary definitions. Subsequently, we present the methods used
in the calculation, and explain the requirements for the construction of a
subtraction scheme in relation to the results for the triple-collinear
splittings. In the third Section, we discuss the evaluation of the occurring
Feynman integrals and provide several new results. After a short outlook on
future work, we provide more infomation on the special functions involved in
this study in two Appendices. The results obtained are lengthy. We believe
that it is crucial to provide them in an easily accessible form, and there is
no better form than a set of ancillary files that can be easily manipulated.
We use the format of the computer algebra system Mathematica Mathematica  and
describe the files in the last Appendix.

##  2 Splitting operators and splitting functions

###  2.1 Definitions and properties

We consider QCD with $n_{l}$ massless quark fields, defined through
perturbation theory supplemented with conventional dimensional
regularisation222Our results assume that gluons have $d-2$ polarisation
states. with spacetime-dimension parameter $d=4-2\epsilon$ and dimension-
setting scale $\mu$. Any parton-scattering amplitude viewed as a vector in
color- and spin-space (see Ref. Catani:1996vz  for a pedagogical introduction
to the formalism) may be expanded in the bare coupling constant, $g_{s}^{B}$:

$\ket{M}\equiv(\mu^{-\epsilon}g_{s}^{B})^{n}\Big{(}\ket{M^{(0)}}+\frac{\mu^{-2\epsilon}\alpha^{B}_{s}}{(4\pi)^{1-\epsilon}}\,\ket{M^{(1)}}+\order{\alpha_{s}^{2}}\Big{)}\;,\qquad\alpha^{B}_{s}\equiv\frac{(g^{B}_{s})^{2}}{4\pi}\;.$
(1)

Renormalisation is not essential in the present context, since splitting
operators that are the topic of this work renormalise as ordinary amplitudes.
Nevertheless, we use $\mu^{-2\epsilon}\alpha^{B}_{s}/(4\pi)^{1-\epsilon}$ as
expansion parameter as it is dimensionless and thus yields fixed-order
amplitudes $\ket{M^{(n)}}$ of integer mass dimension, because $\ket{M}$ has
this property. Furthermore, there are no $\ln(4\pi)$ terms in the results as
in the $\overline{\rm{MS}}$ renormalisation scheme.

Consider now an amplitude’s triple-collinear limit for partons
$f_{1},f_{2},f_{3}$, where $f_{i}$ may denote a gluon, a quark or an anti-
quark, with outgoing momenta $p_{i}$, $i=1,2,3$. The case of some momenta
incoming can be inferred by crossing and analytic continuation of the Feynman
integrals. Our results are provided for a final state splitting for
definiteness. In the present context, it is sufficient to consider two quark
flavors denoted generically by $q$ and $q^{\prime}$. The amplitude is known to
factorize at leading power in $s_{123}\equiv
p_{123}^{2}\equiv(p_{1}+p_{2}+p_{3})^{2}$ as discussed for the particular case
considered here in Ref. Catani:2003vu  (see also Ref. Feige:2014wja  for a
pedagogical general discussion)333We use the standard notation for asymptotic
expansions throughout this publication.:

$\ket{M_{f_{1}f_{2}f_{3}\dots}(p_{1},p_{2},p_{3},\dots)}\quad\sim\quad\textbf{Split}_{f_{1}f_{2}f_{3}}(p_{1},p_{2},p_{3})\,\ket{M_{f\dots}(p,\dots)}\qquad\big{(}s_{123}\to
0)\;.$ (2)

$\textbf{Split}_{f_{1}f_{2}f_{3}}(p_{1},p_{2},p_{3})$ is the splitting
operator that we wish to evaluate in a perturbative expansion up to
$\order{\alpha_{s}^{B}}$:

$\textbf{Split}_{f_{1}f_{2}f_{3}}(p_{1},p_{2},p_{3})=\big{(}\mu^{-\epsilon}g^{B}_{s}\big{)}^{2}\,\Big{(}\textbf{Split}^{(0)}_{f_{1}f_{2}f_{3}}(p_{1},p_{2},p_{3})+\frac{\mu^{-2\epsilon}\alpha^{B}_{s}}{(4\pi)^{1-\epsilon}}\,\textbf{Split}^{(1)}_{f_{1}f_{2}f_{3}}(p_{1},p_{2},p_{3})+\order{(\alpha_{s}^{B})^{2}}\Big{)}\;.$
(3)

The splitting operator acts on the color and polarisation indices of the
parton $f$ in $\ket{M_{f\dots}(p,\dots)}$ and yields a vector in color- and
spin-space of the nearly collinear partons $f_{1,2,3}$. The splitting parton
$f$ is determined from $f_{1,2,3}$ by flavor conservation, while the
factorization is only valid if $f$ is well defined. The momentum, $p$, of the
splitting parton $f$ is defined with the help of an auxiliary light-like
vector $q$ as follows:

$p_{123}\equiv p+\frac{s_{123}}{2\,(p_{123}\cdot q)}\,q\;,\qquad
p^{2}=q^{2}=0\;,\quad p\cdot q\neq 0\;,\quad p^{0},q^{0}>0\;.$ (4)

Suitable $p$ and $q$ may always be uniquely constructed up to normalisation of
$q$, assuming that the three-vectors, $\bm{q}$ and $\bm{p}$, are parallel,
$\bm{q}\parallel\bm{p}\parallel\bm{p}_{123}$. Since the singularity structure
of one-loop QCD amplitudes is completely known Catani:2000ef , it is possible
to derive the singularity structure of the one-loop splitting operators
Catani:2003vu :

$\begin{split}\textbf{Split}&{}^{(1)}_{f_{1}f_{2}f_{3}}(p_{1},p_{2},p_{3})=\\\\[5.69046pt]
&\frac{\Gamma(1+\epsilon)\Gamma^{2}(1-\epsilon)}{\Gamma(1-2\epsilon)}\Bigg{[}\frac{2}{\epsilon^{2}}\sum_{1\leq
i<j\leq 3}\textbf{T}_{i}\cdot\textbf{T}_{j}\,\Bigg{(}\frac{-2(p_{i}\cdot
p_{j})-i0^{+}}{\mu^{2}}\Bigg{)}^{-\epsilon}\\\\[5.69046pt]
&\qquad\qquad\qquad\qquad+\Bigg{(}\frac{-s_{123}-i0^{+}}{\mu^{2}}\Bigg{)}^{-\epsilon}\Bigg{(}\frac{2}{\epsilon^{2}}\sum_{1\leq
i\leq j\leq
3}\textbf{T}_{i}\cdot\textbf{T}_{j}\,\Big{(}2-\Big{(}\frac{p_{i}\cdot
q}{p_{123}\cdot q}\Big{)}^{-\epsilon}-\Big{(}\frac{p_{j}\cdot q}{p_{123}\cdot
q}\Big{)}^{-\epsilon}\Big{)}\\\\[5.69046pt]
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\frac{1}{\epsilon}\Big{(}b_{0}+\gamma_{f}-\sum_{1\leq
i\leq
3}\gamma_{i}\Big{)}\Bigg{)}\Bigg{]}\,\textbf{Split}^{(0)}_{f_{1}f_{2}f_{3}}(p_{1},p_{2},p_{3})\\\\[5.69046pt]
&+\order{\epsilon^{0}}\;,\end{split}$ (5)

where $\textbf{T}_{i}$ are color-space operators corresponding to the
generators, $T^{a}_{c^{\prime}c}$, of the fundamental representation for
quarks, the generators, $-T^{a\,\star}_{c^{\prime}c}$, of the anti-fundamental
representation for anti-quarks, and the generators, $if^{c^{\prime}ac}$, of
the adjoint representation for gluons. Furthermore,
$\gamma_{q}=\gamma_{\bar{q}}=3/2C_{F}$, $\gamma_{g}=b_{0}/2$,
$b_{0}=(11/3C_{A}-4/3T_{F}n_{l})$, $C_{F}=(N_{c}^{2}-1)/2N_{c}=4/3$,
$C_{A}=N_{c}=3$. The presence of the leading coefficient of the QCD
$\beta$-function, $b_{0}$, is due to the lack of renormalisation. We stress
that, in the above formula, all terms of $\order{\epsilon^{0}}$ are arbitrary.
For convenience of the reader, we provide the expressions obtained from Eq.
(5) for each of the splitting operators in an ancillary file attached to this
publication. There, we only keep the pure-pole contributions, proportional to
$1/\epsilon^{2}$ and $1/\epsilon$, in the pre-factor in Eq. (5) acting on the
tree-level splitting operators
$\textbf{Split}^{(0)}_{f_{1}f_{2}f_{3}}(p_{1},p_{2},p_{3})$.

Figure 1: 9 diagrams contributing to the splitting operator $q\to
qq^{\prime}\bar{q}^{\prime}$. The thick incoming lines on the left side of
each diagram correspond to the off-shell splitting parton. The outgoing lines
on the right side of each diagram are on shell. The anti-quark line carries an
arrow of direction opposite to that of the energy flow. Figure 2: 18 diagrams
contributing to the splitting operator $q\to qq\bar{q}$. Description as in
Fig. 1. Figure 3: 30 diagrams contributing to the splitting operator $q\to
qgg$ Description as in Fig. 1. Figure 4: 33 diagrams contributing to the
splitting operator $g\to gq\bar{q}$. Description as in Fig. 1. Figure 5: Part
of the 68 diagrams contributing to the splitting operator $g\to ggg$.
Description as in Fig. 1. Figure 6: Part of the 68 diagrams contributing to
the splitting operator $g\to ggg$. Description as in Fig. 1.

The splitting operators for the five possible splittings444Anti-quark
splitting operators are related to those of quarks by charge-conjugation
symmetry of QCD., $q\to qq^{\prime}\bar{q}^{\prime}$, $q\to qq\bar{q}$, $q\to
qgg$, $g\to gq\bar{q}$ and $g\to ggg$ are obtained from diagrams depicted in
Figs. 1, 2, 3, 4 and 5, 6 respectively. The incoming off-shell line of the
splitting parton is contracted with a massless Dirac spinor $u_{0}\equiv
u(p,\lambda_{0})$ for a quark splitting, and with a massless transverse
polarisation vector $\epsilon_{0}\equiv\epsilon(p,\lambda_{0})$ for a gluon
splitting, with $\lambda_{0}$ the helicity of the splitting parton. This
prescription obviously realises Eq. (2) by approximating the numerator of the
splitting parton’s propagator by a polarisation sum of a product of on-shell
spinor/polarisation vectors. The only non-trivial issue is the necessity of
using a light-cone-gauge gluon propagator:

$\expectationvalue{\tilde{A}_{\mu}^{a}(k)A_{\nu}^{b}(0)}{0}=\frac{i\delta^{ab}}{k^{2}+i0^{+}}\Big{(}-g_{\mu\nu}+\frac{q^{\mu}k^{\nu}+q^{\nu}k^{\mu}}{q\cdot
k}\Big{)}\;,$ (6)

where $A_{\mu}^{a}$ is the free gluon field, tilde denotes Fourier
transformation, and $q$, $q^{2}=0$ is an auxilliary light-like vector. This
is, of course, part of the proof of the factorisation theorem. For more
details we refer the reader again to Ref. Feige:2014wja .

Besides splitting operators, one also defines dimensionless splitting
functions and averaged splitting functions:

$\bm{\hat{P}}_{f_{1}f_{2}f_{3}}\equiv\Big{(}\tfrac{1}{2}s_{123}\Big{)}^{2}\,\textbf{Split}^{\dagger}_{f_{1}f_{2}f_{3}}\textbf{Split}_{f_{1}f_{2}f_{3}}\;,\qquad\expectationvalue{\hat{P}_{f_{1}f_{2}f_{3}}}\equiv\frac{1}{n_{f}^{c}n^{s}_{f}}\Tr(\bm{\hat{P}}_{f_{1}f_{2}f_{3}})\;,$
(7)

where $n_{f}^{c}$ and $n_{f}^{s}$ are the numbers of color and spin degrees-
of-freedom of the splitting parton, $n^{c}_{q}=3$, $n^{c}_{g}=8$,
$n^{s}_{q}=2$ and $n^{s}_{g}=d-2$. The color- and spin-space operators
$\bm{\hat{P}}_{f_{1}f_{2}f_{3}}$ are relevant for the study of the triple-
collinear limit of amplitudes squared and summed over colors and polarisations
of the external states. By color conservation, they are proportional to the
identity operator in color space. Due to helicity conservation in massless
QCD, the same holds for spin space in the case of quarks. The corresponding
averages $\expectationvalue{\hat{P}_{f_{1}f_{2}f_{3}}}$ are necessary for the
determination of the integrated subtraction terms in the construction of
subtraction schemes, see Section 2.3. As a result of our work, we provide:

$\bm{\hat{P}}^{(1)}_{f_{1}f_{2}f_{3}}\equiv\Big{(}\tfrac{1}{2}s_{123}\Big{)}^{2}\,\Big{(}\textbf{Split}^{(0)\,\dagger}_{f_{1}f_{2}f_{3}}\textbf{Split}^{(1)}_{f_{1}f_{2}f_{3}}+\text{hermitian
conjugate}\Big{)}\;,$ (8)

as well as its averaged version.

###  2.2 Evaluation of the diagrams

The calculation proceeds in standard fashion. The diagrams are generated with
private software (although the figures have been produced with FeynArts
Hahn:2000kx ) and simplified with the program FORM Vermaseren:2000nd .
Subsequently, we perform a Passarino-Veltman Passarino:1978jh  reduction of
the tensor integrals, consisting of expressing integrals with numerators
proportional to $l^{\mu_{1}}\cdots l^{\mu_{n}}$, with $l$ the loop momentum
and $n$ the tensor rank, as sums of symmetric tensors built of products of
$p_{1,2,3}$, $q$ and the metric tensor $g$, multiplied with scalar integrals
with each $l^{\mu_{i}}$ contracted with another momentum. The reduction
formulae are efficiently derived with the help of the program Fermat fermat .
In the case of Feynman integrals whose denominators contain up to three
different momenta from the set $\\{p_{1},p_{2},p_{3},q\\}$, we encounter
tensors of rank up to four. In the case of integrals whose denominators
contain all four momenta, the maximal tensor rank is two. Since four linearly-
independent vectors form a basis in four-dimensional space, the metric tensor
itself can be expressed as a linear combination of products of two out of the
four momenta. Hence, the procedure is singular in four dimensions. As we are
working in dimensional regularisation, this leads to the appearance of a pole,
$1/\epsilon$, in the coefficients of the Passarino-Veltman reduction. It has
been noticed in the study of ordinary (i.e. without linear propagators) tensor
pentagon integrals in Ref. Bern:1993kr  that these poles are removed by
switching from four-dimensional to six-dimensional pentagon integrals. The
same phenomenon occurs here as well, see Section 3.2. After Passarino-Veltman
reduction, we perform a further integration-by-parts reduction using the
software package Kira Maierhofer:2017gsa ; Klappert:2020nbg  for convenience,
since it allows for linear propagators. The number of occurring integrals is
reduced to thirty four, see Section 3. The reduction introduces spurious
$1/\epsilon$ poles in coefficients of the bubble integrals Eq. (52) and (53).
This is the price for the removal of triangle integrals, i.e. integrals with
three ordinary Feynman propagators. Since the bubble integrals are known
exactly, these spurious poles do not constitute a difficulty.

At this point, the results for the splitting operators still have a
complicated color and spin structure. Both of these structures can be
simplified by the following two algorithms. As far as color factors are
concerned, we use the formulae:

$if^{abc}=\frac{1}{T_{F}}\Tr(\commutator{T^{a}}{T^{b}}T^{c})\;,\qquad
T^{a}_{ij}T^{a}_{kl}=T_{F}\Big{(}\delta_{il}\delta_{kj}-\frac{1}{N_{c}}\delta_{ij}\delta_{kl}\Big{)}\;,\qquad
T_{F}=\frac{1}{2}\;,\quad N_{c}=3\;.$ (9)

This procedure is commonly referred to as the Cvitanovic algorithm
Cvitanovic:1976am . For definiteness, we assign the following color indices:
$c_{0}$ to the splitting parton, and $c_{i}$, $i=1,2,3$ to the outgoing
partons. This algorithm provides a unique basis of color structures. Below, we
list the occuring ones for each of the splitting functions.

The spin structure is simplified in several steps. First, $\not{q}$ is shifted
to the right in any product of Dirac $\gamma$-matrices using the Dirac
algebra. If possible, this step is followed by the application of:

$\not{q}\,u_{0}=\frac{2(p_{123}\cdot
q)}{p_{123}^{2}}\,\not{p}_{123}\,u_{0}\;.$ (10)

Subsequently, equations of motion for spinors
$\bar{u}_{i}\equiv\bar{u}(p_{i},\lambda_{i})$, $i=1,2$, $v_{3}\equiv
v(p_{3},\lambda_{3})$:

$\bar{u}_{1}\,\not{p}_{1}=\bar{u}_{2}\,\not{p}_{2}=\not{p}_{3}v_{3}=0\;,$ (11)

are used, once $\not{p}_{1,2}$ has been shifted to the left, and $\not{p}_{3}$
to the right in any product of $\gamma$-matrices. Polarisation vectors
$\not{\epsilon}_{0}$ and $\not{\epsilon}^{*}_{1,2}$ are shifted to the left.
The assumed transversality of the polarisation vectors to not only the
respective momentum, but also to $q$:

$p\cdot\epsilon_{0}=q\cdot\epsilon_{0}=p_{i}\cdot\epsilon^{*}_{i}=q\cdot\epsilon^{*}_{i}=0\;,\qquad
i=1,2,3\;,$ (12)

is used at each step together with the implied relation:

$p_{3}\cdot\epsilon_{0}=-(p_{1}+p_{2})\cdot\epsilon_{0}\;.$ (13)

This last relation breaks the explicit permutation symmetry in the case of the
$g\to ggg$ splitting, but is necessary to obtain a basis of spin-space
structures without redundancy.

With the above, the expressions for the splitting operators attain their final
unique representation in dimensional regularisation. The spin structures are,
however, not linearly independent in four dimensions. In Ref. Bern:1993kr , it
was pointed out that six-dimensional scalar pentagon integrals disappear from
the finite $\order{\epsilon^{0}}$ terms of the $\epsilon$-expansion of
scattering amplitudes, if four-dimensional relations between spin structures
are exploited. The same happens to hold in our case for six-dimensional
integrals containing four ordinary and one linear propagator discussed in
Section 3.2. In order to make this fact explicit, we represent any vector,
including $\gamma$-matrices in a basis of the four vectors $k_{i}\equiv
p_{i}$, $i=1,2,3$, $k_{4}\equiv q$:

$v^{\mu}=\sum_{i,j=1}^{4}k_{i}^{\mu}\,\big{(}K^{-1}\big{)}_{ij}\,k_{j}\cdot
v\;,\qquad K_{ij}\equiv k_{i}\cdot k_{j}\;.$ (14)

The resulting splitting operators are only valid in four dimensions. On the
other hand, the remaining spinor chains in the case of splittings involving
quarks are vastly simplified and can be trivially expressed in the popular
spinor-helicity formalism using $v_{3\pm}\;\propto\;u_{3\mp}$ (up to a phase
factor) and with $i,k\in\\{1,2,3\\}$, $j\in\\{0,\dots,3\\}$:

$\begin{gathered}\bar{u}_{i\pm}\,\not{p}_{k}\,u_{j\pm}=\big{(}\bar{u}_{i\pm}u_{k\mp}\big{)}\,\big{(}\bar{u}_{k\mp}u_{j\pm}\big{)}\;,\qquad\bar{u}_{i\pm}\,\not{q}\,u_{j\pm}=2\sqrt{(p_{i}\cdot
q)(p_{j}\cdot q)}\;,\\\\[5.69046pt]
p_{i}\cdot\epsilon^{*}_{j\pm}=\pm\frac{\big{(}\bar{u}_{j\pm}u_{i\mp}\big{)}\,\big{(}\bar{u}_{i\mp}u_{\pm}(q)\big{)}}{\sqrt{2}\bar{u}_{\mp}(q)u_{j\pm}}\;,\end{gathered}$
(15)

where the subscripts $\pm$ denote helicity, while the unlisted cases vanish.
As long as there are spinor chains, this procedure does not yield larger
expressions than those valid in general $d$ dimensions. As splitting operators
expanded up to $\order{\epsilon^{0}}$ are useful to subtract the singularities
of exact amplitudes in practical calculations, we also provide them as
ancillary files attached to this publication, just as we do with for the exact
results which are far too lengthy to be useful in printed text.

In order to illustrate the structure of the splitting operators, we now list
the occurring color and spin structures for each of them.

#### $q\to q_{1}q^{\prime}_{2}\bar{q}^{\prime}_{3}$

In this case, there are two color structures both at tree- and one-loop level:

$\delta_{c_{1}c_{0}}\delta_{c_{2}c_{3}}\;,\qquad\delta_{c_{2}c_{0}}\delta_{c_{1}c_{3}}\;.$
(16)

At tree level, nevertheless, they are actually generated from a single diagram
proportional to $T^{a}_{c_{1}c_{0}}\,T^{a}_{c_{2}c_{3}}$, which is split into
two terms by the Cvitanovic algorithm according to Eqs. (9). Furthermore,
there are only three distinct spin structures at tree level:

$\,\bar{u}_{1}\,\gamma^{\mu}\,u_{0}\,\bar{u}_{2}\,\gamma_{\mu}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,\bar{u}_{2}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\not{q}\,v_{3}\,\;,$
(17)

while there appear the following additional structures at one-loop level:

$\displaystyle\,\bar{u}_{1}\,\gamma^{\mu_{1}}\,\gamma^{\mu_{2}}\,\gamma^{\mu_{3}}\,u_{0}\,\bar{u}_{2}\,\gamma_{\mu_{1}}\,\gamma_{\mu_{2}}\,\gamma_{\mu_{3}}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\gamma^{\mu}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,\gamma_{\mu}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\gamma^{\nu_{1}}\,\gamma^{\nu_{2}}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\gamma_{\nu_{1}}\,\gamma_{\nu_{2}}\,\not{q}\,v_{3}\,\;,$
$\displaystyle\,\bar{u}_{1}\,\gamma^{\nu_{1}}\,\gamma^{\nu_{2}}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,\gamma_{\nu_{1}}\,\gamma_{\nu_{2}}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,\gamma^{\mu}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\gamma_{\mu}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,\gamma^{\mu}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,\gamma_{\mu}\,\not{q}\,v_{3}\,\;,$
$\displaystyle\,\bar{u}_{1}\,\not{p}_{2}\,\gamma^{\nu_{1}}\,\gamma^{\nu_{2}}\,u_{0}\,\bar{u}_{2}\,\gamma_{\nu_{1}}\,\gamma_{\nu_{2}}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,\gamma^{\nu_{1}}\,\gamma^{\nu_{2}}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,\gamma_{\nu_{1}}\,\gamma_{\nu_{2}}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,\;,$
$\displaystyle\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,\;.$
(18)

Since the singularities of the one-loop splitting operators given in Eq. (5)
have the same spin structure as the tree-level splitting operators, the
coefficients of the structures (18) are regular in $\epsilon$. Finally, the
four-dimensional projection Eq. (14) transforms the linear combination of (17)
and (18) into a linear combination of just:

$\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,\bar{u}_{2}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\not{q}\,v_{3}\,\;.$
(19)

#### $q\to q_{1}q_{2}\bar{q}_{3}$

Up to the additional permutation $1\leftrightarrow 2$, the color and spin
structures in this case are very similar to those of the $q\to
q_{1}q^{\prime}_{2}\bar{q}^{\prime}_{3}$ case. The occurring color structures
are, in fact, the same:

$\delta_{c_{1}c_{0}}\delta_{c_{2}c_{3}}\;,\qquad\delta_{c_{2}c_{0}}\delta_{c_{1}c_{3}}\;.$
(20)

The number of spin structures at tree-level is doubled:

$\displaystyle\,\bar{u}_{1}\,\gamma^{\mu}\,u_{0}\,\bar{u}_{2}\,\gamma_{\mu}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,\bar{u}_{2}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\not{q}\,v_{3}\,\;,$
$\displaystyle\,\bar{u}_{2}\,\gamma^{\mu}\,u_{0}\,\bar{u}_{1}\,\gamma_{\mu}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\not{p}_{1}\,u_{0}\,\bar{u}_{1}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\not{p}_{3}\,u_{0}\,\bar{u}_{1}\,\not{q}\,v_{3}\,\;,$
(21)

while the additional spin structures with finite coefficients at one-loop
level are:

$\displaystyle\,\bar{u}_{1}\,\gamma^{\mu_{1}}\,\gamma^{\mu_{2}}\,\gamma^{\mu_{3}}\,u_{0}\,\bar{u}_{2}\,\gamma_{\mu_{1}}\,\gamma_{\mu_{2}}\,\gamma_{\mu_{3}}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\gamma^{\mu}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,\gamma_{\mu}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\gamma^{\nu_{1}}\,\gamma^{\nu_{2}}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\gamma_{\nu_{1}}\,\gamma_{\nu_{2}}\,\not{q}\,v_{3}\,\;,$
$\displaystyle\,\bar{u}_{1}\,\gamma^{\nu_{1}}\,\gamma^{\nu_{2}}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,\gamma_{\nu_{1}}\,\gamma_{\nu_{2}}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,\gamma^{\mu}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\gamma_{\mu}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,\gamma^{\mu}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,\gamma_{\mu}\,\not{q}\,v_{3}\,\;,$
$\displaystyle\,\bar{u}_{1}\,\not{p}_{2}\,\gamma^{\nu_{1}}\,\gamma^{\nu_{2}}\,u_{0}\,\bar{u}_{2}\,\gamma_{\nu_{1}}\,\gamma_{\nu_{2}}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,\gamma^{\nu_{1}}\,\gamma^{\nu_{2}}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,\gamma_{\nu_{1}}\,\gamma_{\nu_{2}}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,\;,$
$\displaystyle\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\gamma^{\mu_{1}}\,\gamma^{\mu_{2}}\,\gamma^{\mu_{3}}\,u_{0}\,\bar{u}_{1}\,\gamma_{\mu_{1}}\,\gamma_{\mu_{2}}\,\gamma_{\mu_{3}}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\gamma^{\mu}\,u_{0}\,\bar{u}_{1}\,\not{p}_{2}\,\gamma_{\mu}\,\not{q}\,v_{3}\,\;,$
$\displaystyle\,\bar{u}_{2}\,\gamma^{\nu_{1}}\,\gamma^{\nu_{2}}\,\not{p}_{3}\,u_{0}\,\bar{u}_{1}\,\gamma_{\nu_{1}}\,\gamma_{\nu_{2}}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\gamma^{\nu_{1}}\,\gamma^{\nu_{2}}\,\not{p}_{3}\,u_{0}\,\bar{u}_{1}\,\not{p}_{2}\,\gamma_{\nu_{1}}\,\gamma_{\nu_{2}}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\not{p}_{1}\,\gamma^{\mu}\,\not{p}_{3}\,u_{0}\,\bar{u}_{1}\,\gamma_{\mu}\,v_{3}\,\;,$
$\displaystyle\,\bar{u}_{2}\,\not{p}_{1}\,\gamma^{\mu}\,\not{p}_{3}\,u_{0}\,\bar{u}_{1}\,\not{p}_{2}\,\gamma_{\mu}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\not{p}_{1}\,\gamma^{\nu_{1}}\,\gamma^{\nu_{2}}\,u_{0}\,\bar{u}_{1}\,\gamma_{\nu_{1}}\,\gamma_{\nu_{2}}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\not{p}_{1}\,\gamma^{\nu_{1}}\,\gamma^{\nu_{2}}\,u_{0}\,\bar{u}_{1}\,\not{p}_{2}\,\gamma_{\nu_{1}}\,\gamma_{\nu_{2}}\,v_{3}\,\;,$
$\displaystyle\,\bar{u}_{2}\,\not{p}_{1}\,u_{0}\,\bar{u}_{1}\,\not{p}_{2}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\not{p}_{3}\,u_{0}\,\bar{u}_{1}\,\not{p}_{2}\,v_{3}\,\;.$
(22)

Finally, the result of the four-dimensional projection Eq. (14) of the spin
structures is given by:

$\displaystyle\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,\bar{u}_{2}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,\bar{u}_{2}\,\not{q}\,v_{3}\,\;,$
$\displaystyle\,\bar{u}_{2}\,\not{p}_{1}\,u_{0}\,\bar{u}_{1}\,\not{p}_{2}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\not{p}_{1}\,u_{0}\,\bar{u}_{1}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\not{p}_{3}\,u_{0}\,\bar{u}_{1}\,\not{p}_{2}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\not{p}_{3}\,u_{0}\,\bar{u}_{1}\,\not{q}\,v_{3}\,\;.$

#### $q\to q_{1}g_{2}g_{3}$

In this case, there are two color structures at tree-level reflecting the
symmetry of the splitting operator with respect to the exchange of the two
gluons:

$\big{(}T^{c_{2}}T^{c_{3}}\big{)}_{c_{1}c_{0}}\;,\qquad\big{(}T^{c_{3}}T^{c_{2}}\big{)}_{c_{1}c_{0}}\;.$
(24)

The one additional structure at one-loop level is already symmetric:

$\delta_{c_{1}c_{0}}\Tr(T^{c_{2}}T^{c_{3}})\;.$ (25)

The spin-structures at tree-level are as follows:

$\displaystyle\,\bar{u}_{1}\,\not{\epsilon}^{*}_{2}\,\not{\epsilon}^{*}_{3}\,\not{p}_{3}\,u_{0}\,\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,\not{\epsilon}^{*}_{2}\,\not{\epsilon}^{*}_{3}\,u_{0}\,\;,\qquad\,\bar{u}_{1}\,\not{\epsilon}^{*}_{2}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{3})\;,\qquad\,\bar{u}_{1}\,\not{\epsilon}^{*}_{2}\,u_{0}\,(p_{2}\cdot\epsilon^{*}_{3})\;,$
$\displaystyle\,\bar{u}_{1}\,\not{\epsilon}^{*}_{3}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{2})\;,\qquad\,\bar{u}_{1}\,\not{\epsilon}^{*}_{3}\,u_{0}\,(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,(\epsilon^{*}_{2}\cdot\epsilon^{*}_{3})\;,\qquad\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,(\epsilon^{*}_{2}\cdot\epsilon^{*}_{3})\;,$
(26)

while the list is extended at one-loop by:

$\displaystyle\,\bar{u}_{1}\,\not{p}_{2}\,\not{\epsilon}^{*}_{2}\,\not{p}_{3}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{3})\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,\not{\epsilon}^{*}_{2}\,\not{p}_{3}\,u_{0}\,(p_{2}\cdot\epsilon^{*}_{3})\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,\not{\epsilon}^{*}_{3}\,\not{p}_{3}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{2})\;,$
$\displaystyle\,\bar{u}_{1}\,\not{p}_{2}\,\not{\epsilon}^{*}_{3}\,\not{p}_{3}\,u_{0}\,(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon^{*}_{3})\;,$
$\displaystyle\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})\;,$
$\displaystyle\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon^{*}_{3})\;,\qquad\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{2})\;.$
(27)

In both cases, Bose symmetry is not explicit due to the simplification
algorithm. While this admittedly makes the expressions less elegant, it has no
influence on actual calculations using the splitting operators. Finally, the
four-dimensional projection Eq. (14) yields:

$\displaystyle\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon^{*}_{3})\;,\qquad\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{2})\;,$
$\displaystyle\,\bar{u}_{1}\,\not{p}_{2}\,u_{0}\,(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})\;,\qquad\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon^{*}_{3})\;,$
$\displaystyle\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,(p_{1}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad\,\bar{u}_{1}\,\not{p}_{3}\,u_{0}\,(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{2})\;.$
(28)

#### $g\to g_{1}q_{2}\bar{q}_{3}$

The external states for this splitting operator are the same as in the
previous case. This translates into the same color structures up to index
permutation. Hence, at tree level, the color structures are:

$\big{(}T^{c_{0}}T^{c_{1}}\big{)}_{c_{3}c_{2}}\;,\qquad\big{(}T^{c_{1}}T^{c_{0}}\big{)}_{c_{3}c_{2}}\;,$
(29)

while at one-loop there again appears:

$\delta_{c_{2}c_{3}}\Tr(T^{c_{0}}T^{c_{1}})\;.$ (30)

The spin structures are not related in this trivial way, because the
simplification algorithm is not symmetric. We observe, for instance, that
there are less spin structures at tree level:

$\displaystyle\,\bar{u}_{2}\,\not{p}_{1}\,\not{\epsilon}^{*}_{1}\,\not{\epsilon}_{0}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\not{\epsilon}_{0}\,v_{3}\,(p_{2}\cdot\epsilon^{*}_{1})\;,\qquad\,\bar{u}_{2}\,\not{\epsilon}_{0}\,v_{3}\,(p_{3}\cdot\epsilon^{*}_{1})\;,\qquad\,\bar{u}_{2}\,\not{\epsilon}^{*}_{1}\,v_{3}\,(p_{1}\cdot\epsilon_{0})\;,$
$\displaystyle\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,(\epsilon_{0}\cdot\epsilon^{*}_{1})\;,\qquad\,\bar{u}_{2}\,\not{q}\,v_{3}\,(\epsilon_{0}\cdot\epsilon^{*}_{1})\;,$
(31)

but more additional spin structures with finite coefficients at one-loop
level:

$\displaystyle\,\bar{u}_{2}\,\not{\epsilon}^{*}_{1}\,\not{\epsilon}_{0}\,\not{q}\,v_{3}\,\;,\qquad\,\bar{u}_{2}\,\not{\epsilon}^{*}_{1}\,v_{3}\,(p_{2}\cdot\epsilon_{0})\;,\qquad\,\bar{u}_{2}\,\not{p}_{1}\,\not{\epsilon}_{0}\,\not{q}\,v_{3}\,(p_{2}\cdot\epsilon^{*}_{1})\;,\qquad\,\bar{u}_{2}\,\not{p}_{1}\,\not{\epsilon}_{0}\,\not{q}\,v_{3}\,(p_{3}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle\,\bar{u}_{2}\,\not{p}_{1}\,\not{\epsilon}^{*}_{1}\,\not{q}\,v_{3}\,(p_{1}\cdot\epsilon_{0})\;,\qquad\,\bar{u}_{2}\,\not{p}_{1}\,\not{\epsilon}^{*}_{1}\,\not{q}\,v_{3}\,(p_{2}\cdot\epsilon_{0})\;,\qquad\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,(p_{1}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,(p_{1}\cdot\epsilon_{0})(p_{3}\cdot\epsilon^{*}_{1})\;,\qquad\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})\;,\qquad\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,(p_{2}\cdot\epsilon_{0})(p_{3}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle\,\bar{u}_{2}\,\not{q}\,v_{3}\,(p_{1}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})\;,\qquad\,\bar{u}_{2}\,\not{q}\,v_{3}\,(p_{1}\cdot\epsilon_{0})(p_{3}\cdot\epsilon^{*}_{1})\;,\qquad\,\bar{u}_{2}\,\not{q}\,v_{3}\,(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle\,\bar{u}_{2}\,\not{q}\,v_{3}\,(p_{2}\cdot\epsilon_{0})(p_{3}\cdot\epsilon^{*}_{1})\;.$
(32)

A simple relation between the spin structures for $q\to qgg$ and $g\to
gq\bar{q}$ is restored after four-dimensional projection, Eq. (14), as seen by
comparing (28) with:

$\displaystyle\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,(p_{1}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})\;,\qquad\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,(p_{1}\cdot\epsilon_{0})(p_{3}\cdot\epsilon^{*}_{1})\;,\qquad\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle\,\bar{u}_{2}\,\not{p}_{1}\,v_{3}\,(p_{2}\cdot\epsilon_{0})(p_{3}\cdot\epsilon^{*}_{1})\;,\qquad\,\bar{u}_{2}\,\not{q}\,v_{3}\,(p_{1}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})\;,\qquad\,\bar{u}_{2}\,\not{q}\,v_{3}\,(p_{1}\cdot\epsilon_{0})(p_{3}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle\,\bar{u}_{2}\,\not{q}\,v_{3}\,(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})\;,\qquad\,\bar{u}_{2}\,\not{q}\,v_{3}\,(p_{2}\cdot\epsilon_{0})(p_{3}\cdot\epsilon^{*}_{1})\;.$
(33)

#### $g\to g_{1}g_{2}g_{3}$

The color structures of the tree-level pure-gluon splitting operator consist,
as for ordinary tree-level gluon-scattering amplitudes, of a trace of a
product of fundamental representation generators in all possible color-index
permutations:

$\displaystyle\Tr(T^{c_{0}}T^{c_{1}}T^{c_{2}}T^{c_{3}})\;,\qquad\Tr(T^{c_{0}}T^{c_{1}}T^{c_{3}}T^{c_{2}})\;,\qquad\Tr(T^{c_{0}}T^{c_{2}}T^{c_{1}}T^{c_{3}})\;,$
$\displaystyle\Tr(T^{c_{0}}T^{c_{2}}T^{c_{3}}T^{c_{1}})\;,\qquad\Tr(T^{c_{0}}T^{c_{3}}T^{c_{1}}T^{c_{2}})\;,\qquad\Tr(T^{c_{0}}T^{c_{3}}T^{c_{2}}T^{c_{1}})\;.$
(34)

Disconnected traces appear at one-loop only:

$\displaystyle\Tr(T^{c_{0}}T^{c_{1}})\Tr(T^{c_{2}}T^{c_{3}})\;,\qquad\Tr(T^{c_{0}}T^{c_{2}})\Tr(T^{c_{1}}T^{c_{3}})\;,\qquad\Tr(T^{c_{0}}T^{c_{3}})\Tr(T^{c_{1}}T^{c_{2}})\;.$
(35)

The spin structures consist of various contractions of polarisation vectors
amongst themselves and with external momenta. With the present algorithm, the
list at tree-level reads:

$\displaystyle(\epsilon_{0}\cdot\epsilon^{*}_{1})(\epsilon^{*}_{2}\cdot\epsilon^{*}_{3})\;,\qquad(\epsilon_{0}\cdot\epsilon^{*}_{2})(\epsilon^{*}_{1}\cdot\epsilon^{*}_{3})\;,\qquad(\epsilon_{0}\cdot\epsilon^{*}_{3})(\epsilon^{*}_{1}\cdot\epsilon^{*}_{2})\;,$
$\displaystyle(p_{1}\cdot\epsilon^{*}_{2})(\epsilon^{*}_{1}\cdot\epsilon^{*}_{3})(p_{1}\cdot\epsilon_{0})\;,\qquad(p_{1}\cdot\epsilon^{*}_{3})(\epsilon_{0}\cdot\epsilon^{*}_{1})(p_{1}\cdot\epsilon^{*}_{2})\;,\qquad(p_{2}\cdot\epsilon_{0})(\epsilon^{*}_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})\;,$
$\displaystyle(p_{2}\cdot\epsilon_{0})(\epsilon^{*}_{1}\cdot\epsilon^{*}_{3})(p_{1}\cdot\epsilon^{*}_{2})\;,\qquad(p_{2}\cdot\epsilon^{*}_{1})(\epsilon_{0}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})\;,\qquad(p_{2}\cdot\epsilon^{*}_{1})(\epsilon^{*}_{2}\cdot\epsilon^{*}_{3})(p_{1}\cdot\epsilon_{0})\;,$
$\displaystyle(p_{2}\cdot\epsilon^{*}_{1})(\epsilon^{*}_{2}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon_{0})\;,\qquad(p_{2}\cdot\epsilon^{*}_{3})(\epsilon_{0}\cdot\epsilon^{*}_{1})(p_{1}\cdot\epsilon^{*}_{2})\;,\qquad(p_{2}\cdot\epsilon^{*}_{3})(\epsilon_{0}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle(p_{2}\cdot\epsilon^{*}_{3})(\epsilon^{*}_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon_{0})\;,\qquad(p_{3}\cdot\epsilon^{*}_{1})(\epsilon_{0}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon^{*}_{3})\;,\qquad(p_{3}\cdot\epsilon^{*}_{1})(\epsilon_{0}\cdot\epsilon^{*}_{3})(p_{1}\cdot\epsilon^{*}_{2})\;,$
$\displaystyle(p_{3}\cdot\epsilon^{*}_{1})(\epsilon^{*}_{2}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon_{0})\;,\qquad(p_{3}\cdot\epsilon^{*}_{2})(\epsilon_{0}\cdot\epsilon^{*}_{1})(p_{1}\cdot\epsilon^{*}_{3})\;,\qquad(p_{3}\cdot\epsilon^{*}_{2})(\epsilon_{0}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle(p_{3}\cdot\epsilon^{*}_{2})(\epsilon_{0}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle(p_{3}\cdot\epsilon^{*}_{2})(\epsilon^{*}_{1}\cdot\epsilon^{*}_{3})(p_{1}\cdot\epsilon_{0})\;,$
(36)

and it is extended at one-loop by:

$\displaystyle(p_{1}\cdot\epsilon_{0})(p_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon^{*}_{1})\;,\qquad(p_{1}\cdot\epsilon_{0})(p_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle(p_{1}\cdot\epsilon_{0})(p_{1}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon^{*}_{1})(p_{2}\cdot\epsilon^{*}_{3})\;,\qquad(p_{1}\cdot\epsilon_{0})(p_{1}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle(p_{1}\cdot\epsilon_{0})(p_{1}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon^{*}_{1})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad(p_{1}\cdot\epsilon_{0})(p_{1}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{1})(p_{3}\cdot\epsilon^{*}_{2})\;,$
$\displaystyle(p_{1}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad(p_{1}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{1})(p_{3}\cdot\epsilon^{*}_{2})\;,$
$\displaystyle(p_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})\;,\qquad(p_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon_{0})(p_{3}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle(p_{1}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})(p_{2}\cdot\epsilon^{*}_{3})\;,\qquad(p_{1}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle(p_{1}\cdot\epsilon^{*}_{3})(\epsilon^{*}_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon_{0})\;,\qquad(p_{1}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad(p_{1}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon_{0})(p_{3}\cdot\epsilon^{*}_{1})(p_{3}\cdot\epsilon^{*}_{2})\;,$
$\displaystyle(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{1})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad(p_{2}\cdot\epsilon^{*}_{1})(\epsilon_{0}\cdot\epsilon^{*}_{3})(p_{1}\cdot\epsilon^{*}_{2})\;,$
$\displaystyle(p_{2}\cdot\epsilon^{*}_{3})(\epsilon^{*}_{1}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon_{0})\;,\qquad(p_{3}\cdot\epsilon^{*}_{1})(\epsilon_{0}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})\;,\qquad(p_{3}\cdot\epsilon^{*}_{1})(\epsilon^{*}_{2}\cdot\epsilon^{*}_{3})(p_{1}\cdot\epsilon_{0})\;,$
$\displaystyle(p_{3}\cdot\epsilon^{*}_{2})(\epsilon_{0}\cdot\epsilon^{*}_{1})(p_{2}\cdot\epsilon^{*}_{3})\;,\qquad(p_{3}\cdot\epsilon^{*}_{2})(\epsilon^{*}_{1}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon_{0})\;.$
(37)

The four-dimensional projection Eq. (14) removes contractions of polarisation
vectors amongst themselves:

$\displaystyle(p_{1}\cdot\epsilon_{0})(p_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon^{*}_{1})\;,\qquad(p_{1}\cdot\epsilon_{0})(p_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle(p_{1}\cdot\epsilon_{0})(p_{1}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon^{*}_{1})(p_{2}\cdot\epsilon^{*}_{3})\;,\qquad(p_{1}\cdot\epsilon_{0})(p_{1}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle(p_{1}\cdot\epsilon_{0})(p_{1}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon^{*}_{1})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad(p_{1}\cdot\epsilon_{0})(p_{1}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{1})(p_{3}\cdot\epsilon^{*}_{2})\;,$
$\displaystyle(p_{1}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad(p_{1}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{1})(p_{3}\cdot\epsilon^{*}_{2})\;,$
$\displaystyle(p_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})\;,\qquad(p_{1}\cdot\epsilon^{*}_{2})(p_{1}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon_{0})(p_{3}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle(p_{1}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})(p_{2}\cdot\epsilon^{*}_{3})\;,\qquad(p_{1}\cdot\epsilon^{*}_{2})(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{1})\;,$
$\displaystyle(p_{1}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad(p_{1}\cdot\epsilon^{*}_{3})(p_{2}\cdot\epsilon_{0})(p_{3}\cdot\epsilon^{*}_{1})(p_{3}\cdot\epsilon^{*}_{2})\;,$
$\displaystyle(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{1})(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{2})\;,\qquad(p_{2}\cdot\epsilon_{0})(p_{2}\cdot\epsilon^{*}_{3})(p_{3}\cdot\epsilon^{*}_{1})(p_{3}\cdot\epsilon^{*}_{2})\;.$
(38)

In this case, however, the expression for the splitting operator is much
larger after projection than before. Hence, we do not provide it in electronic
form. The only advantage of the projection is the removal of terms
proportional to the six-dimensional integral of Section 3.2 for six possible
permutations of the external momenta. On the other hand, since we know (and
have checked) that these integrals disappear in four dimensions, one can
directly set them to zero in the expanded expression up to
$\order{\epsilon^{0}}$.

###  2.3 Requirements posed by the construction of a subtraction scheme

If one could obtain the exact $\epsilon$-dependence of the Feynman integrals
occurring in the expressions of the splitting operators derived with the
methods of the previous section, then this paragraph would not be required.
Unfortunately, as we will see in Section 3.2, there is one integral that
cannot, at least at present, be obtained in this generality. Hence, we must
understand what is actually required of the results for them to be useful in
the construction of a subtraction scheme. The conclusions of this discussion
depend on the perturbation-theory order at which the said subtraction scheme
is to be valid. Here, we restrict ourselves to N3LO.

To make the argument clear, we begin with the simplest possible subtraction
scheme, namely a scheme at NLO (more details in the classic Refs.
Catani:1996vz ; Frixione:1995ms ). Here, the relevant factor in the
dimensionally-regulated phase-space integration measure for a selected
massless parton takes the form:

$\iint_{0}^{1}\frac{\differential{\eta}}{\eta^{\epsilon}}\frac{\differential{\xi}}{\xi^{-1+2\epsilon}}\;,$
(39)

where $\eta\equiv\tfrac{1}{2}(1-\cos\theta)$ and $\theta$ is the angle between
the three-momentum of the chosen parton and that of another massless parton,
while $\xi$ is the normalised energy of the chosen parton. The squared tree-
level matrix element has singular asymptotics at vanishing $\eta$, if the two
partons are either a quark-anti-quark pair of the same flavor, or one of them
is a gluon:

$\innerproduct{M^{(0)}}{M^{(0)}}\equiv\frac{f(\eta)}{\eta}\quad\sim\quad\frac{f(0)}{\eta}\qquad(\eta\to
0)\;.$ (40)

We have suppressed the dependence on the remaining variables parameterising
the final state momenta that are irrelevant to the problem. There may also be
a singularity at vanishing $\xi$ if the energy is that of a gluon, but it is
sufficient to restrict to just the collinear singularity in order to
demonstrate the main issue. A subtraction scheme allows to evaluate the phase
space integral of the squared matrix element with the help of:

$\begin{split}\int_{0}^{1}\frac{\differential{\eta}}{\eta^{1+\epsilon}}f(\eta)&=\Big{[}\int_{0}^{1}\frac{\differential{\eta}}{\eta^{1+\epsilon}}\big{(}f(\eta)-f(0)\big{)}\Big{]}+\Big{[}f(0)\int_{0}^{1}\frac{\differential{\eta}}{\eta^{1+\epsilon}}\Big{]}\\\\[5.69046pt]
&=\Big{[}\int_{0}^{1}\frac{\differential{\eta}}{\eta^{1+\epsilon}}\big{(}f(\eta)-f(0)\big{)}\Big{]}+\Big{[}-\frac{1}{\epsilon}f(0)\Big{]}\;.\end{split}$
(41)

The integrand of the integral in the first square bracket contains a
subtraction term, $-f(0)/\eta^{1+\epsilon}$ that makes it integrable even
after $\epsilon$-expansion, while the integral in the second square bracket is
called the integrated subtraction term.

Suppose now, that we would like to evaluate double-real contributions to a
cross section at N3LO. These are cross-section contributions from processes
that have two additional massless partons in the final state with respect to
the Born-approximation process. Furthermore, the required matrix elements are
evaluated at one-loop order. For simplicity, let us assume that the process is
unpolarised. With these assumptions, the following expression:

$-\Big{(}\frac{2}{s_{123}}\Big{)}^{2}\Big{[}\expectationvalue{\bm{\hat{P}}^{(1)}_{f_{1}f_{2}f_{3}}}{M_{f\dots}^{(0)}}+2\real\matrixelement{M_{f\dots}^{(0)}}{\bm{\hat{P}}^{(0)}_{f_{1}f_{2}f_{3}}}{M_{f\dots}^{(1)}}\Big{]}\;,$
(42)

provides a subtraction term for triple-collinear singularities. Indeed,
according to Eqs. (2), (7) and (8), it correctly removes these singularities
from:

$2\real\innerproduct{M_{f_{1}f_{2}f_{3}\dots}^{(0)}}{M_{f_{1}f_{2}f_{3}\dots}^{(1)}}\;.$
(43)

By extension of Eq. (41), we also need an integrated subtraction term, which
contains at least one explicit $1/\epsilon$ pole. Since the integrated
subtraction term is evaluated in the collinear limit (in the example:
$f(\eta=0)$), there is no dependence on the transverse direction, and we can
use an averaged splitting function. This is, however, secondary. More
importantly, due to the presence of a $1/\epsilon$ pole generated by phase-
space integration, we need the splitting functions to at least
$\order{\epsilon}$, in order for the cross-section contribution to be correct
at $\order{\epsilon^{0}}$.

There is yet another issue to take into account when evaluating integrated
subtraction terms. Let us, for a moment, return to the example of one parton
splitting into two, but this time at the one-loop level. Whereas the tree-
level matrix-element squared is a rational function of scalar products of
momenta, and has the asymptotic (40), one-loop splitting functions for a
splitting into two partons scale as $s_{12}^{-\epsilon}$ when $s_{12}\to 0$
or, when expressed through $\eta$, as $1/\eta^{\epsilon}$. This must be taken
into account when evaluating the one-loop integrated subtraction term with Eq.
(41), otherwise the coefficient of the $1/\epsilon$ pole will be incorrect. In
our case of main interest, the subtraction term (42) scales as
$s_{123}^{-\epsilon}/s_{123}^{2}$ on purely dimensional grounds. However, the
triple-collinear splitting functions have further singularities themselves.
Indeed, let $\theta_{ij}$ be the angle between the momenta of partons $i$ and
$j$, and let $E_{k}$ be the energy of parton $k$, with $i,j,k\in\\{1,2,3\\}$.
Then, in the worst case scenario of a gluon splitting into three gluons, the
following limits are singular additionally to the triple-collinear singularity
already present: 1) $\theta_{ij}\to 0$ (iterated single-collinear limit); 2)
$E_{i}\to 0$ (single-soft limit); 3) $\theta_{ij}\to 0$, $E_{k}\to 0$; 4)
$E_{i},E_{j}\to 0$ (double-soft limit); 5) $E_{i},E_{j}\to 0$, $E_{i}/E_{j}\to
0$ (iterated single-soft limit); 6) $\theta_{ij}\to 0$, $E_{k},E_{l}\to 0$; 7)
$\theta_{ij}\to 0$, $E_{k},E_{l}\to 0$, $E_{k}/E_{l}\to 0$. Each of these
singular configurations requires subtraction and integrated subtraction terms
with the correct scaling in the relevant variables. Up to four poles are thus
generated: one for the triple-collinear limit, and one for each of the three
limits in case 7). Thus, the splitting functions must be known to
$\order{\epsilon^{4}}$. To be more precise, beyond $\order{\epsilon}$ we only
need the approximations to the splitting functions valid in the respective
limit.

Let us conclude this discussion by taking a first look at the results for the
Feynman integrals presented in full detail in Section 3. The results that are
available with exact dependence on $\epsilon$ consist of terms of the form:

$\displaystyle(\dots)^{-\epsilon}\,\frac{1}{\epsilon}\,{}_{2}F_{1}(\dots,x)\;,\qquad(\dots)^{-\epsilon}\,\frac{1}{\epsilon^{2}}\,{}_{2}F_{1}(\dots,x)\;,$
$\displaystyle(\dots)^{-\epsilon}\,(\dots)^{\epsilon}\,\frac{1}{\epsilon^{2}}\,{}_{2}F_{1}(\dots,x)\;,\qquad(\dots)^{-\epsilon}\,\frac{1}{\epsilon}\,F_{1}(\dots,x,y)\;,$
(44)

where ${}_{2}F_{1}$ is the hypergeometric function, see Appendix A, and
$F_{1}$ is the Appell function, see Appendix B. The arguments $x$ and $y$ of
these functions belong to the unit interval, and the functions are regular at
$x\to 0$ and $y\to 0$. The dots in the exponential functions represent
kinematic-dependent variables providing the scaling in the singular limits
discussed above. From the structure follows, that the results for the
splitting functions will be of sufficient quality for N3LO for a given limit
corresponding to $x\to 0$, if the hypergeometric functions are expanded up to
$\order{\epsilon^{5}}$ for a single pole $1/\epsilon$, or up to
$\order{\epsilon^{6}}$ for a double pole $1/\epsilon^{2}$ in the coefficient
of the ${}_{2}F_{1}$ function. These expansions are provided in Appendix A. In
case a singular limit corresponds to $x\to 1$, the scaling might be changed
due to the branch-point singularity of the ${}_{2}F_{1}$ function. In order to
capture the correct behaviour, we provide expansions about unit argument in
Appendix A. The situation is more involved for integrals containing the Appell
$F_{1}$ function. In any of the limits containing a collinear singularity, the
triple-collinear splitting functions must factorize into a product of
splitting functions for one-to-two splittings. These splitting functions may
be expressed through hypergeometric functions as worst, see Ref.
Kosower:1999rx . Hence, the Appell function should reduce to hypergeometric
functions in these cases. We have verified the correctness of this conclusion.
A similar simplification should take place for limits with a single-soft
singularity. The latter has an even simpler functional form, see Ref.
Catani:2000pi . In the remaining, double-soft-singular case, the Appell
$F_{1}$ functions present in the integrals reduce to hypergeometric functions
as well, see Eqs. (3.1), (59) and (83). In the end, the $F_{1}$ function with
arbitrary arguments is only needed to $\order{\epsilon}$. This expansion is
given in Eq. (81).

There still remains the highly-nontrivial six-dimensional integral of Section
3.2. We have verified that it only contributes to the pure triple-collinear
case, as well as to the triple-collinear/double-soft case. In the latter, one
could have expected a non-trivial integral with a square root of a Gram
determinant, just as we have obtained, in view of the general results for the
one-loop double-soft limit presented in Ref. Zhu:2020ftr 555The results of
Ref. Zhu:2020ftr  are restricted to expansions up to $\order{\epsilon^{0}}$.
Nevertheless, the discussion of the calculation of the integrals points to the
presence of a square root of a Gram determinant at higher orders in the
$\epsilon$-expansion.

###  2.4 A short list of checks

The complexity of the calculation warrants extensive checking. We have
successfully performed the following checks:

  1. 1.

comparison of the predicted singularity structure of the splitting operators,
Eq. (5), with that obtained from the direct calculation;

  2. 2.

comparison of the anti-symmetric part of the splitting function for $q\to
qq^{\prime}\bar{q}^{\prime}$ with the result given in Ref. Catani:2003vu ;

  3. 3.

comparison of the splitting functions for $q\to qq^{\prime}\bar{q}^{\prime}$
and $q\to qq\bar{q}$ expanded to $\order{\epsilon^{0}}$ obtained using the
methods of Section 2.2, with the result of an expansion in the triple-
collinear limit of one-loop matrix-elements squared for the processes $V\to
q\bar{q}q^{\prime}\bar{q}^{\prime}$ and $V\to q\bar{q}q\bar{q}$, with $V$ a
massive off-shell vector boson;

  4. 4.

numerical comparison of the triple-collinear limits of one-loop matrix-
elements squared at $\order{\epsilon^{0}}$ for six-parton processes with the
predicted asymptotics, Eq. (42), using the software library NJet Badger:2012pg
;

  5. 5.

comparison of the values of the master integrals obtained from analytic
formulae and from Mellin-Barnes representations up to the provided orders of
$\epsilon$-expansion.

##  3 Master integrals

Figure 7: Diagrammatic representation of the master integrals. Lines 1,2 and 3
correspond to outgoing momenta $p_{1}$, $p_{2}$ and $p_{3}$ respectively. The
thick incoming line is off shell with momentum $p_{123}=p_{1}+p_{2}+p_{3}$.
The red double-line corresponds to a linear (eikonal) propagator introduced
through the use of the lightcone gauge.

After Passarino-Veltman reduction and integration-by-parts reduction using the
software package Kira Maierhofer:2017gsa ; Klappert:2020nbg , the splitting
operators are expressed in terms of 34 master integrals. However, most of
these integrals are related by permutation of the external momenta. Taking
this into account, there remain only 9 master integrals depicted in Fig. 7.
All of them are generated from the following formula:

$\begin{split}&I^{(d)}_{a_{1}a_{2}a_{3}a_{4}a_{5}}\equiv\mu^{2\epsilon}\int\frac{\differential[d]{l}}{i\pi^{d/2}}\frac{1}{\big{(}l^{2}\big{)}^{a_{1}}\big{(}(l+p_{1})^{2}\big{)}^{a_{2}}\big{(}(l+p_{1}+p_{2})^{2}\big{)}^{a_{3}}\big{(}(l+p_{1}+p_{2}+p_{3})^{2}\big{)}^{a_{4}}\big{(}l\cdot
q\big{)}^{a_{5}}}\\\\[11.38092pt]
&\quad=(-1)^{a+a_{5}}\Gamma\Big{(}a-\frac{d}{2}\Big{)}\int_{\mathbb{R}_{+}^{4}}\Big{(}\prod_{i=1}^{4}\frac{\alpha_{i}^{a_{i}-1}\differential{\alpha_{i}}}{\Gamma(a_{i})}\Big{)}\frac{\delta(1-\sum_{i=1}^{4}\alpha_{i})(\alpha_{1}+\alpha_{2}+\alpha_{3}+\alpha_{4})^{a+a_{5}-d}}{\big{(}-s_{12}\;\alpha_{1}\alpha_{3}-s_{123}\;\alpha_{1}\alpha_{4}-s_{23}\;\alpha_{2}\alpha_{4}-i0^{+}\big{)}^{a-d/2}}\\\\[5.69046pt]
&\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\times\frac{1}{\big{(}(p_{1}\cdot
q)\;\alpha_{2}+(p_{12}\cdot q)\;\alpha_{3}+(p_{123}\cdot
q)\;\alpha_{4}-i0^{+}\big{)}^{a_{5}}}\;,\end{split}$ (45)

where $a\equiv a_{1}+a_{2}+a_{3}+a_{4}$, $p_{ij}\equiv p_{i}+p_{j}$,
$s_{ij}\equiv p_{ij}^{2}=2(p_{i}\cdot p_{j})$, $p_{123}\equiv
p_{1}+p_{2}+p_{3}$ and $s_{123}\equiv p_{123}^{2}=s_{12}+s_{13}+s_{23}$.
Expressions for missing propagators, i.e. at $a_{i}\to 0$ are obtained by the
simple replacement:

$\frac{\alpha_{i}^{a_{i}-1}}{\Gamma(a_{i})}\quad\xrightarrow[a_{i}\;\to\;0]{}\quad\delta(\alpha_{i})\;.$
(46)

The results presented in the next subsections are normalised with a universal
function of $\epsilon$:

$r_{\Gamma}\equiv\frac{\Gamma(1-\epsilon)^{2}\Gamma(1+\epsilon)}{\Gamma(1-2\epsilon)}\;.$
(47)

With $r_{\Gamma}$, the one-loop triangle integral in $d=4-2\epsilon$ is a pure
double pole in $\epsilon$ up to dependence on the momentum-transfer squared
dictated by dimensional arguments:

$I^{(4-2\epsilon)}_{11100}=r_{\Gamma}\;\Big{[}\Big{(}\frac{-s_{12}}{\mu^{2}}\Big{)}^{-\epsilon}\frac{1}{s_{12}}\Big{]}\;\frac{1}{\epsilon^{2}}\;.$
(48)

The master integrals depend on a number of kinematics-based variables. It is
useful to define the following dimensionless variables:

$\begin{gathered}x_{1}\equiv\frac{s_{23}}{s_{123}}\;,\quad
x_{2}\equiv\frac{s_{13}}{s_{123}}\;,\quad
x_{3}\equiv\frac{s_{12}}{s_{123}}\;,\qquad z_{i}\equiv\frac{p_{i}\cdot
q}{p_{123}\cdot q}\;,\\\\[5.69046pt] \quad x_{i}\in(0,1)\;,\qquad
z_{i}\in(0,1)\;,\qquad\sum_{i=1}^{3}x_{i}=\sum_{i=1}^{3}z_{i}=1\;.\end{gathered}$
(49)

There are further restrictions on the possible values of $x_{i}$ and $z_{i}$.
Indeed, consider for example the restframe of $p_{23}$ with $m_{23}^{2}\equiv
s_{23}$, and let $c_{\theta}$ be the cosine of the angle between the three-
vectors $\bm{p}_{1}$ and $\bm{q}$. Then:

$\begin{split}x_{1}&=\frac{s_{23}}{s_{123}}=\frac{m_{23}^{2}}{m_{23}^{2}+2p_{1}^{0}\;m_{23}}=\frac{1}{1+\frac{2p_{1}^{0}}{m_{23}}}\;,\\\\[11.38092pt]
1-z_{1}&=\frac{p_{23}\cdot q}{p_{123}\cdot
q}=\frac{m_{23}\;q^{0}}{p_{1}^{0}q^{0}\left(1-c_{\theta}\right)+m_{23}\;q^{0}}=\frac{1}{1+\frac{2p_{1}^{0}}{m_{23}}\left(\frac{1-c_{\theta}}{2}\right)}\geq\frac{1}{1+\frac{2p_{1}^{0}}{m_{23}}}=x_{1}\;.\end{split}$
(50)

Hence, in general:

$0<x_{i}+z_{i}\leq 1\;,\qquad\frac{z_{i}}{1-x_{i}}\leq
1\;,\qquad\frac{x_{i}}{1-z_{i}}\leq 1\;.$ (51)

###  3.1 Results with exact dependence on the spacetime dimension

The integrals listed in this subsection have been evaluated in $d=4-2\epsilon$
dimensions. Their dependence on the dimensional regularisation parameter
$\epsilon$ is exact, while the results are expressed in terms of
hypergeometric ${}_{2}F_{1}$ functions, see Appendix A, and Appell $F_{1}$
functions, see Appendix B. The arguments of these functions belong to the unit
interval in every case due to (49) and (51). We begin with ordinary Feynman
integrals:

$\displaystyle I^{(4-2\epsilon)}_{10010}$
$\displaystyle=r_{\Gamma}\Big{(}\frac{-s_{123}}{\mu^{2}}\Big{)}^{-\epsilon}\;\frac{1}{\epsilon(1-2\epsilon)}\;,$
(52) $\displaystyle I^{(4-2\epsilon)}_{10100}$
$\displaystyle=r_{\Gamma}\Big{(}\frac{-s_{12}}{\mu^{2}}\Big{)}^{-\epsilon}\;\frac{1}{\epsilon(1-2\epsilon)}\;,$
(53) $\displaystyle I^{(4-2\epsilon)}_{11110}$
$\displaystyle=r_{\Gamma}\Big{(}\frac{-s_{12}s_{23}}{\mu^{2}s_{123}}\Big{)}^{-\epsilon}\;\frac{2}{s_{12}s_{23}\;\epsilon^{2}}$
$\displaystyle\quad\times\Big{(}(x_{1}+x_{2})^{\epsilon}\;{}_{2}F_{1}\Big{(}-\epsilon,-\epsilon,1-\epsilon,1-\frac{x_{1}}{x_{1}+x_{2}}\Big{)}+(x_{2}+x_{3})^{\epsilon}\;{}_{2}F_{1}\Big{(}-\epsilon,-\epsilon,1-\epsilon,1-\frac{x_{3}}{x_{2}+x_{3}}\Big{)}$
$\displaystyle\qquad\quad-\big{(}(x_{1}+x_{2})(x_{2}+x_{3})\big{)}^{\epsilon}\;{}_{2}F_{1}\Big{(}-\epsilon,-\epsilon,1-\epsilon,1-\frac{x_{1}}{x_{1}+x_{2}}\frac{x_{3}}{x_{2}+x_{3}}\Big{)}\Big{)}\;.$
(54)

The integrals Eqs. (52) and (53) are two cases of the one-loop bubble integral
also known as the $B_{0}$ function. As such, they are textbook material. The
one-loop off-shell box integral Eq. (3.1) has been evaluated in Ref.
Bern:1993kr . The following integrals with a linear propagator can be found in
Ref. Sborlini:2014eib :

$\displaystyle I^{(4-2\epsilon)}_{01011}$
$\displaystyle=-r_{\Gamma}\Big{(}\frac{-s_{23}}{\mu^{2}}\Big{)}^{-\epsilon}\;\frac{1}{\big{(}p_{123}\cdot
q\big{)}\;\epsilon\;(1-2\epsilon)}\;_{2}F_{1}\big{(}1,1-\epsilon,2-2\epsilon,1-z_{1}\big{)}\;,$
(55) $\displaystyle I^{(4-2\epsilon)}_{11101}$
$\displaystyle=-r_{\Gamma}\Big{(}\frac{-s_{12}}{\mu^{2}}\Big{)}^{-\epsilon}\;\frac{2}{s_{12}\;\big{(}p_{12}\cdot
q\big{)}\;\epsilon^{2}}\;{}_{2}F_{1}\Big{(}1,1,1-\epsilon,1-\frac{z_{1}}{z_{1}+z_{2}}\Big{)}\;,$
(56) $\displaystyle I^{(4-2\epsilon)}_{10111}$
$\displaystyle=-r_{\Gamma}\Big{(}\frac{-s_{123}}{\mu^{2}}\Big{)}^{-\epsilon}\frac{2}{(s_{13}+s_{23})\;\big{(}p_{123}\cdot
q\big{)}\;\epsilon\;(1+\epsilon)}$
$\displaystyle\quad\times\Bigg{(}\,_{2}F_{1}\Big{(}1,1,2+\epsilon,1-\frac{z_{3}}{1-x_{3}}\Big{)}-\frac{x_{3}^{-\epsilon}}{1-z_{3}}\;_{2}F_{1}\Big{(}1,1,2+\epsilon,1-\frac{z_{3}}{1-x_{3}}\frac{x_{3}}{1-z_{3}}\Big{)}\Bigg{)}\;.$
(57)

These results may be verified by directly integrating the corresponding
Feynman-parameter integrals using the general expression Eq. (45). In the
present work, we have obtained two more exact results:

$\displaystyle I^{(4-2\epsilon)}_{11011}$
$\displaystyle=-r_{\Gamma}\Big{(}\frac{-s_{123}}{\mu^{2}}\Big{)}^{-\epsilon}\frac{1}{(s_{12}+s_{13})\;\big{(}p_{123}\cdot
q\big{)}\;\epsilon^{2}}$
$\displaystyle\quad\times\Bigg{(}2\,_{2}F_{1}\Big{(}1,1,1-\epsilon,1-\frac{z_{1}}{1-x_{1}}\Big{)}-x_{1}^{-\epsilon}F_{1}\Big{(}1,-\epsilon,1,1-2\epsilon,1-z_{1},1-\frac{z_{1}}{1-x_{1}}\Big{)}\Bigg{)}\;,$
(58) $\displaystyle I^{(4-2\epsilon)}_{01111}$
$\displaystyle=-r_{\Gamma}\Big{(}\frac{-s_{23}}{\mu^{2}}\Big{)}^{-\epsilon}\frac{1}{s_{23}\;\big{(}p_{123}\cdot
q\big{)}\;\epsilon^{2}}\;F_{1}\Big{(}1,-\epsilon,1,1-2\epsilon,1-z_{1},1-z_{1}-z_{2}\Big{)}\;.$
(59)

The values of these integrals in an $\epsilon$-expansion to
$\order{\epsilon^{0}}$ can be found in Ref. Sborlini:2014eib . These
expansions are, however, insufficient for applications at N3LO. The first
integral, Eq. (3.1), can be evaluated starting from the following Mellin-
Barnes representation, which can be derived with standard methods described
for example in Ref. Smirnov:2006ry :

$\begin{split}I^{(4-2\epsilon)}_{11011}&=-\Big{(}\frac{-s_{123}}{\mu^{2}}\Big{)}^{-\epsilon}\frac{1}{s_{123}\;\big{(}p_{123}\cdot
q\big{)}\;\Gamma(-2\epsilon)}\\\\[5.69046pt] &\quad\times\frac{1}{(2\pi
i)^{2}}\iint_{\mathcal{C}}\differential{z_{1}}\differential{z_{2}}\Gamma\big{(}1+\epsilon+z_{1}\big{)}\Gamma\big{(}1+z_{2}\big{)}\Gamma\big{(}1+z_{1}+z_{2}\big{)}\Gamma\big{(}-\epsilon-
z_{1}\big{)}\Gamma\big{(}-1-\epsilon-z_{2}\big{)}\\\\[5.69046pt]
&\qquad\qquad\qquad\qquad\qquad\times\Gamma\big{(}-z_{1}\big{)}\Gamma\big{(}-z_{2}\big{)}\;\Big{(}\frac{s_{23}}{s_{123}}\Big{)}^{z_{1}}\Big{(}\frac{p_{1}\cdot
q}{p_{123}\cdot q}\Big{)}^{z_{2}}\;.\end{split}$ (60)

Using the method proposed in Appendix C of Ref. Anastasiou:2013srw , one
obtains a one-fold integral, which straightforwardly yields Eq. (3.1). The
second integral, Eq. (59), has the Mellin-Barnes representation:

$\begin{split}I^{(4-2\epsilon)}_{01111}&=-\Big{(}\frac{-s_{23}}{\mu^{2}}\Big{)}^{-\epsilon}\frac{\Gamma(1+\epsilon)}{s_{23}\;\big{(}p_{123}\cdot
q\big{)}\;\Gamma(-2\epsilon)}\\\\[5.69046pt] &\quad\times\frac{1}{(2\pi
i)^{2}}\iint_{\mathcal{C}}\differential{z_{1}}\differential{z_{2}}\Gamma\big{(}-\epsilon+z_{1}\big{)}\Gamma\big{(}1+z_{2}\big{)}\Gamma\big{(}1+z_{1}+z_{2}\big{)}\Gamma\big{(}-1-\epsilon-
z_{1}-z_{2}\big{)}\\\\[5.69046pt]
&\qquad\qquad\qquad\qquad\qquad\times\Gamma\big{(}-z_{1}\big{)}\Gamma\big{(}-z_{2}\big{)}\;\Big{(}\frac{p_{1}\cdot
q}{p_{123}\cdot q}\Big{)}^{z_{1}}\Big{(}\frac{p_{12}\cdot q}{p_{123}\cdot
q}\Big{)}^{z_{2}}\;.\end{split}$ (61)

The result follows by comparing with the Mellin-Barnes representation of the
Appell $F_{1}$ function that appears last in Eq. (80). An analysis of the
soft/collinear limits of the splitting functions shows that the function
$F_{1}(1,-\epsilon,1,1-2\epsilon,1-x,1-y)$ with general arguments $x\neq y$
away from the endpoints 0 and 1, is only needed in an $\epsilon$-expansion to
$\order{\epsilon}$. The respective result can be found in Eq. (81). Indeed, in
any soft/collinear limit, the $F_{1}$ functions with general arguments cancel,
while they are expressible through ${}_{2}F_{1}$ functions in the remaining
cases. The appropriate expressions are given in Appendix B.

###  3.2 The light-cone-gauge box integral with one leg off shell

It turns out that the light-cone-gauge box integral $I^{(d)}_{11111}$ is
convergent in $d=6$ dimensions. This property mostly follows from the well-
known fact that there are no soft/collinear divergences in six dimensions, and
from ultraviolet power counting. In the case of ordinary Feynman integrals,
these arguments are sufficient to prove convergence. Due to the presence of
the linear propagator, it is necessary to consider possible rapidity
divergences. The latter make the six-dimensional integrals
$I^{(6)}_{01111},I^{(6)}_{10111},I^{(6)}_{11011}$ and $I^{(6)}_{11101}$
diverge, although the aforementioned “soft/collinear/ultraviolet argument”
applies in their case as well. Since these rapidity divergences are due to
large longitudinal, $l_{+}$, and transverse, $\bm{l}_{T}$, components of the
loop momentum defined with respect to $q$, with $|l_{+}|\sim\bm{l}_{T}^{2}$,
they only disappear for integrals with at least four ordinary propagators,
e.g. for $I^{(6)}_{11111}$. With the help of the integration-by-parts
reduction relations and using the methods of Refs. Tarasov:1996br ;
Tarasov:1997kx , we obtain the following dimensional-shift relation:

$\displaystyle 2s_{123}\;{I}^{{(d)}}_{11111}$
$\displaystyle=\frac{\big{(}x_{1}z_{1}-x_{2}z_{2}+x_{3}z_{3}\big{)}^{2}-4x_{1}x_{3}z_{1}z_{3}}{x_{1}x_{3}z_{1}\big{(}1-x_{3}-z_{3}\big{)}}\;\big{(}d-4\big{)}\;{I}^{{(d+2)}}_{11111}$
$\displaystyle\quad+\frac{\big{(}x_{1}z_{1}-x_{2}z_{2}+x_{3}z_{3}\big{)}\big{(}z_{1}+z_{2}\big{)}-2x_{3}z_{1}z_{3}}{x_{3}z_{1}\big{(}1-x_{3}-z_{3}\big{)}}\;{I}^{{(d)}}_{01111}-\frac{x_{1}z_{1}-x_{2}z_{2}+x_{3}z_{3}}{x_{1}x_{3}z_{1}}\;{I}^{{(d)}}_{10111}$
$\displaystyle\quad+\frac{\big{(}x_{1}z_{1}-x_{2}z_{2}+x_{3}z_{3}\big{)}-2x_{1}x_{3}}{x_{1}x_{3}\big{(}1-x_{3}-z_{3}\big{)}}\;{I}^{{(d)}}_{11011}-\frac{\big{(}x_{1}z_{1}-x_{2}z_{2}+x_{3}z_{3}\big{)}-2x_{1}(z_{1}+z_{2})}{x_{1}\big{(}1-x_{3}-z_{3}\big{)}}\;{I}^{{(d)}}_{11101}$
$\displaystyle\quad+\frac{\big{(}x_{1}z_{1}-x_{2}z_{2}+x_{3}z_{3}\big{)}-2z_{1}(x_{1}+x_{2})}{z_{1}\big{(}1-x_{3}-z_{3}\big{)}}\;\Big{(}\frac{s_{123}}{p_{123}\cdot
q}\Big{)}\;{I}^{{(d)}}_{11110}\;.$ (62)

The presence of the factor $(d-4)$ on the right-hand-side coefficient of
$I^{(d+2)}_{11111}$ implies that the integral $I^{(4-2\epsilon)}_{11111}$
expanded up to $\order{\epsilon^{0}}$ is entirely given by the integrals of
the previous subsection. In fact, this expansion is given in Ref.
Sborlini:2014eib . On the other hand, since the results of previous studies of
one-loop triple-collinear splittings have been restricted to expansions up to
$\order{\epsilon^{0}}$, this integral has never been really needed in view of
Eq. (3.2). Let us note that the coefficient of $I^{(d+2)}_{11111}$ is
proportional to the Gram determinant of the momenta $k_{i}=p_{i}$, $i=1,2,3$,
$k_{4}=q$. Indeed:

$\begin{split}\Delta_{4}\equiv\det(k_{i}\cdot
k_{j})&=\frac{1}{4}s^{2}_{123}\,(p_{123}\cdot
q)^{2}\big{(}\big{(}x_{1}z_{1}-x_{2}z_{2}+x_{3}z_{3}\big{)}^{2}-4x_{1}x_{3}z_{1}z_{3}\big{)}\\\\[11.38092pt]
&\equiv\Big{(}\frac{1}{2}s_{123}\,(p_{123}\cdot
q)\;d_{4}\Big{)}^{2}\;.\end{split}$ (63)

An analysis of the soft/collinear limits of the splitting functions shows that
only the value of $I^{(6)}_{11111}$ is required for non-singular
configurations, while the double-soft limit $p_{1},p_{2}\to 0$, with
$z_{1}/(z_{1}+z_{2})\in(0,1)$, i.e. excluding the strongly ordered limits
$z_{1}\ll z_{2}$ and $z_{2}\ll z_{1}$ further requires
$I^{(6-2\epsilon)}_{11111}$ up to $\order{\epsilon}$. In order to evaluate
$I^{(6-2\epsilon)}_{11111}$, we turn to the Feynman-parameter representation
Eq. (45):

$\begin{split}I^{(6-2\epsilon)}_{11111}&=\frac{\Big{(}\frac{-s_{123}}{\mu^{2}}\Big{)}^{-\epsilon}\Gamma(1+\epsilon)}{s_{123}\,\big{(}p_{123}\cdot
q\big{)}}\int_{\mathbb{R}_{+}^{4}}\prod_{i=1}^{4}\differential{\alpha_{i}}\frac{\delta(1-\alpha_{1})\big{(}\alpha_{1}+\alpha_{2}+\alpha_{3}+\alpha_{4}\big{)}^{2\epsilon-1}}{\big{(}\alpha_{2}\,z_{1}+\alpha_{3}(z_{1}+z_{2})+\alpha_{4}\big{)}\big{(}\alpha_{1}(\alpha_{3}\,x_{3}+\alpha_{4})+\alpha_{2}\alpha_{4}\,x_{1}\big{)}^{1+\epsilon}}\\\\[11.38092pt]
&=\frac{\Big{(}\frac{-s_{123}\,(p_{1}\cdot q)}{\mu^{2}\,(p_{123}\cdot
q)}\Big{)}^{-\epsilon}\Gamma(1+\epsilon)}{s_{123}\,\big{(}p_{12}\cdot
q\big{)}}\iiint_{0}^{\infty}\differential{\alpha_{2}}\differential{\alpha_{3}}\differential{\alpha_{4}}\frac{\big{(}1+\alpha_{2}+\alpha_{3}\,y_{1}+\alpha_{4}\,z_{1}\big{)}{}^{2\epsilon-1}}{\big{(}\alpha_{2}+\alpha_{3}+\alpha_{4}\big{)}\big{(}\alpha_{3}\,u_{3}+\alpha_{4}(1+\alpha_{2}\,x_{1})\big{)}^{1+\epsilon}}\;.\end{split}$
(64)

In the first line, we have used the Cheng-Wu theorem to change the
$\delta$-function to only restrict $\alpha_{1}$. In the second line, we have
rescaled the integration variables as follows:

$\alpha_{3}\;\to\;\alpha_{3}\,y_{1}\;,\qquad\alpha_{4}\;\to\;\alpha_{4}\,z_{1}\;,$
(65)

and defined:

$y_{1}\equiv\frac{z_{1}}{z_{1}+z_{2}}\in(0,1)\;,\qquad
u_{3}\equiv\frac{x_{3}}{1-z_{3}}\in(0,1)\;.$ (66)

The integral thus depends on the variables $x_{1}$, $y_{1}$, $z_{1}$ and
$u_{3}$. The purpose of the rescaling Eq. (65) is to yield an integrable
integrand in the double-soft limit. In this limit, out of the four variables
$x_{1},y_{1},z_{1},u_{3}$ only $z_{1}$ vanishes. Setting $p_{1}=p_{2}=0$ in
the first line of Eq. (64) yields an integral in $\alpha_{4}$ that is not
integrable at the lower limit $\alpha_{4}=0$. The integral on the second line
does not suffer from this pathological behaviour. At the same time, we factor-
out the non-trivial scaling in the limit in the pre-factor.

The integrations may now be performed with the software package PolyLogTools
Duhr:2019tlz  in the order $\alpha_{3}$, $\alpha_{2}$, $\alpha_{4}$, using a
fibration basis for generalized polylogarithms corresponding to the ordering
of the variables ($\alpha_{3}$, $\alpha_{2}$, $\alpha_{4}$, $z_{1}$, $x_{1}$,
$y_{1}$, $u_{3}$, $r_{1}$, $r_{2}$). The last two variables are roots of a
quadratic polynomial in $\alpha_{4}$ that must be factorized in order to
obtain an alphabet linear in the integration variables:

$r_{1,2}\equiv\frac{x_{1}z_{1}-x_{2}z_{2}+x_{3}z_{3}\mp
d_{4}}{2x_{1}z_{1}z_{3}}\;.$ (67)

Here, $d_{4}$ is related to the Gram determinant as defined in Eq. (63). The
final result is:

$\displaystyle 2\sqrt{\Delta_{4}}\;\Bigg{(}\Big{(}\frac{-s_{123}\,(p_{1}\cdot
q)}{\mu^{2}\,(p_{123}\cdot
q)}\Big{)}^{-\epsilon}\Gamma(1+\epsilon)\Bigg{)}^{-1}\;I^{(6-2\epsilon)}_{11111}=$
$\displaystyle\quad
G\left(0,-r_{1}\right)G\left(0,u_{3}\right)G\left(\frac{u_{3}}{r_{1}},x_{1}\right)+G\left(0,-r_{1}\right)G\left(0,u_{3}\right)G\left(\frac{u_{3}}{r_{1}},y_{1}\right)$
$\displaystyle-G\left(0,-r_{2}\right)G\left(0,u_{3}\right)G\left(\frac{u_{3}}{r_{2}},x_{1}\right)-G\left(0,-r_{2}\right)G\left(0,u_{3}\right)G\left(\frac{u_{3}}{r_{2}},y_{1}\right)$
$\displaystyle+G\left(0,-r_{2}\right)G\left(0,x_{1}\right)G\left(\frac{1-x_{1}}{r_{2}x_{1}},z_{1}\right)-G\left(0,-r_{1}\right)G\left(0,x_{1}\right)G\left(\frac{1-x_{1}}{r_{1}x_{1}},z_{1}\right)$
$\displaystyle+G\left(0,-r_{1}\right)G\left(0,y_{1}\right)G\left(y_{1}-\frac{1}{r_{1}},z_{1}\right)-G\left(0,-r_{2}\right)G\left(0,y_{1}\right)G\left(y_{1}-\frac{1}{r_{2}},z_{1}\right)$
$\displaystyle+G\left(0,-r_{1}\right)G\left(0,y_{1}\right)G\left(\frac{u_{3}}{r_{1}y_{1}},x_{1}\right)-G\left(0,-r_{1}\right)G\left(0,u_{3}\right)G\left(\frac{u_{3}}{r_{1}y_{1}},x_{1}\right)$
$\displaystyle+G\left(0,-r_{2}\right)G\left(0,u_{3}\right)G\left(\frac{u_{3}}{r_{2}y_{1}},x_{1}\right)-G\left(0,u_{3}\right)G\left(0,y_{1}\right)G\left(\frac{u_{3}}{r_{1}y_{1}},x_{1}\right)$
$\displaystyle+G\left(0,u_{3}\right)G\left(0,y_{1}\right)G\left(\frac{u_{3}}{r_{2}y_{1}},x_{1}\right)-G\left(0,-r_{2}\right)G\left(0,y_{1}\right)G\left(\frac{u_{3}}{r_{2}y_{1}},x_{1}\right)$
$\displaystyle+G\left(0,-r_{1}\right)G\left(0,u_{3}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{1}},z_{1}\right)-G\left(0,-r_{1}\right)G\left(0,y_{1}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{1}},z_{1}\right)$
$\displaystyle+G\left(0,u_{3}\right)G\left(0,y_{1}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{1}},z_{1}\right)-G\left(0,-r_{2}\right)G\left(0,u_{3}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{2}},z_{1}\right)$
$\displaystyle+G\left(0,-r_{2}\right)G\left(0,y_{1}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{2}},z_{1}\right)-G\left(0,u_{3}\right)G\left(0,y_{1}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{2}},z_{1}\right)$
$\displaystyle-G\left(0,0,-r_{1}\right)G\left(0,u_{3}\right)+G\left(0,0,-r_{1}\right)G\left(\frac{1}{r_{1}},x_{1}\right)-G\left(0,0,-r_{1}\right)G\left(1-\frac{1}{r_{1}},z_{1}\right)$
$\displaystyle-G\left(0,0,-r_{1}\right)G\left(\frac{u_{3}}{r_{1}},x_{1}\right)-G\left(0,0,-r_{1}\right)G\left(\frac{u_{3}}{r_{1}},y_{1}\right)+G\left(0,0,-r_{1}\right)G\left(\frac{1}{r_{1}},y_{1}\right)$
$\displaystyle+G\left(0,0,-r_{1}\right)G\left(\frac{u_{3}}{r_{1}y_{1}},x_{1}\right)-G\left(0,0,-r_{1}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{1}},z_{1}\right)+G\left(0,0,-r_{1}\right)G\left(y_{1}-\frac{1}{r_{1}},z_{1}\right)$
$\displaystyle+G\left(0,0,-r_{2}\right)G\left(0,u_{3}\right)-G\left(0,0,-r_{2}\right)G\left(\frac{1}{r_{2}},x_{1}\right)+G\left(0,0,-r_{2}\right)G\left(1-\frac{1}{r_{2}},z_{1}\right)$
$\displaystyle+G\left(0,0,-r_{2}\right)G\left(\frac{u_{3}}{r_{2}},x_{1}\right)+G\left(0,0,-r_{2}\right)G\left(\frac{u_{3}}{r_{2}},y_{1}\right)-G\left(0,0,-r_{2}\right)G\left(\frac{1}{r_{2}},y_{1}\right)$
$\displaystyle-G\left(0,0,-r_{2}\right)G\left(\frac{u_{3}}{r_{2}y_{1}},x_{1}\right)+G\left(0,0,-r_{2}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{2}},z_{1}\right)-G\left(0,0,-r_{2}\right)G\left(y_{1}-\frac{1}{r_{2}},z_{1}\right)$
$\displaystyle-G\left(0,0,u_{3}\right)G\left(\frac{u_{3}}{r_{1}},x_{1}\right)+G\left(0,-r_{1}\right)G\left(0,0,u_{3}\right)-G\left(0,-r_{2}\right)G\left(0,0,u_{3}\right)$
$\displaystyle+G\left(0,0,u_{3}\right)G\left(\frac{u_{3}}{r_{2}},x_{1}\right)-G\left(0,0,u_{3}\right)G\left(\frac{u_{3}}{r_{1}},y_{1}\right)+G\left(0,0,u_{3}\right)G\left(\frac{u_{3}}{r_{2}},y_{1}\right)$
$\displaystyle+G\left(0,0,u_{3}\right)G\left(\frac{u_{3}}{r_{1}y_{1}},x_{1}\right)-G\left(0,0,u_{3}\right)G\left(\frac{u_{3}}{r_{2}y_{1}},x_{1}\right)-G\left(0,0,u_{3}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{1}},z_{1}\right)$
$\displaystyle+G\left(0,0,u_{3}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{2}},z_{1}\right)-G\left(0,0,x_{1}\right)G\left(\frac{1-x_{1}}{r_{1}x_{1}},z_{1}\right)+G\left(0,0,x_{1}\right)G\left(\frac{1-x_{1}}{r_{2}x_{1}},z_{1}\right)$
$\displaystyle+G\left(0,0,y_{1}\right)G\left(\frac{u_{3}}{r_{1}y_{1}},x_{1}\right)+G\left(0,0,y_{1}\right)G\left(y_{1}-\frac{1}{r_{1}},z_{1}\right)-G\left(0,0,y_{1}\right)G\left(y_{1}-\frac{1}{r_{2}},z_{1}\right)$
$\displaystyle-G\left(0,0,y_{1}\right)G\left(\frac{u_{3}}{r_{2}y_{1}},x_{1}\right)-G\left(0,0,y_{1}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{1}},z_{1}\right)+G\left(0,0,y_{1}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{2}},z_{1}\right)$
$\displaystyle-G\left(0,-r_{1}\right)G\left(1,0,u_{3}\right)+G\left(0,-r_{2}\right)G\left(1,0,u_{3}\right)-G\left(0,-r_{1}\right)G\left(1,0,x_{1}\right)$
$\displaystyle+G\left(1,0,x_{1}\right)G\left(\frac{1-x_{1}}{r_{1}x_{1}},z_{1}\right)-G\left(1,0,x_{1}\right)G\left(\frac{1-x_{1}}{r_{2}x_{1}},z_{1}\right)+G\left(0,-r_{2}\right)G\left(1,0,x_{1}\right)$
$\displaystyle-G\left(0,-r_{1}\right)G\left(1,0,y_{1}\right)+G\left(0,-r_{2}\right)G\left(1,0,y_{1}\right)+G\left(0,-r_{1}\right)G\left(1-\frac{1}{r_{1}},-\frac{1}{r_{1}},z_{1}\right)$
$\displaystyle+G\left(0,-r_{1}\right)G\left(\frac{1}{r_{1}},0,x_{1}\right)+G\left(0,-r_{1}\right)G\left(\frac{1}{r_{1}},0,y_{1}\right)-G\left(0,-r_{2}\right)G\left(1-\frac{1}{r_{2}},-\frac{1}{r_{2}},z_{1}\right)$
$\displaystyle-G\left(0,-r_{1}\right)G\left(\frac{u_{3}}{r_{1}},0,x_{1}\right)-G\left(0,-r_{2}\right)G\left(\frac{1}{r_{2}},0,x_{1}\right)-G\left(0,-r_{2}\right)G\left(\frac{1}{r_{2}},0,y_{1}\right)$
$\displaystyle+G\left(0,u_{3}\right)G\left(\frac{u_{3}}{r_{1}},0,x_{1}\right)-G\left(0,-r_{1}\right)G\left(\frac{u_{3}}{r_{1}},0,y_{1}\right)+G\left(0,u_{3}\right)G\left(\frac{u_{3}}{r_{1}},0,y_{1}\right)$
$\displaystyle+G\left(0,-r_{2}\right)G\left(\frac{u_{3}}{r_{2}},0,x_{1}\right)-G\left(0,u_{3}\right)G\left(\frac{u_{3}}{r_{2}},0,x_{1}\right)+G\left(0,-r_{2}\right)G\left(\frac{u_{3}}{r_{2}},0,y_{1}\right)$
$\displaystyle-G\left(0,u_{3}\right)G\left(\frac{u_{3}}{r_{2}},0,y_{1}\right)-G\left(0,x_{1}\right)G\left(\frac{1-x_{1}}{r_{1}x_{1}},0,z_{1}\right)-G\left(0,-r_{1}\right)G\left(\frac{1-x_{1}}{r_{1}x_{1}},-\frac{1}{r_{1}},z_{1}\right)$
$\displaystyle+G\left(0,x_{1}\right)G\left(\frac{1-x_{1}}{r_{2}x_{1}},0,z_{1}\right)+G\left(0,-r_{2}\right)G\left(\frac{1-x_{1}}{r_{2}x_{1}},-\frac{1}{r_{2}},z_{1}\right)-G\left(0,-r_{1}\right)G\left(y_{1}-\frac{1}{r_{1}},-\frac{1}{r_{1}},z_{1}\right)$
$\displaystyle+G\left(0,y_{1}\right)G\left(y_{1}-\frac{1}{r_{1}},y_{1},z_{1}\right)+G\left(0,-r_{2}\right)G\left(y_{1}-\frac{1}{r_{2}},-\frac{1}{r_{2}},z_{1}\right)-G\left(0,y_{1}\right)G\left(y_{1}-\frac{1}{r_{2}},y_{1},z_{1}\right)$
$\displaystyle-G\left(0,u_{3}\right)G\left(\frac{u_{3}}{r_{1}y_{1}},0,x_{1}\right)+G\left(0,y_{1}\right)G\left(\frac{u_{3}}{r_{1}y_{1}},0,x_{1}\right)+G\left(0,-r_{1}\right)G\left(\frac{u_{3}}{r_{1}y_{1}},0,x_{1}\right)$
$\displaystyle+G\left(0,u_{3}\right)G\left(\frac{u_{3}}{r_{2}y_{1}},0,x_{1}\right)-G\left(0,y_{1}\right)G\left(\frac{u_{3}}{r_{2}y_{1}},0,x_{1}\right)-G\left(0,-r_{2}\right)G\left(\frac{u_{3}}{r_{2}y_{1}},0,x_{1}\right)$
$\displaystyle+G\left(0,u_{3}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{1}},\frac{y_{1}}{u_{3}},z_{1}\right)-G\left(0,y_{1}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{1}},\frac{y_{1}}{u_{3}},z_{1}\right)+G\left(0,-r_{1}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{1}},-\frac{1}{r_{1}},z_{1}\right)$
$\displaystyle-G\left(0,u_{3}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{2}},\frac{y_{1}}{u_{3}},z_{1}\right)+G\left(0,y_{1}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{2}},\frac{y_{1}}{u_{3}},z_{1}\right)-G\left(0,-r_{2}\right)G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{2}},-\frac{1}{r_{2}},z_{1}\right)$
$\displaystyle-G\left(-1,0,0,-r_{1}\right)+G\left(-1,0,0,-r_{2}\right)+G\left(0,0,0,-r_{1}\right)-G\left(0,0,0,-r_{2}\right)+G\left(1-\frac{1}{r_{1}},1,0,z_{1}\right)$
$\displaystyle+G\left(1-\frac{1}{r_{1}},-\frac{1}{r_{1}},0,z_{1}\right)+G\left(\frac{1}{r_{1}},0,0,x_{1}\right)+G\left(\frac{1}{r_{1}},0,0,y_{1}\right)-G\left(1-\frac{1}{r_{2}},1,0,z_{1}\right)$
$\displaystyle-G\left(\frac{1}{r_{2}},0,0,x_{1}\right)-G\left(\frac{1}{r_{2}},0,0,y_{1}\right)-G\left(1-\frac{1}{r_{2}},-\frac{1}{r_{2}},0,z_{1}\right)-G\left(\frac{u_{3}}{r_{1}},0,0,x_{1}\right)+G\left(\frac{u_{3}}{r_{2}},0,0,x_{1}\right)$
$\displaystyle-G\left(\frac{u_{3}}{r_{1}},0,0,y_{1}\right)+G\left(\frac{u_{3}}{r_{2}},0,0,y_{1}\right)-G\left(\frac{1-x_{1}}{r_{1}x_{1}},-\frac{1}{r_{1}},0,z_{1}\right)+G\left(\frac{1-x_{1}}{r_{2}x_{1}},-\frac{1}{r_{2}},0,z_{1}\right)$
$\displaystyle-G\left(y_{1}-\frac{1}{r_{1}},-\frac{1}{r_{1}},0,z_{1}\right)-G\left(y_{1}-\frac{1}{r_{1}},y_{1},0,z_{1}\right)+G\left(y_{1}-\frac{1}{r_{2}},-\frac{1}{r_{2}},0,z_{1}\right)+G\left(\frac{u_{3}}{r_{1}y_{1}},0,0,x_{1}\right)$
$\displaystyle-G\left(\frac{u_{3}}{r_{2}y_{1}},0,0,x_{1}\right)+G\left(y_{1}-\frac{1}{r_{2}},y_{1},0,z_{1}\right)+G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{1}},-\frac{1}{r_{1}},0,z_{1}\right)+G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{1}},\frac{y_{1}}{u_{3}},0,z_{1}\right)$
$\displaystyle-G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{2}},-\frac{1}{r_{2}},0,z_{1}\right)-G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{2}},\frac{y_{1}}{u_{3}},0,z_{1}\right)+\zeta_{2}\Bigg{(}-3G\left(\frac{u_{3}}{r_{1}},x_{1}\right)+3G\left(\frac{u_{3}}{r_{2}},x_{1}\right)$
$\displaystyle-3G\left(\frac{u_{3}}{r_{1}},y_{1}\right)+3G\left(\frac{u_{3}}{r_{2}},y_{1}\right)-G\left(\frac{1-x_{1}}{r_{1}x_{1}},z_{1}\right)+G\left(\frac{1-x_{1}}{r_{2}x_{1}},z_{1}\right)+3G\left(\frac{u_{3}}{r_{1}y_{1}},x_{1}\right)$
$\displaystyle-3G\left(\frac{u_{3}}{r_{2}y_{1}},x_{1}\right)-3G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{1}},z_{1}\right)+3G\left(\frac{y_{1}}{u_{3}}-\frac{1}{r_{2}},z_{1}\right)+3G\left(\frac{1}{r_{1}},x_{1}\right)-3G\left(\frac{1}{r_{2}},x_{1}\right)$
$\displaystyle+3G\left(y_{1}-\frac{1}{r_{1}},z_{1}\right)-3G\left(y_{1}-\frac{1}{r_{2}},z_{1}\right)+3G\left(\frac{1}{r_{1}},y_{1}\right)-3G\left(\frac{1}{r_{2}},y_{1}\right)-3G\left(1-\frac{1}{r_{1}},z_{1}\right)$
$\displaystyle+3G\left(1-\frac{1}{r_{2}},z_{1}\right)-3G\left(-1,-r_{1}\right)+3G\left(-1,-r_{2}\right)+4G\left(0,-r_{1}\right)-4G\left(0,-r_{2}\right)\Bigg{)}+\order{\epsilon}\;,$
(68)

where the multiple polylogarithms are defined recursively as follows:

$G(a_{1},\dots,a_{n},z)\equiv\int_{0}^{z}\frac{\differential{t}}{t-a_{1}}\,G(a_{2},\dots,a_{n},t)\;,\qquad
G(\underbrace{0,\dots,0}_{n},z)\equiv\frac{1}{n!}\ln^{n}(z)\;,$ (69)

and can be numerically evaluated with PolyLogTools through its interface to
the software package GiNaC Bauer:2000cp ; Vollinga:2004sn . The right-hand-
side of Eq. (68) is regular in $z_{1}$ by construction. The double-soft limit
may be easily obtained with PolyLogTools using the function ExpandPolyLogs to
expand in $z_{1}$ to $\order{z_{1}^{0}}$. The result may then be simplified
using the shuffle algebra of the multiple polylogarithms with the function
ShuffleG. It turns out that the occurring multiple polylogarithms may be
rewritten in terms of classical polylogarithms up to weight three only by
simply transforming the expression with the function GToLi. Of course, the
result Eq. (68) may also be rewritten in terms of classical polylogarithms up
to weight three, since arbitrary multiple polylogarithms up to this weight
have this property. Apart from the double-soft limit, however, the result is
not simpler, but involves even more polylogarithms than there are multiple
polylogarithms in Eq. (68). Finally, we relegate the lengthy expression for
the $\order{\epsilon}$ term of the expansion of the left-hand-side of Eq. (68)
in the double-soft limit to an ancillary file attached to this publication,
see Appendix C.

##  4 Outlook

The present publication completes the study of one-loop triple-collinear
splitting operators and the related splitting functions. The results are
sufficient for the construction of an N3LO subtraction scheme. In fact, the
only missing input for complete generality is the exact $\epsilon$-dependence
of the six-dimensional light-cone-gauge box integral discussed in Section 3.2.
Should this result be derived in the future, it can be readily substituted in
the provided expressions.

For applications at N3LO, there still remains the need to provide the double-
soft asymptotics of one-loop matrix-elements squared to sufficient order in
the $\epsilon$-expansion. Of course, even with the knowledge of all the
limits, the construction of an N3LO subtraction scheme remains a daunting
task.

###### Acknowledgements.

The work of M.C. was supported by the Deutsche Forschungsgemeinschaft (DFG)
under grant 396021762 - TRR 257: Particle Physics Phenomenology after the
Higgs Discovery. The work of S.S. was supported by the Polish National Science
Centre grant no. 2017/27/B/ST2/02004.

##  Appendix A Definition, representations and expansions of hypergeometric
${}_{2}F_{1}$ functions

The hypergeometric functions occurring in the expressions of the splitting
operators have been chosen to have arguments $0<x<1$, which ensures that their
defining series converge:

$\begin{split}{}_{2}F_{1}(a,b,c,x)&\equiv\sum_{n=0}^{\infty}\frac{\Gamma(a+n)\Gamma(b+n)\Gamma(c)}{\Gamma(a)\Gamma(b)\Gamma(c+n)}\frac{x^{n}}{n!}\\\\[5.69046pt]
&=\frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)}\int_{0}^{1}\differential{t}t^{b-1}(1-t)^{c-b-1}(1-tx)^{-a}\\\\[5.69046pt]
&=\frac{1}{2\pi
i}\int_{\mathcal{C}}\differential{z}\frac{\Gamma(a+z)\Gamma(b+z)\Gamma(c)\Gamma(-z)}{\Gamma(a)\Gamma(b)\Gamma(c+z)}(-x)^{z}\\\\[5.69046pt]
&=\frac{1}{2\pi
i}\int_{\mathcal{C}}\differential{z}\frac{\Gamma(a+z)\Gamma(b+z)\Gamma(c)\Gamma(c-a-
b-z)\Gamma(-z)}{\Gamma(a)\Gamma(b)\Gamma(c-a)\Gamma(c-b)}(1-x)^{z}\;.\end{split}$
(70)

For completeness, we have also recalled the standard ordinary integral
representation of the ${}_{2}F_{1}$ function and two of its Mellin-Barnes
representations. The integration contours, $\mathcal{C}$, in the latter case
should be chosen along the imaginary axis and separate the poles of the
$\Gamma$-functions $\Gamma(\dots+z)$ and $\Gamma(\dots-z)$. If this cannot be
achieved for the desired value of the dimensional regularisation parameter
$\epsilon$, then the correct representation may be obtained with the help of
an analytic continuation, for example using the software package MB
Czakon:2005rk . However, this is not needed for any of the functions listed
below. The behaviour of the functions at the endpoint $x=1$ requires an
asymptotic expansion which may be derived with the help of a direct
consequence of the integral representation on the second line of Eq. (70):

$\begin{split}{}_{2}F_{1}(a,b,c,x)&=\frac{\Gamma(c-a-b)\Gamma(c)}{\Gamma(c-a)\Gamma(c-b)}\,_{2}F_{1}(a,b,a+b-c+1,1-x)\\\\[5.69046pt]
&\quad+\frac{\Gamma(a+b-c)\Gamma(c)}{\Gamma(a)\Gamma(b)}(1-x)^{c-a-b}\,_{2}F_{1}(c-a,c-b,c-a-b+1,1-x)\;.\end{split}$
(71)

The $\epsilon$-expansions of the hypergeometric functions required by the
present work may be obtained with the help of the software package HypExp
Maitre:2005uu ; Maitre:2007kp ; Huber:2005yg . Below, we reproduce them
retaining terms up to transcendental weight six as necessary for applications
at N3LO. The results are given in terms of Riemann’s $\zeta_{n}$, classical
polylogarithms $\text{Li}_{n}(x)$ and Nielsen polylogarithms
$\text{S}_{n,p}(x)$:

$\text{S}_{n,p}(x)\equiv\frac{(-1)^{n+p-1}}{(n-1)!p!}\int_{0}^{1}\differential{t}\frac{\ln^{n-1}(t)\ln^{p}(1-tx)}{t}\;,\qquad\text{Li}_{n}(x)=\text{S}_{n-1,1}(x)\;,\qquad\zeta_{n}=\text{Li}_{n}(1)\;.$
(72)

We also provide expansions at $x=1$ and, in one case where it is needed, at
$x=0$. The following functions are present in the expressions of master
integrals:

${}_{2}F_{1}$ $\displaystyle(-\epsilon,-\epsilon,1-\epsilon,x)=$
$\displaystyle
1+\epsilon^{2}\text{Li}_{2}(x)+\epsilon^{3}\Big{(}-\zeta_{3}+\text{Li}_{3}(1-x)+\text{Li}_{3}(x)-\text{Li}_{2}(1-x)\ln(1-x)-\frac{1}{2}\ln(x)\ln^{2}(1-x)\Big{)}$
$\displaystyle+\epsilon^{4}\Big{(}-S_{2,2}(x)+\zeta_{4}-\text{Li}_{4}(1-x)+\text{Li}_{4}(x)-\frac{1}{2}\text{Li}_{2}(1-x)\ln^{2}(1-x)+\text{Li}_{3}(1-x)\ln(1-x)$
$\displaystyle-\frac{1}{6}\ln(x)\ln^{3}(1-x)\Big{)}+\epsilon^{5}\Big{(}S_{2,3}(x)-S_{3,2}(x)-\zeta_{5}+\text{Li}_{5}(1-x)+\text{Li}_{5}(x)-\frac{1}{6}\text{Li}_{2}(1-x)\ln^{3}(1-x)$
$\displaystyle+\frac{1}{2}\text{Li}_{3}(1-x)\ln^{2}(1-x)-\text{Li}_{4}(1-x)\ln(1-x)-\frac{1}{24}\ln(x)\ln^{4}(1-x)\Big{)}+\epsilon^{6}\Big{(}-S_{2,4}(x)+S_{3,3}(x)$
$\displaystyle-
S_{4,2}(x)+\zeta_{6}-\text{Li}_{6}(1-x)+\text{Li}_{6}(x)-\frac{1}{24}\text{Li}_{2}(1-x)\ln^{4}(1-x)+\frac{1}{6}\text{Li}_{3}(1-x)\ln^{3}(1-x)$
$\displaystyle-\frac{1}{2}\text{Li}_{4}(1-x)\ln^{2}(1-x)+\text{Li}_{5}(1-x)\ln(1-x)-\frac{1}{120}\ln(x)\ln^{5}(1-x)\Big{)}+\mathcal{O}\big{(}\epsilon^{7}\big{)}$
$\displaystyle=1+\frac{\epsilon^{2}}{1-\epsilon}x+\mathcal{O}\big{(}x^{2}\big{)}$
$\displaystyle=\Big{(}\Gamma(1-\epsilon)\Gamma(1+\epsilon)+\mathcal{O}\big{(}1-x\big{)}\Big{)}+(1-x)^{\epsilon}\mathcal{O}\big{(}1-x\big{)}\;,$
(73) ${}_{2}F_{1}$ $\displaystyle(1,1-\epsilon,2-2\epsilon,x)=$
$\displaystyle\frac{1}{x}\Big{(}-\ln(1-x)+\epsilon\Big{(}2\text{Li}_{2}(x)+\frac{1}{2}\ln(1-x)(\ln(1-x)+4)\Big{)}+\epsilon^{2}\Big{(}-2\zeta_{3}+2\text{Li}_{3}(1-x)+4\text{Li}_{3}(x)$
$\displaystyle\quad\;-2\text{Li}_{2}(1-x)\ln(1-x)-2\text{Li}_{2}(x)(\ln(1-x)+2)-\frac{1}{6}\ln^{3}(1-x)-\ln(x)\ln^{2}(1-x)-\ln^{2}(1-x)\Big{)}$
$\displaystyle\quad\;+\epsilon^{3}\Big{(}-4S_{2,2}(x)+4\zeta_{3}+2\zeta_{4}-4\text{Li}_{3}(1-x)-8\text{Li}_{3}(x)-2\text{Li}_{4}(1-x)+8\text{Li}_{4}(x)$
$\displaystyle\quad\;+\text{Li}_{2}(1-x)\ln(1-x)(\ln(1-x)+4)+\text{Li}_{2}(x)\ln(1-x)(\ln(1-x)+4)-4\text{Li}_{3}(x)\ln(1-x)$
$\displaystyle\quad\;+2\zeta_{3}\ln(1-x)+\frac{1}{24}\ln^{4}(1-x)+\frac{2}{3}\ln(x)\ln^{3}(1-x)+\frac{1}{3}\ln^{3}(1-x)+2\ln(x)\ln^{2}(1-x)\Big{)}$
$\displaystyle\quad\;+\epsilon^{4}\Big{(}8S_{2,2}(x)+4S_{2,3}(x)-8S_{3,2}(x)+4\ln(1-x)S_{2,2}(x)-4\zeta_{4}-2\zeta_{5}+4\text{Li}_{4}(1-x)-16\text{Li}_{4}(x)$
$\displaystyle\quad\;+2\text{Li}_{5}(1-x)+16\text{Li}_{5}(x)-\frac{1}{3}\text{Li}_{2}(1-x)(\ln(1-x)+6)\ln^{2}(1-x)$
$\displaystyle\quad\;-\frac{1}{3}\text{Li}_{2}(x)(\ln(1-x)+6)\ln^{2}(1-x)+2\text{Li}_{3}(x)\ln^{2}(1-x)+8\text{Li}_{3}(x)\ln(1-x)-8\text{Li}_{4}(x)\ln(1-x)$
$\displaystyle\quad\;-\zeta_{3}\ln^{2}(1-x)-4\zeta_{3}\ln(1-x)-2\zeta_{4}\ln(1-x)-\frac{1}{120}\ln^{5}(1-x)-\frac{1}{4}\ln(x)\ln^{4}(1-x)$
$\displaystyle\quad\;-\frac{1}{12}\ln^{4}(1-x)-\frac{4}{3}\ln(x)\ln^{3}(1-x)\Big{)}+\epsilon^{5}\Big{(}-8S_{2,3}(x)-4S_{2,4}(x)+16S_{3,2}(x)+8S_{3,3}(x)$
$\displaystyle\quad\;-16S_{4,2}(x)-2\ln^{2}(1-x)S_{2,2}(x)-8\ln(1-x)S_{2,2}(x)-4\ln(1-x)S_{2,3}(x)+8\ln(1-x)S_{3,2}(x)$
$\displaystyle\quad\;+4\zeta_{5}+2\zeta_{6}-4\text{Li}_{5}(1-x)-32\text{Li}_{5}(x)-2\text{Li}_{6}(1-x)+32\text{Li}_{6}(x)$
$\displaystyle\quad\;+\frac{1}{12}\text{Li}_{2}(1-x)(\ln(1-x)+8)\ln^{3}(1-x)+\frac{1}{12}\text{Li}_{2}(x)(\ln(1-x)+8)\ln^{3}(1-x)$
$\displaystyle\quad\;-\frac{2}{3}\text{Li}_{3}(x)\ln^{3}(1-x)-4\text{Li}_{3}(x)\ln^{2}(1-x)+4\text{Li}_{4}(x)\ln^{2}(1-x)+16\text{Li}_{4}(x)\ln(1-x)$
$\displaystyle\quad\;-16\text{Li}_{5}(x)\ln(1-x)+\frac{1}{3}\zeta_{3}\ln^{3}(1-x)+2\zeta_{3}\ln^{2}(1-x)+\zeta_{4}\ln^{2}(1-x)+4\zeta_{4}\ln(1-x)$
$\displaystyle\quad\;+2\zeta_{5}\ln(1-x)+\frac{1}{720}\ln^{6}(1-x)+\frac{1}{15}\ln(x)\ln^{5}(1-x)+\frac{1}{60}\ln^{5}(1-x)+\frac{1}{2}\ln(x)\ln^{4}(1-x)\Big{)}$
$\displaystyle\quad\;+\mathcal{O}\big{(}\epsilon^{6}\big{)}\Big{)}$
$\displaystyle=\Big{(}\Big{(}-\frac{1}{\epsilon}+2\Big{)}+\mathcal{O}\big{(}1-x\big{)}\Big{)}+(1-x)^{-\epsilon}\Big{(}\frac{(1-2\epsilon)\Gamma(1-2\epsilon)\Gamma(1+\epsilon)}{\epsilon\Gamma(1-\epsilon)}+\mathcal{O}\big{(}1-x\big{)}\Big{)}\;,$
(74) ${}_{2}F_{1}$ $\displaystyle(1,1,1-\epsilon,x)=$
$\displaystyle\frac{1}{1-x}\Big{(}1-\epsilon\ln(1-x)+\epsilon^{2}\Big{(}\text{Li}_{2}(x)+\frac{1}{2}\ln^{2}(1-x)\Big{)}+\epsilon^{3}\Big{(}-\zeta_{3}+\text{Li}_{3}(1-x)+\text{Li}_{3}(x)$
$\displaystyle\qquad\quad-\text{Li}_{2}(1-x)\ln(1-x)-\text{Li}_{2}(x)\ln(1-x)-\frac{1}{6}\ln^{3}(1-x)-\frac{1}{2}\ln(x)\ln^{2}(1-x)\Big{)}$
$\displaystyle\qquad\quad+\epsilon^{4}\Big{(}-S_{2,2}(x)+\zeta_{4}-\text{Li}_{4}(1-x)+\text{Li}_{4}(x)+\frac{1}{2}\text{Li}_{2}(1-x)\ln^{2}(1-x)+\frac{1}{2}\text{Li}_{2}(x)\ln^{2}(1-x)$
$\displaystyle\qquad\quad-\text{Li}_{3}(x)\ln(1-x)+\zeta_{3}\ln(1-x)+\frac{1}{24}\ln^{4}(1-x)+\frac{1}{3}\ln(x)\ln^{3}(1-x)\Big{)}+\epsilon^{5}\Big{(}S_{2,3}(x)$
$\displaystyle\qquad\quad-
S_{3,2}(x)+\ln(1-x)S_{2,2}(x)-\zeta_{5}+\text{Li}_{5}(1-x)+\text{Li}_{5}(x)-\frac{1}{6}\text{Li}_{2}(1-x)\ln^{3}(1-x)$
$\displaystyle\qquad\quad-\frac{1}{6}\text{Li}_{2}(x)\ln^{3}(1-x)+\frac{1}{2}\text{Li}_{3}(x)\ln^{2}(1-x)-\text{Li}_{4}(x)\ln(1-x)-\frac{1}{2}\zeta_{3}\ln^{2}(1-x)$
$\displaystyle\qquad\quad-\zeta_{4}\ln(1-x)-\frac{1}{120}\ln^{5}(1-x)-\frac{1}{8}\ln(x)\ln^{4}(1-x)\Big{)}+\epsilon^{6}\Big{(}-S_{2,4}(x)+S_{3,3}(x)-S_{4,2}(x)$
$\displaystyle\qquad\quad-\frac{1}{2}\ln^{2}(1-x)S_{2,2}(x)-\ln(1-x)S_{2,3}(x)+\ln(1-x)S_{3,2}(x)+\zeta_{6}-\text{Li}_{6}(1-x)+\text{Li}_{6}(x)$
$\displaystyle\qquad\quad+\frac{1}{24}\text{Li}_{2}(1-x)\ln^{4}(1-x)+\frac{1}{24}\text{Li}_{2}(x)\ln^{4}(1-x)-\frac{1}{6}\text{Li}_{3}(x)\ln^{3}(1-x)+\frac{1}{2}\text{Li}_{4}(x)\ln^{2}(1-x)$
$\displaystyle\qquad\quad-\text{Li}_{5}(x)\ln(1-x)+\frac{1}{6}\zeta_{3}\ln^{3}(1-x)+\frac{1}{2}\zeta_{4}\ln^{2}(1-x)+\zeta_{5}\ln(1-x)+\frac{1}{720}\ln^{6}(1-x)$
$\displaystyle\qquad\quad+\frac{1}{30}\ln(x)\ln^{5}(1-x)\Big{)}+\mathcal{O}\big{(}\epsilon^{7}\big{)}\Big{)}$
$\displaystyle=(1-x)^{-1-\epsilon}\Big{(}\Gamma(1-\epsilon)\Gamma(1+\epsilon)+\mathcal{O}\big{(}1-x\big{)}\Big{)}+\mathcal{O}\big{(}(1-x)^{0}\big{)}\;,$
(75) ${}_{2}F_{1}$ $\displaystyle(1,1,2+\epsilon,x)=$
$\displaystyle\frac{1}{x}\Big{(}-\ln(1-x)+\epsilon\Big{(}-\text{Li}_{2}(x)-\frac{1}{2}\ln(1-x)(\ln(1-x)+2)\Big{)}+\epsilon^{2}\Big{(}-\zeta_{3}+\text{Li}_{3}(1-x)+\text{Li}_{3}(x)$
$\displaystyle\quad\;-\text{Li}_{2}(1-x)\ln(1-x)-\text{Li}_{2}(x)(\ln(1-x)+1)-\frac{1}{6}\ln^{3}(1-x)-\frac{1}{2}\ln(x)\ln^{2}(1-x)$
$\displaystyle\quad\;-\frac{1}{2}\ln^{2}(1-x)\Big{)}+\epsilon^{3}\Big{(}S_{2,2}(x)-\zeta_{3}-\zeta_{4}+\text{Li}_{3}(1-x)+\text{Li}_{3}(x)+\text{Li}_{4}(1-x)-\text{Li}_{4}(x)$
$\displaystyle\quad\;-\frac{1}{2}\text{Li}_{2}(1-x)\ln(1-x)(\ln(1-x)+2)-\frac{1}{2}\text{Li}_{2}(x)\ln(1-x)(\ln(1-x)+2)+\text{Li}_{3}(x)\ln(1-x)$
$\displaystyle\quad\;-\zeta_{3}\ln(1-x)-\frac{1}{24}\ln^{4}(1-x)-\frac{1}{3}\ln(x)\ln^{3}(1-x)-\frac{1}{6}\ln^{3}(1-x)-\frac{1}{2}\ln(x)\ln^{2}(1-x)\Big{)}$
$\displaystyle\quad\;+\epsilon^{4}\Big{(}S_{2,2}(x)+S_{2,3}(x)-S_{3,2}(x)+\ln(1-x)S_{2,2}(x)-\zeta_{4}-\zeta_{5}+\text{Li}_{4}(1-x)-\text{Li}_{4}(x)+\text{Li}_{5}(1-x)$
$\displaystyle\quad\;+\text{Li}_{5}(x)-\frac{1}{6}\text{Li}_{2}(1-x)(\ln(1-x)+3)\ln^{2}(1-x)-\frac{1}{6}\text{Li}_{2}(x)(\ln(1-x)+3)\ln^{2}(1-x)$
$\displaystyle\quad\;+\frac{1}{2}\text{Li}_{3}(x)\ln^{2}(1-x)+\text{Li}_{3}(x)\ln(1-x)-\text{Li}_{4}(x)\ln(1-x)-\frac{1}{2}\zeta_{3}\ln^{2}(1-x)-\zeta_{3}\ln(1-x)$
$\displaystyle\quad\;-\zeta_{4}\ln(1-x)-\frac{1}{120}\ln^{5}(1-x)-\frac{1}{8}\ln(x)\ln^{4}(1-x)-\frac{1}{24}\ln^{4}(1-x)-\frac{1}{3}\ln(x)\ln^{3}(1-x)\Big{)}$
$\displaystyle\quad\;+\epsilon^{5}\Big{(}S_{2,3}(x)+S_{2,4}(x)-S_{3,2}(x)-S_{3,3}(x)+S_{4,2}(x)+\frac{1}{2}\ln^{2}(1-x)S_{2,2}(x)+\ln(1-x)S_{2,2}(x)$
$\displaystyle\quad\;+\ln(1-x)S_{2,3}(x)-\ln(1-x)S_{3,2}(x)-\zeta_{5}-\zeta_{6}+\text{Li}_{5}(1-x)+\text{Li}_{5}(x)+\text{Li}_{6}(1-x)-\text{Li}_{6}(x)$
$\displaystyle\quad\;-\frac{1}{24}\text{Li}_{2}(1-x)(\ln(1-x)+4)\ln^{3}(1-x)-\frac{1}{24}\text{Li}_{2}(x)(\ln(1-x)+4)\ln^{3}(1-x)$
$\displaystyle\quad\;+\frac{1}{6}\text{Li}_{3}(x)\ln^{3}(1-x)+\frac{1}{2}\text{Li}_{3}(x)\ln^{2}(1-x)-\frac{1}{2}\text{Li}_{4}(x)\ln^{2}(1-x)-\text{Li}_{4}(x)\ln(1-x)$
$\displaystyle\quad\;+\text{Li}_{5}(x)\ln(1-x)-\frac{1}{6}\zeta_{3}\ln^{3}(1-x)-\frac{1}{2}\zeta_{3}\ln^{2}(1-x)-\frac{1}{2}\zeta_{4}\ln^{2}(1-x)-\zeta_{4}\ln(1-x)$
$\displaystyle\quad\;-\zeta_{5}\ln(1-x)-\frac{1}{720}\ln^{6}(1-x)-\frac{1}{30}\ln(x)\ln^{5}(1-x)-\frac{1}{120}\ln^{5}(1-x)-\frac{1}{8}\ln(x)\ln^{4}(1-x)\Big{)}$
$\displaystyle\quad\;+\mathcal{O}\big{(}\epsilon^{6}\big{)}\Big{)}$
$\displaystyle=\Big{(}\Big{(}\frac{1}{\epsilon}+1\Big{)}+\mathcal{O}\big{(}1-x\big{)}\Big{)}+(1-x)^{\epsilon}\Big{(}-\frac{(1+\epsilon)\Gamma(1-\epsilon)\Gamma(1+\epsilon)}{\epsilon}+\mathcal{O}\big{(}1-x\big{)}\Big{)}\;.$
(76)

The following functions are either special cases or occur in the asymptotic
behaviour of the Appell $F_{1}$ function discussed in Appendix B:

${}_{2}F_{1}$ $\displaystyle(1,1-\epsilon,1-2\epsilon,x)=$
$\displaystyle\frac{1}{1-x}\Big{(}1-\epsilon\ln(1-x)+\epsilon^{2}\Big{(}2\text{Li}_{2}(x)+\frac{1}{2}\ln^{2}(1-x)\Big{)}+\epsilon^{3}\Big{(}-2\zeta_{3}+2\text{Li}_{3}(1-x)+4\text{Li}_{3}(x)$
$\displaystyle\qquad\quad-2\text{Li}_{2}(1-x)\ln(1-x)-2\text{Li}_{2}(x)\ln(1-x)-\frac{1}{6}\ln^{3}(1-x)-\ln(x)\ln^{2}(1-x)\Big{)}$
$\displaystyle\qquad\quad+\epsilon^{4}\Big{(}-4S_{2,2}(x)+2\zeta_{4}-2\text{Li}_{4}(1-x)+8\text{Li}_{4}(x)+\text{Li}_{2}(1-x)\ln^{2}(1-x)+\text{Li}_{2}(x)\ln^{2}(1-x)$
$\displaystyle\qquad\quad-4\text{Li}_{3}(x)\ln(1-x)+2\zeta_{3}\ln(1-x)+\frac{1}{24}\ln^{4}(1-x)+\frac{2}{3}\ln(x)\ln^{3}(1-x)\Big{)}+\epsilon^{5}\Big{(}4S_{2,3}(x)$
$\displaystyle\qquad\quad-8S_{3,2}(x)+4\ln(1-x)S_{2,2}(x)-2\zeta_{5}+2\text{Li}_{5}(1-x)+16\text{Li}_{5}(x)$
$\displaystyle\qquad\quad-\frac{1}{3}\text{Li}_{2}(1-x)\ln^{3}(1-x)-\frac{1}{3}\text{Li}_{2}(x)\ln^{3}(1-x)+2\text{Li}_{3}(x)\ln^{2}(1-x)-8\text{Li}_{4}(x)\ln(1-x)$
$\displaystyle\qquad\quad-\zeta_{3}\ln^{2}(1-x)-2\zeta_{4}\ln(1-x)-\frac{1}{120}\ln^{5}(1-x)-\frac{1}{4}\ln(x)\ln^{4}(1-x)\Big{)}+\epsilon^{6}\Big{(}-4S_{2,4}(x)$
$\displaystyle\qquad\quad+8S_{3,3}(x)-16S_{4,2}(x)-2\ln^{2}(1-x)S_{2,2}(x)-4\ln(1-x)S_{2,3}(x)+8\ln(1-x)S_{3,2}(x)$
$\displaystyle\qquad\quad+2\zeta_{6}-2\text{Li}_{6}(1-x)+32\text{Li}_{6}(x)+\frac{1}{12}\text{Li}_{2}(1-x)\ln^{4}(1-x)+\frac{1}{12}\text{Li}_{2}(x)\ln^{4}(1-x)$
$\displaystyle\qquad\quad-\frac{2}{3}\text{Li}_{3}(x)\ln^{3}(1-x)+4\text{Li}_{4}(x)\ln^{2}(1-x)-16\text{Li}_{5}(x)\ln(1-x)+\frac{1}{3}\zeta_{3}\ln^{3}(1-x)$
$\displaystyle\qquad\quad+\zeta_{4}\ln^{2}(1-x)+2\zeta_{5}\ln(1-x)+\frac{1}{720}\ln^{6}(1-x)+\frac{1}{15}\ln(x)\ln^{5}(1-x)\Big{)}+\mathcal{O}\big{(}\epsilon^{7}\big{)}\Big{)}$
$\displaystyle=(1-x)^{-1-\epsilon}\Big{(}\frac{\Gamma(1-2\epsilon)\Gamma(1+\epsilon)}{\Gamma(1-\epsilon)}+\mathcal{O}\big{(}1-x\big{)}\Big{)}+\mathcal{O}\big{(}(1-x)^{0}\big{)}\;,$
(77) ${}_{2}F_{1}$ $\displaystyle(1,-\epsilon,1-2\epsilon,x)=$ $\displaystyle
1+\epsilon\ln(1-x)+\epsilon^{2}\Big{(}-2\text{Li}_{2}(x)-\frac{1}{2}\ln^{2}(1-x)\Big{)}+\epsilon^{3}\Big{(}2\zeta_{3}-2\text{Li}_{3}(1-x)-4\text{Li}_{3}(x)$
$\displaystyle+2\text{Li}_{2}(1-x)\ln(1-x)+2\text{Li}_{2}(x)\ln(1-x)+\frac{1}{6}\ln^{3}(1-x)+\ln(x)\ln^{2}(1-x)\Big{)}+\epsilon^{4}\Big{(}4S_{2,2}(x)$
$\displaystyle-2\zeta_{4}+2\text{Li}_{4}(1-x)-8\text{Li}_{4}(x)-\text{Li}_{2}(1-x)\ln^{2}(1-x)-\text{Li}_{2}(x)\ln^{2}(1-x)+4\text{Li}_{3}(x)\ln(1-x)$
$\displaystyle-2\zeta_{3}\ln(1-x)-\frac{1}{24}\ln^{4}(1-x)-\frac{2}{3}\ln(x)\ln^{3}(1-x)\Big{)}+\epsilon^{5}\Big{(}-4S_{2,3}(x)+8S_{3,2}(x)$
$\displaystyle-4\ln(1-x)S_{2,2}(x)+2\zeta_{5}-2\text{Li}_{5}(1-x)-16\text{Li}_{5}(x)+\frac{1}{3}\text{Li}_{2}(1-x)\ln^{3}(1-x)+\frac{1}{3}\text{Li}_{2}(x)\ln^{3}(1-x)$
$\displaystyle-2\text{Li}_{3}(x)\ln^{2}(1-x)+8\text{Li}_{4}(x)\ln(1-x)+\zeta_{3}\ln^{2}(1-x)+2\zeta_{4}\ln(1-x)+\frac{1}{120}\ln^{5}(1-x)$
$\displaystyle+\frac{1}{4}\ln(x)\ln^{4}(1-x)\Big{)}+\epsilon^{6}\Big{(}4S_{2,4}(x)-8S_{3,3}(x)+16S_{4,2}(x)+2\ln^{2}(1-x)S_{2,2}(x)$
$\displaystyle+4\ln(1-x)S_{2,3}(x)-8\ln(1-x)S_{3,2}(x)-2\zeta_{6}+2\text{Li}_{6}(1-x)-32\text{Li}_{6}(x)-\frac{1}{12}\text{Li}_{2}(1-x)\ln^{4}(1-x)$
$\displaystyle-\frac{1}{12}\text{Li}_{2}(x)\ln^{4}(1-x)+\frac{2}{3}\text{Li}_{3}(x)\ln^{3}(1-x)-4\text{Li}_{4}(x)\ln^{2}(1-x)+16\text{Li}_{5}(x)\ln(1-x)$
$\displaystyle-\frac{1}{3}\zeta_{3}\ln^{3}(1-x)-\zeta_{4}\ln^{2}(1-x)-2\zeta_{5}\ln(1-x)-\frac{1}{720}\ln^{6}(1-x)-\frac{1}{15}\ln(x)\ln^{5}(1-x)\Big{)}+\mathcal{O}\big{(}\epsilon^{7}\big{)}$
$\displaystyle=\Big{(}2+\mathcal{O}\big{(}1-x\big{)}\Big{)}+(1-x)^{-\epsilon}\Big{(}-\frac{\Gamma(1-2\epsilon)\Gamma(1+\epsilon)}{\Gamma(1-\epsilon)}+\mathcal{O}\big{(}1-x\big{)}\Big{)}\;,$
(78) ${}_{2}F_{1}$ $\displaystyle(-\epsilon,1+\epsilon,1-\epsilon,x)=$
$\displaystyle
1+\epsilon\ln(1-x)+\epsilon^{2}\Big{(}-2\text{Li}_{2}(x)-\frac{1}{2}\ln^{2}(1-x)\Big{)}+\epsilon^{3}\Big{(}-2\zeta_{3}+2\text{Li}_{3}(1-x)-2\text{Li}_{3}(x)$
$\displaystyle-2\text{Li}_{2}(1-x)\ln(1-x)+\frac{1}{6}\ln^{3}(1-x)-\ln(x)\ln^{2}(1-x)\Big{)}+\epsilon^{4}\Big{(}-2S_{2,2}(x)-2\zeta_{4}+2\text{Li}_{4}(1-x)$
$\displaystyle-2\text{Li}_{4}(x)+\text{Li}_{2}(1-x)\ln^{2}(1-x)-2\text{Li}_{3}(1-x)\ln(1-x)-\frac{1}{24}\ln^{4}(1-x)+\frac{1}{3}\ln(x)\ln^{3}(1-x)\Big{)}$
$\displaystyle+\epsilon^{5}\Big{(}-2S_{2,3}(x)-2S_{3,2}(x)-2\zeta_{5}+2\text{Li}_{5}(1-x)-2\text{Li}_{5}(x)-\frac{1}{3}\text{Li}_{2}(1-x)\ln^{3}(1-x)$
$\displaystyle+\text{Li}_{3}(1-x)\ln^{2}(1-x)-2\text{Li}_{4}(1-x)\ln(1-x)+\frac{1}{120}\ln^{5}(1-x)-\frac{1}{12}\ln(x)\ln^{4}(1-x)\Big{)}$
$\displaystyle+\epsilon^{6}\Big{(}-2S_{2,4}(x)-2S_{3,3}(x)-2S_{4,2}(x)-2\zeta_{6}+2\text{Li}_{6}(1-x)-2\text{Li}_{6}(x)+\frac{1}{12}\text{Li}_{2}(1-x)\ln^{4}(1-x)$
$\displaystyle-\frac{1}{3}\text{Li}_{3}(1-x)\ln^{3}(1-x)+\text{Li}_{4}(1-x)\ln^{2}(1-x)-2\text{Li}_{5}(1-x)\ln(1-x)-\frac{1}{720}\ln^{6}(1-x)$
$\displaystyle+\frac{1}{60}\ln(x)\ln^{5}(1-x)\Big{)}+\mathcal{O}\big{(}\epsilon^{7}\big{)}$
$\displaystyle=\Big{(}\frac{2\Gamma(1-\epsilon)^{2}}{\Gamma(1-2\epsilon)}+\mathcal{O}\big{(}1-x\big{)}\Big{)}+(1-x)^{-\epsilon}\Big{(}-1+\mathcal{O}\big{(}1-x\big{)}\Big{)}\;.$
(79)

##  Appendix B Definition, representations and expansions of the Appell
$F_{1}$ function

The Appell $F_{1}$ function is defined as follows:

$\begin{split}F_{1}&(a,b_{1},b_{2},c,x,y)\equiv\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\frac{\Gamma(a+m+n)\Gamma(b_{1}+m)\Gamma(b_{2}+n)\Gamma(c)}{\Gamma(a)\Gamma(b_{1})\Gamma(b_{2})\Gamma(c+m+n)}\frac{x^{m}}{m!}\frac{y^{n}}{n!}\\\\[11.38092pt]
&=\frac{\Gamma(c)}{\Gamma(a)\Gamma(c-a)}\int_{0}^{1}\differential{t}t^{a-1}(1-t)^{c-a-1}(1-tx)^{-b_{1}}(1-ty)^{-b_{2}}\\\\[11.38092pt]
&=\frac{1}{(2\pi
i)^{2}}\iint_{\mathcal{C}}\differential{z_{1}}\differential{z_{2}}\frac{\Gamma(a+z_{1}+z_{2})\Gamma(b_{1}+z_{1})\Gamma(b_{2}+z_{2})\Gamma(c)\Gamma(-z_{1})\Gamma(-z_{2})}{\Gamma(a)\Gamma(b_{1})\Gamma(b_{2})\Gamma(c+z_{1}+z_{2})}(-x)^{z_{1}}(-y)^{z_{2}}\\\\[11.38092pt]
&=\frac{1}{(2\pi
i)^{2}}\iint_{\mathcal{C}}\differential{z_{1}}\differential{z_{2}}\frac{\Gamma(a+z_{1}+z_{2})\Gamma(b_{1}+z_{1})\Gamma(b_{2}+z_{2})\Gamma(c-a-
b_{1}-z_{1})\Gamma(c)}{\Gamma(a)\Gamma(b_{1})\Gamma(b_{2})\Gamma(c-a)\Gamma(c-b_{1}+z_{2})}\\\\[5.69046pt]
&\qquad\qquad\qquad\qquad\times\Gamma(-z_{1})\Gamma(-z_{2})(1-x)^{z_{1}}(-y)^{z_{2}}\\\\[11.38092pt]
&=\frac{1}{(2\pi
i)^{2}}\iint_{\mathcal{C}}\differential{z_{1}}\differential{z_{2}}\frac{\Gamma(a+z_{1}+z_{2})\Gamma(b_{1}+z_{1})\Gamma(b_{2}+z_{2})\Gamma(c-a-
b_{1}-b_{2}-z_{1}-z_{2})\Gamma(c)}{\Gamma(a)\Gamma(b_{1})\Gamma(b_{2})\Gamma(c-a)\Gamma(c-b_{1}-b_{2})}\\\\[5.69046pt]
&\qquad\qquad\qquad\qquad\times\Gamma(-z_{1})\Gamma(-z_{2})(1-x)^{z_{1}}(1-y)^{z_{2}}\;.\end{split}$
(80)

As in the case of the hypergeometric function Appendix A, the integration
contours of the Mellin-Barnes representations starting with the third
equality, $\mathcal{C}$, should be chosen along the imaginary axis and
separate the poles of the $\Gamma$-functions $\Gamma(\dots+z)$ and
$\Gamma(\dots-z)$. It turns out that only one Appell function is actually
needed in the present case. We choose it to have arguments $0<x\leq y<1$
ensuring the convergence of the series representation Eq. (80):

$\begin{split}F_{1}(1&,-\epsilon,1,1-2\epsilon,1-x,1-y)=\\\\[5.69046pt]
&\frac{1}{y}\Big{(}1+\epsilon(\ln(x)-2\ln(y))+\frac{1}{2}\epsilon^{2}\Big{(}4\text{Li}_{2}\Big{(}\frac{x}{y}\Big{)}-4\text{Li}_{2}(y)-4\ln(x)\ln(y)-4\ln(y)\ln(y-x)\\\\[5.69046pt]
&\quad+4\ln(x)\ln(y-x)-\ln^{2}(x)+6\ln^{2}(y)-4\ln(1-y)\ln(y)\Big{)}+\frac{1}{6}\epsilon^{3}\Big{(}-12\zeta_{3}-12\text{Li}_{3}\Big{(}1-\frac{1-y}{1-x}\Big{)}\\\\[5.69046pt]
&\quad-24\text{Li}_{3}\Big{(}1-\frac{x}{y}\Big{)}+12\text{Li}_{3}\Big{(}1-\frac{x(1-y)}{(1-x)y}\Big{)}-12\text{Li}_{3}\Big{(}\frac{1-y}{1-x}\Big{)}+12\text{Li}_{3}\Big{(}\frac{x(1-y)}{(1-x)y}\Big{)}\\\\[5.69046pt]
&\quad-12\ln(x)\text{Li}_{2}\Big{(}\frac{1-y}{1-x}\Big{)}-12\ln(x)\text{Li}_{2}\Big{(}\frac{x}{y}\Big{)}-12\text{Li}_{2}(y)\ln(x)-12\text{Li}_{3}(x)+12\text{Li}_{2}(x)\ln(x)\\\\[5.69046pt]
&\quad+24\text{Li}_{3}(1-y)+24\text{Li}_{3}(y)+24\zeta_{2}\ln(x)-12\ln^{2}(1-x)\ln(y)+6\ln^{2}(x)\ln(y)\\\\[5.69046pt]
&\quad-6\ln^{2}(x)\ln(y-x)-12\ln(1-x)\ln^{2}(y)+6\ln^{2}(y)\ln(y-x)\\\\[5.69046pt]
&\quad+12\ln(1-x)\ln(x)\ln(y)+12\ln(1-x)\ln(1-y)\ln(y)-12\ln(x)\ln(1-y)\ln(y)\\\\[5.69046pt]
&\quad+12\ln(1-x)\ln(y)\ln(y-x)-12\ln(1-y)\ln(y)\ln(y-x)+\ln^{3}(x)-12\zeta_{2}\ln(y)-8\ln^{3}(y)\\\\[5.69046pt]
&\quad+18\ln(1-y)\ln^{2}(y)\Big{)}+\mathcal{O}\big{(}\epsilon^{4}\big{)}\Big{)}\;.\end{split}$
(81)

The above result has been obtained with the software package PolyLogTools
Duhr:2019tlz  by direct integration of the integral representation after
remapping the integration region from $[0,1]$ to $[0,\infty)$. The series
expansion has been truncated after enough terms to match the requirements of
N3LO applications as discussed in the main part of the text. Care has been
taken, so that the arguments of all functions belong to the unit interval. The
following special cases are also required with deeper $\epsilon$-expansions:

$\displaystyle\text{F}_{1}(1,-\epsilon,1,1-2\epsilon,1-x,1-x)$
$\displaystyle={}_{2}F_{1}(1,1-\epsilon,1-2\epsilon,1-x)\;,$ (82)
$\displaystyle\text{F}_{1}(1,-\epsilon,1,1-2\epsilon,1-x,0)$
$\displaystyle={}_{2}F_{1}(1,-\epsilon,1-2\epsilon,1-x)\;.$

The corresponding $\epsilon$-expansions of the hypergeometric functions are
given in Eqs. (77) and (78), while the above expressions are a direct
consequence of the integral representations in Eqs. (80) and (70). Finally,
the following asymptotics of the Appell $F_{1}$ function are needed as well:

$\begin{split}\text{F}_{1}(1,-\epsilon,1,1-2\epsilon,1-x,1-y)\;&\sim\;2\,_{2}F_{1}(1,1,1-\epsilon,1-y)-\frac{x^{-\epsilon}\Gamma(1-2\epsilon)\Gamma(1+\epsilon)}{y\Gamma(1-\epsilon)}\quad\big{(}x\to
0\big{)}\;,\\\\[11.38092pt]
\text{F}_{1}(1,-\epsilon,1,1-2\epsilon,1-x,1-y)\;&\sim\;\frac{y^{-1-\epsilon}\Gamma(1-2\epsilon)\Gamma(1+\epsilon)}{\Gamma(1-\epsilon)}\,_{2}F_{1}\Big{(}-\epsilon,1+\epsilon,1-\epsilon,1-\frac{x}{y}\Big{)}\quad\big{(}y\to
0\big{)}\;.\end{split}$ (83)

These results can be obtained by an asymptotic expansion of the Mellin-Barnes
representations Eq. (80), while the two hypergeometric functions are given in
Eqs. (75) and (79).

##  Appendix C Description of the results provided in ancillary files

The following files are provided:

  * •

integrals.m, expandedIntegrals.m \- substitution lists for the exact
(integrals) and expanded up to $\order{\epsilon^{0}}$ (expandedIntegrals)
master integrals, Eqs. (52), (53), (3.1), (55), (3.1), (3.1) and (59), divided
by $r_{\Gamma}$, Eq. (47), including the results with required permutations of
the external momenta; the notation for the integrals in the substitutions
matches the occurrences in the files containing the splitting operators and
splitting functions listed below; the result for the light-cone-gauge six-
dimensional box integral of Section 3.2 is not provided here but rather
substituted as LCGBox[pi,pj,pk,d+2], where pi,pj,pk correspond to a
permutation of $p_{1},p_{2}$ and $p_{3}$ for which the result is provided in
Eq. (68) and in the file LCGBox.m;

  * •

AppellF1.m \- substitution giving the $\epsilon$-expansion of the Appell
function Eq. (81);

  * •

LCGBox.m, LCGBoxDoubleSoft.m \- right-hand side of Eq. (68) named LCGBoxEp0,
and its double-soft limit at $\order{\epsilon^{0}}$ (LCGBoxDoubleSoftEp0) and
$\order{\epsilon}$ (LCGBoxDoubleSoftEp1);

  * •

Pggg.m, Pgqqbar.m, Pqgg.m, Pqqqbar.m, Pqqpqpbar.m \- exact results at tree and
one-loop level, the latter expressed through the master integrals, for the
splitting (P0Pol and P1Pol) and averaged splitting (P0Avg and P1Avg)
functions, Eq. (7), for the processes $g\to ggg$, $g\to gq\bar{q}$, $q\to
qgg$, $q\to qq\bar{q}$ and $q\to qq^{\prime}\bar{q}^{\prime}$ respectively;
the tree-level results, P0Pol and P0Avg, additionally contain the
singularities, proportional to the tag Ioperator, of the one-loop splitting
and averaged splitting functions obtained using Eq. (5) that has been
partially expanded in $\epsilon$ as described in the text after the equation;
the one-loop level results, P1Pol and P1Avg, defined in Eq. (8), are given
without the hermitian conjugate;

  * •

P1gggExp.m, P1gqqbarExp.m, P1qggExp.m, P1qqqbarExp.m, P1qqpqpbarExp.m \- one-
loop splitting (P1PolExp) and averaged splitting (P1AvgExp) functions from the
files listed in the previous item, expanded up to $\order{\epsilon^{0}}$ after
substitution of the master integrals from expandedIntegrals, which implies
that they are divided by $r_{\Gamma}$ as well;

  * •

ggg0l.m, gqqbar0l.m, qgg0l.m, qqqbar0l.m, qqpqpbar0l.m, ggg1l.m, gqqbar1l.m,
qgg1l.m, qqqbar1l.m, qqpqpbar1l.m \- exact splitting operators at tree (*0l.m)
and one-loop (*1l.m) level, the latter expressed through the master integrals;
each operator has been multiplied with $1/2\,s_{123}$ to make it dimensionless
and easier to relate to the respective splitting function; the operators are
named correspondingly to the file name (without the extension .m of course);
the tree-level results contain the singularities, tagged with the factor
Ioperator, of the one-loop splitting operators obtained using Eq. (5) that has
partially been expanded in $\epsilon$ as described in the text after the
equation;

  * •

gqqbar0l4D.m, qgg0l4D.m, qqqbar0l4D.m, qqpqpbar0l4D.m, gqqbar1l4D.m,
qgg1l4D.m,
qqqbar1l4D.m, qqpqpbar1l4D.m \- four-dimensional projections, obtained using
Eq. (14), of the results listed in the previous item at tree and one-loop
level; the projections of the pure-gluon splitting operators are not provided,
since they are lengthier than the original expressions and lack usefulness;

  * •

README \- notation used in the files.

Although the splitting operators and splitting functions are not renormalised,
subtracting from them the singularities given in Eq. (5) as provided in the
above files, gives the expressions needed to correctly subtract the triple-
collinear limit of one-loop amplitudes according to Eq. (42), in the ‘t Hooft-
Veltman scheme (four-dimensional external gluon polarisation vectors, but
$d$-dimensional internal virtual gluon fields). These expressions thus match
the default conventions of the majority of software one-loop amplitude
providers, e.g. NJet Badger:2012pg .

## References

  * (1) S. Catani and M.H. Seymour, _A General algorithm for calculating jet cross-sections in NLO QCD_ , _Nucl. Phys. B_ 485 (1997) 291 [hep-ph/9605323].
  * (2) S. Frixione, Z. Kunszt and A. Signer, _Three jet cross-sections to next-to-leading order_ , _Nucl. Phys. B_ 467 (1996) 399 [hep-ph/9512328].
  * (3) F. Caola, W. Chen, C. Duhr, X. Liu, B. Mistlberger, F. Petriello et al., _The Path forward to N 3LO_, in _2022 Snowmass Summer Study_ , Seattle, WA, U.S.A., 17–26 July 2022 [2203.06730].
  * (4) J.M. Campbell and E.W.N. Glover, _Double unresolved approximations to multiparton scattering amplitudes_ , _Nucl. Phys. B_ 527 (1998) 264 [hep-ph/9710255].
  * (5) S. Catani and M. Grazzini, _Collinear factorization and splitting functions for next-to-next-to-leading order QCD calculations_ , _Phys. Lett. B_ 446 (1999) 143 [hep-ph/9810389].
  * (6) S. Catani and M. Grazzini, _Infrared factorization of tree level QCD amplitudes at the next-to-next-to-leading order and beyond_ , _Nucl. Phys. B_ 570 (2000) 287 [hep-ph/9908523].
  * (7) V. Del Duca, A. Frizzo and F. Maltoni, _Factorization of tree QCD amplitudes in the high-energy limit and in the collinear limit_ , _Nucl. Phys. B_ 568 (2000) 211 [hep-ph/9909464].
  * (8) D.A. Kosower, _Multiple singular emission in gauge theories_ , _Phys. Rev. D_ 67 (2003) 116003 [hep-ph/0212097].
  * (9) S. Catani, D. Colferai and A. Torrini, _Triple (and quadruple) soft-gluon radiation in QCD hard scattering_ , _JHEP_ 01 (2020) 118 [1908.01616].
  * (10) V. Del Duca, C. Duhr, R. Haindl, A. Lazopoulos and M. Michel, _Tree-level splitting amplitudes for a quark into four collinear partons_ , _JHEP_ 02 (2020) 189 [1912.06425].
  * (11) V. Del Duca, C. Duhr, R. Haindl, A. Lazopoulos and M. Michel, _Tree-level splitting amplitudes for a gluon into four collinear partons_ , _JHEP_ 10 (2020) 093 [2007.05345].
  * (12) O. Braun-White and N. Glover, _Decomposition of Triple Collinear Splitting Functions_ , 2204.10755.
  * (13) Z. Bern, L.J. Dixon, D.C. Dunbar and D.A. Kosower, _One loop n point gauge theory amplitudes, unitarity and collinear limits_ , _Nucl. Phys. B_ 425 (1994) 217 [hep-ph/9403226].
  * (14) Z. Bern, V. Del Duca and C.R. Schmidt, _The Infrared behavior of one loop gluon amplitudes at next-to-next-to-leading order_ , _Phys. Lett. B_ 445 (1998) 168 [hep-ph/9810409].
  * (15) D.A. Kosower and P. Uwer, _One loop splitting amplitudes in gauge theory_ , _Nucl. Phys. B_ 563 (1999) 477 [hep-ph/9903515].
  * (16) Z. Bern, V. Del Duca, W.B. Kilgore and C.R. Schmidt, _The infrared behavior of one loop QCD amplitudes at next-to-next-to leading order_ , _Phys. Rev. D_ 60 (1999) 116001 [hep-ph/9903516].
  * (17) S. Catani and M. Grazzini, _The soft gluon current at one loop order_ , _Nucl. Phys. B_ 591 (2000) 435 [hep-ph/0007142].
  * (18) S. Catani, D. de Florian and G. Rodrigo, _The Triple collinear limit of one loop QCD amplitudes_ , _Phys. Lett. B_ 586 (2004) 323 [hep-ph/0312067].
  * (19) G.F.R. Sborlini, D. de Florian and G. Rodrigo, _Triple collinear splitting functions at NLO for scattering processes with photons_ , _JHEP_ 10 (2014) 161 [1408.4821].
  * (20) G.F.R. Sborlini, D. de Florian and G. Rodrigo, _Polarized triple-collinear splitting functions at NLO for processes with photons_ , _JHEP_ 03 (2015) 021 [1409.6137].
  * (21) S. Badger, F. Buciuni and T. Peraro, _One-loop triple collinear splitting amplitudes in QCD_ , _JHEP_ 09 (2015) 188 [1507.05070].
  * (22) Y.J. Zhu, _Double soft current at one-loop in QCD_ , 2009.08919.
  * (23) S. Catani and L. Cieri, _Multiple soft radiation at one-loop order and the emission of a soft quark–antiquark pair_ , _Eur. Phys. J. C_ 82 (2022) 97 [2108.13309].
  * (24) Z. Bern, L.J. Dixon and D.A. Kosower, _Two-loop g — $>$ gg splitting amplitudes in QCD_, _JHEP_ 08 (2004) 012 [hep-ph/0404293].
  * (25) S.D. Badger and E.W.N. Glover, _Two loop splitting functions in QCD_ , _JHEP_ 07 (2004) 040 [hep-ph/0405236].
  * (26) C. Duhr, T. Gehrmann and M. Jaquier, _Two-loop splitting amplitudes and the single-real contribution to inclusive Higgs production at N 3LO_, _JHEP_ 02 (2015) 077 [1411.3587].
  * (27) Y. Li and H.X. Zhu, _Single soft gluon emission at two loops_ , _JHEP_ 11 (2013) 080 [1309.4391].
  * (28) C. Duhr and T. Gehrmann, _The two-loop soft current in dimensional regularization_ , _Phys. Lett. B_ 727 (2013) 452 [1309.4393].
  * (29) Wolfram Research Inc., “Mathematica, Version 13.0.0.”
  * (30) I. Feige and M.D. Schwartz, _Hard-Soft-Collinear Factorization to All Orders_ , _Phys. Rev. D_ 90 (2014) 105020 [1403.6472].
  * (31) S. Catani, S. Dittmaier and Z. Trocsanyi, _One loop singular behavior of QCD and SUSY QCD amplitudes with massive partons_ , _Phys. Lett. B_ 500 (2001) 149 [hep-ph/0011222].
  * (32) T. Hahn, _Generating Feynman diagrams and amplitudes with FeynArts 3_ , _Comput. Phys. Commun._ 140 (2001) 418 [hep-ph/0012260].
  * (33) J.A.M. Vermaseren, _New features of FORM_ , math-ph/0010025.
  * (34) G. Passarino and M.J.G. Veltman, _One Loop Corrections for e+ e- Annihilation Into mu+ mu- in the Weinberg Model_ , _Nucl. Phys. B_ 160 (1979) 151.
  * (35) R.H. Lewis, “Computer algebra system fermat.”
  * (36) Z. Bern, L.J. Dixon and D.A. Kosower, _Dimensionally regulated pentagon integrals_ , _Nucl. Phys. B_ 412 (1994) 751 [hep-ph/9306240].
  * (37) P. Maierhöfer, J. Usovitsch and P. Uwer, _Kira—A Feynman integral reduction program_ , _Comput. Phys. Commun._ 230 (2018) 99 [1705.05610].
  * (38) J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch, _Integral reduction with Kira 2.0 and finite field methods_ , _Comput. Phys. Commun._ 266 (2021) 108024 [2008.06494].
  * (39) P. Cvitanovic, _Group theory for Feynman diagrams in non-Abelian gauge theories_ , _Phys. Rev. D_ 14 (1976) 1536.
  * (40) S. Badger, B. Biedermann, P. Uwer and V. Yundin, _Numerical evaluation of virtual corrections to multi-jet production in massless QCD_ , _Comput. Phys. Commun._ 184 (2013) 1981 [1209.0100].
  * (41) G.F.R. Sborlini, _Comportamiento singular de las amplitudes de scattering en el límite colineal_ , Ph.D. thesis, University of Buenos Aires, Buenos Aires, Argentina, 2014.
  * (42) V.A. Smirnov, _Feynman integral calculus_ (2006).
  * (43) C. Anastasiou, C. Duhr, F. Dulat and B. Mistlberger, _Soft triple-real radiation for Higgs production at N3LO_ , _JHEP_ 07 (2013) 003 [1302.4379].
  * (44) O.V. Tarasov, _Connection between Feynman integrals having different values of the space-time dimension_ , _Phys. Rev. D_ 54 (1996) 6479 [hep-th/9606018].
  * (45) O.V. Tarasov, _Generalized recurrence relations for two loop propagator integrals with arbitrary masses_ , _Nucl. Phys. B_ 502 (1997) 455 [hep-ph/9703319].
  * (46) C. Duhr and F. Dulat, _PolyLogTools — polylogs for the masses_ , _JHEP_ 08 (2019) 135 [1904.07279].
  * (47) C.W. Bauer, A. Frink and R. Kreckel, _Introduction to the GiNaC framework for symbolic computation within the C++ programming language_ , _J. Symb. Comput._ 33 (2002) 1 [cs/0004015].
  * (48) J. Vollinga and S. Weinzierl, _Numerical evaluation of multiple polylogarithms_ , _Comput. Phys. Commun._ 167 (2005) 177 [hep-ph/0410259].
  * (49) M. Czakon, _Automatized analytic continuation of Mellin-Barnes integrals_ , _Comput. Phys. Commun._ 175 (2006) 559 [hep-ph/0511200].
  * (50) D. Maitre, _HPL, a mathematica implementation of the harmonic polylogarithms_ , _Comput. Phys. Commun._ 174 (2006) 222 [hep-ph/0507152].
  * (51) D. Maitre, _Extension of HPL to complex arguments_ , _Comput. Phys. Commun._ 183 (2012) 846 [hep-ph/0703052].
  * (52) T. Huber and D. Maitre, _HypExp: A Mathematica package for expanding hypergeometric functions around integer-valued parameters_ , _Comput. Phys. Commun._ 175 (2006) 122 [hep-ph/0507094].