# Same-hemisphere three-gluon-emission contribution to the zero-jettiness soft
function at N3LO QCD

Daniel Baranowski<EMAIL_ADDRESS>Institut für Theoretische
Teilchenphysik, Karlsruher Institut für Technologie (KIT), 76128 Karlsruhe,
Germany  Maximilian Delto<EMAIL_ADDRESS>Physik Department,
Technische Universität München, James-Franck-Strasse 1, 85748, Germany  Kirill
Melnikov<EMAIL_ADDRESS>Institut für Theoretische Teilchenphysik,
Karlsruher Institut für Technologie (KIT), 76128 Karlsruhe, Germany  Chen-Yu
Wang<EMAIL_ADDRESS>Institut für Theoretische Teilchenphysik,
Karlsruher Institut für Technologie (KIT), 76128 Karlsruhe, Germany

###### Abstract

We complete the calculation of the three-gluon-emission contribution to the
same-hemisphere part of the zero-jettiness soft function at next-to-next-to-
next-to-leading order in perturbative QCD.

††preprint: P3H-22-039,TTP22-025,TUM-HEP-1396/22

##  I Introduction

The goal of this paper is to present the result for the same-hemisphere three-
gluon-emission contribution to the zero-jettiness soft function at next-to-
next-to-next-to-leading order in perturbative QCD. Computation of such
corrections needs to overcome several technical challenges that were discussed
in Ref. [1]. These challenges stem from the fact that an observable that
defines the soft function – the so-called jettiness – involves Heaviside
functions. These Heaviside functions are needed to distinguish between
emissions of soft gluons into two hemispheres, defined relative to directions
of incoming hard radiators.

Presence of Heaviside functions complicates application of generalized
unitarity [2] and integration-by-parts identities [3] to phase-space
integrals. We have discussed in Ref. [1] how to overcome this problem and
explained how to derive useful integration-by-parts relations for integrals
with Heaviside functions. To show the efficacy of this method, we employed an
eikonal function derived in Ref. [4], which describes emissions of three soft
gluons, and integrated it over soft-gluon phase space subject to zero-
jettiness constraints. We restricted ourselves to contributions where all
gluons are emitted into the same hemisphere.

Unfortunately, in Ref. [1], we have not completed the computation of this
“same-hemisphere” contribution. Indeed, the representation of the eikonal
function derived in Ref. [4] involves four terms and in Ref. [1] we have fully
integrated three of them. The fourth contribution can be written in the
following way

$\displaystyle S_{d}=\int{\rm
d}\Phi^{nnn}_{\theta\theta\theta}\,\omega^{(3),d}_{n\bar{n}}(k_{1},k_{2},k_{3})\,,$
(1)

where $k_{1,2,3}$ are four-momenta of final-state gluons,
$\Phi^{nnn}_{\theta\theta\theta}$ is the phase space subject to zero-jettiness
conditions, cf. Eq. (3), and $\omega^{(3),d}_{n\bar{n}}(k_{1},k_{2},k_{3})$ is
the eikonal function defined in Eq. (5.49) in Ref. [1]. We note that the
integral in Eq. (1) is the the most complicated one, as it contains a
propagator that depends on the relative orientation of four-momenta of all
three soft gluons.

Although we described a possible way to calculate this contribution in Ref.
[1], we did not complete its computation there. The goal of this paper is to
compute the missing piece and to present the result for the same-hemisphere
three-gluon-emission contribution to the zero-jettiness soft function.

The rest of the paper is organized as follows. In Section II we explain the
general strategy that we used to integrate the function
$\omega_{n\bar{n}}^{3,(d)}$. In Section III we describe the computation of the
boundary condition for one of the master integrals that provides a
contribution that has unusual sensitivity to an analytic regulator. In Section
IV we discuss checks that we performed to ensure the correctness of the
computation. In Section V we present the result for the same-hemisphere three-
gluon-emission contribution to the zero-jettiness soft function and conclude
in Section VI.

##  II Integrating $\omega_{n\bar{n}}^{3,(d)}$

The missing part of the same-hemisphere contribution to the zero-jettiness
soft function, displayed in Eq. (1), requires integration of the function
$\omega_{n\bar{n}}^{3,(d)}$, which contains terms with a propagator
$1/k_{123}^{2}$, where $k_{123}=k_{1}+k_{2}+k_{3}$.

To integrate the function $\omega_{n\bar{n}}^{3,(d)}$, we consider a class of
integrals

$\displaystyle I_{\theta\theta\theta}$ $\displaystyle=\int\;{\rm
d}\Phi^{nnn}_{\theta\theta\theta}\;\frac{(k_{1}n)^{\nu}(k_{2}n)^{\nu}(k_{3}n)^{\nu}}{k_{123}^{2}\;(k_{1}k_{2})(k_{1}k_{3})\cdots},$
(2)

where $n$ and $\bar{n}$ are two light-like four-vectors pointing in the
direction of incoming partons, ellipses stand for eikonal propagators,111These
are all possible scalar products $q\\!\cdot\\!v$, where $q$ is a linear
combination of four-vectors of soft gluons and $v$ is one of the two light-
cone vectors $n$ or $\bar{n}$. and ${\rm d}\Phi^{nnn}_{\theta\theta\theta}$ is
the normalized phase-space measure. It is defined as follows [1]

${\rm
d}\Phi^{nnn}_{f_{1}f_{2}f_{3}}=N_{\epsilon}^{-3}\left(\prod_{i=1}^{3}[{\rm
d}k_{i}]f_{i}(k_{i})\right)\delta\left(1-\sum_{i=1}^{3}k_{i}n\right),$ (3)

where the normalization factor is

$N_{\epsilon}=\frac{\Omega^{(d-2)}}{4(2\pi)^{d-1}}=\frac{(4\pi)^{\epsilon}}{16\pi^{2}\Gamma(1-\epsilon)},$
(4)

the function $f_{i}$ is either $\theta(k_{i}\bar{n}-k_{i}n)$ or
$\delta(k_{i}\bar{n}-k_{i}n)$, and

$[{\rm d}k_{i}]=\frac{{\rm
d}^{d}k}{(2\pi)^{d}}2\pi\delta(k_{i}^{2})\theta(k_{i}^{0}).$ (5)

We note that we set the jettiness variable $\tau$ to one since the final
result is a uniform function of $\tau$. We also note that in the spirit of
integration-by-parts, all propagators that appear in the integrand in Eq. (2)
can be raised to arbitrary integer powers.

Finally, as can be seen from Eq. (2), we introduced scalar products of the
gluon four-momenta $k_{1,2,3}$ with the light-cone vector $n$ raised to power
$\nu$ into the integrand of $I_{\theta\theta\theta}$. As we explained in Ref.
[1] they are required because some integrals, that appear in the course of the
IBP reduction, contain divergences that are not regulated dimensionally; the
analytic regulator $\nu$ is introduced to regulate them. To simplify the
notation, we define the regulated phase-space measure to be

${\rm d}\Phi^{\nu}_{f_{1}f_{2}f_{3}}={\rm
d}\Phi^{nnn}_{f_{1}f_{2}f_{3}}(k_{1}n)^{\nu}(k_{2}n)^{\nu}(k_{3}n)^{\nu}.$ (6)

Unfortunately, even after integrals defined in Eq. (2) are reduced to master
integrals, the master integrals with $1/k_{123}^{2}$ propagators appear to be
too complicated for direct analytic integration. For this reason, as explained
in Ref. [1], we derive differential equations satisfied by these integrals,
and solve them numerically. To do that, we introduce a mass parameter into the
propagator that contains the momenta of all three gluons222Similar ideas that
use auxiliary mass-scales have been presented, for example, in Refs. [5, 6, 7,
8].

$\frac{1}{k_{123}^{2}}\to\frac{1}{k_{123}^{2}+m^{2}}.$ (7)

The appearance of the mass parameter $m$ allows us to differentiate with
respect to it and use integration-by-parts to derive differential equations
for relevant integrals. We then fix boundary conditions at $m\to\infty$, solve
differential equations numerically and determine relevant integrals at $m=0$.
This can be done by matching the numerical solution with the formal solution
at $m=0$ and then taking the $m\to 0$ limit in an appropriate manner. Since,
as we explained in Ref. [1], the differential equations can be solved to an
arbitrary precision as a matter of principle, and to more than 2000 digits in
practice, we have used the high-precision numerical results for master
integrals to find the analytic form of the solution by fitting them to a basis
of transcendental constants and rational numbers.

To reiterate, once it is understood how to use generalized unitarity to write
down integration-by-parts identities for integrals with $\theta$-functions, it
becomes a fairly standard problem to derive differential equations for
relevant integrals. However, great care is needed when choosing the basis of
master integrals because of the analytic regulator $\nu$; in essence, we need
to find a basis that admits a simple $\nu\to 0$ limit. Ideally, this should
happen for master integrals that appear in the soft function $S_{d}$, as well
as in $m$-dependent differential equations that we need to solve.

To find a suitable basis, we use the following consideration. In spite of the
fact that the reduction to master integrals performed after setting $\nu=0$ in
Eq. (2) is incorrect, it gives us a good idea about integrals that are
independent of each other in the $\nu\to 0$ limit. Hence, to find a suitable
basis for master integrals, we start by performing a reduction of integrals
shown in Eq. (2) at $\nu=0$. We then insist that the master integrals found in
the course of such a reduction should be also chosen as master integrals for
the reduction of integrals with $\nu\neq 0$, to an extent possible. We do this
for integrals with and without auxiliary mass parameter $m$. With this choice
of master integrals we find that the integral of the soft function
$\omega_{n\bar{n}}^{3,(d)}$ written in terms of master integrals, as well as
the differential equations that these integrals obey, admit a simple $\nu\to
0$ limit.

For example the result of the reduction of $S_{d}$ to master integrals can be
written in the following way

$S_{d}=\sum\limits_{\alpha}c_{\alpha}(\nu)I_{\alpha}^{\nu}+\nu\sum\limits_{\alpha}\tilde{c}_{\alpha}(\nu){\bar{I}}_{\alpha}^{\nu}\,,$
(8)

where the coefficients $c_{\alpha}$ and $\tilde{c}_{\alpha}$ are regular in
the $\nu\to 0$ limit. The list of integrals $\\{I_{\alpha}^{\nu}\\}$ coincides
with the list of master integrals for $S_{d}$ that one obtains performing a
reduction at $\nu=0$. New integrals that appear in the reduction, which we
denote as ${\bar{I}}_{\alpha}^{\nu}$ in the above equation, are multiplied
with the parameter $\nu$ and, therefore, disappear if the naive $\nu\to 0$
limit is taken.

However, the naive $\nu\to 0$ limit in Eq. (8) cannot be taken because some of
${\bar{I}}_{\alpha}^{\nu}$ integrals are $1/\nu$ divergent and, therefore,
need to be retained. Examples of such integrals without the $1/k_{123}^{2}$
propagator can be found in Ref. [1].

Unfortunately, many integrals among $I_{\alpha}^{\nu}$ and
${\bar{I}}_{\alpha}^{\nu}$ contain the $1/k_{123}^{2}$ propagator; to study
these integrals, we use Eq. (7) and turn them into integrals with the mass
parameter $m$. We construct a system of differential equations w.r.t. the mass
parameter $m$ by including every integral with the $1/k^{2}_{123}$ propagator
that appears in Eq. (8) and other integrals that are needed to close it. Such
a system takes the following form

$\begin{split}&\frac{\partial}{\partial
m^{2}}\boldsymbol{J}^{\nu}=\boldsymbol{M}_{1}(\nu)\boldsymbol{J}^{\nu}+\nu\boldsymbol{N}_{1}(\nu)\boldsymbol{\bar{J}}^{\nu},\\\
&\frac{\partial}{\partial
m^{2}}\boldsymbol{\bar{J}}^{\nu}=\boldsymbol{M}_{2}(\nu)\boldsymbol{\bar{J}}^{\nu}+\boldsymbol{N}_{2}(\nu)\boldsymbol{J}^{\nu}.\end{split}$
(9)

In Eq. (9) matrices $\boldsymbol{M}_{1,2}$ and $\boldsymbol{N}_{1,2}$ are
regular in the $\nu\to 0$ limit, $\boldsymbol{J}^{\nu}$ integrals are the
master integrals that need to be considered for computing
$\\{I_{\alpha}^{\nu}\\}$ integrals and $\boldsymbol{\bar{J}}^{\nu}$ are
integrals that are needed for computing $\\{{\bar{I}_{\alpha}}^{\nu}\\}$
integrals.

From the structure of differential equations in Eq. (9) it is clear that
taking the limit $\nu\to 0$ at the early stages of the computation is
beneficial since this makes the system of differential equations significantly
simpler. However, this is only possible if we know which master integrals are
singular in the $\nu\to 0$ limit and which master integrals are not.
Unfortunately, it is not trivial to answer this question and we use several
approaches to clarify it.

First, from the computation of integrals of the eikonal functions
$\omega^{a,b,c}$ reported in Ref. [1], we know $1/\nu$-divergent integrals
that appear in cases when the propagator $1/k_{123}^{2}$ is absent. Upon
inspection, we find that the only $m$-independent $1/\nu$-divergent integral
that appears in the amplitude $S_{d}$ and in the differential equations reads

${\bar{I}}_{1}^{\nu}={\bar{J}}_{1}^{\nu}=\int\frac{\mathrm{d}\Phi_{\theta\delta\theta}^{\nu}}{(k_{1}\cdot
k_{3})(k_{1}\cdot n)(k_{12}\cdot\overline{n})(k_{3}\cdot\overline{n})}\,.$
(10)

This integral was computed in Ref. [1] and for this reason we do not discuss
it here.

Second, to determine which of the integrals with the propagator
$1/k_{123}^{2}$ are singular in the $\nu\to 0$ limit, we can study these
integrals at finite values of $m$ and then employ differential equations to
determine the $1/\nu$ behavior of the corresponding ${\bar{I}}_{\alpha}^{\nu}$
integrals. To this end, we employed Mellin-Barnes representation of the
relevant integrals and used public programs MB [9] and MBresolve [10] to
numerically compute all $m$-dependent integrals that appear in the
differential equations at finite values of $m$. We also used the program
pySecDec [11, 12] for cross-checks of the numerical computation. Upon doing
that, we discovered yet another integral that is singular in the $\nu\to 0$
limit. It reads

$\displaystyle{\bar{J}}_{2}^{\nu}=\int\frac{\mathrm{d}\Phi_{\theta\delta\theta}^{\nu}}{(k_{123}^{2}+m^{2})(k_{1}\cdot
k_{3})(k_{1}\cdot n)(k_{12}\cdot\overline{n})}\,.$ (11)

This integral is quite peculiar and we explain in the next section why this is
the case and how to compute it.

To proceed further, we re-scale $\boldsymbol{\bar{J}}^{\nu}$ integrals that
appear in Eq. (9) by a factor $\nu$. Since, as we just explained, only two
integrals ${\bar{J}}_{1}^{\nu}$ and ${\bar{J}}_{2}^{\nu}$ diverge in the
$\nu\to 0$ limit, we need to keep them in
$\mathop{}\\!\mathcal{O}\\!\left(\nu^{0}\right)$ part of the differential
equations. Therefore, we combine $\boldsymbol{J}^{\nu}$ integrals together
with $\nu{\bar{J}}_{1}^{\nu}$ and $\nu{\bar{J}}_{2}^{\nu}$ integrals into a
vector $\boldsymbol{\mathcal{J}}^{\nu}$ and use remaining ${\bar{J}}^{\nu}$
integrals, rescaled by a parameter $\nu$, to define a new vector
$\boldsymbol{\mathcal{\bar{J}}}^{\nu}$

$\displaystyle\begin{split}\boldsymbol{\mathcal{J}}^{\nu}&=\left(\boldsymbol{J}^{\nu},\nu\bar{J}_{1}^{\nu},\nu\bar{J}_{2}^{\nu}\right),\\\
\boldsymbol{\mathcal{\bar{J}}}^{\nu}&=\left(\nu\bar{J}_{3}^{\nu},\cdots\right).\end{split}$
(12)

We then obtain a new system of differential equations

$\displaystyle\begin{split}&\frac{\partial}{\partial
m^{2}}\boldsymbol{\mathcal{J}}^{\nu}=\boldsymbol{\cal
M}_{1}(\nu)\boldsymbol{\cal J}^{\nu}+\boldsymbol{{\cal
N}}_{1}(\nu)\boldsymbol{\bar{\cal J}}^{\nu},\\\ &\frac{\partial}{\partial
m^{2}}\boldsymbol{\bar{\cal J}}^{\nu}=\boldsymbol{{\cal
M}}_{2}(\nu)\boldsymbol{\bar{\cal J}}^{\nu}+\nu\boldsymbol{\cal
N}_{2}(\nu)\boldsymbol{{\cal J}}^{\nu}.\end{split}$ (13)

It is straightforward to solve the above equation expanding in $\nu\to 0$
because all the matrices that appear there have smooth $\nu\to 0$ limits and
because the integral vectors satisfy
$\boldsymbol{\mathcal{J}}^{\nu}\sim\mathop{}\\!\mathcal{O}\\!\left(\nu^{0}\right)$,
$\boldsymbol{\mathcal{\bar{J}}}^{\nu}\sim\mathop{}\\!\mathcal{O}\\!\left(\nu\right)$.

Working to order $\mathop{}\\!\mathcal{O}\\!\left(\nu^{0}\right)$, we can drop
$\boldsymbol{\mathcal{\bar{J}}}^{\nu}$ integrals and set $\nu$ to zero in Eq.
(13). This leads to a significant reduction in the number of integrals that
appear in Eq. (13) and allows for solving the system of differential equations
in a more efficient way.

As follows from the differential equation satisfied by
$\nu{\bar{J}}_{1}^{\nu}$ and $\nu{\bar{J}}_{2}^{\nu}$, these integrals are
independent of $m$ through $\mathop{}\\!\mathcal{O}\\!\left(\nu^{0}\right)$.
In addition, they also have to obey the following relation

$\lim_{\nu\to
0}\nu{\bar{J}}_{1}^{\nu}=\frac{1+6\epsilon}{1+4\epsilon}\;\lim_{\nu\to
0}\nu{\bar{J}}_{2}^{\nu},$ (14)

to cancel the $1/\nu$ singularity which is naively present in the rescaled
differential equation. We will show in the following section that this
condition is indeed satisfied.

We note that the above discussion applies to differential equations at
_finite_ values of $m$ whereas, eventually, we are interested in the solutions
at $m=0$. This limit is, potentially, non-trivial. Indeed, to find required
values of integrals at $m=0$, we need to compute master integrals using the
following sequence of limits: $m\to 0,\nu\to 0,\epsilon\to 0$.

However, as explained above, we would like to simplify differential equations
by taking the $\nu\to 0$ limit first and there are two problems that may arise
if the order of limits is changed. First, since the mass parameter $m$ can
_also_ serve as a regulator of collinear and soft singularities, $m\to 0$ and
$\nu\to 0$ limits should not necessarily commute. Second, additional
contributions can mix into the Taylor $m^{0}$-branch of the solution that we
require, if the $\nu\to 0$ limit is taken before the $m\to 0$ limit. We will
discuss these two problems now.

Suppose we take the $\nu\to 0$ limit for finite $m$-integrals but the
resulting integrals are still not regulated dimensionally at $m=0$. This
feature can be detected in the following way. The dependence of the integral
on $m$ and $\epsilon$ at small $m$ has the following form

$\mathcal{J}\sim\sum\limits_{n_{1},n_{2},n_{3}}c_{n_{1}n_{2}n_{3}}m^{n_{1}+n_{2}\epsilon}\ln^{n_{3}}(m).$
(15)

We are interested in taking the $m\to 0$ limit of this solution at fixed
$\epsilon$. However, this is only possible if the coefficients
$c_{n_{1}0n_{3}}$ with $n_{1}<0$ and $n_{1}=0,n_{3}>0$ vanish so that there
are no $1/m$ and $\log m$ terms that are not multiplied by additional powers
of $m^{\epsilon}$ or sufficiently high powers of $m$. We have checked that
this condition is satisfied for all $\nu=0$ integrals that we considered; this
implies that the $m\to 0$ limit does not lead to divergencies in integrals
that are not regulated dimensionally. It follows that indeed
${\bar{I}}_{1,2}^{\nu}$ are the only integrals in ${\bar{I}}_{\alpha}^{\nu}$
that contribute to the amplitude, all other integrals can be safely discarded.

The second problem concerns possible mixing between different branches of
integrals if the $\nu\to 0$ limit is taken too early. To see how this comes
about, consider a general solution in the limit $m\to 0,\nu\to 0,\epsilon\to
0$

$\mathcal{J}\sim\sum\limits_{n_{1},n_{2},n_{3},n_{4}}c_{n_{1}n_{2}n_{3}n_{4}}m^{n_{1}+n_{2}\epsilon+n_{4}\nu}\ln^{n_{3}}(m).$
(16)

If there are terms that correspond to $n_{1}=0,n_{2}=0,n_{3}=0$, $n_{4}\neq
0$, they will mix with the contribution $n_{1}=0,n_{2}=0,n_{3}=0,n_{4}=0$,
i.e. the $m^{0}$-branch that we are interested in.

Hence, we need to understand if such contributions exist and, if they are
there, how to isolate and remove them. This can be done by studying exact
differential equation Eq. (13) at small values of $m$ but with full $\nu$\-
and $\epsilon$-dependence, and checking if $m^{n_{4}\nu}$ solutions without
additional dependencies of exponents on $\epsilon$ and additional powers of
$m$ are possible. We find that this does not happen at $m=0$ and, therefore,
the Taylor branch does not receive any unwanted contributions.

In summary, we can solve the differential equations Eq. (13) as an expansion
in $\nu$. As discussed in Ref. [1], we can compute the boundary conditions at
$m=\infty$ and then find values of integrals $\boldsymbol{\mathcal{J}}$ at
$m=0$ by discarding all terms that have non-analytic dependencies on $m$ at
$m=0$. We present the results of such a computation in Section V; in the next
section we describe computation of a peculiar boundary condition which
illustrates that our worry about potential mixing of a Taylor- and
$m^{n_{4}\nu}$-branches is not unfounded.

##  III Integral ${\bar{J}}_{2}^{\nu}$ with $1/\nu$ divergence

We can illustrate some points discussed in the previous section by considering
the integral ${\bar{J}}_{2}^{\nu}$ and its contribution to differential
equations. As mentioned there, this integral is singular in the limit $\nu\to
0$. Multiplied by a factor $\nu$, it appears in the differential equations for
two sets of $\nu$-regular integrals. We will consider one of them for the sake
of example

$\displaystyle J_{a_{1}}^{\nu}$
$\displaystyle=\int\frac{\mathrm{d}\Phi_{\theta\delta\theta}^{\nu}}{(k_{123}^{2}+m^{2})(k_{1}\cdot
k_{3})(k_{23}\cdot\bar{n})}\,,$ $\displaystyle J_{a_{2}}^{\nu}$
$\displaystyle=\int\frac{\mathrm{d}\Phi_{\theta\delta\theta}^{\nu}}{(k_{123}^{2}+m^{2})(k_{1}\cdot
k_{3})(k_{23}\cdot\bar{n})^{2}}\,.$ (17)

To compute these integrals from differential equations, we require a boundary
condition for $\nu{\bar{J}}_{2}^{\nu}$, which we compute at $m=\infty$. It
follows from the analysis of the differential equation for
$\nu{\bar{J}}_{2}^{\nu}$ at $m=\infty$ that the required boundary condition
can be determined from the leading, mass-independent term in the
$1/m$-expanison. Upon inspecting the various contributions to the asymptotic
behaviour of the integral ${\bar{J}}_{2}^{\nu}$ discussed in Ref. [1], we find
that they do not produce mass-independent terms that are $1/\nu$-divergent.

It turns out that the integral ${\bar{J}}_{2}^{\nu}$ provides an example of a
situation where the analysis of different regions that contribute to
$m\to\infty$ asymptotic behaviour of integrals, performed in Ref. [1], is
incomplete and that there is another region that needs to be considered. In
fact we have found that the following scaling of integration variables

$\displaystyle k_{3}\cdot\bar{n}=\alpha_{3}\sim m^{2}\gg 1,$ $\displaystyle
k_{1}\cdot\bar{n}=\alpha_{1}\sim 1,$ (18) $\displaystyle k_{1}\cdot
n=\beta_{1}\sim m^{-2}\ll 1,$

leads to a $1/\nu$-divergent $\mathop{}\\!\mathcal{O}\\!\left(1/m^{0}\right)$
contribution to ${\bar{J}}_{2}^{\nu}$.

To show this, we write an approximation to the integrand of
${\bar{J}}_{2}^{\nu}$-integral in the region defined by Eq. (18)

$\displaystyle{\bar{J}}_{2}^{\nu}\sim\widetilde{J_{2}^{\nu}}=\int\\!\frac{\mathrm{d}\Phi_{\theta\delta\theta}^{\nu}}{((k_{3}\\!\cdot\\!\bar{n})\,(k_{2}\\!\cdot\\!n)+m^{2})(k_{1}\\!\cdot\\!k_{3})(k_{1}\\!\cdot\\!n)(k_{12}\\!\cdot\\!\bar{n})}\,.$
(19)

We use the Sudakov variables $\alpha_{i}$ and $\beta_{i}$ (see Ref. [1] for
details) and take into account the condition $\beta_{1}\ll 1$ to remove
$\beta_{1}$ from the “jettiness” delta-function
$\delta(1-\beta_{1}-\beta_{2}-\beta_{3})\to\delta(1-\beta_{2}-\beta_{3})$. We
then extend the integration over $\beta_{1}$ to infinity. We find

$\displaystyle\widetilde{J_{2}^{\nu}}$
$\displaystyle=2\int\limits_{0}^{\infty}{\rm
d}\beta_{1}\int\mathrm{d}\beta_{2}\mathrm{d}\beta_{3}\;\mathrm{d}\alpha_{1}\mathrm{d}\alpha_{3}\;\beta_{1}^{-\epsilon+\nu}\beta_{2}^{-2\epsilon+\nu}\beta_{3}^{-\epsilon+\nu}\alpha_{1}^{-\epsilon}\alpha_{3}^{-\epsilon}\;\frac{\delta(1-\beta_{23})\theta(\alpha_{1}-\beta_{1})\theta(\alpha_{3}-\beta_{3})}{(\alpha_{3}\beta_{2}+m^{2})\beta_{1}(\alpha_{1}+\beta_{2})}$
$\displaystyle\phantom{={}}\times\Bigg{[}\frac{\theta(\beta_{1}/\alpha_{1}-\beta_{3}/\alpha_{3})}{\beta_{1}\alpha_{3}}{}_{2}F_{1}\left(1,1+\epsilon;1-\epsilon;\frac{\alpha_{1}\beta_{3}}{\alpha_{3}\beta_{1}}\right)+\frac{\theta(\beta_{3}/\alpha_{3}-\beta_{1}/\alpha_{1})}{\beta_{3}\alpha_{1}}{}_{2}F_{1}\left(1,1+\epsilon;1-\epsilon;\frac{\alpha_{3}\beta_{1}}{\alpha_{1}\beta_{3}}\right)\Bigg{]}\,.$
(20)

We change integration variables $\alpha_{1}=\beta_{1}/\xi_{1}$ and
${\alpha_{3}=\beta_{3}/\xi_{3}}$ and obtain

$\displaystyle\widetilde{J_{2}^{\nu}}$
$\displaystyle=2\int\mathrm{d}\beta_{2}\mathrm{d}\beta_{3}\mathrm{d}\xi_{1}\mathrm{d}\xi_{3}\beta_{2}^{-2\epsilon+\nu}\beta_{3}^{-2\epsilon+\nu}\xi_{1}^{\epsilon-1}\xi_{3}^{\epsilon-1}\,\frac{\delta(1-\beta_{23})\theta(1-\xi_{1})\theta(1-\xi_{3})}{\beta_{3}\beta_{2}+m^{2}\xi_{3}}$
$\displaystyle\phantom{={}}\times\left(\xi_{3}\theta(\xi_{1}-\xi_{3})\mathop{}\\!_{2}F_{1}\\!\left[\\{1,1+\epsilon\\},\\{1-\epsilon\\};\frac{\xi_{3}}{\xi_{1}}\right]+\xi_{1}\theta(\xi_{3}-\xi_{1})\mathop{}\\!_{2}F_{1}\\!\left[\\{1,1+\epsilon\\},\\{1-\epsilon\\};\frac{\xi_{1}}{\xi_{3}}\right]\right)\;$
(21)
$\displaystyle\phantom{={}}\times\int\limits_{0}^{\infty}\mathrm{d}\beta_{1}\frac{\beta_{1}^{-2\epsilon+\nu-1}}{\beta_{1}+\beta_{2}\xi_{1}}\,.$

Integrating over $\beta_{1}$, we find

$\displaystyle\int\limits_{0}^{\infty}\mathrm{d}\beta_{1}\frac{\beta_{1}^{-2\epsilon+\nu-1}}{\beta_{1}+\beta_{2}\xi_{1}}$
$\displaystyle=(\beta_{2}\xi_{1})^{-2\epsilon+\nu-1}\Gamma(-2\epsilon+\nu)$
$\displaystyle\phantom{={}}\times\Gamma(2\epsilon-\nu+1).$ (22)

We use this result in Eq. (21), change variables $\xi_{1}=r\xi_{3}$, and
arrive at

$\displaystyle\widetilde{J_{2}^{\nu}}$
$\displaystyle=2\Gamma(-2\epsilon+\nu)\Gamma(2\epsilon-\nu+1)$
$\displaystyle\phantom{={}}\times\int\mathrm{d}\beta_{2}\mathrm{d}\beta_{3}\beta_{2}^{-4\epsilon+2\nu-1}\beta_{3}^{-2\epsilon+\nu}\delta(1-\beta_{23})$
(23)
$\displaystyle\phantom{={}}\times\int_{0}^{1}\frac{\mathrm{d}\xi_{3}\xi_{3}^{\nu-1}J(\nu,\xi_{3})}{\beta_{3}\beta_{2}+m^{2}\xi_{3}}\,,$

where

$\displaystyle J(\nu,\xi_{3})$
$\displaystyle=\int_{0}^{1}\mathrm{d}r\left[\theta(1-\xi_{3}/r)r^{\epsilon-\nu}+r^{-\epsilon+\nu-1}\right]$
$\displaystyle\phantom{={}}\times\mathop{}\\!_{2}F_{1}\\!\left[\\{1,1+\epsilon\\},\\{1-\epsilon\\};r\right]\,.$
(24)

It follows from the above expression that the $\nu$-pole originates from a
singularity at $\xi_{3}=0$. It is straightforward to compute it since
$J(\nu,\xi_{3})$ is regular at $\xi_{3}=0$.333As with any analytic regulator,
the $\nu\to 0$ limit should be computed keeping $\epsilon$ fixed. We find

$\displaystyle\int_{0}^{1}\mathrm{d}\xi_{3}\frac{\xi_{3}^{\nu-1}J(\nu,\xi_{3})}{\beta_{3}\beta_{2}+m^{2}\xi_{3}}=\frac{J(0,0)}{\nu}(\beta_{3}\beta_{2})^{-1}+\mathcal{O}(\nu^{0})\,.$
(25)

Upon further integration, we obtain a $1/\nu$-divergent contribution to
$\tilde{J}_{2}^{\nu}$. It reads

$\displaystyle\widetilde{J_{2}^{\nu}}=\frac{C_{2}}{\nu},$ (26)

where

$\displaystyle C_{2}={}$
$\displaystyle\frac{2\Gamma^{2}(-2\epsilon)\Gamma(-4\epsilon-1)\Gamma(2\epsilon+1)}{\Gamma(-6\epsilon-1)}$
$\displaystyle\times\left(\frac{\mathop{}\\!_{3}F_{2}\\!\left[\\{1,1+\epsilon,1+\epsilon\\},\\{1-\epsilon,2+\epsilon\\};1\right]}{1+\epsilon}-\frac{\mathop{}\\!_{3}F_{2}\\!\left[\\{1,-\epsilon,1+\epsilon\\},\\{1-\epsilon,1-\epsilon\\};1\right]}{\epsilon}\right).$
(27)

Using the result of the explicit computation of the $1/\nu$ pole of
${\bar{J}}_{2}^{\nu}$ and the result for ${\bar{J}}_{1}^{\nu}$ reported in
Ref. [1], we find that they satisfy the relation shown Eq. (14). As we pointed
out earlier, this relation is needed to ensure the smooth $\nu\to 0$ limit of
the differential equations.

The above result provides the required boundary condition at $m=\infty$ and
allows us to start solving differential equations numerically. However, it is
interesting to point out that, from the perspective of the differential
equations, integral ${\bar{J}}_{2}^{\nu}$ provides an example of a
contribution proportional to $m^{-2\nu}$ which, therefore, can mix with the
Taylor branch of the required integrals if the $\nu\to 0$ limit is taken
first.

Indeed, if we first apply the scaling defined in Eq. (18) to the integration
variables in Eq. (20), we find that this integration region leads to an
overall factor $(1/m)^{2\nu}$

$\displaystyle\nu{\bar{J}}_{2}^{\nu}\sim\left(1/m\right)^{2\nu}\left[C_{2}+\mathcal{O}(\nu^{1})\right]\,.$
(28)

We note that in the limit $\nu\to 0$, this region looks like a “normal” Taylor
region $1/m^{0}$

$\displaystyle\nu{\bar{J}}_{2}^{\nu}\sim\left(1/m\right)^{0}C_{2}+\mathcal{O}(\nu^{1})\,,$
(29)

and, once the $\nu\to 0$ limit is taken, the two regions cannot be
distinguished. This is an illustration of the problem of mixing between
different $m$-branches that we discussed in the previous section; however, in
contrast to the discussion there, our example in this section refers to
$m\to\infty$ limit.

The existence of $(1/m)^{2\nu}$ branch in the integral $\nu J_{2}^{\nu}$ leads
to the appearance of similar contributions in other integrals that $\nu
J_{2}^{\nu}$ couples to through differential equations. For the sake of
example, we consider the differential equation for the so-called $J_{a}^{\nu}$
sector, cf. Eq. (17). The differential equations take the following form

$\displaystyle\frac{\partial}{\partial m^{2}}\begin{pmatrix}J_{a_{1}}^{\nu}\\\
J_{a_{2}}^{\nu}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}-\frac{2\epsilon}{m^{2}}-\frac{2(1+2\epsilon)}{4m^{2}+1}&\frac{1}{4m^{2}+1}\\\
\frac{2\epsilon(1+2\epsilon)}{m^{2}}-\frac{8\epsilon(1+2\epsilon)}{4m^{2}+1}&\frac{1+4\epsilon}{m^{2}}-\frac{28\epsilon}{4m^{2}+1}\end{pmatrix}\begin{pmatrix}J_{a_{1}}^{\nu}\\\
J_{a_{2}}^{\nu}\end{pmatrix}-\begin{pmatrix}\frac{(1+4\epsilon)}{(1+6\epsilon)}\left(\frac{1}{m^{2}}-\frac{4}{4m^{2}+1}\right)\\\\[6.00006pt]
\frac{(1+4\epsilon)}{(1+6\epsilon)}\left(\frac{4\epsilon}{m^{2}}-\frac{16\epsilon}{4m^{2}+1}\right)\end{pmatrix}\nu{\bar{J}}_{2}^{\nu}$
(30) $\displaystyle\phantom{={}}+\text{contributions from other $J^{\nu}$
integrals.}$

It follows from the above equation that the integral $\nu{\bar{J}}_{2}^{\nu}$
plays a role of an inhomogeneous contribution to the differential equation
that $J_{a_{1}}^{\nu}$ and $J_{a_{2}}^{\nu}$ satisfy. In fact, analyzing the
homogeneous terms of the above equation in the $m\to\infty$ limit, we find
that we do not need to compute the boundary conditions for these integrals and
that the solution of the differential equation in this limit are obtained by
integrating the inhomogeneous part. We obtain

$\displaystyle J_{a_{1}}^{\nu}$
$\displaystyle=m^{-2}\left(\frac{4\epsilon+1}{4(6\epsilon+1)}C_{2}\right)+\dots,$
(31) $\displaystyle J_{a_{2}}^{\nu}$
$\displaystyle=m^{-2}\left(\frac{\epsilon(4\epsilon+1)}{6\epsilon+1}C_{2}\right)+\dots,$
(32)

where $\nu\to 0$ limit has already been taken where appropriate and dots stand
for other contributions, including homogeneous and inhomogeneous ones.

We emphasize one more time that from the point of view of the differential
equation, the contribution shown in the above equation scales as
$\sim(1/m)^{2\nu}$ but, once the $\nu\to 0$ limit is taken it becomes
indistinguishable from a regular Taylor part of the integral. If a similar
situation occurs at $m=0$, we should have identified and removed all the
$m^{n_{4}\nu}$ regions in all integrals since the correct sequence of limits
that is needed is $m\to 0,\nu\to 0,{\epsilon\to 0}$. However, as we mentioned
earlier, an analysis of the differential equation at $m\to 0$ leads to the
conclusion that there are no eigenvalues that vanish if $\nu\to 0$ limit is
taken so that the problem described above does not occur.

##  IV Checks

Given the highly unusual nature of the integrals that need to be computed to
obtain the zero-jettiness soft function and the complex interplay of the
various infra-red regulators it is important to perform as many checks as
possible to ensure correctness of the result.

The most comprehensive check that can be performed is the numerical
computation of all $m$-dependent integrals which appear in the differential
equations, as well as their derivatives. We constructed Mellin-Barnes
representation of the relevant integrals using public programs MB [9] and
MBresolve [10] for this purpose. We also used the program pySecDec [11, 12] as
an alternative for the numerical computation. Using these programs, we have
computed all integrals that appeared in the differential equations at finite
values of $m$ and checked them against numerical solutions of the differential
equations.

Next, we compared the solutions of the differential equations at $m=0$ and
$\nu=0$ with the results of the direct numerical computation. Unfortunately,
although this can be done for some integrals that contribute to $S_{d}$, there
are many integrals for which the numerical integration becomes next to
impossible. To enlarge the set of integrals at $m=0$ that can be checked, we
have derived linear relations between various integrals at $m=0$ using the
integration-by-parts identities and checked that integrals obtained from
$m$-dependent differential equations and extrapolated to $m=0$, satisfy them.

##  V Results

Solving differential equations and separating the Taylor branch at $m=0$, we
obtain numerical result for the integral of the function
$\omega_{n\bar{n}}^{3,(d)}$. Since we can determine the solution of the
differential equation to, essentially, arbitrary precision, we can try to
obtain the analytic result for $S_{d}$ by fitting the numerical results to a
linear combination of various transcendental and rational numbers. By making
use of the PSLQ [13] and LLL [14] algorithms, and choosing appropriate basis
of transcendental numbers [15], we find the following result

$\displaystyle S_{d}=\int{\rm
d}\Phi^{nnn}_{\theta\theta\theta}\,\omega^{(3),d}_{n\bar{n}}(k_{1},k_{2},k_{3})=\frac{12}{\epsilon^{5}}+\frac{142}{3\epsilon^{4}}+\frac{1}{\epsilon^{3}}\left(\frac{46\pi^{2}}{3}+\frac{628}{3}\right)+\frac{1}{\epsilon^{2}}\left(196\zeta_{3}+\frac{650\pi^{2}}{9}+\frac{18161}{27}\right)$
$\displaystyle+\frac{1}{\epsilon}\left(\frac{397\pi^{4}}{45}+1380\zeta_{3}+\frac{6808\pi^{2}}{27}+\frac{165323}{81}\right)+\left(8982\zeta_{5}-\frac{2146\zeta_{3}\pi^{2}}{3}+\frac{191\pi^{4}}{9}+4224\mathop{}\\!{\mathrm{Li}}_{4}\\!\left(\frac{1}{2}\right)\right.$
$\displaystyle\left.{}+3696\zeta_{3}\mathop{}\\!{\mathrm{ln}}(2)-176\pi^{2}\mathop{}\\!{\mathrm{ln}}^{2}(2)+176\mathop{}\\!{\mathrm{ln}}^{4}(2)+\frac{46184\zeta_{3}}{9}+\frac{66614\pi^{2}}{81}+96\mathop{}\\!{\mathrm{ln}}(2)+\frac{413971}{81}\right)$
$\displaystyle+\epsilon\left(2304\mathop{}\zeta_{-5,-1}-4464\zeta_{5}\mathop{}\\!{\mathrm{ln}}(2)-8380\zeta_{3}^{2}+\frac{46934\pi^{6}}{2835}-6336G_{R}(0,0,r_{2},1,-1)\right.$
$\displaystyle\phantom{={}}-6336G_{R}(0,0,1,r_{2},-1)-3168G_{R}(0,0,1,r_{2},r_{4})-6336G_{R}(0,0,r_{2},-1)\mathop{}\\!{\mathrm{ln}}(2)+\frac{324215\zeta_{5}}{3}$
(33)
$\displaystyle\phantom{={}}-45056\mathop{}\\!{\mathrm{Li}}_{5}\\!\left(\frac{1}{2}\right)-45056\mathop{}\\!{\mathrm{Li}}_{4}\\!\left(\frac{1}{2}\right)\mathop{}\\!{\mathrm{ln}}(2)+176\mathop{}\\!{\mathrm{Cl}}_{4}\\!\left(\frac{\pi}{3}\right)\pi-1056\zeta_{3}\mathop{}\\!{\mathrm{Li}}_{2}\\!\left(\frac{1}{4}\right)-\frac{9634\zeta_{3}\pi^{2}}{3}$
$\displaystyle\phantom{={}}-21824\zeta_{3}\mathop{}\\!{\mathrm{ln}}^{2}(2)+2112\zeta_{3}\mathop{}\\!{\mathrm{ln}}(2)\mathop{}\\!{\mathrm{ln}}(3)-1584\mathop{}\\!{\mathrm{Cl}}_{2}^{2}\\!\left(\frac{\pi}{3}\right)\mathop{}\\!{\mathrm{ln}}(3)-\frac{4400\mathop{}\\!{\mathrm{Cl}}_{2}\\!\left(\frac{\pi}{3}\right)\pi^{3}}{27}+\frac{88\pi^{4}\mathop{}\\!{\mathrm{ln}}(2)}{45}$
$\displaystyle\phantom{={}}-\frac{616\pi^{4}\mathop{}\\!{\mathrm{ln}}(3)}{27}+\frac{11264\pi^{2}\mathop{}\\!{\mathrm{ln}}^{3}(2)}{9}-\frac{22528\mathop{}\\!{\mathrm{ln}}^{5}(2)}{15}+8576\mathop{}\\!{\mathrm{Li}}_{4}\\!\left(\frac{1}{2}\right)+7504\zeta_{3}\mathop{}\\!{\mathrm{ln}}(2)+\frac{4646\pi^{4}}{27}$
$\displaystyle\phantom{={}}-\frac{1072\pi^{2}\mathop{}\\!{\mathrm{ln}}^{2}(2)}{3}+\frac{1072\mathop{}\\!{\mathrm{ln}}^{4}(2)}{3}+\frac{496592\zeta_{3}}{27}-32\pi^{2}\mathop{}\\!{\mathrm{ln}}(2)+\frac{587380\pi^{2}}{243}-384\mathop{}\\!{\mathrm{ln}}^{2}(2)+832\mathop{}\\!{\mathrm{ln}}(2)$
$\displaystyle\phantom{={}}+\frac{7857076}{729}+\left.\sqrt{3}\left(192\,\Im\Big{\\{}\\!\mathop{}\\!{\mathrm{Li}}_{3}\\!\left(\frac{\exp(\mathrm{i}\pi/3)}{2}\right)\\!\Big{\\}}+160\mathop{}\\!{\mathrm{Cl}}_{2}\\!\left(\frac{\pi}{3}\right)\mathop{}\\!{\mathrm{ln}}(2)-16\pi\mathop{}\\!{\mathrm{ln}}^{2}(2)-\frac{560\pi^{3}}{81}\right)\right)$
$\displaystyle\phantom{={}}+\mathop{}\\!\mathcal{O}\\!\left(\epsilon^{2}\right)\,,$

where $\zeta_{-5,-1}\approx-0.029902$ is a multiple zeta value, and
$\text{Cl}_{n}(x)$ are Clausen functions. $G_{R}(a_{1},\ldots,a_{w})$ is the
real part of the multiple polylogarithm $G(a_{1},\ldots,a_{w};z)$ evaluated at
$z=1$ [15]

$\displaystyle G_{R}(a_{1},\ldots,a_{w})=\Re\\{G(a_{1},\ldots,a_{w};1)\\}\,.$
(34)

Finally, $r_{2}=\exp(-\mathrm{i}\pi/3)$ and $r_{4}=\exp(-\mathrm{i}2\pi/3)$.
We note that we have computed the master integrals to more than two thousand
digits to check the validity of the analytic result.

Having obtained the result for the integral of the function
$\omega_{n\bar{n}}^{3,(d)}$ we are now in position to present the complete
result for the same-hemisphere _three-gluon_ -emission contribution to the
N3LO soft function. To this end, we write

$\displaystyle S^{nnn}$ $\displaystyle=\int{\rm
d}\Phi^{nnn}_{\theta\theta\theta}\left|\boldsymbol{J}(k_{1},k_{2},k_{3})\right|^{2}=\tau^{-1-6\epsilon}\frac{N_{\epsilon}^{3}}{3!}\left[C_{a}^{3}S^{nnn}_{1+1+1}+C_{a}^{2}C_{A}S^{nnn}_{1+2}+C_{a}C_{A}^{2}S^{nnn}_{3}\right],$
(35)

where we re-introduced the dependence on $\tau$, recovered the symmetry factor
$1/3!$ and the normalization factor $N_{\epsilon}^{3}$ (cf. Eq. (3)), and
split the integral into three color factors following Eq. (7.10) in Ref. [4].
We also note that $C_{a}=C_{F,A}$ for the quark (gluon) soft function,
respectively.

The computation of the Abelian contributions $S^{nnn}_{1+1+1}$ and
$S^{nnn}_{1+2}$ is described in Appendix A. We obtain the maximally non-
Abelian contribution $S^{nnn}_{3}$ by adding results obtained in Ref. [1] and
the result of this paper given in Eq. (33). We find

$\displaystyle S^{nnn}_{1+1+1}$
$\displaystyle=\frac{48\,\Gamma^{3}(1-2\epsilon)}{\epsilon^{5}\Gamma(1-6\epsilon)},$
(36) $\displaystyle S^{nnn}_{1+2}$
$\displaystyle=-\frac{9\,\Gamma(1-4\epsilon)\,\Gamma(1-2\epsilon)}{\epsilon^{2}\Gamma(1-6\epsilon)}\times\bigg{[}\frac{8}{\epsilon^{3}}+\frac{44}{3\epsilon^{2}}+\frac{1}{\epsilon}\left(\frac{268}{9}-8\zeta_{2}\right)+\left(\frac{1544}{27}+\frac{88}{3}\zeta_{2}-72\zeta_{3}\right)$
$\displaystyle\phantom{={}}+\epsilon\left(\frac{9568}{81}+\frac{536\zeta_{2}}{9}+\frac{352}{3}\zeta_{3}-300\zeta_{4}\right)+\epsilon^{2}\left(\frac{55424}{243}+\frac{3520\zeta_{2}}{27}+\frac{2144\zeta_{3}}{9}+352\zeta_{4}+96\zeta_{2}\zeta_{3}-1208\zeta_{5}\right)$
$\displaystyle\phantom{={}}+\epsilon^{3}\left(\frac{297472}{729}+\frac{22592\zeta_{2}}{81}+\frac{14080\zeta_{3}}{27}+\frac{2144}{3}\zeta_{4}-\frac{4576}{3}\zeta_{2}\zeta_{3}+3696\zeta_{5}+424\zeta_{3}^{2}-3596\zeta_{6}\right)+\mathop{}\\!\mathcal{O}\\!\left(\epsilon^{4}\right)\bigg{]},$
(37) $\displaystyle S^{nnn}_{3}$
$\displaystyle=\frac{24}{\epsilon^{5}}+\frac{308}{3\epsilon^{4}}+\frac{1}{\epsilon^{3}}\left(-12\pi^{2}+\frac{3380}{9}\right)+\frac{1}{\epsilon^{2}}\left(-1000\zeta_{3}+\frac{440\pi^{2}}{9}+\frac{10048}{9}\right)$
$\displaystyle\phantom{={}}+\frac{1}{\epsilon}\left(-\frac{2377\pi^{4}}{45}+\frac{440\zeta_{3}}{3}+\frac{7192\pi^{2}}{27}+\frac{253252}{81}\right)$
$\displaystyle\phantom{={}}+\left(-28064\zeta_{5}+\frac{1972\zeta_{3}\pi^{2}}{3}-\frac{638\pi^{4}}{15}+4224\mathop{}\\!{\mathrm{Li}}_{4}\\!\left(\frac{1}{2}\right)+3696\zeta_{3}\mathop{}\\!{\mathrm{ln}}(2)-176\pi^{2}\mathop{}\\!{\mathrm{ln}}^{2}(2)+176\mathop{}\\!{\mathrm{ln}}^{4}(2)\right.$
$\displaystyle\phantom{={}}+\left.\frac{13208\zeta_{3}}{3}+\frac{78848\pi^{2}}{81}+96\mathop{}\\!{\mathrm{ln}}(2)+\frac{1925074}{243}\right)$
$\displaystyle\phantom{={}}+\epsilon\left(2304\,\mathop{}\zeta_{-5,-1}-4464\zeta_{5}\mathop{}\\!{\mathrm{ln}}(2)+25784\zeta_{3}^{2}-\frac{67351\pi^{6}}{567}-6336G_{R}(0,0,r_{2},1,-1)\right.$
$\displaystyle\phantom{={}}-6336G_{R}(0,0,1,r_{2},-1)-3168G_{R}(0,0,1,r_{2},r_{4})-6336G_{R}(0,0,r_{2},-1)\mathop{}\\!{\mathrm{ln}}(2)+\frac{268895\zeta_{5}}{3}$
$\displaystyle\phantom{={}}-45056\mathop{}\\!{\mathrm{Li}}_{5}\\!\left(\frac{1}{2}\right)-45056\mathop{}\\!{\mathrm{Li}}_{4}\\!\left(\frac{1}{2}\right)\mathop{}\\!{\mathrm{ln}}(2)+176\mathop{}\\!{\mathrm{Cl}}_{4}\\!\left(\frac{\pi}{3}\right)\pi-1056\zeta_{3}\mathop{}\\!{\mathrm{Li}}_{2}\\!\left(\frac{1}{4}\right)-3982\zeta_{3}\pi^{2}$
$\displaystyle\phantom{={}}-21824\zeta_{3}\mathop{}\\!{\mathrm{ln}}^{2}(2)+2112\zeta_{3}\mathop{}\\!{\mathrm{ln}}(2)\mathop{}\\!{\mathrm{ln}}(3)-1584\mathop{}\\!{\mathrm{Cl}}_{2}^{2}\\!\left(\frac{\pi}{3}\right)\mathop{}\\!{\mathrm{ln}}(3)-\frac{4400\mathop{}\\!{\mathrm{Cl}}_{2}\\!\left(\frac{\pi}{3}\right)\pi^{3}}{27}+\frac{88\pi^{4}\mathop{}\\!{\mathrm{ln}}(2)}{45}$
$\displaystyle\phantom{={}}-\frac{616\pi^{4}\mathop{}\\!{\mathrm{ln}}(3)}{27}+\frac{11264\pi^{2}\mathop{}\\!{\mathrm{ln}}^{3}(2)}{9}-\frac{22528\mathop{}\\!{\mathrm{ln}}^{5}(2)}{15}+8576\mathop{}\\!{\mathrm{Li}}_{4}\\!\left(\frac{1}{2}\right)+7504\zeta_{3}\mathop{}\\!{\mathrm{ln}}(2)+\frac{4174\pi^{4}}{27}$
$\displaystyle\phantom{={}}-\frac{1072\pi^{2}\mathop{}\\!{\mathrm{ln}}^{2}(2)}{3}+\frac{1072\mathop{}\\!{\mathrm{ln}}^{4}(2)}{3}+\frac{554032\zeta_{3}}{27}-32\pi^{2}\mathop{}\\!{\mathrm{ln}}(2)+\frac{730378\pi^{2}}{243}-384\mathop{}\\!{\mathrm{ln}}^{2}(2)+832\mathop{}\\!{\mathrm{ln}}(2)$
$\displaystyle\phantom{={}}+\frac{1408681}{81}+\left.\sqrt{3}\left(192\,\Im\Big{\\{}\\!\mathop{}\\!{\mathrm{Li}}_{3}\\!\left(\frac{\exp(\mathrm{i}\pi/3)}{2}\right)\\!\Big{\\}}+160\mathop{}\\!{\mathrm{Cl}}_{2}\\!\left(\frac{\pi}{3}\right)\mathop{}\\!{\mathrm{ln}}(2)-16\pi\mathop{}\\!{\mathrm{ln}}^{2}(2)-\frac{560\pi^{3}}{81}\right)\right)+\mathop{}\\!\mathcal{O}\\!\left(\epsilon^{2}\right).$
(38)

##  VI Conclusions

In this paper, we have discussed the computation of the same-hemisphere three-
gluon-emission contribution to the zero-jettiness soft function at N3LO in
perturbative QCD. We have used the approach of Ref. [1], which allows us to
apply integration-by-parts technology and the method of differential equations
to phase-space integrals that contain Heaviside functions. While the
appearance of integrals that are not regulated dimensionally requires an
analytic regulator and thus complicates the use of differential equations, we
have described a way to bypass this problem in an efficient way.

Finally, we note that the missing kinematic configuration, in which one of the
three gluons is emitted into the opposite hemisphere can be computed in a
similar fashion. Once a complete result for the three-gluon-emission
contribution is known, the contribution that arises from the emission of a
soft $q\bar{q}$-pair and a soft gluon can be computed in a straightforward
way. Similarly, we expect that virtual corrections to the double real
emissions can be dealt with using the same method. We leave both problems to
future investigations.

Acknowledgments This research was partially supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under grant
396021762-TRR 257 and by Karlsruhe School of Particle and Astroparticle
physics (KSETA). MD is supported by the Excellence Cluster ORIGINS funded by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC-2094 - 390783311 and by the ERC Starting
Grant 949279 HighPHun.

##  Appendix A Computation of abelian contributions

In this appendix we describe how to compute Abelian contributions to the zero-
jettiness soft function, i.e. the first two terms in Eq. (35). The first, so-
called fully-Abelian contribution $S_{1+1+1}^{nnn}$ reads

$\displaystyle S_{1+1+1}^{nnn}=\int{\rm
d}\Phi^{nnn}_{\theta\theta\theta}\,\omega_{n\bar{n}}^{(1)}(k_{1})\,\omega_{n\bar{n}}^{(1)}(k_{2})\,\omega_{n\bar{n}}^{(1)}(k_{3})\,,$
(39)

where [4]

$\displaystyle\omega_{n\bar{n}}^{(1)}(q)=\frac{4}{(nq)\,(\bar{n}q)}\,.$ (40)

Thanks to its fully-factorized structure, the integral in Eq. (39) is
straightforward to compute. Using Sudakov variables $\alpha_{i}$ and
$\beta_{i}$, we find

$\displaystyle S_{1+1+1}^{nnn}$
$\displaystyle=64\int\limits_{0}^{\infty}\left(\prod\limits_{i=1}^{3}\mathop{}\\!\mathrm{d}\alpha_{i}\mathop{}\\!\mathrm{d}\beta_{i}\,\left(\alpha_{i}\beta_{i}\right)^{-1-\epsilon}\mathop{}\\!\theta\\!\left(\alpha_{i}-\beta_{i}\right)\right)$
$\displaystyle\times\mathop{}\\!\delta\\!\left(1-\beta_{123}\right)=\frac{64}{\epsilon^{3}}\,\frac{\mathop{}\\!\Gamma^{3}\\!\left(-2\epsilon\right)}{\mathop{}\\!\Gamma\\!\left(-6\epsilon\right)}\,.$
(41)

Making the $1/\epsilon$-poles in the above equation explicit yields Eq. (36).

The second Abelian contribution reads

$\displaystyle S_{1+2}^{nnn}={}$ $\displaystyle\int{\rm
d}\Phi^{nnn}_{\theta\theta\theta}\Bigg{[}\omega_{n\bar{n}}^{(1)}(k_{1})\,\omega_{n\bar{n}}^{(2)}(k_{2},k_{3})$
$\displaystyle+(k_{1}\leftrightarrow k_{2})+(k_{1}\leftrightarrow
k_{3})\Bigg{]}$ $\displaystyle={}$ $\displaystyle 3\int{\rm
d}\Phi^{nnn}_{\theta\theta\theta}\left[\omega_{n\bar{n}}^{(1)}(k_{1})\,\omega_{n\bar{n}}^{(2)}(k_{2},k_{3})\right]\,.$
(42)

We write it as

$\displaystyle S_{1+2}^{nnn}=3N_{\epsilon}^{-1}\int[{\rm
d}k_{1}]\mathop{}\\!\theta\\!\left(k_{1}\bar{n}-k_{1}n\right)\,\omega_{n\bar{n}}^{(1)}(k_{1})$
$\displaystyle\times N_{\epsilon}^{-2}\int\left(\prod\limits_{i=2}^{3}[{\rm
d}k_{i}]\mathop{}\\!\theta\\!\left(k_{i}\bar{n}-k_{i}n\right)\right)$
$\displaystyle\times\mathop{}\\!\delta\\!\left(1-k_{123}n\right)\omega_{n\bar{n}}^{(2)}(k_{2},k_{3})\,.$
(43)

The inner integral in Eq. (A) over $[{\rm d}k_{2}][{\rm d}k_{3}]$ can be
obtained from the same-hemisphere double-real gluon emission contribution to
the NNLO soft function. We find

$\displaystyle
S_{1+2}^{nnn}=3\,\,\int_{0}^{\infty}\mathop{}\\!\mathrm{d}\alpha_{1}\mathop{}\\!\mathrm{d}\beta_{1}\mathop{}\\!\theta\\!\left(\alpha_{1}-\beta_{1}\right)\frac{4(\alpha_{1}\beta_{1})^{-\epsilon}}{\alpha_{1}\beta_{1}}$
$\displaystyle\times
C^{nn}_{2}(1-\beta_{1})^{-1-4\epsilon}=\frac{12}{\epsilon}\frac{\mathop{}\\!\Gamma\\!\left(-4\epsilon\right)\mathop{}\\!\Gamma\\!\left(-2\epsilon\right)}{\mathop{}\\!\Gamma\\!\left(-6\epsilon\right)}\,C^{nn}_{2}\,,$
(44)

where the factor $C^{nn}_{2}$ can be extracted from Refs. [16, 1]. Upon doing
so, we obtain the result displayed in Eq. (37).

## References

  * Baranowski  _et al._ [2022] D. Baranowski, M. Delto, K. Melnikov, and C.-Y. Wang, On phase-space integrals with Heaviside functions, JHEP 02, 081, arXiv:2111.13594 [hep-ph] .
  * Anastasiou and Melnikov [2002] C. Anastasiou and K. Melnikov, Higgs boson production at hadron colliders in NNLO QCD, Nucl. Phys. B 646, 220 (2002), arXiv:hep-ph/0207004 .
  * Chetyrkin and Tkachov [1981] K. G. Chetyrkin and F. V. Tkachov, Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops, Nucl. Phys. B 192, 159 (1981).
  * Catani  _et al._ [2020] S. Catani, D. Colferai, and A. Torrini, Triple (and quadruple) soft-gluon radiation in QCD hard scattering, JHEP 01, 118, arXiv:1908.01616 [hep-ph] .
  * Henn  _et al._ [2014] J. M. Henn, A. V. Smirnov, and V. A. Smirnov, Evaluating single-scale and/or non-planar diagrams by differential equations, JHEP 03, 088, arXiv:1312.2588 [hep-th] .
  * Papadopoulos [2014] C. G. Papadopoulos, Simplified differential equations approach for Master Integrals, JHEP 07, 088, arXiv:1401.6057 [hep-ph] .
  * Liu  _et al._ [2018] X. Liu, Y.-Q. Ma, and C.-Y. Wang, A Systematic and Efficient Method to Compute Multi-loop Master Integrals, Phys. Lett. B 779, 353 (2018), arXiv:1711.09572 [hep-ph] .
  * Lee  _et al._ [2022] R. N. Lee, A. von Manteuffel, R. M. Schabinger, A. V. Smirnov, V. A. Smirnov, and M. Steinhauser, Quark and gluon form factors in four-loop QCD, (2022), arXiv:2202.04660 [hep-ph] .
  * Czakon [2006] M. Czakon, Automatized analytic continuation of Mellin-Barnes integrals, Comput. Phys. Commun. 175, 559 (2006), arXiv:hep-ph/0511200 .
  * Smirnov and Smirnov [2009] A. V. Smirnov and V. A. Smirnov, On the Resolution of Singularities of Multiple Mellin-Barnes Integrals, Eur. Phys. J. C 62, 445 (2009), arXiv:0901.0386 [hep-ph] .
  * Borowka  _et al._ [2018] S. Borowka, G. Heinrich, S. Jahn, S. Jones, M. Kerner, J. Schlenk, and T. Zirke, pySecDec: a toolbox for the numerical evaluation of multi-scale integrals, Comput. Phys. Commun. 222, 313–326 (2018), arXiv:1703.09692 [hep-ph] .
  * Borowka  _et al._ [2019] S. Borowka, G. Heinrich, S. Jahn, S. Jones, M. Kerner, and J. Schlenk, A GPU compatible quasi-Monte Carlo integrator interfaced to pySecDec, Comput. Phys. Commun. 240, 120–137 (2019), arXiv:1811.11720 [physics.comp-ph] .
  * Ferguson and Bailey [1992] H. R. P. Ferguson and D. H. Bailey,  _A polynomial time, numerically stable integer relation algorithm_ , Tech. Rep. (1992).
  * Lenstra  _et al._ [1982] A. K. Lenstra, H. W. Lenstra, and L. Lovász, Factoring polynomials with rational coefficients, Mathematische annalen 261, 515 (1982).
  * Henn  _et al._ [2017] J. M. Henn, A. V. Smirnov, and V. A. Smirnov, Evaluating Multiple Polylogarithm Values at Sixth Roots of Unity up to Weight Six, Nucl. Phys. B 919, 315 (2017), arXiv:1512.08389 [hep-th] .
  * Baranowski [2020] D. Baranowski, NNLO zero-jettiness beam and soft functions to higher orders in the dimensional-regularization parameter $\epsilon$, Eur. Phys. J. C 80, 523 (2020), arXiv:2004.03285 [hep-ph] .