# A Toolbox for $q_{T}$ and 0-Jettiness Subtractions at N3LO

Georgios Billis  Markus A. Ebert  Johannes K. L. Michel  and Frank J. Tackmann

(September 2, 2019)

###### Abstract

We derive the leading-power singular terms at three loops for both $q_{T}$ and
0-jettiness, $\mathcal{T}_{0}$, for generic color-singlet processes. Our
results provide the complete set of differential subtraction terms for $q_{T}$
and $\mathcal{T}_{0}$ subtractions at N3LO, which are an important ingredient
for matching N3LO calculations with parton showers. We obtain the full three-
loop structure of the relevant beam and soft functions, which are necessary
ingredients for the resummation of $q_{T}$ and $\mathcal{T}_{0}$ at N3LL′ and
N4LL order, and which constitute important building blocks in other contexts
as well. The nonlogarithmic boundary coefficients of the beam functions, which
contribute to the integrated subtraction terms, are not yet fully known at
three loops. By exploiting consistency relations between different
factorization limits, we derive results for the $q_{T}$ and $\mathcal{T}_{0}$
beam function coefficients at N3LO in the $z\to 1$ threshold limit, and we
also estimate the size of the unknown terms beyond threshold.

##  1 Introduction

The ever increasing precision of experimental measurements at the LHC requires
equally precise theoretical predictions in order to be fully exploited. Color-
singlet processes play a central role in the LHC physics program. The $pp\to
Z,W$ Drell-Yan processes are key benchmark processes that have been measured
at the percent level and below [1, 2, 3, 4]. Precise measurements of Higgs and
diboson processes provide strong sensitivity to possible contributions beyond
the Standard Model [5, 6, 7, 8, 9, 10]. They are also important irreducible
backgrounds in direct searches for dark-matter production at the LHC.

The inclusion of higher-order QCD corrections is crucial to obtain reliable
predictions. Depending on the specific process and phase-space region,
reducing the current theoretical uncertainties requires the calculation of the
full corrections at the next order in $\alpha_{s}$ and/or the resummation of
the dominant higher-order terms to all orders in $\alpha_{s}$. For color-
singlet processes, theory predictions are being pushed to the third order in
the fixed-order expansion [11, 12, 13, 14, 15, 16, 17, 18, 19, 20] as well as
in resummed perturbation theory [21, 22, 23, 24, 25, 26, 27, 28, 29].

A key requirement in both cases is to understand the infrared singular
structure of QCD at N3LO. For fixed-order calculations this is crucial for the
cancellation of infrared divergences between real and virtual emissions, as
evidenced by the variety of methods that have been developed by now at NNLO
[30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45]. Resummed
predictions are intimately linked to the singular limit, and the N3LO singular
structure is a key ingredient to extend the resummation to the full three-loop
level.

One way to study the infrared singular limit of QCD is to consider a suitable
resolution variable $\tau$, whose differential cross section can be factorized
in the singular limit $\tau\to 0$. In this paper, we consider two such
variables, $0$-jettiness $\mathcal{T}_{0}$ and the total color-singlet
transverse momentum $q_{T}$, and derive their singular structure to N3LO. Our
results are necessary ingredients for carrying out the resummation for $q_{T}$
and $\mathcal{T}_{0}$ at N3LL′ and N4LL order. For the associated $q_{T}$ and
$\mathcal{T}_{0}$ subtraction methods [35, 38, 39], we provide the complete
differential N3LO subtraction terms in analytic form. The structure and
required ingredients for the $\mathcal{T}_{0}$ subtractions at N3LO were
already discussed in ref. [39]. The $q_{T}$ slicing method at N3LO was also
considered in ref. [18]. Differential $\mathcal{T}_{0}$ subtractions are the
basis of the NNLO$+$PS (parton shower) matching in Geneva [46, 47], and our
results are an important ingredient for its extension to N3LO$+$PS.

To continue our discussion we need to set up some basic notation. We consider
the production of a generic color-singlet final state $L$ in hadronic
collisions. In the singular limit, the only hard interaction process that
contributes is the Born process, which we denote as

$\kappa_{a}(q_{a})\,\kappa_{b}(q_{b})\to L(q)\quad\text{with}\quad
q_{a}^{\mu}+q_{b}^{\mu}=q^{\mu}\,.$ (1.1)

We always use the indices $a$ and $b$ to label the initial states, and
$\kappa_{i}\in\\{g,u,\bar{u},d,\bar{d},s,\ldots\\}$ denotes the parton type
and flavor. When there is no ambiguity we simply identify $\kappa_{i}\equiv
i$, e.g., we just write $ab\to L$. The $q_{a,b}^{\mu}$ are lightlike Born
reference (label) momenta given by

$q_{a,b}^{\mu}=\omega_{a,b}\frac{n_{a,b}^{\mu}}{2}\,,\qquad n_{a}^{\mu}\equiv
n^{\mu}=(1,\hat{z})\,,\quad
n_{b}^{\mu}\equiv{\bar{n}}^{\mu}=(1,-\hat{z})\,,\quad
x_{a,b}=\frac{\omega_{a,b}}{E_{\mathrm{cm}}}=\frac{Q}{E_{\mathrm{cm}}}e^{\pm
Y}\,,$ (1.2)

where $\hat{z}$ is the beam direction and $E_{\rm cm}$ is the hadronic center-
of-mass energy. The precise definition of the Born momentum fractions
$x_{a,b}$ and the associated $\omega_{a,b}=x_{a,b}E_{\mathrm{cm}}$ depends on
how we choose to parametrize the Born phase space in terms of physical
observables. For definiteness, in eq. (1.2) we have chosen the total invariant
mass $Q=\sqrt{q^{2}}$ and rapidity $Y$ of the color singlet. Other choices are
possible as well, e.g.,
$\omega_{a,b}=q^{\mp}\equiv\sqrt{Q^{2}+q_{T}^{2}}e^{\pm Y}$. In the singular
limit, all possible choices are equivalent and yield the same factorized cross
section. The specific choice affects the nonsingular power-suppressed
corrections.

The cross section for color-singlet production for a (suitably factorizable)
resolution variable $\tau$ factorizes for $\tau\to 0$ as [48, 49]

$\frac{\mathrm{d}\sigma}{\mathrm{d}Q^{2}\mathrm{d}Y\mathrm{d}\tau}=\sum_{a,b}H_{ab}(\omega_{a}\omega_{b})\,\bigl{[}B_{a}(x_{a})\otimes
B_{b}(x_{b})\otimes
S_{c}\bigr{]}(\tau)\times\bigl{[}1+\mathcal{O}(\tau)\bigr{]}\,.$ (1.3)

The convolution structure denoted by $\otimes$ depends on the precise
definition of $\tau$. The key properties of eq. (1.3) are that it captures all
QCD singularities in $\tau$ and that it factorizes the dependence on the
underlying process from the dependence on $\tau$.

The process dependence is carried by the hard function $H_{ab}$, which
describes the Born process $ab\to L$, with the sum running over all relevant
parton channels. At lowest order, $H_{ab}^{(0)}$ is equivalent to the partonic
Born cross section for $ab\to L$. At higher orders, it encodes the finite
virtual corrections to the Born process and thus can be obtained from the
corresponding quark or gluon form factors. Results at three loops are known
for $gg\to H$ in the $m_{t}\to\infty$ limit, $b\bar{b}\to H$, and Drell-Yan
production [50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63]. Explicit
expressions for the hard functions in our notation can be found in ref. [27].
The hard function also encodes any additional cuts or measurements on the
constituents of $L$, which we keep implicit.

The entire $\tau$ dependence in eq. (1.3) is encoded by the beam and soft
functions. The beam functions $B_{a,b}$ describe collinear emissions from the
incoming partons $a$ and $b$, while the soft function $S_{c}$ encodes soft
radiation between them. They are universal objects that do not depend on the
details of the hard process. Namely, $B_{i}$ only depends on the type of its
incoming parton $i\equiv\kappa_{i}$, while $S_{c}$ only depends on the color
channel of the Born process. In our case, the only possible color channels are
$c=\\{q\bar{q},gg\\}$, which are equivalent to the color representation of the
incoming partons, so we simply label it by $c\equiv i=\\{q,g\\}$.

The beam and soft functions do depend on the definition of $\tau$, which also
determines their convolution structure in eq. (1.3). They can be formally
defined as renormalized operator matrix elements in soft-collinear effective
theory (SCET) [64, 65, 66, 67, 68]. The beam and soft functions relevant for
$\mathcal{T}_{0}$ and $q_{T}$ are the most basic of their type, measuring the
small light-cone momentum or the total transverse momentum of the inclusive
sum of all collinear and soft emissions, respectively. For this reason, they
are important objects in their own right, encoding fundamental properties of
the singular structure of QCD, and also appear in a variety of other contexts.
In particular, they often serve as building blocks for constructing the beam
and soft functions necessary for more complicated scenarios or observables,
see e.g. refs. [69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81].

In this paper, we derive the analytic structure of the $\mathcal{T}_{0}$ and
$q_{T}$ beam and soft functions at three loops from their known
renormalization group equations (RGEs). The nonlogarithmic boundary
coefficients are not predicted by the RGE and require an explicit three-loop
calculation. While they are not required for the differential subtraction
terms, they are the essential ingredient required for the integrated
subtraction terms. So far, they are known at three loops for the $q_{T}$ soft
function [82].

The most complicated are the beam function boundary coefficients, because they
are nontrivial functions of a partonic momentum fraction $z$. However, they
drastically simplify in the limit $z\to 1$. In this limit, the energy of
collinear emissions is constrained to be small which means their interactions
with the primary collinear parton can be described in the eikonal
approximation where they only resolve its color charge and direction. This was
already pointed out and exploited at NNLO in refs. [85, 86]. Here we exploit
this to obtain for the first time the three-loop beam function coefficients in
the $z\to 1$ limit for both $\mathcal{T}_{0}$ and $q_{T}$ by relating them via
appropriate consistency relations to known soft matrix elements. The required
consistency relations were only partially known so far. We give their detailed
derivation and show explicitly that they hold to all orders, allowing one to
obtain the $z\to 1$ limits of the beam functions also to higher orders once
the relevant soft matrix elements are known. In case of $q_{T}$, we provide
their general structure for illustration up to six loops. We find that a
previous conjecture for this limit [137] only holds up to N3LO but fails
starting at N4LO. Since our results capture the complete singular structure
for $z\to 1$, they can also simplify the full calculation because it can be
carried out strictly for $z<1$ which reduces the degree of divergences. We
also employ the obtained eikonal terms of the beam function coefficients to
construct an ansatz for the missing next-to-eikonal coefficients to estimate
their numerical size.

When this paper first appeared, the $\mathcal{T}_{0}$ beam function was only
known at NNLO [87, 88, 85, 86], while only partial results were known at N3LO
[83, 84]. By now, the complete results for the three-loop beam functions have
become available [91, 94], for which our results provided important cross
checks. In particular, the predicted $z\to 1$ limit was the only available
check of the genuine three-loop contribution. Similarly, the complete results
for the three-loop beam functions for $q_{T}$ have become available [92, 93],
for which our results in the $z\to 1$ limit again provided important checks.

For $q_{T}$, a similar study of the logarithmic structure at N3LO was
performed in ref. [18] to construct an approximate $q_{T}$ subtraction at this
order. Here, we present a more detailed derivation of its fixed-order
structure and the ensuing $q_{T}$ subtraction, which differs from the method
employed in ref. [18]. While the RGEs necessary to derive the three-loop
differential subtractions are in principle known in the literature, we provide
here for the first time a comprehensive account of the complete structure for
both $q_{T}$ and $\mathcal{T}_{0}$. All required perturbative ingredients are
collected in the Appendix, while the results for the three-loop beam functions
can be directly used together with our results. This provides the complete
results for $q_{T}$ and $\mathcal{T}_{0}$ subtractions for $q\bar{q}$ and $gg$
processes at N3LO.

In the remainder of this section, we summarize important conventions used
throughout this paper. The three-loop structure of the beam and soft functions
and the eikonal limit of the beam functions are derived for $\mathcal{T}_{0}$
in section 2 and for $q_{T}$ in section 3. The application to
$\mathcal{T}_{0}$ and $q_{T}$ subtractions at N3LO is discussed in section 4.
Readers primarily interested in this application may directly skip ahead to
section 4. We conclude in section 5. In appendix A we collect the needed
definitions and relations for plus distributions. In appendices B and C we
discuss in more detail the soft matrix elements that are involved in
extracting the eikonal limits of the beam functions. Explicit expressions for
required perturbative ingredients are collected in appendix D.

###  1.1 Notation and conventions

Throughout the paper, we denote the perturbative expansion of any function
$F(\mu)$ as

$F(\mu)=\sum_{n=0}^{\infty}F^{(n)}(\mu)\,\Bigl{[}\frac{\alpha_{s}(\mu)}{4\pi}\Bigr{]}^{n}\,.$
(1.4)

All anomalous dimensions $\gamma^{i}_{x}(\alpha_{s})$ and the QCD splitting
functions are expanded as

$\gamma^{i}_{x}(\alpha_{s})=\sum_{n=0}^{\infty}\gamma_{x\,n}^{i}\,\Bigl{(}\frac{\alpha_{s}}{4\pi}\Bigr{)}^{n+1}\,,\qquad
P_{ij}(z,\mu)=\sum_{n=0}^{\infty}P_{ij}^{(n)}(z)\biggl{[}\frac{\alpha_{s}(\mu)}{4\pi}\biggr{]}^{n+1}\,.$
(1.5)

We use the following notation to abbreviate Mellin convolutions and flavor
sums

$\displaystyle(I^{(m)}P^{(n)})_{ij}(z)$
$\displaystyle\equiv\sum_{k}\int\\!\frac{\mathrm{d}z^{\prime}}{z^{\prime}}\,I^{(m)}_{ik}\Bigl{(}\frac{z}{z^{\prime}}\Bigr{)}P^{(n)}_{kj}(z^{\prime})\,,$
$\displaystyle[\mathcal{I}(t,\mu)P^{(n)}]_{ij}(z)$
$\displaystyle\equiv\sum_{k}\int\\!\frac{\mathrm{d}z^{\prime}}{z^{\prime}}\,\mathcal{I}_{ik}\Bigl{(}t,\frac{z}{z^{\prime}},\mu\Bigr{)}P^{(n)}_{kj}(z^{\prime})\,,$
(1.6)

where $i,j,k\in\\{g,u,\bar{u},d,\bar{d},s,\ldots\\}$ label parton type and
flavor. We also define a corresponding identity operator as

$\mathbf{1}_{ij}(z)\equiv\delta_{ij}\,\delta(1-z)\quad\text{with}\quad(\mathbf{1}P)_{ij}(z)=P_{ij}(z)\,.$
(1.7)

For Fourier-type convolutions, we use the notation

$\displaystyle(FG)(k,\mu)$
$\displaystyle\equiv\int\\!\mathrm{d}k^{\prime}\,F(k-k^{\prime},\mu)\,G(k^{\prime},\mu)\,,$
$\displaystyle\bigl{[}F\,\mathcal{I}_{ij}(z)\bigr{]}(t,\mu^{2})$
$\displaystyle\equiv\int\\!\mathrm{d}t^{\prime}\,F(t-t^{\prime},\mu^{2})\,\mathcal{I}_{ij}(t^{\prime},z,\mu^{2})\,.$
(1.8)

Here, the corresponding identity elements are simply $\delta(k)$ or
$\delta(t)$.

We denote logarithmic plus distributions as

$\mathcal{L}_{n}(x)=\biggl{[}\frac{\theta(x)\ln^{n}x}{x}\biggr{]}_{+}\qquad\text{with}\qquad\int_{0}^{1}\\!\mathrm{d}x\,\mathcal{L}_{n}(x)=0\,.$
(1.9)

For dimensionful arguments, we define

$\mathcal{L}_{n}(k,\mu)=\frac{1}{\mu}\mathcal{L}_{n}\Bigl{(}\frac{k}{\mu}\Bigr{)}\,,\quad\mathcal{L}_{n}(t,\mu^{2})=\frac{1}{\mu^{2}}\mathcal{L}_{n}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}\,,\quad\mathcal{L}_{n}({\vec{q}}_{T},\mu)=\frac{1}{\pi\mu^{2}}\mathcal{L}_{n}\biggl{(}\frac{q_{T}^{2}}{\mu^{2}}\biggr{)}\,.$
(1.10)

More details are given in appendix A.

##  2 $\mathcal{T}_{0}$ factorization to three loops

###  2.1 Factorization

The factorization for $N$-jettiness, $\mathcal{T}_{N}$, has been derived using
SCET in refs. [49, 95, 69]. Here we focus on $0$-jettiness, $\mathcal{T}_{0}$,
which is relevant for color-singlet production and coincides with beam thrust
[49, 88]. It can be defined in terms of generic measures as [95, 69]

$\displaystyle\mathcal{T}_{0}=\sum_{i}{\rm min}\Bigl{\\{}\frac{2q_{a}\cdot
k_{i}}{Q_{a}},\,\frac{2q_{b}\cdot k_{i}}{Q_{b}}\Bigr{\\}}\,,$ (2.1)

where the sum runs over the momenta $k_{i}$ of all final-state particles
excluding $L$ and any of its constituents. The measures $Q_{a,b}$ determine
different definitions of 0-jettiness. Two possible choices, corresponding to
the original definitions in refs. [49, 88], are

leptonic: $\displaystyle Q_{a}$
$\displaystyle=Q_{b}=\sqrt{\omega_{a}\omega_{b}}=Q\,,\qquad$
$\displaystyle\mathcal{T}_{0}^{\rm lep}$
$\displaystyle=\sum_{i}\min\Bigl{\\{}e^{Y}n_{a}\cdot
k_{i}\,,\,e^{-Y}n_{b}\cdot k_{i}\Bigr{\\}}$ hadronic: $\displaystyle Q_{a,b}$
$\displaystyle=\omega_{a,b}=Q\,e^{\pm Y}\,,$
$\displaystyle\mathcal{T}_{0}^{\rm cm}$
$\displaystyle=\sum_{i}\min\Bigl{\\{}n_{a}\cdot k_{i}\,,\,n_{b}\cdot
k_{i}\Bigr{\\}}\,.$ (2.2)

For our present purposes, the precise choice of the $Q_{i}$ is not important,
so we will simply use the symbol $\mathcal{T}_{0}$.

The factorization for $\mathcal{T}_{0}$ is given by [49]

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}Q^{2}\mathrm{d}Y\mathrm{d}\mathcal{T}_{0}}$
$\displaystyle=\sum_{a,b}H_{ab}(Q^{2},\mu)\int\\!\mathrm{d}t_{a}\,\mathrm{d}t_{b}\,B_{a}(t_{a},x_{a},\mu)\,B_{b}(t_{b},x_{b},\mu)\,S_{i}\Bigl{(}\mathcal{T}_{0}-\frac{t_{a}}{Q_{a}}-\frac{t_{b}}{Q_{b}},\mu\Bigr{)}$
$\displaystyle\quad\times\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{\mathcal{T}_{0}}{Q}\Bigr{)}\Bigr{]}\,.$
(2.3)

Explicit definitions of the beam and soft functions for $\mathcal{T}_{0}$ in
terms of operator matrix elements in SCET can be found in refs. [49, 87].

The beam function appearing in eq. (2.1) is the inclusive virtuality-dependent
(SCETI) beam function. It appears in the $N$-jettiness factorization for any
$N$ [95], including deep-inelastic scattering [96]. Recently, it was shown
that it also arises in generalized threshold factorization theorems for
inclusive color-singlet production in hadronic collisions [97]. The
virtuality-dependent quark and gluon beam functions are known to NNLO [87, 88,
85, 86], and they are being calculated at N3LO [83, 84].

The soft function in eq. (2.1) is the hemisphere soft function for incoming
Wilson lines. It is closely related to the hemisphere soft function for
$e^{+}e^{-}\to$ jets, which is known to NNLO [98, 99, 100, 101, 102]. They
have the same anomalous dimensions to all orders [49, 87], and are equal to
NNLO [49, 103]. It is an open question whether they remain equivalent at
higher fixed orders.

The factorization in eq. (2.1) receives power corrections suppressed by
$\mathcal{T}_{0}/Q$, as indicated. In addition, starting at N4LO it also
receives contributions from perturbative Glauber-gluon exchanges, which are
not captured by eq. (2.1) [104, 105].

###  2.2 $\mathcal{T}_{0}$ soft function

The beam thrust soft function satisfies the all-order RGE [49, 87]

$\displaystyle\mu\frac{\mathrm{d}}{\mathrm{d}\mu}S_{i}(k,\mu)$
$\displaystyle=\int\mathrm{d}k^{\prime}\gamma_{S}^{i}(k-k^{\prime},\mu)\,S_{i}(k^{\prime},\mu)\equiv(\gamma_{S}^{i}\,S_{i})(k,\mu)\,,$
$\displaystyle\gamma_{S}^{i}(k,\mu)$
$\displaystyle=4\Gamma_{\mathrm{cusp}}^{i}[\alpha_{s}(\mu)]\,\mathcal{L}_{0}(k,\mu)+\gamma_{S}^{i}[\alpha_{s}(\mu)]\delta(k)\,,$
(2.4)

where $\Gamma_{\mathrm{cusp}}^{i}(\alpha_{s})$ and
$\gamma_{S}^{i}(\alpha_{s})$ are the cusp and soft noncusp anomalous
dimensions.

The RGE in eq. (2.2) fully predicts the structure of $S_{i}(k,\mu)$ in $k$ and
$\mu$ to all orders in perturbation theory. By solving it recursively order by
order in $\alpha_{s}$, we can derive this structure at any given fixed order.
Expanding both sides of eq. (2.2) to fixed order in $\alpha_{s}(\mu)$ and
accounting for the running of $\alpha_{s}(\mu)$, we obtain a relation for the
$(n+1)$-loop term in terms of the terms up to $n$ loops,

$\mu\frac{\mathrm{d}}{\mathrm{d}\mu}S_{i}^{(n+1)}(k,\mu)=\sum_{m=0}^{n}\Bigl{[}4\Gamma^{i}_{n-m}\bigl{(}\mathcal{L}_{0}S_{i}^{(m)}\bigr{)}(k,\mu)+\bigl{(}2m\beta_{n-m}+\gamma_{S\,n-m}^{i}\bigr{)}S_{i}^{(m)}(k,\mu)\Bigr{]},$
(2.5)

where we used the short-hand notation in eq. (1.1) for the convolution in $k$.
This can be integrated to give

$\displaystyle S_{i}^{(n+1)}(k,\mu)$
$\displaystyle=\int_{k|_{+}}^{\mu}\\!\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\sum_{m=0}^{n}\Bigl{[}4\Gamma^{i}_{n-m}\bigl{(}\mathcal{L}_{0}S_{i}^{(m)}\bigr{)}(k,\mu)+\bigl{(}2m\beta_{n-m}+\gamma_{S\,n-m}^{i}\bigr{)}S_{i}^{(m)}(k,\mu)\Bigr{]}$
$\displaystyle\quad+\delta(k)\,s_{i}^{(n+1)}\,,$ (2.6)

where the soft function boundary coefficients are defined by

$S_{i}^{(n)}(k,\mu=k|_{+})=\delta(k)\,s_{i}^{(n)}\qquad\text{with}\qquad
s_{i}^{(0)}=1\,.$ (2.7)

Here, we have used distributional scale setting $\mu_{0}=k|_{+}$ [106], which
is defined such that it effectively allows us to treat the $\mu$ dependence of
the logarithmic distributions like ordinary logarithms. In particular, it
satisfies [106]

$\displaystyle\mathcal{L}_{n}(k,\mu_{0}=k|_{+})$ $\displaystyle=0\,,$
$\displaystyle\delta(k)\ln^{n+1}\frac{\mu_{0}}{\mu}\bigg{|}_{\mu_{0}=k|_{+}}$
$\displaystyle=(n+1)\mathcal{L}_{n}(k,\mu)\,,$
$\displaystyle\int_{\mu_{0}=k|_{+}}^{\mu}\\!\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\,\mathcal{L}_{n}(k,\mu^{\prime})$
$\displaystyle=-\frac{1}{n+1}\mathcal{L}_{n+1}(k,\mu)\,.$ (2.8)

The first relation is used in eq. (2.7) to define the boundary coefficients as
the coefficients of the $\delta(k)$ by setting all logarithmic distributions
in $S_{i}^{(n)}(k,\mu)$ to zero. The other two relations allow us to easily
perform the $\mu^{\prime}$ integral in eq. (2.2), essentially turning a
$\delta(k)$ into a $\mathcal{L}_{0}(k)$ and a $\mathcal{L}_{n}(k)$ into a
$\mathcal{L}_{n+1}(k)$. In addition, to evaluate the cross terms for $m\geq 1$
in eq. (2.2), we need the convolutions [107]

$\displaystyle(\mathcal{L}_{0}\mathcal{L}_{0})(k,\mu)$
$\displaystyle=2\mathcal{L}_{1}(k,\mu)-\zeta_{2}\,\delta(k)\,,$
$\displaystyle(\mathcal{L}_{0}\mathcal{L}_{1})(k,\mu)$
$\displaystyle=\frac{3}{2}\mathcal{L}_{2}(k,\mu)-\zeta_{2}\mathcal{L}_{0}(k,\mu)+\zeta_{3}\,\delta(k)\,,$
$\displaystyle(\mathcal{L}_{0}\mathcal{L}_{2})(k,\mu)$
$\displaystyle=\frac{4}{3}\mathcal{L}_{3}(k,\mu)-2\zeta_{2}\mathcal{L}_{1}(k,\mu)+2\zeta_{3}\mathcal{L}_{0}(k,\mu)-2\zeta_{4}\,\delta(k)\,,$
$\displaystyle(\mathcal{L}_{0}\mathcal{L}_{3})(k,\mu)$
$\displaystyle=\frac{5}{4}\mathcal{L}_{4}(k,\mu)-3\zeta_{2}\mathcal{L}_{2}(k,\mu)+6\zeta_{3}\mathcal{L}_{1}(k,\mu)-6\zeta_{4}\mathcal{L}_{0}(k,\mu)+6\zeta_{5}\,\delta(k)\,.$
(2.9)

Starting from the LO result, $s_{i}^{(0)}=1$, eq. (2.2) yields up to two loops

$\displaystyle S_{i}^{(0)}(k,\mu)$ $\displaystyle=\delta(k)\,,$ $\displaystyle
S_{i}^{(1)}(k,\mu)$
$\displaystyle=-\mathcal{L}_{1}(k,\mu)\,4\Gamma_{0}^{i}-\mathcal{L}_{0}(k,\mu)\,\gamma^{i}_{S\,0}+\delta(k)\,s_{i}^{(1)}\,,$
$\displaystyle S_{i}^{(2)}(k,\mu)$
$\displaystyle=\mathcal{L}_{3}(k,\mu)\,8(\Gamma^{i}_{0})^{2}+\mathcal{L}_{2}(k,\mu)\,2\Gamma^{i}_{0}(2\beta_{0}+3\gamma_{S\,0}^{i})$
$\displaystyle\quad+\mathcal{L}_{1}(k,\mu)\Bigl{[}-16\zeta_{2}(\Gamma^{i}_{0})^{2}+(2\beta_{0}+\gamma_{S\,0}^{i})\gamma_{S\,0}^{i}-4\Gamma^{i}_{1}-4\Gamma^{i}_{0}\,s_{i}^{(1)}\Bigr{]}$
$\displaystyle\quad+\mathcal{L}_{0}(k,\mu)\Bigl{[}4\Gamma^{i}_{0}(4\zeta_{3}\Gamma^{i}_{0}-\zeta_{2}\gamma_{S\,0}^{i})-\gamma_{S\,1}^{i}-(2\beta_{0}+\gamma_{S\,0}^{i})\,s_{i}^{(1)}\Bigr{]}+\delta(k)\,s_{i}^{(2)}\,,$
(2.10)

which agrees with ref. [39]. Evaluating eq. (2.2) at the next order, we obtain
the three-loop result,

$S^{(3)}(k,\mu)=\delta(k)\,s_{i}^{(3)}+\sum_{\ell=0}^{5}S_{i,\ell}^{(3)}\,\mathcal{L}_{\ell}(k,\mu)$
(2.11)

with the coefficients of the logarithmic distributions given by

$\displaystyle S_{i,5}^{(3)}$ $\displaystyle=-8(\Gamma^{i}_{0})^{3}$
$\displaystyle S_{i,4}^{(3)}$
$\displaystyle=-\frac{10}{3}(\Gamma^{i}_{0})^{2}\,(4\beta_{0}+3\gamma_{S\,0}^{i})$
$\displaystyle S_{i,3}^{(3)}$
$\displaystyle=4\Gamma^{i}_{0}\Bigl{[}16\zeta_{2}(\Gamma^{i}_{0})^{2}-\frac{4}{3}\beta_{0}^{2}-\Bigl{(}\frac{10}{3}\beta_{0}+\gamma_{S\,0}^{i}\Bigr{)}\gamma_{S\,0}^{i}+4\Gamma_{1}^{i}+2\Gamma_{0}^{i}\,s_{i}^{(1)}\Bigr{]}$
$\displaystyle S_{i,2}^{(3)}$
$\displaystyle=16(\Gamma^{i}_{0})^{2}\bigl{[}-10\zeta_{3}\Gamma^{i}_{0}+3\zeta_{2}(\beta_{0}+\gamma^{i}_{S\,0})\bigr{]}+(4\beta_{0}+3\gamma^{i}_{S\,0})(-\beta_{0}\gamma^{i}_{S\,0}+2\Gamma^{i}_{1})-\frac{(\gamma^{i}_{S\,0})^{3}}{2}$
$\displaystyle\quad+2\Gamma^{i}_{0}\Bigl{[}2\beta_{1}+3\gamma^{i}_{S\,1}+(8\beta_{0}+3\gamma^{i}_{S\,0})s_{i}^{(1)}\Bigr{]}$
$\displaystyle S_{i,1}^{(3)}$
$\displaystyle=32(\Gamma^{i}_{0})^{2}\bigl{[}\zeta_{4}\Gamma^{i}_{0}-\zeta_{3}(3\beta_{0}+2\gamma^{i}_{S\,0})\bigr{]}+8\zeta_{2}\Gamma^{i}_{0}\bigl{[}(3\beta_{0}+\gamma^{i}_{S\,0})\gamma^{i}_{S\,0}-4\Gamma^{i}_{1}\bigr{]}$
$\displaystyle\quad+4\beta_{0}\gamma^{i}_{S\,1}+2\gamma^{i}_{S\,0}(\beta_{1}+\gamma^{i}_{S\,1})+\bigl{[}-16\zeta_{2}(\Gamma^{i}_{0})^{2}+8\beta_{0}^{2}+(6\beta_{0}+\gamma^{i}_{S\,0})\gamma^{i}_{S\,0}-4\Gamma^{i}_{1}\bigr{]}s_{i}^{(1)}$
$\displaystyle\quad-4\Gamma^{i}_{2}-4\Gamma^{i}_{0}\,s_{i}^{(2)}$
$\displaystyle S_{i,0}^{(3)}$
$\displaystyle=16(\Gamma^{i}_{0})^{2}\Bigl{[}4\Gamma^{i}_{0}(2\zeta_{2}\zeta_{3}-3\zeta_{5})+\zeta_{4}\Bigl{(}2\beta_{0}+\frac{\gamma^{i}_{S\,0}}{2}\Bigr{)}\Bigr{]}-4\zeta_{3}\Gamma^{i}_{0}\bigl{[}(2\beta_{0}+\gamma^{i}_{S\,0})\gamma^{i}_{S\,0}-8\Gamma^{i}_{1}\bigr{]}$
$\displaystyle\quad-4\zeta_{2}\bigl{(}\gamma^{i}_{S\,0}\Gamma^{i}_{1}+\Gamma^{i}_{0}\gamma^{i}_{S\,1}\bigr{)}+\Bigl{\\{}4\Gamma^{i}_{0}\bigl{[}4\zeta_{3}\Gamma^{i}_{0}-\zeta_{2}(2\beta_{0}+\gamma^{i}_{S\,0})\bigr{]}-(2\beta_{1}+\gamma^{i}_{S\,1})\Bigr{\\}}s_{i}^{(1)}$
$\displaystyle\quad-\gamma^{i}_{S\,2}-(4\beta_{0}+\gamma^{i}_{S\,0})s_{i}^{(2)}\,.$
(2.12)

This agrees with a corresponding numerical expression in ref. [22]. The
required anomalous dimension coefficients up to three loops and boundary
coefficients up to two loops are given in appendix D.

#### Numerical impact

The soft function $S_{i}(k,\mu)$ has an explicit dependence on $\mu$, which
cancels against that of the hard and beam functions in eq. (2.1). Therefore,
simply varying the scale $\mu$ is not very meaningful for illustrating the
numerical impact of the $\mu$-dependent three-loop terms. Instead, we consider
the resummed soft function,

$S_{i}(k,\mu)=\int\\!\mathrm{d}k^{\prime}\,S_{i}(k^{\prime},\mu_{S})\,U_{S}^{i}(k-k^{\prime},\mu_{S},\mu)\,,$
(2.13)

where the evolution factor $U_{S}^{i}(k,\mu_{S},\mu)$ encodes the solution of
eq. (2.2), with $U_{S}^{i}(k,\mu_{S},\mu_{S})=\delta(k)$. It can be found e.g.
in refs. [87, 88]. Formally, the $\mu_{S}$ dependence on the right-hand side
cancels, but when the starting condition $S_{i}(k,\mu_{S})$ is evaluated at
fixed order, it only cancels up to higher-order terms. For ease of
presentation, we consider the cumulant of the soft function

$(S_{i}\otimes
U_{S}^{i})_{\mathrm{cut}}(\mathcal{T}_{\mathrm{cut}},\mu)=\int^{\mathcal{T}_{\mathrm{cut}}}\\!\mathrm{d}k\,\int\\!\mathrm{d}k^{\prime}\,S_{i}(k^{\prime},\mu_{S})\,U_{S}^{i}(k-k^{\prime},\mu_{S},\mu)\,,$
(2.14)

for which the distributions turn into ordinary logarithms of
$\mathcal{T}_{\mathrm{cut}}$.

Figure -1139: Residual scale dependence of the integrated resummed
$\mathcal{T}_{0}$ soft function for $i=q$ (left) and $i=g$ (right). Shown are
the relative deviations from the NNLL′ result
$S_{i,\mathrm{cut}}^{\mathrm{central}}$ at the central scale
$\mu_{S}=\mathcal{T}_{\mathrm{cut}}$.

In figure -1139, we take as an example
$\mathcal{T}_{\mathrm{cut}}=\mu=20\,\mathrm{GeV}$ and show the residual
$\mu_{S}$ dependence of the resummed soft function when varying $\mu_{S}$
around the canonical central value $\mu_{S}=\mathcal{T}_{\mathrm{cut}}$ at
NLL′ (dotted green), NNLL′ (dashed blue), and N3LL′ (solid orange). For the
latter we set the currently unknown three-loop constant term $s_{i}^{(3)}=0$.
In all cases, we show the relative difference to the central value at NNLL′.
For simplicity we always use the same four-loop (N3LL) running for
$\alpha_{s}$, which formally amounts to a higher-order effect at (N)NLL′. For
the quark soft function (left panel), the $\mu_{S}$ dependence is more than
halved going from NLL′ to NNLL′, and roughly halved again at N3LL′. In the
gluon case (right panel), the $\mu_{S}$ dependence is noticeably larger, but
also reduces significantly at each order as it should. Note that the missing
three-loop constant term will add an additional source of $\mu_{S}$ dependence
due to its $\alpha_{s}^{3}(\mu_{S})$ prefactor, which however should not
change the general picture.

We stress that the residual $\mu_{S}$ dependence in the resummed soft function
by itself is not necessarily a good indicator of the perturbative uncertainty.
Nevertheless, the reduction in the scale dependence still provides a useful
cross check and an indication of the typical reduction of perturbative
uncertainties one might expect at each order. We also emphasize that the size
of the variations in figure -1139 does not necessarily reflect the variations
one should expect in the resummed cross section, where the evolution of the
soft function happens in conjunction with the beam and hard functions.

###  2.3 $\mathcal{T}_{0}$ beam function

The beam function $B_{i}(t,x,\mu)$ obeys the all-order RGE [49, 87]

$\displaystyle\mu\frac{\mathrm{d}}{\mathrm{d}\mu}B_{i}(t,x,\mu)$
$\displaystyle=\int\\!\mathrm{d}t^{\prime}\,\gamma_{B}^{i}(t-t^{\prime},\mu)\,B_{i}(t^{\prime},x,\mu)\,,$
$\displaystyle\gamma_{B}^{i}(t,\mu)$
$\displaystyle=-2\Gamma_{\mathrm{cusp}}^{i}[\alpha_{s}(\mu)]\,\mathcal{L}_{0}(t,\mu^{2})+\gamma_{B}^{i}[\alpha_{s}(\mu)]\,\delta(t)\,,$
(2.15)

where $\Gamma_{\mathrm{cusp}}^{i}(\alpha_{s})$ and
$\gamma^{i}_{B}(\alpha_{s})$ are the cusp and beam noncusp anomalous
dimensions. For $t\gg\Lambda_{\mathrm{QCD}}$, the beam function satisfies an
OPE in terms of standard PDFs [49, 87]

$B_{i}(t,x,\mu)=\sum_{j}\int\\!\frac{\mathrm{d}z}{z}\,\mathcal{I}_{ij}(t,z,\mu)\,f_{j}\Bigl{(}\frac{x}{z},\mu\Bigr{)}\,\biggl{[}1+\mathcal{O}\biggl{(}\frac{\Lambda_{\mathrm{QCD}}^{2}}{t}\biggr{)}\biggr{]}\,,$
(2.16)

where the $\mathcal{I}_{ij}(t,z,\mu)$ are perturbatively calculable matching
coefficients. Taking into account the evolution of the PDFs, they obey the RGE
[87]

$\mu\frac{\mathrm{d}}{\mathrm{d}\mu}\mathcal{I}_{ij}(t,z,\mu)=\sum_{k}\int\\!\mathrm{d}t^{\prime}\,\frac{\mathrm{d}z^{\prime}}{z^{\prime}}\,\mathcal{I}_{ik}\Bigl{(}t-t^{\prime},\frac{z}{z^{\prime}},\mu\Bigr{)}\Bigl{[}\gamma_{B}^{i}(t^{\prime},\mu)\,\mathbf{1}_{kj}(z^{\prime})-\delta(t^{\prime})\,2P_{kj}(z^{\prime},\mu)\Bigr{]},$
(2.17)

where $\mathbf{1}_{ij}(z)\equiv\delta_{ij}\,\delta(1-z)$ and $2P_{ij}(z,\mu)$
are the PDF anomalous dimensions.

By solving the RGE in eq. (2.17) recursively order by order, we can derive the
complete structure of $\mathcal{I}_{ij}(t,z,\mu)$ at any given fixed order, as
was done in refs. [88, 85] to NNLO. Following the same procedure as in section
2.2, keeping track of the flavor indices and Mellin convolutions, the
$(n+1)$-loop term is determined from the up to $n$-loop terms as

$\displaystyle\mathcal{I}_{ij}^{(n+1)}(t,z,\mu^{2})$
$\displaystyle=\delta(t)\,I_{ij}^{(n+1)}(z)+\int_{t|_{+}}^{\mu^{2}}\\!\frac{\mathrm{d}\tilde{\mu}^{2}}{\tilde{\mu}^{2}}\sum_{m=0}^{n}\biggl{\\{}-\Gamma^{i}_{n-m}\bigl{[}\mathcal{L}_{0}\mathcal{I}^{(m)}_{ij}(z)\bigr{]}(t,\tilde{\mu}^{2})$
(2.18)
$\displaystyle\quad+\Bigl{(}m\beta_{n-m}+\frac{\gamma_{B\,n-m}^{i}}{2}\Bigr{)}\mathcal{I}^{(m)}_{ij}(t,z,\tilde{\mu}^{2})-\bigl{[}\mathcal{I}^{(m)}(t,\tilde{\mu}^{2})P^{(n-m)}\bigr{]}_{ij}(z)\biggr{\\}}\,,$

where the $\mu$-independent boundary coefficients are defined as

$\mathcal{I}^{(n)}_{ij}(t,z,\mu^{2}=t|_{+})=\delta(t)\,I^{(n)}_{ij}(z)\,.$
(2.19)

Starting from the LO result,
$I^{(0)}_{ij}(z)=\mathbf{1}_{ij}(z)\equiv\delta_{ij}\,\delta(1-z)$, we obtain
up to two loops

$\displaystyle\mathcal{I}_{ij}^{(0)}(t,z,\mu^{2})$
$\displaystyle=\delta(t)\,\mathbf{1}_{ij}(z)\,,$
$\displaystyle\mathcal{I}_{ij}^{(1)}(t,z,\mu^{2})$
$\displaystyle=\mathcal{L}_{1}(t,\mu^{2})\,\Gamma_{0}^{i}\,\mathbf{1}_{ij}(z)+\mathcal{L}_{0}(t,\mu^{2})\Bigl{[}-\frac{\gamma_{B\,0}^{i}}{2}\,\mathbf{1}_{ij}(z)+P^{(0)}_{ij}(z)\Bigr{]}+\delta(t)\,I^{(1)}_{ij}(z)\,,$
$\displaystyle\mathcal{I}_{ij}^{(2)}(t,z,\mu^{2})$
$\displaystyle=\mathcal{L}_{3}(t,\mu^{2})\,\frac{(\Gamma_{0}^{i})^{2}}{2}\,\mathbf{1}_{ij}(z)$
$\displaystyle\quad+\mathcal{L}_{2}(t,\mu^{2})\,\frac{\Gamma_{0}^{i}}{2}\Bigl{[}-\Bigl{(}\beta_{0}+\frac{3}{2}\gamma_{B\,0}^{i}\Bigr{)}\mathbf{1}_{ij}(z)+3P^{(0)}_{ij}(z)\Bigr{]}$
$\displaystyle\quad+\mathcal{L}_{1}(t,\mu^{2})\biggl{\\{}\Bigl{[}-\zeta_{2}(\Gamma_{0}^{i})^{2}+\Bigl{(}\beta_{0}+\frac{\gamma_{B\,0}^{i}}{2}\Bigr{)}\frac{\gamma_{B\,0}^{i}}{2}+\Gamma_{1}^{i}\Bigr{]}\mathbf{1}_{ij}(z)$
$\displaystyle\qquad\qquad\qquad\quad-(\beta_{0}+\gamma_{B\,0}^{i})P^{(0)}_{ij}(z)+(P^{(0)}\\!P^{(0)})_{ij}(z)+\Gamma_{0}^{i}\,I^{(1)}_{ij}(z)\biggr{\\}}$
$\displaystyle\quad+\mathcal{L}_{0}(t,\mu^{2})\biggl{\\{}\Bigl{[}\Gamma_{0}^{i}\Bigl{(}\zeta_{3}\Gamma_{0}^{i}+\zeta_{2}\frac{\gamma_{B\,0}^{i}}{2}\Bigr{)}-\frac{\gamma_{B\,1}^{i}}{2}\Bigr{]}\mathbf{1}_{ij}(z)-\zeta_{2}\Gamma_{0}^{i}\,P^{(0)}_{ij}(z)$
$\displaystyle\qquad\qquad\qquad\quad+P^{(1)}_{ij}(z)-\Bigl{(}\beta_{0}+\frac{\gamma_{B\,0}^{i}}{2}\Bigr{)}I^{(1)}_{ij}(z)+(I^{(1)}\\!P^{(0)})_{ij}(z)\biggr{\\}}$
$\displaystyle\quad+\delta(t)\,I^{(2)}_{ij}(z)\,,$ (2.20)

which agrees with refs. [88, 85, 39]. The NLO and NNLO boundary coefficients
$I_{ij}^{(1,2)}(z)$ together with the required Mellin convolutions
$(P^{(0)}P^{(0)})_{ij}(z)$ and $(I^{(1)}P^{(0)})_{ij}(z)$ can be found in
refs. [85, 86].111We caution that the functions $P_{ij}(z)$ and $I_{ij}(z)$ in
refs. [85, 86, 39] are expanded in powers of $\alpha_{s}/(2\pi)$ while here we
expand them in powers of $\alpha_{s}/(4\pi)$.

Plugging eq. (2.3) back into eq. (2.18), we obtain the N3LO result

$\mathcal{I}_{ij}^{(3)}(t,z,\mu^{2})=\delta(t)\,I_{ij}^{(3)}(z)+\sum_{\ell=0}^{5}\mathcal{I}_{ij,\ell}^{(3)}(z)\,\mathcal{L}_{\ell}(t,\mu^{2})$
(2.21)

with the coefficients

$\displaystyle\mathcal{I}_{ij,5}^{(3)}(z)$
$\displaystyle=\frac{(\Gamma^{i}_{0})^{3}}{8}\mathbf{1}_{ij}(z)$
$\displaystyle\mathcal{I}_{ij,4}^{(3)}(z)$
$\displaystyle=\frac{5}{8}(\Gamma^{i}_{0})^{2}\Bigl{[}-\Bigl{(}\frac{2}{3}\beta_{0}+\frac{\gamma_{B\,0}^{i}}{2}\Bigr{)}\mathbf{1}_{ij}(z)+P_{ij}^{(0)}\Bigr{]}$
$\displaystyle\mathcal{I}_{ij,3}^{(3)}(z)$
$\displaystyle=\Gamma^{i}_{0}\biggl{\\{}\Bigl{[}-\zeta_{2}(\Gamma^{i}_{0})^{2}+\frac{\beta_{0}^{2}}{3}+\Bigl{(}\frac{5}{3}\beta_{0}+\frac{\gamma_{B\,0}^{i}}{2}\Bigr{)}\frac{\gamma_{B\,0}^{i}}{2}+\Gamma_{1}^{i}\Bigr{]}\mathbf{1}_{ij}(z)$
$\displaystyle\quad-\Bigl{(}\frac{5}{3}\beta_{0}+\gamma_{B\,0}^{i}\Bigr{)}P_{ij}^{(0)}(z)+(P^{(0)}\\!P^{(0)})_{ij}(z)+\frac{\Gamma_{0}^{i}}{2}I_{ij}^{(1)}(z)\biggr{\\}}$
$\displaystyle\mathcal{I}_{ij,2}^{(3)}(z)$
$\displaystyle=\biggl{\\{}(\Gamma^{i}_{0})^{2}\Bigl{[}\frac{5}{2}\zeta_{3}\Gamma^{i}_{0}+\frac{3}{2}\zeta_{2}(\beta_{0}+\gamma^{i}_{B\,0})\Bigr{]}-\Bigl{(}\beta_{0}+\frac{3}{4}\gamma^{i}_{B\,0}\Bigr{)}\Bigl{(}\beta_{0}\frac{\gamma^{i}_{B\,0}}{2}+\Gamma^{i}_{1}\Bigr{)}-\frac{(\gamma^{i}_{B\,0})^{3}}{16}$
$\displaystyle\quad-\frac{\Gamma^{i}_{0}}{2}\Bigl{(}\beta_{1}+\frac{3}{2}\gamma^{i}_{B\,1}\Bigr{)}\biggr{\\}}\mathbf{1}_{ij}(z)+3\Bigl{[}-\zeta_{2}(\Gamma^{i}_{0})^{2}+\frac{\beta_{0}^{2}}{3}+\Bigl{(}\beta_{0}+\frac{\gamma^{i}_{B\,0}}{4}\Bigr{)}\frac{\gamma^{i}_{B\,0}}{2}+\frac{\Gamma^{i}_{1}}{2}\Bigr{]}P_{ij}^{(0)}(z)$
$\displaystyle\quad-\frac{3}{2}\Bigl{(}\beta_{0}+\frac{\gamma^{i}_{B\,0}}{2}\Bigr{)}(P^{(0)}\\!P^{(0)})_{ij}(z)+\frac{1}{2}(P^{(0)}\\!P^{(0)}\\!P^{(0)})_{ij}(z)$
$\displaystyle\quad+\frac{3}{2}\Gamma^{i}_{0}\Bigl{[}P_{ij}^{(1)}(z)-\Bigl{(}\frac{4}{3}\beta_{0}+\frac{\gamma^{i}_{B\,0}}{2}\Bigr{)}I_{ij}^{(1)}(z)+(I^{(1)}\\!P^{(0)})_{ij}(z)\Bigr{]}$
$\displaystyle\mathcal{I}_{ij,1}^{(3)}(z)$
$\displaystyle=\biggl{\\{}-(\Gamma^{i}_{0})^{2}\Bigl{[}\frac{\zeta_{4}}{2}\Gamma^{i}_{0}+\zeta_{3}(3\beta_{0}+2\gamma^{i}_{B\,0})\Bigr{]}-\zeta_{2}\Gamma^{i}_{0}\Bigl{[}(3\beta_{0}+\gamma^{i}_{B\,0})\frac{\gamma^{i}_{B\,0}}{2}+2\Gamma^{i}_{1}\Bigr{]}$
$\displaystyle\quad+\beta_{0}\gamma^{i}_{B\,1}+\frac{\gamma^{i}_{B\,0}}{2}(\beta_{1}+\gamma^{i}_{B\,1})+\Gamma^{i}_{2}\biggr{\\}}\mathbf{1}_{ij}(z)+\Bigl{\\{}\Gamma^{i}_{0}\bigl{[}4\zeta_{3}\Gamma^{i}_{0}+\zeta_{2}(3\beta_{0}+2\gamma^{i}_{B\,0})\bigr{]}$
$\displaystyle\quad-(\beta_{1}+\gamma^{i}_{B\,1})\Bigr{\\}}P^{(0)}_{ij}(z)-2\zeta_{2}\Gamma^{i}_{0}(P^{(0)}\\!P^{(0)})_{ij}(z)-(2\beta_{0}+\gamma^{i}_{B\,0})P^{(1)}_{ij}(z)$
$\displaystyle\quad+(P^{(0)}\\!P^{(1)}\\!+\\!P^{(1)}\\!P^{(0)})_{ij}(z)+\Bigl{[}-\zeta_{2}(\Gamma^{i}_{0})^{2}+2\beta_{0}^{2}+\Bigl{(}3\beta_{0}+\frac{\gamma^{i}_{B\,0}}{2}\Bigr{)}\frac{\gamma^{i}_{B\,0}}{2}+\Gamma^{i}_{1}\Bigr{]}I^{(1)}_{ij}(z)$
$\displaystyle\quad-(3\beta_{0}+\gamma^{i}_{B\,0})(I^{(1)}\\!P^{(0)})_{ij}(z)+(I^{(1)}\\!P^{(0)}\\!P^{(0)})_{ij}(z)+\Gamma^{i}_{0}\,I^{(2)}_{ij}(z)$
$\displaystyle\mathcal{I}_{ij,0}^{(3)}(z)$
$\displaystyle=\biggl{\\{}(\Gamma^{i}_{0})^{2}\Bigl{[}-\Gamma^{i}_{0}(2\zeta_{2}\zeta_{3}-3\zeta_{5})+\zeta_{4}\Bigl{(}\beta_{0}+\frac{\gamma^{i}_{B\,0}}{4}\Bigr{)}\Bigr{]}+\zeta_{3}\Gamma^{i}_{0}\Bigl{[}\Bigl{(}\beta_{0}+\frac{\gamma^{i}_{B\,0}}{2}\Bigr{)}\frac{\gamma^{i}_{B\,0}}{2}+2\Gamma^{i}_{1}\Bigr{]}$
$\displaystyle\quad+\frac{\zeta_{2}}{2}\bigl{(}\gamma^{i}_{B\,0}\Gamma^{i}_{1}+\Gamma^{i}_{0}\gamma^{i}_{B\,1}\bigr{)}-\frac{\gamma^{i}_{B\,2}}{2}\biggr{\\}}\mathbf{1}_{ij}(z)-\biggl{\\{}\Gamma^{i}_{0}\Bigl{[}\frac{\zeta_{4}}{2}\Gamma^{i}_{0}+\zeta_{3}(\beta_{0}+\gamma^{i}_{B\,0})\Bigr{]}$
$\displaystyle\quad+\zeta_{2}\Gamma^{i}_{1}\biggr{\\}}P^{(0)}_{ij}(z)+\Gamma^{i}_{0}\bigl{[}\zeta_{3}(P^{(0)}\\!P^{(0)})_{ij}(z)-\zeta_{2}P^{(1)}_{ij}(z)\bigr{]}+P^{(2)}_{ij}(z)$
$\displaystyle\quad+\biggl{\\{}\Gamma^{i}_{0}\Bigl{[}\zeta_{3}\Gamma^{i}_{0}+\zeta_{2}\Bigl{(}\beta_{0}+\frac{\gamma^{i}_{B\,0}}{2}\Bigr{)}\Bigr{]}-\Bigl{(}\beta_{1}+\frac{\gamma^{i}_{B\,1}}{2}\Bigr{)}\biggr{\\}}I^{(1)}_{ij}(z)-\Gamma^{i}_{0}\zeta_{2}(I^{(1)}\\!P^{(0)})_{ij}(z)$
$\displaystyle\quad+(I^{(1)}\\!P^{(1)})_{ij}(z)-\Bigl{(}2\beta_{0}+\frac{\gamma^{i}_{B\,0}}{2}\Bigr{)}I^{(2)}_{ij}(z)+(I^{(2)}\\!P^{(0)})_{ij}(z)\,.$
(2.22)

The required anomalous dimensions and splitting functions up to three loops
are given in appendix D. We have evaluated all Mellin convolutions appearing
in eq. (2.3) using the MT package [108]. To calculate
$(I^{(2)}P^{(0)})_{ij}(z)$ this required employing the identity

$\mathrm{Li}_{3}\Bigl{(}\frac{1}{1+z}\Bigr{)}+\mathrm{Li}_{3}(-z)+\mathrm{Li}_{3}\Bigl{(}\frac{z}{1+z}\Bigr{)}=\zeta_{3}-\zeta_{2}\ln(1+z)-\frac{1}{2}\ln^{2}(1+z)\ln
z+\frac{1}{3}\ln^{3}(1+z)\,.$ (2.23)

#### Numerical impact

As for the soft function above, to illustrate the numerical impact of the
three-loop corrections, we consider the integrated resummed beam function

$\displaystyle(B_{i}\otimes U_{B}^{i})_{\mathrm{cut}}(t_{\mathrm{cut}},x,\mu)$
$\displaystyle=\int^{t_{\mathrm{cut}}}\\!\mathrm{d}t\,\int\\!\mathrm{d}t^{\prime}\,B_{i}(t,x,\mu_{B})\,U_{B}^{i}(t-t^{\prime},\mu_{B},\mu)\,.$
(2.24)

The explicit expression for the beam function evolution kernel
$U_{B}^{i}(t,\mu_{B},\mu)$ can be found in refs. [87, 88].

Figure -1138: Residual scale dependence of the resummed integrated
$\mathcal{T}_{N}$ beam function for $i=d$ (top left), $u$ (top right),
$\bar{d}$ (bottom left), and $g$ (bottom right). Shown are the relative
deviations from the NNLL′ result $B_{i,\mathrm{cut}}^{\mathrm{central}}$ at
the central scale $\mu_{B}=\sqrt{t_{\mathrm{cut}}}$.

In figure -1138, we show the residual $\mu_{B}$ dependence of the resummed
integrated beam function at fixed representative values of $x=10^{-2}$ and
$\sqrt{t_{\mathrm{cut}}}=\mu=30\,\mathrm{GeV}$. We again show the relative
difference to the NNLL′ central result at $\mu_{B}=\sqrt{t_{\mathrm{cut}}}$ at
NLL′ (dotted green), NNLL′ (dashed blue), and N3LL′ with the unknown three-
loop $I_{ij}^{(3)}(z)=0$ (solid orange). We use the MMHT2014nnlo68cl [109]
NNLO PDFs and four-loop running of $\alpha_{s}$ everywhere. These evolution
orders are sufficient to ensure the formal cancellation of the $\mu_{B}$
dependence at N3LL′, while at lower orders they amount to a higher-order
effect. Numerical results to N3LL with PDFs and $\alpha_{s}$ evolution at
corresponding lower orders can be found in ref. [86]. The residual dependence
on $\mu_{B}$ is noticeably reduced by about a factor of two at N3LL′ compared
to NNLL′. The missing three-loop constant terms will again add an additional
source of $\mu_{B}$ dependence due to its $\alpha_{s}^{3}(\mu_{B})$ prefactor
and also the scale dependence of the PDFs, which however should not change the
qualitative picture.

###  2.4 Beam function coefficients in the eikonal limit

We now obtain the beam function coefficients $I_{ij}^{(n)}(z)$ in the $z\to 1$
limit. As was already pointed out and exploited in the NNLO calculation in
refs. [85, 86], the beam function in this limit is effectively determined by a
matrix element of eikonal Wilson lines. Here, we exploit a recently derived
consistency relation [97] that explicitly relates the $I_{ij}^{(n)}(z\to 1)$
to the threshold soft function to all orders in $\alpha_{s}$. Consistency
relations of this kind generically arise from different factorization theorems
that apply in different limits of the same multi-differential cross section.
In particular, a soft or collinear matrix element of several arguments will
refactorize into a product (or convolution) of simpler pieces of fewer
arguments by taking a stronger limit.

We start by defining the color-singlet lightcone momenta $q^{\mp}$ and
corresponding momentum fractions $x_{\mp}$,

$q^{-}\equiv{\bar{n}}\cdot q=\sqrt{Q^{2}+q_{T}^{2}}\,e^{+Y}\,,\qquad
q^{+}\equiv n\cdot q=\sqrt{Q^{2}+q_{T}^{2}}\,e^{-Y}\,,\qquad
x_{\mp}=\frac{q^{\mp}}{E_{\mathrm{cm}}}\,.$ (2.25)

As recently shown in ref. [97], in the generalized threshold limit $x_{-}\to
1$ but generic $x_{+}$, the inclusive color-singlet cross section differential
in $q^{\pm}$ factorizes as

$\frac{\mathrm{d}\sigma}{\mathrm{d}q^{-}\mathrm{d}q^{+}}=\sum_{a,b}H_{ab}(q^{+}q^{-},\mu)\int\\!\mathrm{d}t\,f_{a}^{\mathrm{thr}}\Bigl{[}x_{-}\Bigl{(}1+\frac{t}{q^{+}q^{-}}\Bigr{)},\mu\Bigr{]}B_{b}(t,x_{+},\mu)\,\bigl{[}1+\mathcal{O}(1-x_{-})\bigr{]}\,.$
(2.26)

Here, $H_{ab}$ is the same hard function as in eq. (2.1), and $B_{b}$ is the
same inclusive beam function as in eq. (2.1). The threshold PDF
$f^{\mathrm{thr}}_{a}(x)$ encodes the extraction of parton $a$ from the proton
for $x\to 1$.

On the other hand, in the well-known and stronger soft threshold limit, where
both $x_{-}\to 1$ and $x_{+}\to 1$, the cross section factorizes as [110, 111,
112, 113, 114]

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}q^{-}\mathrm{d}q^{+}}$
$\displaystyle=\sum_{a,b}H_{ab}(q^{+}q^{-},\mu)\int\\!\mathrm{d}k^{-}\mathrm{d}k^{+}\,f^{\mathrm{thr}}_{a}\Bigl{[}x_{-}\Bigl{(}1+\frac{k^{-}}{q^{-}}\Bigr{)},\mu\Bigr{]}\,f^{\mathrm{thr}}_{b}\Bigl{[}x_{+}\Bigl{(}1+\frac{k^{+}}{q^{+}}\Bigr{)},\mu\Bigr{]}$
$\displaystyle\quad\times
S^{\mathrm{thr}}_{i}(k^{-},k^{+},\mu)\,\bigl{[}1+\mathcal{O}(1-x_{-},1-x_{+})\bigr{]}\,.$
(2.27)

The new ingredient is the threshold soft function
$S^{\mathrm{thr}}_{i}(k^{-},k^{+},\mu)$. It describes the process-independent
contribution of soft emissions with total lightcone momenta $k^{+}=n\cdot k$
and $k^{-}={\bar{n}}\cdot k$. It also only depends on the color representation
$c\equiv i=\\{q,g\\}$ of the incoming partons.

The threshold soft function is defined as a vacuum matrix element of Wilson
lines that are invariant under longitudinal boosts, and therefore satisfies
the rescaling property

$\displaystyle
S^{\mathrm{thr}}_{i}(k^{-},k^{+},\mu)=S^{\mathrm{thr}}_{i}(e^{+y}k^{-},e^{-y}k^{+},\mu)\,.$
(2.28)

More specifically, in the context of SCET, the soft function is invariant
under RPI-III transformations [115, 116]. Exploiting this property, the soft
function can be extracted [111, 14, 117, 82, 19] from the soft-virtual limit
of the total color-singlet production cross section
$\mathrm{d}\sigma/\mathrm{d}Q^{2}$, which is known to
$\mathcal{O}(\alpha_{s}^{3})$ [11, 12]. In appendix B, we review this
procedure and give explicit results for
$S^{\mathrm{thr}}_{i}(k^{-},k^{+},\mu)$ to three loops.

The factorization theorems eqs. (2.26) and (2.4) describe the same cross
section and share a number of common ingredients. In particular, only the beam
function depends on $x_{+}$ in eq. (2.26). Further expanding eq. (2.26) in the
limit $x=x_{+}\to 1$, it must reproduce eq. (2.4). As a result, the eikonal
$x\to 1$ limit of the beam function must coincide with the threshold soft
function [97],

$B_{i}(t,x,\mu)=\int\\!\frac{\mathrm{d}k}{\omega}\,S^{\mathrm{thr}}_{i}\Bigl{(}\frac{t}{\omega},k,\mu\Bigr{)}\,f_{i}^{\mathrm{thr}}\Bigl{[}x\Bigl{(}1+\frac{k}{\omega}\Bigr{)},\mu\Bigr{]}\,\bigl{[}1+\mathcal{O}(1-x)\bigr{]}\,.$
(2.29)

Replacing $f^{\mathrm{thr}}_{i}[x(1+1-z)]$ by $f_{i}(x/z)/z$, which is
justified at leading power in $1-z$, yields the corresponding relation for the
matching coefficients [97],

$\mathcal{I}_{ij}(t,z,\mu)=\delta_{ij}\,S_{i}^{\mathrm{thr}}\Bigl{[}\frac{t}{\omega},\omega(1-z),\mu\Bigr{]}\,\bigl{[}1+\mathcal{O}(1-z)\bigr{]}\,.$
(2.30)

This relation captures all terms in $\mathcal{I}_{ij}(t,z,\mu)$ that are
singular for $z\to 1$, while power corrections have at most an integrable
singularity for $z\to 1$. Notably, the beam function becomes flavor diagonal
as $z\to 1$, while offdiagonal channels are $\mathcal{O}(1-z)$ suppressed. By
eq. (2.30), the matching coefficients also inherit the rescaling property in
eq. (2.28), i.e., in the limit $z\to 1$, they become invariant under a
simultaneous rescaling $t\mapsto e^{+y}t$ and $1-z\mapsto e^{-y}(1-z)$. In
other words, they are symmetric in $t/\omega$ and $\omega(1-z)$ such that the
dependence on $\omega$ cancels on the right-hand side.

In ref. [97], eq. (2.30) was explicitly confirmed at two loops by comparison
to refs. [85, 86]. We now use it to predict the beam function coefficients in
the eikonal limit at three loops. They are given by the coefficient of
$\delta(k^{-})$ in the threshold soft function upon identifying
$\delta(k^{+})\mapsto\delta(1-z)$ and
$\mathcal{L}_{n}(k^{+},\mu)\mapsto\mathcal{L}_{n}(1-z)$. Including the one-
loop and two-loop results for reference, we find

$\displaystyle I_{ij}^{(1)}(z)$
$\displaystyle=\delta_{ij}\Bigl{[}\mathcal{L}_{1}(1-z)\,\Gamma_{0}^{i}+\delta(1-z)\,s^{\mathrm{thr}{(1)}}_{i}\Bigr{]}+\mathcal{O}\bigl{[}(1-z)^{0}\bigr{]}\,,$
$\displaystyle I_{ij}^{(2)}(z)$
$\displaystyle=\delta_{ij}\biggl{\\{}\mathcal{L}_{3}(1-z)\,\frac{(\Gamma_{0}^{i})^{2}}{2}-\mathcal{L}_{2}(1-z)\,\frac{\Gamma_{0}^{i}}{2}\beta_{0}+\mathcal{L}_{1}(1-z)\Bigl{[}-2\zeta_{2}(\Gamma_{0}^{i})^{2}+\Gamma_{1}^{i}+\Gamma_{0}^{i}\,s^{\mathrm{thr}{(1)}}_{i}\Bigr{]}$
$\displaystyle\quad\qquad+\mathcal{L}_{0}(1-z)\Bigl{[}2\zeta_{3}(\Gamma_{0}^{i})^{2}+\frac{\gamma_{S\,1}^{i}}{2}-\beta_{0}\,s^{\mathrm{thr}{(1)}}_{i}\Bigr{]}+\delta(1-z)\,s^{\mathrm{thr}{(2)}}_{i}\biggr{\\}}+\mathcal{O}\bigl{[}(1-z)^{0}\bigr{]}\,,$
$\displaystyle I_{ij}^{(3)}(z)$
$\displaystyle=\delta_{ij}\biggl{\\{}\mathcal{L}_{5}(1-z)\,\frac{(\Gamma^{i}_{0})^{3}}{8}-\mathcal{L}_{4}(1-z)\,\frac{5}{12}(\Gamma^{i}_{0})^{2}\beta_{0}$
$\displaystyle\quad\qquad+\mathcal{L}_{3}(1-z)\,\Gamma^{i}_{0}\Bigl{[}-2\zeta_{2}(\Gamma^{i}_{0})^{2}+\frac{\beta_{0}^{2}}{3}+\Gamma_{1}^{i}+\frac{\Gamma_{0}^{i}}{2}s_{i}^{\mathrm{thr}{(1)}}\Bigr{]}$
$\displaystyle\quad\qquad+\mathcal{L}_{2}(1-z)\,\Bigl{[}(\Gamma^{i}_{0})^{2}(5\zeta_{3}\Gamma^{i}_{0}+3\zeta_{2}\beta_{0})-\beta_{0}\Gamma^{i}_{1}-\frac{\Gamma^{i}_{0}}{2}\Bigl{(}\beta_{1}-\frac{3}{2}\gamma^{i}_{S\,1}+4\beta_{0}\,s_{i}^{\mathrm{thr}{(1)}}\Bigr{)}\Bigr{]}$
$\displaystyle\quad\qquad+\mathcal{L}_{1}(1-z)\,\Bigl{[}(\Gamma^{i}_{0})^{2}(4\zeta_{4}\Gamma^{i}_{0}-6\zeta_{3}\beta_{0})-4\zeta_{2}\Gamma^{i}_{0}\Gamma^{i}_{1}-\beta_{0}\gamma^{i}_{S\,1}+\Gamma^{i}_{2}$
$\displaystyle\quad\qquad\hskip
55.97205pt+\bigl{(}-2\zeta_{2}(\Gamma^{i}_{0})^{2}+2\beta_{0}^{2}+\Gamma^{i}_{1}\bigr{)}s_{i}^{\mathrm{thr}{(1)}}+\Gamma^{i}_{0}\,s_{i}^{\mathrm{thr}{(2)}}\Bigr{]}$
$\displaystyle\quad\qquad+\mathcal{L}_{0}(1-z)\,\Bigl{[}(\Gamma^{i}_{0})^{2}\bigl{(}-\Gamma^{i}_{0}(8\zeta_{2}\zeta_{3}-6\zeta_{5})+2\zeta_{4}\beta_{0}\bigr{)}+\Gamma^{i}_{0}(4\zeta_{3}\Gamma^{i}_{1}-\zeta_{2}\gamma^{i}_{S\,1})+\frac{\gamma^{i}_{S\,2}}{2}$
$\displaystyle\quad\qquad\hskip
55.97205pt+\Bigl{(}(\Gamma^{i}_{0})^{2}2\zeta_{3}+\Gamma^{i}_{0}2\zeta_{2}\beta_{0}-\beta_{1}+\frac{\gamma^{i}_{S\,1}}{2}\Bigr{)}s_{i}^{\mathrm{thr}{(1)}}-2\beta_{0}\,s_{i}^{\mathrm{thr}{(2)}}\Bigr{]}$
$\displaystyle\quad\qquad+\delta(1-z)\,s^{\mathrm{thr}{(3)}}_{i}\biggr{\\}}+\mathcal{O}\bigl{[}(1-z)^{0}\bigr{]}\,.$
(2.31)

The boundary coefficients $s^{\mathrm{thr}(n)}_{i}$ of the threshold soft
function are given in eq. (B.2). We have exploited that the noncusp anomalous
dimension of the threshold soft function is given by
$-\gamma_{S}^{i}(\alpha_{s})$, see appendix B.2. For brevity, we also used
that $\gamma_{S\,0}^{i}=0$. The result for generic $\gamma_{S\,0}^{i}$ can be
read off from the full expression for the threshold soft function in eq.
(B.2).

The three-loop result in eq. (2.4) is new and a genuine prediction of the
consistency relation in eq. (2.30). We stress that the information provided by
it goes beyond the RGE predicted three-loop structure in eq. (2.3). The fact
that the leading $z\to 1$ terms must be symmetric in $t/\omega$ and
$\omega(1-z)$ allows one to directly determine (or check) the
$\delta(t)\mathcal{L}_{n}(1-z)$ terms from the RGE-predicted
$\mathcal{L}_{n}(t)\delta(1-z)$ terms, which was already noted in refs. [118,
85]. However, the $\delta(t)\delta(1-z)$ coefficient cannot be predicted in
this way, and eq. (2.30) explicitly identifies it with the threshold soft
function coefficients $s^{\mathrm{thr}{(3)}}_{i}$.

As was shown in ref. [97], a factorization theorem analogous to eq. (2.26)
also holds for the inclusive cross section differential in $Q$ and $Y$, with
$B_{i}$ replaced by a closely related, modified beam function
$\tilde{B}_{i}(t,x,\mu)$.222Not to be confused with the $q_{T}$ beam function
$\tilde{B}_{i}(x,{\vec{b}}_{T},\mu,\nu)$ in the following section. Note that
in contrast to eqs. (1.3) and (2.1), here the difference between $q^{\pm}$ and
$(Q,Y)$ matters. The RGE for $\tilde{B}_{i}(t,x,\mu)$ is the same as for
$B_{i}(t,x,\mu)$ in eq. (2.3), and hence eq. (2.3) also holds for
$\tilde{B}_{i}$ just with different boundary coefficients
$\tilde{I}_{ij}^{(n)}(z)$. In the limit $z\to 1$, the difference between
$B_{i}$ and $\tilde{B}_{i}$ becomes power suppressed in $1-z$. As a result,
the $z\to 1$ limit of the modified $\tilde{I}_{i}^{(n)}$ is also given by eq.
(2.4).

###  2.5 Estimating beam function coefficients beyond the eikonal limit

Having the eikonal limit of the beam function coefficients at hand, we can
study to what extent it can be used to approximate the full result and/or
estimate the uncertainty due to the missing terms beyond the eikonal limit.

In figure -1137, we compare the full $\mathcal{T}_{0}$ beam function
coefficient (solid) to its eikonal (LP dotted green) and next-to-eikonal (NLP
dashed blue) expansions at NLO and NNLO for the $u$ quark and gluon channels.
We always show the convolution $(I_{ij}\otimes f_{j})(x)/f_{i}(x)$ with the
appropriate PDF $f_{j}$ and normalize to the PDF $f_{i}(x)$, corresponding to
the LO result, where $i=u$ for the $u$-quark case and $i=g$ for the gluon
case. With this normalization, the shape gives an indication of the rapidity
dependence of the beam function coefficient relative to the LO rapidity
dependence induced by the shape of the PDFs. We also include the appropriate
powers of $\alpha_{s}/(4\pi)$ at each order, so the overall normalization
shows the percent impact relative to the LO result. For definiteness, we
choose $\mu=30\,\,\mathrm{GeV}$ for the scale entering the PDFs and
$\alpha_{s}$.

The eikonal approximation reproduces the correct divergent behavior of the
full flavor-diagonal contributions, denoted as $qqV$ and $gg$, toward large
$x$ but is off away from large $x$. On the other hand, including the next-to-
eikonal terms yields an excellent approximation for all $x$, particularly for
the quark beam function. The rise at very small $x$ for the gluon, which is
not reproduced at NLP, is due to the $z\to 0$ divergent behavior in the gluon
coefficient, which is not reproduced by its $z\to 1$ expansion. If desired, it
can be captured by including the leading $z\to 0$ behavior of the
coefficients, which for simplicity we refrain from doing here. For
illustration, we also show the total contribution from all other corresponding
nondiagonal channels (gray dot-dashed). In each case, they are numerically
subdominant to the flavor-diagonal channel and also much flatter in $x$, since
they only start at NLP.

Figure -1137: Comparison of the full beam function coefficients to their
leading eikonal (LP) and next-to-eikonal (NLP) expansion at NLO (top) and NNLO
(bottom). The $u$-quark channel is shown on the left and the gluon channel on
the right. In each case we also show the sum of all nondiagonal partonic
channels for comparison.

The fact that the NLP result reproduces the full result very well, motivates
us to construct an approximate ansatz for it, which we can then use at three
loops to get a good estimate of the size of the unknown three-loop beam
function coefficient beyond the eikonal limit.

We consider the following ansatz to approximate the coefficient,

$\displaystyle I_{ij,\mathrm{approx}}^{(n)}(z)$
$\displaystyle=I_{ij}^{\mathrm{LP}(n)}(z)+I_{ij,\mathrm{approx}}^{\mathrm{NLP}(n)}(z)+X_{2}\,(1-z)I_{ij,\mathrm{approx}}^{\mathrm{NLP}(n)}(z)\,,$
(2.32)

where the NLP coefficient itself is approximated as

$\displaystyle I_{ij,\mathrm{approx}}^{\mathrm{NLP}(n)}(z)$
$\displaystyle=-(1-z)I_{ij}^{\mathrm{LP}(n)}(z)+X_{1}\,(1-z)\frac{\mathrm{d}}{\mathrm{d}(1-z)}\Bigl{[}(1-z)I_{ij}^{\mathrm{LP}(n)}(z)\Bigr{]}\,,$
(2.33)

and $X_{1}$ and $X_{2}$ are free parameters that can be varied to estimate
residual uncertainties. This ansatz is motivated by the known general
logarithmic structure at NLP

$\displaystyle
I_{ij}^{\mathrm{NLP}(n)}(z)=\sum_{k=0}^{2n-1}c_{n,k}^{\mathrm{NLP}}\ln^{k}(1-z)\,.$
(2.34)

By multiplying the LP term by $(1-z)$ in eq. (2.33), we generate the
appropriate logarithmic structure at NLP. The first term in eq. (2.33)
reproduces the correct NLO and NNLO coefficients for the leading logarithm at
NLP $c_{1,1}^{\mathrm{NLP}}=-4C_{i}$ and $c_{2,3}^{\mathrm{NLP}}=-8C_{i}^{2}$
for both quarks and gluons. Here, the additional double logarithm is
determined by the same power of $(\Gamma_{0}^{i})^{n}$ as at LP, and this
pattern can be expected to hold at higher orders. The second term in eq.
(2.33) generates a next-to-leading logarithmic NLP series. We fix the central
value for $X_{1}=1$ to reproduce the NLL constant term at NLO
$c_{1,0}=4C_{i}$. Interestingly, we find that this choice also reproduces all
NNLO coefficients $c_{2,k}$ very well, typically to within 10%, for both
quarks and gluons and also independently of the choice of $n_{f}$. This
provides a very nontrivial check and so we expect that eq. (2.33) provides a
very good model of the true NLP structure also at higher orders. To estimate
the uncertainties, we vary $X_{1}$ by $\pm 0.5$, which effectively varies the
coefficients of the subleading terms. At NNLO, this variation covers the exact
value for all coefficients. In addition, the last term in eq. (2.32) estimates
the possible effect of terms beyond NLP. Here, we simply take the central
value $X_{2}=0$ and vary it by $\pm 1$.

Since $X_{i}$ probe independent structures, we can consider them as
uncorrelated. Hence, we add the impacts $\Delta_{i}$ on the final result of
their variation in quadrature

$\displaystyle\Delta=\Delta_{1}\oplus\Delta_{2}=\sqrt{\Delta_{1}^{2}+\Delta_{2}^{2}}\,.$
(2.35)

Figure -1136: Approximate ansatzes for the NNLO (top) and N3LO (bottom)
kernels, in the $u$-quark (left) and gluon (right) channels.

In figure -1136, we show the approximate kernel at NNLO (top) and N3LO
(bottom) for the $u$-quark (left) and gluon (right) channels. The dashed
orange line shows the central result from our ansatz and the yellow band its
estimated uncertainty. The gray lines show the impact of the individual
variations of the $X_{i}$ as indicated. In the top panel (NNLO), we also show
the known full two-loop results (red solid). It shows that the ansatz
including uncertainties approximates the true result very well, except for the
gluon at very small $x$ where we do not expect it to hold.

At N3LO, we see that the approximate result gives rise to a sizable shift from
the pure eikonal limit, which we believe to be genuine. Hence, we expect the
full three-loop coefficients to have a nontrivial impact in the one to few
percent range. The uncertainties at N3LO are reduced compared to NNLO as
expected, but are still sizable, which adds motivation for the exact
calculation of the full three-loop coefficients.

##  3 $q_{T}$ factorization to three loops

###  3.1 Factorization theorem

The factorization of the ${\vec{q}}_{T}$ distribution in the limit
$q_{T}\equiv\lvert{\vec{q}}_{T}\rvert\ll Q$ was first established by Collins,
Soper, and Sterman (CSS) [119, 120, 48], and was later elaborated on in refs.
[121, 122, 123, 124]. The factorization for $q_{T}$ was also shown within the
framework of SCET in refs. [125, 126, 127, 117]. Sometimes it is also referred
to as transverse-momentum dependent (TMD) factorization. We write the
factorized cross section as

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}Q^{2}\mathrm{d}Y\mathrm{d}^{2}{\vec{q}}_{T}}$
$\displaystyle=\int\\!\mathrm{d}^{2}{\vec{b}}_{T}\,e^{\mathrm{i}\,{\vec{q}}_{T}\cdot{\vec{b}}_{T}}\,\frac{\mathrm{d}\tilde{\sigma}^{\mathrm{sing}}(\vec{b}_{T})}{\mathrm{d}Q^{2}\mathrm{d}Y}\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{q_{T}^{2}}{Q^{2}}\Bigr{)}\Bigr{]}\,,$
$\displaystyle\frac{\mathrm{d}\tilde{\sigma}^{\mathrm{sing}}(\vec{b}_{T})}{\mathrm{d}Q^{2}\mathrm{d}Y}$
$\displaystyle=\sum_{a,b}H_{ab}(Q^{2},\mu)\tilde{B}_{a}(x_{a},{\vec{b}}_{T},\mu,\nu)\,\tilde{B}_{b}(x_{b},{\vec{b}}_{T},\mu,\nu)\,\tilde{S}_{i}(b_{T},\mu,\nu)\,.$
(3.1)

It receives power corrections suppressed by $q_{T}^{2}/Q^{2}$ as indicated. As
is common, we consider the factorized singular cross section in Fourier-
conjugate ${\vec{b}}_{T}$ space, where convolutions in ${\vec{q}}_{T}$ space
turn into simple products. In particular, Fourier transforming the
${\vec{q}}_{T}$-dependent plus distributions
$\mathcal{L}_{n}({\vec{q}}_{T},\mu)$ turns them into powers of the canonical
${\vec{b}}_{T}$-space logarithms, which we denote as

$L_{b}=\ln\frac{b_{T}^{2}\mu^{2}}{b_{0}^{2}}\,,\qquad
b_{0}=2e^{-\gamma_{E}}\,.$ (3.2)

More details on their Fourier transformation are given in appendix A.2.

The $q_{T}$ factorization is affected by rapidity divergences that must be
regulated by a dedicated rapidity regulator. This gives rise to an additional
rapidity scale, denoted as $\nu$ in eq. (3.1). We use the exponential
regulator of ref. [117], which up to two loops gives results equivalent to the
$\eta$ regulator of refs. [128, 127].

The beam function appearing in eq. (3.1) is the inclusive transverse-momentum
dependent (SCETII) beam function, which also appears in the $q_{T}$
factorization of $Z+j$ and $\gamma+j$ [129, 130]. The $q_{T}$-dependent soft
function in eq. (3.1) is the renormalized vacuum matrix element of two
incoming soft Wilson lines. Note that for simplicity, we generically refer to
them as $q_{T}$ beam and soft functions, even though we mostly consider their
$b_{T}$-dependent Fourier conjugates. The $q_{T}$ beam and soft functions are
known at two loops for several regulators [131, 132, 133, 134, 135, 136, 137,
138]. The soft function is known at three loops using the exponential
regulator [82].

We also note that one can equivalently define $\nu$-independent TMDPDFs as

$\tilde{f}_{i}(x,{\vec{b}}_{T},\mu,\zeta)=\tilde{B}_{i}(x,{\vec{b}}_{T},\mu,\nu)\sqrt{\tilde{S}_{i}(b_{T},\mu,\nu)}\,,\qquad\zeta=\omega^{2}=(xE_{\mathrm{cm}})^{2}\,,$
(3.3)

as is done e.g. in refs. [48, 121, 122, 123, 124, 125, 126]. Here, the
Collins-Soper scale [119, 120] $\zeta=\omega^{2}$ is given in terms of the
lightcone momentum $\omega=xP^{-}$ carried by the struck parton.

###  3.2 Rapidity anomalous dimension

The $\nu$ dependence of the beam and soft functions is encoded in their
rapidity RGEs [127],

$\displaystyle\nu\frac{\mathrm{d}}{\mathrm{d}\nu}\tilde{B}_{i}(x,{\vec{b}}_{T},\mu,\nu)$
$\displaystyle=\tilde{\gamma}_{\nu,B}^{i}(b_{T},\mu)\,\tilde{B}_{i}(x,{\vec{b}}_{T},\mu,\nu)\,,$
$\displaystyle\nu\frac{\mathrm{d}}{\mathrm{d}\nu}\tilde{S}_{i}(b_{T},\mu,\nu)$
$\displaystyle=\tilde{\gamma}_{\nu,S}^{i}(b_{T},\mu)\,\tilde{S}_{i}(b_{T},\mu,\nu)\,,$
(3.4)

where $\tilde{\gamma}_{\nu,B}^{i}$ and $\tilde{\gamma}_{\nu,S}^{i}$ are the
beam and soft rapidity anomalous dimensions, which are closely related to the
Collins-Soper kernel [119, 120]. Because the cross section in eq. (3.1) is
independent of $\nu$, they are related by

$\displaystyle\tilde{\gamma}_{\nu}^{i}(b_{T},\mu)\equiv\tilde{\gamma}_{\nu,S}^{i}(b_{T},\mu)=-2\tilde{\gamma}_{\nu,B}^{i}(b_{T},\mu)\,,$
(3.5)

and we will simply refer to $\tilde{\gamma}_{\nu}^{i}(b_{T},\mu)$ as the
rapidity anomalous dimension.

An important property of $\tilde{\gamma}_{\nu}^{i}(b_{T},\mu)$ is that like
the soft function it only depends on the color representation $i=\\{q,g\\}$
but not on the specific massless quark flavor. While we only need its fixed-
order expansion, we note that it becomes genuinely nonperturbative for
$b_{T}^{-1}\lesssim\Lambda_{\rm QCD}$, and recently a proposal was made to
calculate it nonperturbatively using lattice QCD [139, 140].

The rapidity anomalous dimension itself satisfies an RGE in $\mu$,

$\displaystyle\mu\frac{\mathrm{d}}{\mathrm{d}\mu}\tilde{\gamma}^{i}_{\nu}(b_{T},\mu)=-4\Gamma_{\mathrm{cusp}}^{i}[\alpha_{s}(\mu)]\,,$
(3.6)

which predicts its all-order structure in $b_{T}$ and $\mu$. Similar to the
$\mathcal{T}_{0}$ soft function in section 2.2, it can be solved recursively
order by order in $\alpha_{s}$. Expanding both sides of eq. (3.6) to fixed
order in $\alpha_{s}(\mu)$ and accounting for the running of
$\alpha_{s}(\mu)$, the $(n+1)$-loop term is related to the lower-order terms
by

$\tilde{\gamma}_{\nu}^{i\,(n+1)}(b_{T},\mu)=-2\Gamma^{i}_{n+1}L_{b}+\sum_{m=0}^{n}2(m+1)\,\beta_{n-m}\int_{b_{0}/b_{T}}^{\mu}\\!\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\,\tilde{\gamma}_{\nu}^{i\,(m)}(b_{T},\mu^{\prime})+\tilde{\gamma}_{\nu\,n+1}^{i}\,,$
(3.7)

where the nonlogarithmic boundary coefficients are defined as

$\tilde{\gamma}_{\nu\,n}^{i}=\tilde{\gamma}_{\nu}^{i\,(n)}(b_{T},\mu=b_{0}/b_{T})\,.$
(3.8)

The result up to three loops is

$\displaystyle\tilde{\gamma}_{\nu}^{i\,(0)}(b_{T},\mu)$
$\displaystyle=-L_{b}\,2\Gamma_{0}^{i}+\tilde{\gamma}_{\nu\,0}^{i}\,,$
$\displaystyle\tilde{\gamma}_{\nu}^{i\,(1)}(b_{T},\mu)$
$\displaystyle=-L_{b}^{2}\,\Gamma_{0}^{i}\beta_{0}+L_{b}\bigl{(}\beta_{0}\tilde{\gamma}_{\nu\,0}^{i}-2\Gamma_{1}^{i}\bigr{)}+\tilde{\gamma}_{\nu\,1}^{i}\,,$
$\displaystyle\tilde{\gamma}_{\nu}^{i\,(2)}(b_{T},\mu)$
$\displaystyle=-L_{b}^{3}\,\frac{2}{3}\Gamma_{0}^{i}\beta_{0}^{2}+L_{b}^{2}\bigl{(}\beta_{0}^{2}\tilde{\gamma}_{\nu\,0}^{i}-2\Gamma_{1}^{i}\beta_{0}-\Gamma_{0}^{i}\beta_{1}\bigr{)}+L_{b}\bigl{(}2\beta_{0}\tilde{\gamma}_{\nu\,1}^{i}+\beta_{1}\tilde{\gamma}_{\nu\,0}^{i}-2\Gamma_{2}^{i}\bigr{)}$
$\displaystyle\quad+\tilde{\gamma}_{\nu\,2}^{i}\,.$ (3.9)

The boundary coefficients $\tilde{\gamma}_{\nu\,n}^{i}$ are known up to three
loops [136, 82, 141] and are summarized in eq. (D.3).

###  3.3 Soft function

The soft function is explicitly known to three loops [82]. For completeness,
we explicitly derive its fixed-order structure to illustrate the joint
solution of its $\mu$ and $\nu$ RGEs,

$\displaystyle\mu\frac{\mathrm{d}}{\mathrm{d}\mu}\tilde{S}_{i}(b_{T},\mu,\nu)$
$\displaystyle=\tilde{\gamma}_{S}^{i}(\mu,\nu)\,\tilde{S}_{i}(b_{T},\mu,\nu)\,,$
$\displaystyle\nu\frac{\mathrm{d}}{\mathrm{d}\nu}\tilde{S}_{i}(b_{T},\mu,\nu)$
$\displaystyle=\tilde{\gamma}_{\nu}^{i}(b_{T},\mu)\,\tilde{S}_{i}(b_{T},\mu,\nu)\,.$
(3.10)

The perturbative structure of $\gamma_{\nu}^{i}$ is discussed in section 3.2
above. The $\mu$ anomalous dimension has the all-order structure

$\tilde{\gamma}_{S}^{i}(\mu,\nu)=4\Gamma_{\mathrm{cusp}}^{i}[\alpha_{s}(\mu)]\ln\frac{\mu}{\nu}+\tilde{\gamma}_{S}^{i}[\alpha_{s}(\mu)]\,,$
(3.11)

where $\Gamma_{\mathrm{cusp}}^{i}(\alpha_{s})$ and
$\tilde{\gamma}_{S}^{i}(\alpha_{s})$ are the cusp and the soft noncusp
anomalous dimensions. Expanding both sides of eq. (3.3) order by order in
$\alpha_{s}$, we obtain the coupled RGEs

$\displaystyle\mu\frac{\mathrm{d}}{\mathrm{d}\mu}{\tilde{S}_{i}^{(n+1)}(b_{T},\mu,\nu)}$
$\displaystyle=\sum_{m=0}^{n}\Bigl{(}4\Gamma^{i}_{n-m}\ln\frac{\mu}{\nu}+2m\beta_{n-m}+\tilde{\gamma}_{S\,n-m}^{i}\Bigr{)}\tilde{S}_{i}^{(m)}(b_{T},\mu,\nu)\,,$
$\displaystyle\nu\frac{\mathrm{d}}{\mathrm{d}\nu}{\tilde{S}_{i}^{(n+1)}(b_{T},\mu,\nu)}$
$\displaystyle=\sum_{m=0}^{n}\tilde{\gamma}_{\nu}^{i\,(n-m)}(b_{T},\mu)\,\tilde{S}_{i}^{(m)}(b_{T},\mu,\nu)\,.$
(3.12)

These are easily integrated to give

$\displaystyle\tilde{S}_{i}^{(n+1)}(b_{T},\mu,\nu)$
$\displaystyle=\sum_{m=0}^{n}\biggl{[}\int_{b_{0}/b_{T}}^{\mu}\\!\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\,\Bigl{(}4\Gamma^{i}_{n-m}\ln\frac{\mu^{\prime}}{\nu}+2m\beta_{n-m}+\tilde{\gamma}_{S\,n-m}^{i}\Bigr{)}\tilde{S}_{i}^{(m)}(b_{T},\mu^{\prime},\nu)$
$\displaystyle\quad+\int_{b_{0}/b_{T}}^{\nu}\\!\frac{\mathrm{d}\nu^{\prime}}{\nu^{\prime}}\,\tilde{\gamma}_{\nu\,n-m}^{i}\,\tilde{S}_{i}^{(m)}(b_{T},b_{0}/b_{T},\nu^{\prime})\biggr{]}+\tilde{s}_{i}^{(n+1)}\,,$
(3.13)

where we first integrated the $\nu$ RGE at fixed $\mu=b_{0}/b_{T}$ and then
the $\mu$ RGE at arbitrary $\nu$. In this way, the rapidity anomalous
dimension reduces to its boundary coefficients $\gamma_{\nu\,n}^{i}$. The soft
boundary coefficients in eq. (3.3) are defined as

$\tilde{s}_{i}^{(n)}=\tilde{S}_{i}^{(n)}(b_{T},\mu=b_{0}/b_{T},\nu=b_{0}/b_{T})\,.$
(3.14)

Starting from the LO result, $\tilde{s}_{i}^{(0)}=1$, and expressing the
results in terms of

$L_{b}=\ln\frac{b_{T}^{2}\mu^{2}}{b_{0}^{2}}\,,\quad
b_{0}=2e^{-\gamma_{E}}\,,\qquad L_{\nu}=\ln\frac{\mu}{\nu}\,,$ (3.15)

eq. (3.3) yields up to two loops

$\displaystyle\tilde{S}_{i}^{(0)}(b_{T},\mu,\nu)$ $\displaystyle=1\,,$
$\displaystyle\tilde{S}_{i}^{(1)}(b_{T},\mu,\nu)$
$\displaystyle=-L_{b}^{2}\,\frac{\Gamma_{0}^{i}}{2}+L_{b}\Bigl{(}L_{\nu}\,2\Gamma_{0}^{i}+\frac{\tilde{\gamma}_{S\,0}^{i}}{2}+\frac{\tilde{\gamma}_{\nu\,0}^{i}}{2}\Bigr{)}-L_{\nu}\tilde{\gamma}_{\nu\,0}^{i}+\tilde{s}_{i}^{(1)}\,,$
$\displaystyle\tilde{S}_{i}^{(2)}(b_{T},\mu,\nu)$
$\displaystyle=L_{b}^{4}\,\frac{(\Gamma_{0}^{i})^{2}}{8}-L_{b}^{3}\,\Gamma_{0}^{i}\Bigl{(}L_{\nu}\Gamma_{0}^{i}+\frac{\beta_{0}}{3}+\frac{\tilde{\gamma}_{S\,0}^{i}}{4}+\frac{\tilde{\gamma}_{\nu\,0}^{i}}{4}\Bigr{)}$
$\displaystyle\quad+L_{b}^{2}\biggl{[}L_{\nu}^{2}\,2(\Gamma_{0}^{i})^{2}+L_{\nu}\Gamma_{0}^{i}\Bigl{(}\beta_{0}+\tilde{\gamma}_{S\,0}^{i}+\frac{3}{2}\tilde{\gamma}_{\nu\,0}^{i}\Bigr{)}$
$\displaystyle\quad\hskip
25.83325pt+\beta_{0}\Bigl{(}\frac{\tilde{\gamma}_{S\,0}^{i}}{4}+\frac{\tilde{\gamma}_{\nu\,0}^{i}}{2}\Bigr{)}+\frac{1}{8}(\tilde{\gamma}_{S\,0}^{i}+\tilde{\gamma}_{\nu\,0}^{i})^{2}-\frac{\Gamma_{1}^{i}}{2}-\frac{\Gamma_{0}^{i}}{2}\,\tilde{s}_{i}^{(1)}\biggr{]}$
$\displaystyle\quad+L_{b}\biggl{\\{}-L_{\nu}^{2}\,2\Gamma_{0}^{i}\tilde{\gamma}_{\nu\,0}^{i}+L_{\nu}\Bigl{[}-\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{S\,0}^{i}}{2}+\frac{\tilde{\gamma}_{\nu\,0}^{i}}{2}\Bigr{)}\tilde{\gamma}_{\nu\,0}^{i}+2\Gamma_{1}^{i}+2\Gamma_{0}^{i}\,\tilde{s}_{i}^{(1)}\Bigr{]}$
$\displaystyle\quad\hskip
25.83325pt+\frac{\tilde{\gamma}_{S\,1}^{i}}{2}+\frac{\tilde{\gamma}_{\nu\,1}^{i}}{2}+\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{S\,0}^{i}}{2}+\frac{\tilde{\gamma}_{\nu\,0}^{i}}{2}\Bigr{)}\tilde{s}_{i}^{(1)}\biggr{\\}}$
$\displaystyle\quad+L_{\nu}^{2}\,\frac{(\tilde{\gamma}_{\nu\,0}^{i})^{2}}{2}-L_{\nu}\bigl{(}\tilde{\gamma}_{\nu\,1}^{i}+\tilde{\gamma}_{\nu\,0}^{i}\,\tilde{s}_{i}^{(1)}\bigr{)}+\tilde{s}_{i}^{(2)}\,.$
(3.16)

At three loops, we write the result as

$\tilde{S}_{i}^{(3)}(b_{T},\mu,\nu)=\sum_{\ell=0}^{6}\tilde{S}_{i,\ell}^{(3)}(L_{\nu})\,L_{b}^{\ell}\,,$
(3.17)

where the $\tilde{S}_{i,k}^{(3)}$ coefficients themselves are polynomials in
$L_{\nu}$. Inserting $\tilde{\gamma}_{S\,0}^{i}=\tilde{\gamma}_{\nu\,0}^{i}=0$
for brevity, they are given by

$\displaystyle\tilde{S}_{i,6}^{(3)}(L_{\nu})$
$\displaystyle=-\frac{(\Gamma^{i}_{0})^{3}}{48}\,,$
$\displaystyle\tilde{S}_{i,5}^{(3)}(L_{\nu})$
$\displaystyle=L_{\nu}\,\frac{(\Gamma^{i}_{0})^{3}}{4}+\frac{(\Gamma^{i}_{0})^{2}}{6}\beta_{0}\,,$
$\displaystyle\tilde{S}_{i,4}^{(3)}(L_{\nu})$
$\displaystyle=-L_{\nu}^{2}\,(\Gamma^{i}_{0})^{3}-L_{\nu}\,\frac{7}{6}(\Gamma^{i}_{0})^{2}\beta_{0}+\frac{\Gamma_{0}^{i}}{4}\Bigl{(}-\beta_{0}^{2}+\Gamma^{i}_{1}+\frac{\Gamma^{i}_{0}}{2}\,\tilde{s}_{i}^{(1)}\Bigr{)}\,,$
$\displaystyle\tilde{S}_{i,3}^{(3)}(L_{\nu})$
$\displaystyle=L_{\nu}^{3}\,\frac{4}{3}(\Gamma^{i}_{0})^{3}+L_{\nu}^{2}\,2(\Gamma^{i}_{0})^{2}\beta_{0}+L_{\nu}\Gamma_{0}^{i}\Bigl{(}\frac{2}{3}\beta_{0}^{2}-2\Gamma^{i}_{1}-\Gamma^{i}_{0}\,\tilde{s}_{i}^{(1)}\Bigr{)}$
$\displaystyle\quad-\frac{2}{3}\Gamma^{i}_{1}\beta_{0}-\Gamma_{0}^{i}\Bigl{(}\frac{\beta_{1}}{3}+\frac{\tilde{\gamma}_{S\,1}^{i}}{4}+\frac{\tilde{\gamma}_{\nu\,1}^{i}}{4}+\frac{5}{6}\beta_{0}\,\tilde{s}_{i}^{(1)}\Bigr{)}\,,$
$\displaystyle\tilde{S}_{i,2}^{(3)}(L_{\nu})$
$\displaystyle=L_{\nu}^{2}\,\Gamma^{i}_{0}\bigl{(}4\Gamma^{i}_{1}+2\Gamma^{i}_{0}\,\tilde{s}_{i}^{(1)}\bigr{)}+L_{\nu}\Bigl{[}2\Gamma^{i}_{1}\beta_{0}+\Gamma^{i}_{0}\Bigl{(}\beta_{1}+\tilde{\gamma}_{S\,1}^{i}+\frac{3}{2}\tilde{\gamma}_{\nu\,1}^{i}+3\beta_{0}\,\tilde{s}_{i}^{(1)}\Bigr{)}\Bigr{]}$
$\displaystyle\quad+\beta_{0}\Bigl{(}\frac{\tilde{\gamma}_{S\,1}^{i}}{2}+\tilde{\gamma}_{\nu\,1}^{i}\Bigr{)}-\frac{\Gamma^{i}_{2}}{2}+\Bigl{(}\beta_{0}^{2}-\frac{\Gamma^{i}_{1}}{2}\Bigr{)}\tilde{s}_{i}^{(1)}-\frac{\Gamma^{i}_{0}}{2}\,\tilde{s}_{i}^{(2)}\,,$
$\displaystyle\tilde{S}_{i,1}^{(3)}(L_{\nu})$
$\displaystyle=-L_{\nu}^{2}\,2\Gamma^{i}_{0}\tilde{\gamma}_{\nu\,1}^{i}+L_{\nu}\,2\bigl{(}-\beta_{0}\tilde{\gamma}_{\nu\,1}^{i}+\Gamma^{i}_{2}+\Gamma^{i}_{1}\,\tilde{s}_{i}^{(1)}+\Gamma^{i}_{0}\,\tilde{s}_{i}^{(2)}\bigr{)}$
$\displaystyle\quad+\frac{\tilde{\gamma}_{S\,2}^{i}}{2}+\frac{\tilde{\gamma}_{\nu\,2}^{i}}{2}+\Bigl{(}\beta_{1}+\frac{\tilde{\gamma}_{S\,1}^{i}}{2}+\frac{\tilde{\gamma}_{\nu\,1}^{i}}{2}\Bigr{)}\tilde{s}_{i}^{(1)}+2\beta_{0}\,\tilde{s}_{i}^{(2)}\,,$
$\displaystyle\tilde{S}_{i,0}^{(3)}(L_{\nu})$
$\displaystyle=-L_{\nu}\bigl{(}\tilde{\gamma}_{\nu\,2}^{i}+\tilde{\gamma}_{\nu\,1}^{i}\,\tilde{s}_{i}^{(1)}\bigr{)}+\tilde{s}_{i}^{(3)}\,.$
(3.18)

Eqs. (3.3) and (3.3) agree with refs. [136, 82]. The required anomalous
dimension and boundary coefficients up to three loops are given in appendix
D.3.

#### Numerical impact

The soft function $\tilde{S}_{i}(b_{T},\mu,\nu)$ has an explicit dependence on
the scales $\mu$ and $\nu$ that cancels against that of the hard and beam
functions in eq. (3.1). Therefore, varying $\mu$ and $\nu$ is not very
meaningful for illustrating the numerical impact of the scale-dependent three-
loop terms. Instead, we consider the resummed soft function,

$\displaystyle\tilde{S}_{i}(b_{T},\mu,\nu)$
$\displaystyle=\tilde{S}_{i}(b_{T},\mu_{S},\nu_{S})\,\tilde{U}_{S}^{i}(b_{T},\mu_{S},\mu,\nu_{S},\nu)\,,$
$\displaystyle\tilde{U}_{S}^{i}(b_{T},\mu_{S},\mu,\nu_{S},\nu)$
$\displaystyle=\exp\biggl{[}\ln\frac{\nu}{\nu_{S}}\,\tilde{\gamma}_{\nu}^{i}(b_{T},\mu_{S})\biggr{]}\exp\biggl{[}\int_{\mu_{S}}^{\mu}\\!\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\,\tilde{\gamma}_{S}^{i}(\mu^{\prime},\nu)\biggr{]}\,,$
(3.19)

where we have chosen to first evolve in $\nu$ and then in $\mu$.

To probe the full set of terms in the fixed-order expansion of
$\tilde{S}_{i}(b_{T},\mu_{S},\nu_{S})$, we consider simultaneous variations of
$(\mu_{S},\nu_{S})$ around the canonical central scales
$\mu_{S}=\nu_{S}=\mu=\nu=b_{0}/b_{T}$. In figure -1135 we show the residual
scale dependence of the resummed soft function at the representative value
$b_{0}/b_{T}=20\,\mathrm{GeV}$ at NLL′ (dotted green), NNLL′ (dashed blue),
and N3LL′ (solid orange), normalized to the NNLL′ result at the central scale.
The three-loop finite term is included in figure -1135, so the NNLL′ and N3LL′
results do not coincide at the central scales. We use four-loop running of
$\alpha_{s}$ throughout, which formally amounts to a higher-order effect at
(N)NLL′. As expected, the scale dependence reduces from NLL′ to NNLL′, where
it is already quite small. At N3LL′, it further stabilizes over a wider range
of scales. As in section 2.2, we stress that the residual scale dependence in
the resummed soft function by itself is not necessarily a good indicator of
the perturbative uncertainty, but gives an indication of the typical reduction
of perturbative uncertainties one might expect at each order.

Figure -1135: Residual scale dependence of the resummed $q_{T}$ soft function
in Fourier space for $i=q$ (left) and $i=g$ (right). Shown are the relative
deviations from the NNLL′ result $\tilde{S}_{i}^{\mathrm{central}}$ at the
central scales $(\mu_{S},\nu_{S})=(\mu,\nu)=(b_{0}/b_{T},b_{0}/b_{T})$.

###  3.4 Beam function

The beam function obeys the coupled RGEs

$\displaystyle\mu\frac{\mathrm{d}}{\mathrm{d}\mu}{\tilde{B}_{i}(x,{\vec{b}}_{T},\mu,\nu)}$
$\displaystyle=\tilde{\gamma}_{B}^{i}(\mu,\nu/\omega)\,\tilde{B}_{i}(x,{\vec{b}}_{T},\mu,\nu)\,,$
$\displaystyle\nu\frac{\mathrm{d}}{\mathrm{d}\nu}{\tilde{B}_{i}(x,{\vec{b}}_{T},\mu,\nu)}$
$\displaystyle=-\frac{1}{2}\tilde{\gamma}_{\nu}^{i}(b_{T},\mu)\,\tilde{B}_{i}(x,{\vec{b}}_{T},\mu,\nu)\,,$
(3.20)

where the $\nu$ anomalous dimension was discussed in section 3.2, and the
$\mu$ anomalous dimension has the all-order form

$\tilde{\gamma}_{B}^{i}(\mu,\nu/\omega)=2\Gamma_{\mathrm{cusp}}^{i}[\alpha_{s}(\mu)]\ln\frac{\nu}{\omega}+\tilde{\gamma}_{B}^{i}[\alpha_{s}(\mu)]\,,$
(3.21)

where $\Gamma_{\mathrm{cusp}}^{i}(\alpha_{s})$ and
$\tilde{\gamma}_{B}^{i}(\alpha_{s})$ are the cusp and the beam noncusp
anomalous dimensions.

For perturbative $b_{0}/b_{T}\gg\Lambda_{\mathrm{QCD}}$, the TMD beam function
satisfies an OPE in terms of standard PDFs [48],

$\displaystyle\tilde{B}_{q}(x,b_{T},\mu,\nu)$
$\displaystyle=\sum_{j}\int\\!\frac{\mathrm{d}z}{z}\,\tilde{\mathcal{I}}_{qj}(z,b_{T},\mu,\nu/\omega)f_{j}\Bigl{(}\frac{x}{z},\mu\Bigr{)}\bigl{[}1+\mathcal{O}(b_{T}\Lambda_{\mathrm{QCD}})\bigr{]}\,,$
$\displaystyle\tilde{B}^{\rho\lambda}_{g}(x,{\vec{b}}_{T},\mu,\nu)$
$\displaystyle=\sum_{j}\int\\!\frac{\mathrm{d}z}{z}\biggl{[}\frac{g_{\perp}^{\rho\lambda}}{2}\,\tilde{\mathcal{I}}_{gj}(z,b_{T},\mu,\nu/\omega)+\biggl{(}\frac{g_{\perp}^{\rho\lambda}}{2}-\frac{b_{\perp}^{\rho}b_{\perp}^{\lambda}}{b_{\perp}^{2}}\biggr{)}\tilde{\mathcal{J}}_{gj}(z,b_{T},\mu,\nu/\omega)\biggr{]}$
$\displaystyle\quad\hskip 43.05542pt\times
f_{j}\Bigl{(}\frac{x}{z},\mu\Bigr{)}\bigl{[}1+\mathcal{O}(b_{T}\Lambda_{\mathrm{QCD}})\bigr{]}\,.$
(3.22)

For the gluon beam function, we have made its dependence on the gluon helicity
explicit, and decomposed it into two orthogonal structures, namely the
polarization-independent piece $\tilde{\mathcal{I}}_{gj}$ and the
polarization-dependent piece $\tilde{\mathcal{J}}_{gj}$, where
$g_{\perp}^{\rho\lambda}=g^{\rho\lambda}-(n^{\rho}{\bar{n}}^{\lambda}+{\bar{n}}^{\rho}n^{\lambda})/2$
is the transverse metric and $b_{\perp}^{\rho}$ is the transverse four vector
with $b_{\perp}^{2}=-{\vec{b}}_{T}^{2}$. Due to the multiplicative structure
of eq. (3.4), both $\tilde{\mathcal{I}}_{gj}$ and $\tilde{\mathcal{J}}_{gj}$
obey the same RGE, and in the following we will only consider the RGEs for
$\tilde{\mathcal{I}}_{ij}$.

The $\tilde{\mathcal{I}}_{ij}$ are perturbatively calculable matching
coefficients, whose RGEs follow from eq. (3.4) by taking the evolution of the
PDFs into account,

$\displaystyle\mu\frac{\mathrm{d}}{\mathrm{d}\mu}{\tilde{\mathcal{I}}_{ij}(z,b_{T},\mu,\nu/\omega)}$
$\displaystyle=\sum_{k}\int\\!\frac{\mathrm{d}z^{\prime}}{z^{\prime}}\,\tilde{\mathcal{I}}_{ik}\Bigl{(}\frac{z}{z^{\prime}},b_{T},\mu,\nu/\omega\Bigr{)}\bigl{[}\tilde{\gamma}_{B}^{i}(\mu,\nu/\omega)\mathbf{1}_{kj}(z^{\prime})-2P_{kj}(z^{\prime},\mu)\bigr{]}\,,$
$\displaystyle\nu\frac{\mathrm{d}}{\mathrm{d}\nu}{\tilde{\mathcal{I}}_{ij}(z,b_{T},\mu,\nu/\omega)}$
$\displaystyle=-\frac{1}{2}\tilde{\gamma}_{\nu}^{i}(b_{T},\mu)\,\tilde{\mathcal{I}}_{ij}(z,b_{T},\mu,\nu/\omega)\,.$
(3.23)

Similar to the soft function, these coupled RGEs can be solved recursively as

$\displaystyle\tilde{\mathcal{I}}_{ij}^{(n+1)}(z,b_{T},\mu,\nu/\omega)$
$\displaystyle=\sum_{m=0}^{n}\biggl{[}\int_{b_{0}/b_{T}}^{\mu}\\!\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\Bigl{(}2\Gamma^{i}_{n-m}\ln\frac{\nu}{\omega}+\tilde{\gamma}_{B\,n-m}^{i}+2m\beta_{n-m}\Bigr{)}\tilde{\mathcal{I}}_{ij}^{(m)}(z,b_{T},\mu^{\prime},\nu/\omega)$
$\displaystyle\quad-2\int_{b_{0}/b_{T}}^{\mu}\\!\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\,\bigl{[}\tilde{\mathcal{I}}^{(m)}(b_{T},\mu^{\prime},\nu/\omega)P^{(n-m)}\bigr{]}_{ij}(z)$
$\displaystyle\quad-\int_{\omega}^{\nu}\\!\frac{\mathrm{d}\nu^{\prime}}{\nu^{\prime}}\,\frac{\tilde{\gamma}_{\nu\,n-m}^{i}}{2}\,\tilde{\mathcal{I}}_{ij}^{(m)}(z,b_{T},b_{0}/b_{T},\nu^{\prime}/\omega)\biggr{]}+\tilde{I}_{ij}^{(n+1)}(z)\,,\\!$
(3.24)

where the nonlogarithmic boundary coefficients are defined as

$\displaystyle\tilde{I}_{ij}^{(n+1)}(z)=\tilde{\mathcal{I}}_{ij}^{(n+1)}(z,b_{T},\mu=b_{0}/b_{T},\nu/\omega=1)\,.$
(3.25)

Starting from the LO result,
$\tilde{I}_{ij}^{(0)}(z)=\mathbf{1}_{ij}(z)\equiv\delta_{ij}\,\delta(1-z)$, we
obtain up to two loops

$\displaystyle\tilde{\mathcal{I}}_{ij}^{(0)}(z,b_{T},\mu,\nu/\omega)$
$\displaystyle=\mathbf{1}_{ij}(z)\,,$
$\displaystyle\tilde{\mathcal{I}}_{ij}^{(1)}(z,b_{T},\mu,\nu/\omega)$
$\displaystyle=L_{b}\Bigl{[}\Bigl{(}L_{\omega}\Gamma^{i}_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{)}\mathbf{1}_{ij}(z)-P^{(0)}_{ij}(z)\Bigr{]}-L_{\omega}\frac{\tilde{\gamma}_{\nu\,0}^{i}}{2}\,\mathbf{1}_{ij}(z)+\tilde{I}_{ij}^{(1)}(z)\,,$
$\displaystyle\tilde{\mathcal{I}}_{ij}^{(2)}(z,b_{T},\mu,\nu/\omega)$
$\displaystyle=L_{b}^{2}\biggl{\\{}\Bigl{[}L_{\omega}^{2}\,\frac{(\Gamma_{0}^{i})^{2}}{2}+L_{\omega}\,\frac{\Gamma_{0}^{i}}{2}(\beta_{0}+\tilde{\gamma}_{B\,0}^{i})+\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{)}\frac{\tilde{\gamma}_{B\,0}^{i}}{4}\Bigr{]}\mathbf{1}_{ij}(z)$
$\displaystyle\qquad\quad-\Bigl{(}L_{\omega}\Gamma_{0}^{i}+\frac{\beta_{0}}{2}+\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{)}P^{(0)}_{ij}(z)+\frac{1}{2}(P^{(0)}\\!P^{(0)})_{ij}(z)\biggr{\\}}$
$\displaystyle\quad+L_{b}\biggl{\\{}\Bigl{[}-L_{\omega}^{2}\,\Gamma_{0}^{i}\frac{\tilde{\gamma}_{\nu\,0}^{i}}{2}+L_{\omega}\Bigl{[}-\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{)}\frac{\tilde{\gamma}_{\nu\,0}^{i}}{2}+\Gamma_{1}^{i}\Bigr{]}+\frac{\tilde{\gamma}_{B\,1}^{i}}{2}\Bigr{]}\mathbf{1}_{ij}(z)$
$\displaystyle\qquad\qquad+L_{\omega}\,\frac{\tilde{\gamma}_{\nu\,0}^{i}}{2}\,P^{(0)}_{ij}(z)-P^{(1)}_{ij}(z)$
$\displaystyle\qquad\qquad+\Bigl{(}L_{\omega}\Gamma_{0}^{i}+\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{)}\tilde{I}^{(1)}_{ij}(z)-(\tilde{I}^{(1)}\\!P^{(0)})_{ij}(z)\biggr{\\}}$
$\displaystyle\quad+\Bigl{[}L_{\omega}^{2}\,\frac{(\tilde{\gamma}_{\nu\,0}^{i})^{2}}{8}-L_{\omega}\,\frac{\tilde{\gamma}_{\nu\,1}^{i}}{2}\Bigr{]}\mathbf{1}_{ij}(z)-L_{\omega}\,\frac{\tilde{\gamma}_{\nu\,0}^{i}}{2}\,\tilde{I}_{ij}^{(1)}(z)+\tilde{I}_{ij}^{(2)}(z)\,,$
(3.26)

where we abbreviated

$L_{b}=\ln\frac{b_{T}^{2}\mu^{2}}{b_{0}^{2}}\,,\quad
b_{0}=2e^{-\gamma_{E}}\,,\qquad L_{\omega}=\ln\frac{\nu}{\omega}\,.$ (3.27)

Note that $L_{\omega}$ differs from the characteristic logarithm of the soft
function in the previous section. The $\tilde{I}_{ij}^{(n)}(z)$ are given in
ref. [136] for quark and gluon beam functions in terms of the results of ref.
[134], and were directly calculated at NNLO using the exponential regulator
for the quark case in ref. [138].

At three loops we write

$\tilde{\mathcal{I}}_{ij}^{(3)}(z,b_{T},\mu,\nu/\omega)=\sum_{\ell=0}^{3}\tilde{\mathcal{I}}_{ij,\ell}^{(3)}(z,L_{\omega})\,L_{b}^{\ell}\,,$
(3.28)

and using $\tilde{\gamma}_{\nu\,0}^{i}=0$ for brevity, the coefficients are

$\displaystyle\tilde{\mathcal{I}}_{ij,3}^{(3)}(z,L_{\omega})$
$\displaystyle=\biggl{\\{}L_{\omega}^{3}\,\frac{(\Gamma^{i}_{0})^{3}}{6}+L_{\omega}^{2}\,\frac{(\Gamma^{i}_{0})^{2}}{2}\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{)}+L_{\omega}\Gamma^{i}_{0}\Bigl{[}\frac{\beta_{0}^{2}}{3}+\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{4}\Bigr{)}\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{]}$
$\displaystyle\qquad+\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{)}\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{4}\Bigr{)}\frac{\tilde{\gamma}_{B\,0}^{i}}{6}\biggr{\\}}\mathbf{1}_{ij}(z)$
$\displaystyle\quad-\biggl{[}L_{\omega}^{2}\,\frac{(\Gamma^{i}_{0})^{2}}{2}+L_{\omega}\Gamma^{i}_{0}\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{)}+\frac{\beta_{0}^{2}}{3}+\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{4}\Bigr{)}\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\biggr{]}P^{(0)}_{ij}(z)$
$\displaystyle\quad+\Bigl{(}L_{\omega}\,\frac{\Gamma^{i}_{0}}{2}+\frac{\beta_{0}}{2}+\frac{\tilde{\gamma}_{B\,0}^{i}}{4}\Bigr{)}(P^{(0)}\\!P^{(0)})_{ij}(z)-\frac{1}{6}(P^{(0)}\\!P^{(0)}\\!P^{(0)})_{ij}(z)\,,$
$\displaystyle\tilde{\mathcal{I}}_{ij,2}^{(3)}(z,L_{\omega})$
$\displaystyle=\biggl{\\{}L_{\omega}^{2}\,\Gamma^{i}_{0}\Gamma^{i}_{1}+L_{\omega}\Bigl{[}\Gamma_{1}^{i}\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{)}+\frac{\Gamma_{0}^{i}}{2}(\beta_{1}+\tilde{\gamma}_{B\,1}^{i})\Bigr{]}$
$\displaystyle\qquad+\beta_{0}\frac{\tilde{\gamma}_{B\,1}^{i}}{2}+\frac{\tilde{\gamma}_{B\,0}^{i}}{4}(\beta_{1}+\tilde{\gamma}_{B\,1}^{i})\biggr{\\}}\mathbf{1}_{ij}(z)-\Bigl{(}L_{\omega}\Gamma^{i}_{1}+\frac{\beta_{1}}{2}+\frac{\tilde{\gamma}_{B\,1}^{i}}{2}\Bigr{)}P_{ij}^{(0)}(z)$
$\displaystyle\quad-\Bigl{(}L_{\omega}\Gamma^{i}_{0}+\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{)}P_{ij}^{(1)}(z)+\frac{1}{2}(P^{(0)}\\!P^{(1)}+P^{(1)}\\!P^{(0)})_{ij}(z)$
$\displaystyle\quad+\Bigl{[}L_{\omega}^{2}\,\frac{(\Gamma^{i}_{0})^{2}}{2}+L_{\omega}\,\frac{\Gamma^{i}_{0}}{2}(3\beta_{0}+\tilde{\gamma}_{B\,0}^{i})+\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{)}\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{4}\Bigr{)}\Bigr{]}\tilde{I}_{ij}^{(1)}(z)$
$\displaystyle\quad-\Bigl{(}L_{\omega}\Gamma^{i}_{0}+\frac{3}{2}\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{)}(\tilde{I}^{(1)}\\!P^{(0)})_{ij}(z)+\frac{1}{2}(\tilde{I}^{(1)}\\!P^{(0)}\\!P^{(0)})_{ij}(z)\,,$
$\displaystyle\tilde{\mathcal{I}}_{ij,1}^{(3)}(z,L_{\omega})$
$\displaystyle=\biggl{\\{}-L_{\omega}^{2}\,\Gamma^{i}_{0}\frac{\tilde{\gamma}_{\nu\,1}^{i}}{2}+L_{\omega}\Bigl{[}-\Bigl{(}\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{4}\Bigr{)}\tilde{\gamma}_{\nu\,1}^{i}+\Gamma^{i}_{2}\Bigr{]}+\frac{\tilde{\gamma}_{B\,2}^{i}}{2}\biggr{\\}}\mathbf{1}_{ij}(z)$
$\displaystyle\quad+L_{\omega}\frac{\tilde{\gamma}_{\nu\,1}^{i}}{2}P_{ij}^{(0)}(z)-P_{ij}^{(2)}(z)+\Bigl{(}L_{\omega}\Gamma^{i}_{1}+\beta_{1}+\frac{\tilde{\gamma}_{B\,1}^{i}}{2}\Bigr{)}\tilde{I}_{ij}^{(1)}(z)-(\tilde{I}^{(1)}\\!P^{(1)})_{ij}(z)$
$\displaystyle\quad+\Bigl{(}L_{\omega}\Gamma^{i}_{0}+2\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{i}}{2}\Bigr{)}\tilde{I}_{ij}^{(2)}(z)-(\tilde{I}^{(2)}\\!P^{(0)})_{ij}(z)\,,$
$\displaystyle\tilde{\mathcal{I}}_{ij,0}^{(3)}(z,L_{\omega})$
$\displaystyle=-L_{\omega}\,\frac{\tilde{\gamma}_{\nu\,2}^{i}}{2}\,\mathbf{1}_{ij}(z)-L_{\omega}\,\frac{\tilde{\gamma}_{\nu\,1}^{i}}{2}\,\tilde{I}_{ij}^{(1)}(z)+\tilde{I}_{ij}^{(3)}(z)\,,$
(3.29)

where the three-loop boundary coefficients $\tilde{I}_{ij}^{(3)}(z)$ are
currently unknown. We have evaluated all Mellin convolutions appearing in eqs.
(3.4) and (3.4) with the help of the MT package [108]. In contrast to ref.
[18], we were able to perform all required convolutions in terms of standard
harmonic polylogarithms without encountering multiple polylogarithms, after
using the identity in eq. (2.23) to simplify some of the inputs.

The polarization-dependent kernels $\tilde{\mathcal{J}}_{gj}$ have a simpler
structure than the $\tilde{\mathcal{I}}_{ij}$ because their LO contribution
vanishes. For unpolarized gluon-fusion processes, the accompanying tensor
structures are only contracted with each other, and hence we only require
their NNLO expressions for the N3LO cross section. They are given by

$\displaystyle\tilde{\mathcal{J}}_{gj}^{(0)}(z,b_{T},\mu,\nu/\omega)$
$\displaystyle=0\,,$
$\displaystyle\tilde{\mathcal{J}}_{gj}^{(1)}(z,b_{T},\mu,\nu/\omega)$
$\displaystyle=\tilde{J}_{gj}^{(1)}(z)=4C_{j}\,\frac{1-z}{z}\,,$
$\displaystyle\tilde{\mathcal{J}}_{gj}^{(2)}(z,b_{T},\mu,\nu/\omega)$
$\displaystyle=L_{b}\Bigl{[}\Bigl{(}L_{\omega}\Gamma^{g}_{0}+\beta_{0}+\frac{\tilde{\gamma}_{B\,0}^{g}}{2}\Bigr{)}\tilde{J}_{gj}^{(1)}(z)-(\tilde{J}^{(1)}\\!P^{(0)})_{gj}(z)\Bigr{]}$
$\displaystyle\quad-
L_{\omega}\,\frac{\tilde{\gamma}_{\nu\,0}^{g}}{2}\,\tilde{J}_{gj}^{(1)}(z)+\tilde{J}_{gj}^{(2)}(z)\,.$
(3.30)

The two-loop terms $\tilde{J}_{gj}^{(2)}$ have recently been calculated in
ref. [142] using the exponential regulator and in ref. [143] using the
$\delta$ regulator. They can be converted to our convention via the relation

$\displaystyle\tilde{J}_{gj}^{(2)}(z)$ $\displaystyle=I^{\prime(2)}_{gi}(z)$
$\displaystyle=-\delta^{L}C_{g\leftarrow
j}^{(2;0,0)}(z)-\frac{1}{2}\tilde{s}_{g}^{{(1)}}\tilde{J}_{gj}^{(1)}(z)\,.$
(3.31)

In the first line of eq. (3.4), $I^{\prime(2)}_{gi}(z)$ is the two-loop
boundary term as given in ref. [142]. In the second line of eq. (3.4),
$\tilde{s}_{g}^{{(1)}}=-2C_{A}\zeta_{2}$ is the soft function constant at one
loop and $\delta^{L}C_{g\leftarrow j}^{(2;,0,0)}(z)$ is the two-loop finite
piece of the TMDPDF given in ref. [143].

#### Numerical impact

As for the soft function above, we illustrate the numerical impact of the
three-loop corrections for the resummed beam function

$\displaystyle\tilde{B}_{i}(x,\vec{b}_{T},\mu,\nu)$
$\displaystyle=\tilde{B}_{i}(x,\vec{b}_{T},\mu_{B},\nu_{B})\,\tilde{U}_{B}^{i}(\omega,b_{T},\mu_{B},\mu,\nu_{B},\nu)\,,$
$\displaystyle\tilde{U}_{B}^{i}(\omega,b_{T},\mu_{B},\mu,\nu_{B},\nu)$
$\displaystyle=\exp\biggl{[}-\frac{1}{2}\ln\frac{\nu}{\nu_{B}}\tilde{\gamma}_{\nu}^{i}(b_{T},\mu_{B})\biggr{]}\exp\biggl{[}\int_{\mu_{B}}^{\mu}\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\tilde{\gamma}_{B}^{i}(\mu^{\prime},\nu/\omega)\biggr{]}\,.$
(3.32)

For $i=g$, we restrict to the polarization-independent piece
$\tilde{\mathcal{I}}_{gj}$ and write $\tilde{B}_{g}\equiv-
g_{\perp,\rho\lambda}\tilde{B}_{g}^{\rho\lambda}$ for short. As for the soft
function, we restrict to simultaneous variations of $\mu_{B}$ and $\nu_{B}$.

Figure -1134: Residual scale dependence of the resummed $q_{T}$ beam function
in Fourier-space for $i=d$ (top left), $u$ (top right), $\bar{d}$ (bottom
left), and $g$ (bottom right). Shown are the relative deviations from the
NNLL′ result $\tilde{B}_{i}^{\mathrm{central}}$ at the central scales
$(\mu_{B},\nu_{B})=(\mu,\nu)=(b_{0}/b_{T},\omega)$.

In figure -1134, we show the residual dependence on $(\mu_{B},\nu_{B})$ at
NLL′ (dotted green), NNLL′ (dashed blue), and N3LL′ with the unknown
$\mathcal{I}^{(3)}_{ij}(z)=0$ (solid orange) as the relative difference to the
central NNLL′ result at $(\mu_{B},\nu_{B})=(\mu,\nu)=(b_{0}/b_{T},\omega)$ for
$b_{0}/b_{T}=20\,\mathrm{GeV}$ and $\omega=100\,\mathrm{GeV}$. As for
$\mathcal{T}_{0}$, we use four-loop running of $\alpha_{s}$ and
MMHT2014nnlo68cl [109] NNLO PDFs throughout. In all cases, the scale
dependence is substantially reduced at each order. We again anticipate that
this qualitative behavior continues to hold when the full result for
$\tilde{I}^{(3)}_{ij}(z)$ is included.

###  3.5 Beam function coefficients in the eikonal limit

We now proceed to extract the three-loop beam function coefficients in the
$z\to 1$ limit from consistency relations with known soft matrix elements. For
the $q_{T}$ beam function, these consistency relations arise from
factorization theorems for the triple-differential cross section
$\mathrm{d}\sigma_{pp\to L}/\mathrm{d}Q^{2}\mathrm{d}Y\mathrm{d}q_{T}$ that
enable the joint $q_{T}$ and soft threshold resummation [144, 145, 146, 147].
In terms of the momentum fractions $x_{a,b}$ defined in eq. (1.2), the soft
threshold limit is equivalent to taking both $x_{a}\to 1$ and $x_{b}\to 1$. As
$x_{a,b}\to 1$, initial state radiation is constrained to have energy
$\lesssim\lambda_{-}\lambda_{+}Q$, where

$\lambda_{-}^{2}\sim 1-x_{a}\qquad\text{and}\qquad\lambda_{+}^{2}\sim 1-x_{b}$
(3.33)

are power-counting parameters that encode the distance from the kinematic
endpoint.

The all-order factorization relevant for different hierarchies in $q_{T}/Q$
and the threshold constraint $\lambda_{-}\lambda_{+}$ was derived in refs.
[117, 148]. Some consequences of the resulting consistency relations have
already been explored in refs. [117, 148]. In fact, the exponential regulator
is _defined_ by its action on the refactorized pieces in these consistency
relations. In the following, we briefly review the relevant factorization
theorems and derive the all-order structure that arises for the $q_{T}$ beam
function in the eikonal limit.

#### $q_{T}/Q\ll\lambda_{-}\lambda_{+}\sim 1$

In this regime, initial-state radiation is not yet subject to the threshold
constraint, and the standard $q_{T}$ factorization theorem eq. (3.1) holds. It
receives power corrections $\mathcal{O}(q_{T}^{2}/Q^{2})$, but captures the
exact dependence on $x_{a,b}$ via the beam functions.

#### $q_{T}/Q\ll\lambda_{-}\lambda_{+}\ll 1$

For this hierarchy, the factorization takes a form similar to eq. (3.1), but
real collinear radiation into the final state is constrained in energy by
$1-x_{a,b}\ll 1$. The leftover radiation in this limit is described by
intermediate collinear-soft modes [149, 74] in terms of $n_{a,b}$-collinear-
soft functions $\tilde{\mathcal{S}}_{i}(k,b_{T},\mu,\nu)$. They are matrix
elements of collinear-soft Wilson lines and depend on the small additional
momentum $k=k^{\mp}$ available from either one of the (threshold) PDFs and on
the color charge of the colliding partons. The factorization theorem in this
regime reads [117, 148]

$\displaystyle\frac{\mathrm{d}\tilde{\sigma}(\vec{b}_{T})}{\mathrm{d}Q^{2}\,\mathrm{d}Y}$
$\displaystyle=\sum_{a,b}H_{ab}(Q^{2},\mu)\int\\!\mathrm{d}k^{-}\,\tilde{\mathcal{S}}_{i}(k^{-},b_{T},\mu,\nu)\,f^{\mathrm{thr}}_{a}\Bigl{[}x_{a}\Bigl{(}1+\frac{k^{-}}{\omega_{a}}\Bigr{)},\mu\Bigr{]}$
$\displaystyle\quad\times\int\\!\mathrm{d}k^{+}\,\tilde{\mathcal{S}}_{i}(k^{+},b_{T},\mu,\nu)\,f^{\mathrm{thr}}_{b}\Bigl{[}x_{b}\Bigl{(}1+\frac{k^{+}}{\omega_{b}}\Bigr{)},\mu\Bigr{]}\,\tilde{S}_{i}(b_{T},\mu,\nu)$
$\displaystyle\quad\times\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{1}{b_{T}^{2}\lambda_{-}^{2}\lambda_{+}^{2}Q^{2}},\lambda_{-}^{2},\lambda_{+}^{2}\Bigr{)}\Bigr{]}\,.$
(3.34)

Collinear-soft emissions do not contribute angular momentum, so the
polarization indices for gluon-induced processes become trivial in this limit
and we suppress them in the following.

#### $q_{T}/Q\sim\lambda_{-}\lambda_{+}\ll 1$

In this regime, the threshold constraint dominates and all radiation is forced
to be soft. The recoil against soft radiation with transverse momentum
$\vec{k}_{T}=-\vec{q}_{T}$ is encoded in the fully-differential threshold soft
function $S^{\mathrm{thr}}_{i}(k^{-},k^{+},\vec{k}_{T})$. In terms of its
Fourier transform with respect to $\vec{k}_{T}$,
$\tilde{S}^{\mathrm{thr}}_{i}(k^{-},k^{+},b_{T})$, the factorization theorem
reads

$\displaystyle\frac{\mathrm{d}\tilde{\sigma}(\vec{b}_{T})}{\mathrm{d}Q^{2}\,\mathrm{d}Y}$
$\displaystyle=\sum_{a,b}H_{ab}(Q^{2},\mu)\int\\!\mathrm{d}k^{-}\mathrm{d}k^{+}f^{\mathrm{thr}}_{a}\Bigl{[}x_{a}\Bigl{(}1+\frac{k^{-}}{\omega_{a}}\Bigr{)},\mu\Bigr{]}f^{\mathrm{thr}}_{b}\Bigl{[}x_{b}\Bigl{(}1+\frac{k^{+}}{\omega_{b}}\Bigr{)},\mu\Bigr{]}$
$\displaystyle\quad\times\tilde{S}^{\mathrm{thr}}_{i}(k^{-},k^{+},b_{T},\mu)\,\bigl{[}1+\mathcal{O}(\lambda_{-}^{2},\lambda_{+}^{2})\bigr{]}\,.$
(3.35)

Notably, the fully-differential threshold soft function is free of rapidity
divergences because they are regulated by the threshold constraint. (This is
the starting point of the exponential regularization procedure.) The fully-
differential soft function was calculated to $\mathcal{O}(\alpha_{s}^{2})$ in
ref. [150], albeit in a different context, and to
$\mathcal{O}(\alpha_{s}^{3})$ in ref. [82]. By construction, it satisfies

$\int\\!\mathrm{d}^{2}\vec{k}_{T}\,S^{\mathrm{thr}}_{i}(k^{-},k^{+},\vec{k}_{T},\mu)=\tilde{S}^{\mathrm{thr}}_{i}(k^{-},k^{+},b_{T}=0,\mu)=S^{\mathrm{thr}}_{i}(k^{-},k^{+},\mu)\,,$
(3.36)

where $S^{\mathrm{thr}}_{i}(k^{-},k^{+},\mu)$ is the double-differential
threshold soft function appearing in eq. (2.4).

#### Consistency relations

Consistency between eqs. (3.1) and (3.5) implies that the $x\to 1$ limit of
the $q_{T}$ beam function is captured by the collinear-soft function [117,
148],

$\tilde{B}_{i}(x,\vec{b}_{T},\mu,\nu)=\int\\!\mathrm{d}k\,\tilde{\mathcal{S}}_{i}(k,b_{T},\mu,\nu)\,f_{i}^{\mathrm{thr}}\Bigl{[}x\Bigl{(}1+\frac{k}{\omega}\Bigr{)},\mu\Bigr{]}\,\bigl{[}1+\mathcal{O}(1-x)\bigr{]}\,.$
(3.37)

This is the analog of eq. (2.29) for $q_{T}$, but this time relates the
eikonal limit of the beam function to an exclusive collinear-soft matrix
element instead of the inclusive threshold soft function. At the partonic
level, eq. (3.37) implies [117, 148]

$\tilde{\mathcal{I}}_{ij}(z,b_{T},\mu,\nu/\omega)=\delta_{ij}\,\omega\,\tilde{\mathcal{S}}_{i}\bigl{[}\omega(1-z),b_{T},\mu,\nu\bigr{]}\,\bigl{[}1+\mathcal{O}(1-z)\bigr{]}\,.$
(3.38)

Note that eq. (3.38) is true for any rapidity regulator as long as the same
regulator is used on both sides. The consistency between eqs. (3.5) and (3.5)
implies [117, 148]

$\tilde{S}_{i}^{\mathrm{thr}}(k^{-},k^{+},b_{T},\mu)=\tilde{\mathcal{S}}_{i}(k^{-},b_{T},\mu,\nu)\,\tilde{\mathcal{S}}_{i}(k^{+},b_{T},\mu,\nu)\,\tilde{S}_{i}(b_{T},\mu,\nu)\,\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{1}{b_{T}^{2}k^{-}k^{+}}\Bigr{)}\Bigr{]}\,,$
(3.39)

which again holds for any choice of rapidity regulator. In particular, the
left-hand side has no rapidity divergences, so the dependence on the rapidity
regulator cancels between the terms on the right-hand side. Together, eqs.
(3.37) and (3.39) uniquely determine the eikonal limit of the beam function in
any given rapidity regulator scheme in terms of the fully-differential soft
function (which is independent of the scheme) and the $q_{T}$ soft function
$\tilde{S}_{i}(b_{T},\mu,\nu)$ (which determines the scheme). Furthermore, the
scheme ambiguity amounts to moving terms from the soft function boundary
coefficients into the coefficient of $\delta(1-z)$ in the beam function
coefficients. Since $\delta(1-z)$ is a leading-power contribution as $z\to 1$,
it follows that up to lower-order cross terms, all scheme-dependent terms in
the beam function are contained in the leading eikonal terms predicted by eq.
(3.38).

#### Extraction of the finite terms

For the exponential regulator, the relation between the fully-differential and
standard TMD soft function is particularly simple, leading to an all-order
result for the collinear-soft function in terms of the rapidity anomalous
dimension, see appendix C. Inserting this result into eq. (3.38), we find for
the eikonal limit of the $b_{T}$-space beam function matching coefficient
$\tilde{\mathcal{I}}_{ij}$ in the exponential regulator scheme,

$\displaystyle\tilde{\mathcal{I}}_{ij}(z,b_{T},\mu,\nu/\omega)=\delta_{ij}\,\frac{\omega}{\nu}\mathcal{V}_{\tilde{\gamma}_{\nu}^{i}(b_{T},\mu)/2}\Bigl{[}\frac{\omega}{\nu}(1-z)\Bigr{]}\,\bigl{[}1+\mathcal{O}(1-z)\bigr{]}\,,$
(3.40)

where the plus distribution $\mathcal{V}_{a}(x)$ is defined in eq. (A.4). The
simplicity of this result is a direct consequence of the specific rapidity
regulator, i.e., one may equally well have imposed this form of the eikonal
limit as the renormalization condition. Nonetheless, when combined with the
soft function to a given order, the scheme dependence cancels and leaves
behind a unique set of terms that capture the threshold limit of the singular
cross section in eq. (3.1). We note that a close relation between the rapidity
anomalous dimension and the eikonal limit of the beam function is a scheme-
independent feature [148], and was also conjectured for the $\delta$-regulator
in ref. [137].

It is straightforward to expand eq. (3.40) in $\alpha_{s}$ to obtain the
finite terms in the matching coefficient at any given fixed order using eqs.
(3.2) and (A.1). Up to two loops we have

$\displaystyle\tilde{I}_{ij}^{(1)}(z)$
$\displaystyle=\mathcal{O}\bigl{[}(1-z)^{0}\bigr{]}\,,$
$\displaystyle\tilde{I}_{ij}^{(2)}(z)$
$\displaystyle=\delta_{ij}\,\frac{\tilde{\gamma}_{\nu\,1}^{i}}{2}\,\mathcal{L}_{0}(1-z)+\mathcal{O}\bigl{[}(1-z)^{0}\bigr{]}\,,$
(3.41)

in agreement with the full two-loop result [136], and where we have used that
$\tilde{\gamma}_{\nu\,0}^{i}=0$. Including terms up to six loops for
illustration, we find

$\displaystyle\tilde{I}_{ij}^{(3)}(z)$
$\displaystyle=\delta_{ij}\,\frac{\tilde{\gamma}_{\nu\,2}^{i}}{2}\,\mathcal{L}_{0}(1-z)+\mathcal{O}\bigl{[}(1-z)^{0}\bigr{]}\,,$
$\displaystyle\tilde{I}_{ij}^{(4)}(z)$
$\displaystyle=\delta_{ij}\,\frac{\tilde{\gamma}_{\nu\,3}^{i}}{2}\,\mathcal{L}_{0}(1-z)+\frac{(\tilde{\gamma}_{\nu\,1}^{i})^{2}}{4}\Bigl{[}\mathcal{L}_{1}(1-z)-\frac{\zeta_{2}}{2}\delta(1-z)\Bigr{]}+\mathcal{O}\bigl{[}(1-z)^{0}\bigr{]}\,,$
$\displaystyle\tilde{I}_{ij}^{(5)}(z)$
$\displaystyle=\delta_{ij}\,\frac{\tilde{\gamma}_{\nu\,4}^{i}}{2}\,\mathcal{L}_{0}(1-z)+\frac{\tilde{\gamma}_{\nu\,1}^{i}\tilde{\gamma}_{\nu\,2}^{i}}{2}\Bigl{[}\mathcal{L}_{1}(1-z)-\frac{\zeta_{2}}{2}\delta(1-z)\Bigr{]}+\mathcal{O}\bigl{[}(1-z)^{0}\bigr{]}\,,$
$\displaystyle\tilde{I}_{ij}^{(6)}(z)$
$\displaystyle=\delta_{ij}\,\frac{\tilde{\gamma}_{\nu\,5}^{i}}{2}\,\mathcal{L}_{0}(1-z)+\frac{(\tilde{\gamma}_{\nu\,2}^{i})^{2}+2\tilde{\gamma}_{\nu\,1}^{i}\tilde{\gamma}_{\nu\,3}^{i}}{4}\Bigl{[}\mathcal{L}_{1}(1-z)-\frac{\zeta_{2}}{2}\delta(1-z)\Bigr{]}$
$\displaystyle\quad+\frac{(\tilde{\gamma}_{\nu\,1}^{i})^{3}}{8}\Bigl{[}\frac{\mathcal{L}_{2}(1-z)}{2}-\frac{\zeta_{2}}{2}\mathcal{L}_{0}(1-z)+\frac{\zeta_{3}}{3}\delta(1-z)\Bigr{]}+\mathcal{O}\bigl{[}(1-z)^{0}\bigr{]}\,.$
(3.42)

We again stress that these expressions are a direct consequence of the
renormalization condition in the exponential regulator scheme and must be
combined with the soft function in the same scheme to obtain a scheme-
independent result. It is interesting to note that starting at four loops, eq.
(3.40) does in fact predict a term proportional to $\delta(1-z)$ in the beam
function matching coefficient due to the inverse Fourier transform to
$k^{\pm}$ back from the conjugate $b^{\pm}$ space, where the regularization
procedure is applied.

###  3.6 Estimating beam function coefficients beyond the eikonal limit

As in section 2.5, we can use the eikonal limit of the beam function
coefficients to study to what extent it can be used to approximate the full
result and/or estimate the uncertainty due to the missing terms beyond the
eikonal limit.

In figure -1133, we compare the full $q_{T}$ beam function coefficient (solid)
to its eikonal (LP dotted green) and next-to-eikonal (NLP dashed blue)
expansions at NNLO for the $u$-quark and gluon channels. Since the NLO
coefficients are not singular, we do not show the corresponding NLO results.
We always show the convolution $(I_{ij}\otimes f_{j})(x)/f_{i}(x)$ with the
appropriate PDF $f_{j}$ and normalize to the PDF $f_{i}(x)$, corresponding to
the LO result, where $i=u$ for the $u$-quark case and $i=g$ for the gluon
case. With this normalization, the shape gives an indication of the rapidity
dependence of the beam function coefficient relative to the LO rapidity
dependence induced by the shape of the PDFs. We also include the appropriate
powers of $\alpha_{s}/(4\pi)$ at each order, so the overall normalization
shows the percent impact relative to the LO result. For definiteness, the
renormalization scale entering the PDFs is chosen as
$\mu=30~{}\,\mathrm{GeV}$.

In both flavor-diagonal contributions, denoted as $qqV$ and $gg$, the eikonal
limit correctly reproduces the divergent behavior as $x\to 1$, but is off away
from very large $x$. Including the next-to-eikonal terms yields a sizable
shift from the eikonal limit, and provides a very good approximation in the
shown $x$ region. In the gluon channel, one can see a rise of the full kernel
towards small $x$, arising from an overall $1/z$ divergence in the coefficient
$I_{gg}^{(2)}(z)$, which is not captured by the expansion around $z\to 1$. If
desired, one could also include the leading $z\to 0$ behavior of the
coefficients, which for simplicity is not done here. For illustration, we also
show the total contribution from all other corresponding nondiagonal channels
(gray dot-dashed). In both cases, they are numerically subdominant to the
flavor-diagonal channel and also much flatter in $x$, since they only start at
NLP.

Figure -1133: Comparison of the full beam function coefficients to their
leading eikonal (LP) and next-to-eikonal (NLP) expansion at NNLO. The
$u$-quark channel is shown on the left and the gluon channel on the right. In
both cases we also show the sum of all nondiagonal partonic channels for
comparison.

Similar to the $\mathcal{T}_{0}$ coefficients in section 2.5, we now wish to
make an ansatz for the unknown three-loop NLP terms to get an estimate of
their size. A peculiar feature of the $q_{T}$ coefficients is that up to three
loops, its eikonal limit contains no logarithmic distributions
$\mathcal{L}_{n}(1-z)$ with $n>0$, but only $\mathcal{L}_{0}(1-z)$. In
contrast, the NLP NNLO coefficient does contain a double logarithm
$\ln^{2}(1-z)$. Based on this observation, we make the following ansatz for
the NnLO beam coefficient,

$\displaystyle\tilde{I}_{ij,\mathrm{approx}}^{(n)}(z)=\tilde{I}_{ij}^{(n)\,\mathrm{LP}}(z)$
$\displaystyle+\Bigl{[}X_{1}\Gamma_{0}^{i}\ln^{2}(1-z)+X_{2}\gamma_{X}^{i}\ln(1-z)\Bigr{]}\tilde{I}_{ij,\mathrm{reg}}^{(n-1)}(z)$
$\displaystyle-X_{3}(1-z)\tilde{I}_{ij}^{(n)\,\mathrm{LP}}(z)\,.$ (3.43)

Here, $\tilde{I}_{ij,\mathrm{reg}}^{(n)}$ refers to the full regular (non-
eikonal) piece of the beam coefficient at $\mathcal{O}(\alpha_{s}^{n})$. At
NLO, there is no NLP term, so at this order we simply define the regular piece
to be the appropriate color factor. More explicitly, we use

$\displaystyle\tilde{I}_{ij,\mathrm{reg}}^{(1)}(z)$
$\displaystyle=-\delta_{ij}C_{i}\,,\qquad\tilde{I}_{ij,\mathrm{reg}}^{(2)}(z)=\tilde{I}_{ij}^{(2)}(z)-\delta_{ij}\frac{\tilde{\gamma}_{\nu\,1}^{i}}{2}\mathcal{L}_{0}(1-z)\,.$
(3.44)

The ansatz in eq. (3.6) dresses the lower-order regular kernel with two
additional logarithms $\ln(1-z)$. The coefficients of these logarithms are
chosen such that at the central choices $X_{1}=X_{2}=1$, they reproduce the
known double and single logarithms at NNLO. The effective noncusp anomalous
dimension $\gamma_{X}^{i}$ needed to achieve this is given by

$\displaystyle\gamma_{X}^{g}=3C_{A}-\beta_{0}\,,\qquad\gamma_{X}^{q}=10(C_{F}-C_{A})\,.$
(3.45)

The size of these additional logarithms can be probed by varying the
coefficients $X_{1,2}$ by $\pm 1$ around the central choice. Furthermore, we
add the eikonal limit $\tilde{I}_{ij}^{(n)\,\mathrm{LP}}$ suppressed by one
power of $(1-z)$ to estimate the pure NLP constant. Its coefficient $X_{3}$ is
varied by $\pm 1$ around the central choice $X_{3}=0$.

Since the $X_{i}$ probe independent structures, we can consider them as
uncorrelated. Hence, we add the impacts $\Delta_{i}$ on the final result of
their variation in quadrature

$\displaystyle\Delta=\Delta_{1}\oplus\Delta_{2}\oplus\Delta_{3}=\sqrt{\Delta_{1}^{2}+\Delta_{2}^{2}+\Delta_{3}^{2}}\,.$
(3.46)

Figure -1132: Approximate ansatzes for the NNLO (top) and N3LO (bottom)
kernels, in the $u$-quark (left) and gluon (right) channels.

In figure -1132, we show the approximate kernel at NNLO (top) and N3LO
(bottom) for the $u$-quark (left) and gluon (right) channels. The dashed
orange line shows the central result from our ansatz and the yellow band its
estimated uncertainty. The gray lines show the impact of the individual
variations of the $X_{i}$ as indicated. In the top panel (NNLO), we also show
the known full two-loop results (red dashed). It shows that the ansatz
including uncertainties approximates the true result relatively well, even for
the gluon case in the shown $x$ region. In particular, the rather large shift
from LP to the approximate NLP result is needed to correctly capture the full
result within uncertainties.

At N3LO, we see again that the approximate result gives rise to a sizable
shift from the pure eikonal limit, which by itself is a very small correction.
This large shift arises on the one hand because the LP limit only contains
$\mathcal{L}_{0}(1-z)$ with a rather small coefficient $\gamma_{\nu\,2}^{i}$,
while the NLP now contains up to $\ln^{4}(1-z)$. The fact that the uncertainty
bands are of similar size at NNLO and N3LO reflects their numerical importance
and that relatively little is known about the NLP structure, which also
motivates an exact calculation of the three-loop coefficients.

Finally, we briefly comment on the treatment of the unknown three-loop beam
function coefficients in ref. [18], where the $q_{T}$ subtraction was first
applied at N3LO for Higgs production. There, the employed approximation was
$\tilde{I}_{gg}^{(3)}(z)=\tilde{C}_{N3}\,\delta(1-z)$, with $\tilde{C}_{N3}$
fixed such that the inclusive cross section is correctly reproduced. This
effectively absorbs the averaged effect of the actual $z$ dependence into an
effective $\delta(1-z)$ coefficient. From our results we know the exact
$\delta(1-z)$ coefficient, and so our approximate results give an independent
estimate of the actual rapidity dependence and total size of these unknown
terms.

##  4 N3LO subtractions

The factorization in eq. (1.3) fully describes the limit $\tau\to 0$ and thus
captures the singular structure of QCD in this limit. Hence, it can be used to
construct a subtraction method for fixed-order calculations. In principle,
this works for any resolution variable $\tau$ and any process for which a
corresponding factorization is known [35, 151, 152, 153, 154, 155, 156, 38,
39]. The subtractions can be formulated differential in $\tau$ or as a global
$\tau$ slicing, which we briefly review in the following. For a more extensive
discussion we refer to ref. [39].

Our starting point is to write the inclusive cross section as the integral
over the differential cross section in $\tau$,

$\sigma(X)=\int\\!\mathrm{d}\tau\,\frac{\mathrm{d}\sigma(X)}{\mathrm{d}\tau}\,,\qquad\sigma(X,\tau_{\mathrm{cut}})=\int^{\tau_{\mathrm{cut}}}\\!\mathrm{d}\tau\,\frac{\mathrm{d}\sigma(X)}{\mathrm{d}\tau}\,,$
(4.1)

where the second relation defines the cumulant in $\tau_{\mathrm{cut}}$. Here,
$X$ denotes any measurements performed, which can include $Q$ and $Y$ of the
color singlet $L$ but also additional measurements or cuts on its
constituents. For $\tau\to 0$, the cross sections scales like $\sim 1/\tau$,
so performing the $\tau$ integral requires knowing the full analytic
distributional structure involving $\delta(\tau)$ and $\mathcal{L}_{n}(\tau)$,
which encodes the cancellation of real and virtual IR divergences. To separate
out the singular structure in $\tau$, we introduce a subtraction term,

$\displaystyle\sigma(X)$
$\displaystyle=\sigma^{\mathrm{sub}}(X,\tau_{\mathrm{off}})+\int\\!\mathrm{d}\tau\,\biggl{[}\frac{\mathrm{d}\sigma(X)}{\mathrm{d}\tau}-\frac{\mathrm{d}\sigma^{\mathrm{sub}}(X)}{\mathrm{d}\tau}\theta(\tau<\tau_{\mathrm{off}})\biggr{]}\,,$
(4.2)

where $\mathrm{d}\sigma^{\mathrm{sub}}(X)/\mathrm{d}\tau$ captures all
singularities for $\tau\to 0$,

$\frac{\mathrm{d}\sigma(X)}{\mathrm{d}\tau}=\frac{\mathrm{d}\sigma^{\mathrm{sing}}(X)}{\mathrm{d}\tau}\,[1+\mathcal{O}(\tau)]\,,\qquad\frac{\mathrm{d}\sigma^{\mathrm{sub}}(X)}{\mathrm{d}\tau}=\frac{\mathrm{d}\sigma^{\mathrm{sing}}(X)}{\mathrm{d}\tau}\,[1+\mathcal{O}(\tau)]\,,$
(4.3)

and $\sigma^{\mathrm{sub}}(X,\tau_{\mathrm{off}})$ is the integrated
subtraction term,

$\sigma^{\mathrm{sub}}(X,\tau_{\mathrm{off}})=\int\\!\mathrm{d}\tau\,\frac{\mathrm{d}\sigma^{\mathrm{sub}}(X)}{\mathrm{d}\tau}\theta(\tau<\tau_{\mathrm{off}})\,.$
(4.4)

By construction, the integrand in square brackets in eq. (4.2) contains at
most integrable singularities for $\tau\to 0$ and so the integral can be
performed numerically. Hence, the full cross section
$\mathrm{d}\sigma(X)/\mathrm{d}\tau$ is only ever evaluated at finite
$\tau>0$, which means it can be obtained from a calculation of the
corresponding $ab\to L+1$-parton process at one lower order. In practice, one
always has a small IR cutoff $\delta$ on the $\tau$ integral,

$\sigma(X)=\sigma^{\mathrm{sub}}(X,\tau_{\mathrm{off}})+\int_{\delta}\\!\mathrm{d}\tau\,\biggl{[}\frac{\mathrm{d}\sigma(X)}{\mathrm{d}\tau}-\frac{\mathrm{d}\sigma^{\mathrm{sub}}(X)}{\mathrm{d}\tau}\theta(\tau<\tau_{\mathrm{off}})\biggr{]}+\Delta\sigma(X,\delta)\,,$
(4.5)

where the last term contains the integral over $\tau\leq\delta$,

$\Delta\sigma(X,\delta)=\sigma(X,\delta)-\sigma^{\mathrm{sub}}(X,\delta)\sim\mathcal{O}(\delta)\,.$
(4.6)

which is neglected for $\delta\to 0$.

The above is a _differential_ $\tau$-subtraction scheme, where the parameter
$\tau_{\mathrm{off}}\sim 1$ determines the range over which the subtraction
acts. The key advantage of formulating the subtractions in terms of a physical
resolution variable $\tau$, is that the subtraction terms are given by the
singular limit of a physical cross section. Hence, they are precisely given by
the factorization formula for $\tau\to 0$, which is also the basis for the
resummation in $\tau$. In fact, this form of the subtraction is routinely used
when the resummed and fixed-order results are combined via an additive
matching. In this case, $\tau_{\mathrm{off}}$ corresponds to the point where
the $\tau$ resummation is turned off, and the term in square brackets in eq.
(4.2) is the nonsingular cross section that is added to the pure resummed
result. Differential $\mathcal{T}_{0}$ subtractions are used in this way in
the Geneva Monte Carlo to combine the fully-differential NNLO calculation
together with the NNLL′ $\mathcal{T}_{0}$ resummation with a parton shower
[46, 157, 47]. The differential subtractions at N3LO are a key ingredient for
using this method to combine N3LO calculations with parton showers.

In contrast to a fully local subtraction scheme, all singularities are
projected onto the resolution variable $\tau$, so the subtractions are local
in $\tau$ but nonlocal in the additional radiation phase space that is
integrated over. As discussed in ref. [39], the subtractions can be made more
local by considering a factorization theorem that is differential in more
variables. For example, the combined $q_{T}$ and $\mathcal{T}_{0}$ resummation
[74, 75] offers the possibility to construct double-differential
$q_{T}-\mathcal{T}_{0}$ subtractions.

The key point of the differential subtraction is that $\delta$ can in
principle be made arbitrarily small, because the integrand of the $\tau$
integral is nonsingular, which also means that the numerically expensive small
$\tau$ region does not need to be sampled with weight $1/\tau$. On the other
hand, by letting $\delta=\tau_{\mathrm{cut}}$ be a small but finite cutoff and
setting $\tau_{\mathrm{off}}=\tau_{\mathrm{cut}}$, eq. (4.5) turns into a
global $\tau$ subtraction or _slicing_ ,

$\sigma(X)=\sigma^{\mathrm{sub}}(X,\tau_{\mathrm{cut}})+\int_{\tau_{\mathrm{cut}}}\\!\mathrm{d}\tau\,\frac{\mathrm{d}\sigma(X)}{\mathrm{d}\tau}+\Delta\sigma(X,\tau_{\mathrm{cut}})\,.$
(4.7)

The practical advantage of the slicing method is that it allows one to readily
turn an existing $L+1$-jet Nn-1LO calculation into a NnLO calculation for $L$,
and so most implementations use this approach [158, 159, 160, 38, 161, 162,
163, 18, 164]. The main disadvantage is that the cancellation of the
divergences now only happens after the integration over $\tau$. This makes the
$L+1$-jet calculation very demanding, both in terms of computational expense
and numerical stability, because the $1/\tau$-divergent integral of
$\mathrm{d}\sigma(X)/\mathrm{d}\tau$ must be computed with sufficient accuracy
down to sufficiently small $\tau_{\mathrm{cut}}$, which in practice limits how
small one can take $\tau_{\mathrm{cut}}$. Since the integral is divergent, one
cannot let $\tau_{\mathrm{cut}}\to 0$ even in principle, so one always has a
leftover systematic uncertainty from the neglected power corrections
$\Delta\sigma(X,\tau_{\mathrm{cut}})$.

The numerical efficiency of the subtractions can be improved by including the
power corrections in the subtractions for both $\mathcal{T}_{0}$ [165, 166,
167, 168, 169, 170] and $q_{T}$ [171, 172]. The size of the missing power
corrections also strongly depends on the precise definition of
$\mathcal{T}_{0}$ [165, 166, 167]. The hadronic definition in eq. (2.1)
exhibits power corrections that grow like $e^{\lvert Y\rvert}$ at large $Y$,
which is not the case for the leptonic definition. The power corrections also
depend on the Born measurement $X$. In particular, additional selection or
isolation cuts on the color-singlet constituents typically enhance the power
corrections from $\mathcal{O}(\tau)$ to $\mathcal{O}(\sqrt{\tau})$ [173].

###  4.1 Subtraction terms

The singular terms needed for the subtractions only depend on the Born phase
space, so we can write them as

$\frac{\mathrm{d}\sigma^{\mathrm{sing}}(X)}{\mathrm{d}\tau}=\int\\!\mathrm{d}\Phi_{0}\,\frac{\mathrm{d}\sigma^{\mathrm{sing}}(\Phi_{0})}{\mathrm{d}\tau}\,X(\Phi_{0})\,,$
(4.8)

where $\Phi_{0}\equiv\Phi_{0}(\kappa_{a},\kappa_{b},\omega_{a},\omega_{b})$
denotes the full Born phase space, including the parton labels $\kappa_{a,b}$,
the total color-singlet momentum $q^{\mu}$ parametrized in terms of
$\omega_{a,b}$ as in eqs. (1.1) and (1.2) as well as the internal phase space
of $L$. The $X(\Phi_{0})$ denotes the measurement function that implements the
measurement $X$ on a Born configuration.

The singular terms are defined such that their $\tau$ dependence is minimal
and given by

$\displaystyle\frac{\mathrm{d}\sigma^{\mathrm{sing}}(\Phi_{0})}{\mathrm{d}\tau}$
$\displaystyle=\mathcal{C}_{-1}(\Phi_{0})\,\delta(\tau)+\sum_{n\geq
0}\mathcal{C}_{n}(\Phi_{0})\,\mathcal{L}_{n}(\tau)$ $\displaystyle=\sum_{m\geq
0}\biggl{[}\mathcal{C}_{-1}^{(m)}(\Phi_{0})\,\delta(\tau)+\sum_{n=0}^{2m-1}\mathcal{C}_{n}^{(m)}(\Phi_{0})\mathcal{L}_{n}(\tau)\biggr{]}\Bigl{(}\frac{\alpha_{s}}{4\pi}\Bigr{)}^{m}\,.$
(4.9)

Their integral over $\tau\leq\tau_{\mathrm{cut}}$ immediately follows as

$\displaystyle\sigma^{\mathrm{sing}}(\Phi_{0},\tau_{\mathrm{cut}})$
$\displaystyle=\mathcal{C}_{-1}(\Phi_{0})+\sum_{n\geq
0}\mathcal{C}_{n}(\Phi_{0})\,\frac{\ln^{n}\tau_{\mathrm{cut}}}{n+1}$
$\displaystyle=\sum_{m\geq
0}\biggl{[}\mathcal{C}_{-1}^{(m)}(\Phi_{0})+\sum_{n=0}^{2m-1}\mathcal{C}_{n}^{(m)}(\Phi_{0})\frac{\ln^{n}\tau_{\mathrm{cut}}}{n+1}\biggr{]}\Bigl{(}\frac{\alpha_{s}}{4\pi}\Bigr{)}^{m}\,.$
(4.10)

The differential subtractions are given by using eq. (4.1) for $\tau>0$, which
amounts to dropping the $\mathcal{C}_{-1}(\Phi_{0})\delta(\tau)$ term and
using $\mathcal{L}_{n}(\tau>0)=\ln^{n-1}(\tau)/\tau$. The integrated
subtractions are directly given by eq. (4.1).

The precise definition of the $\mathcal{C}_{n}(\Phi_{0})$ coefficients depends
on the normalization of the dimensionless variable $\tau$ or equivalently on
the boundary condition of the $\mathcal{L}_{n}(\tau)$. Rescaling
$\tau\to\lambda\tau$ moves contributions from $\mathcal{C}_{n}(\Phi_{0})$ to
$\mathcal{C}_{m<n}(\Phi_{0})$. This freedom was used in ref. [39] to absorb
all terms with $n\geq 0$ in eq. (4.1) into a
$\mathcal{C}_{-1}(\Phi_{0},\mathcal{T}_{\mathrm{off}})$ by taking
$\tau\equiv\mathcal{T}_{0}/\mathcal{T}_{\mathrm{off}}$. Here, we prefer to
keep the cutoff dependence explicit as in eq. (4.1) and take

$\displaystyle\tau$
$\displaystyle\equiv\frac{\mathcal{T}_{0}}{Q}\qquad\text{(for
$\mathcal{T}_{0}$)}\,,\qquad\tau\equiv\frac{q_{T}^{2}}{Q^{2}}\qquad\text{(for
$q_{T}$)}\,.$ (4.11)

The $m$-loop subtraction coefficients $\mathcal{C}_{n}^{(m)}(\Phi_{0})$
directly follow from expanding eq. (2.1) for $\mathcal{T}_{0}$ or eq. (3.1)
for $q_{T}$ to $m$th order in $\alpha_{s}$. For the three-loop coefficients
this yields

$\displaystyle\mathcal{C}_{-1}^{(3)}(\Phi_{0})$
$\displaystyle=H^{(3)}(\Phi_{0})\,f_{a}(x_{a})\,f_{b}(x_{b})+\sum_{m=1}^{3}H^{(3-m)}(\Phi_{0})\,\bigl{[}B_{a}(x_{a})B_{b}(x_{b})S\bigr{]}_{-1}^{(m)}\,,$
$\displaystyle\mathcal{C}_{n\geq 0}^{(3)}(\Phi_{0})$
$\displaystyle=\sum_{k=1}^{3}H^{(3-k)}(\Phi_{0})\,\bigl{[}B_{a}(x_{a})B_{b}(x_{b})S\bigr{]}_{n}^{(k)}\,,$
(4.12)

where for simplicity we have suppressed the dependence on $\mu$ and the
distinction of the $\mathcal{T}_{0}$ vs. $q_{T}$ beam and soft functions. The
virtual three-loop corrections to the Born process are contained in
$H^{(3)}(\Phi_{0})$, which only enters in $\mathcal{C}_{-1}^{(3)}$. The
$m$-loop soft/collinear contribution $[BBS]_{n}^{(m)}$ follows from inserting
the fixed-order expansions of the respective beam and soft function,
reexpanding their product to $m$th order and picking out the coefficients of
$\delta(\tau)$ and $\mathcal{L}_{n}(\tau)$. The three-loop boundary
coefficients of the beam and soft functions only enter in
$\mathcal{C}_{-1}^{(3)}$ and thus are needed for the integrated subtraction
terms but not the differential ones. Note also that most of the process and
$\Phi_{0}$ dependence resides in the hard coefficients, while the
soft/collinear contributions only depend on $x_{a,b}$ and the parton types,

$\bigl{[}B_{a}(x_{a})B_{b}(x_{b})S\bigr{]}_{n}^{(m)}=\int\\!\frac{\mathrm{d}z_{a}}{z_{a}}\frac{\mathrm{d}z_{b}}{z_{b}}\,\sum_{i,j}\bigl{[}\mathcal{I}_{ai}(z_{a})\mathcal{I}_{bi}(z_{b})S\bigr{]}_{n}^{(m)}f_{i}\Bigl{(}\frac{x_{a}}{z_{a}}\Bigr{)}f_{j}\Bigl{(}\frac{x_{b}}{z_{b}}\Bigr{)}\,.$
(4.13)

The results for the subtraction coefficients $\mathcal{C}_{n}(\Phi_{0})$ in
eq. (4.1) up to three loops for both $\mathcal{T}_{0}$ and $q_{T}$ have been
implemented in the C++ library SCETlib [174] and will be made publicly
available.

Note that evaluating eq. (4.1) for $\mathcal{T}_{0}$ requires rescaling and
convolving the plus distributions in the beam and soft functions, as discussed
in ref. [39]. For $q_{T}$, expanding the ${\vec{b}}_{T}$-space result
$\mathrm{d}\tilde{\sigma}^{\mathrm{sing}}({\vec{b}}_{T})$ yields powers of the
${\vec{b}}_{T}$-space logarithm $L_{b}^{n}$ up to $n\leq 6$. Their Fourier
transform, given in table 1 in appendix A.2, yields simple
$\delta({\vec{q}}_{T})$ and $\mathcal{L}_{n}({\vec{q}}_{T},\mu)$, which are
easily rescaled to $\delta(\tau)$ and $\mathcal{L}_{n}(\tau)$.

Note also that the original $q_{T}$ subtraction method in ref. [35] was based
on the $q_{T}$ resummation framework of ref. [175], where the canonical
${\vec{b}}_{T}$-space logarithms are replaced by

$L_{b}\to\tilde{L}_{b}\equiv\ln\Bigl{(}\frac{b_{T}^{2}\mu^{2}}{b_{0}^{2}}+1\Bigr{)}\,.$
(4.14)

This form is also used e.g. in refs. [160, 18]. While using $\tilde{L}_{b}$
has certain advantages in the context of $q_{T}$ resummation, it is
unnecessary for the purpose of $q_{T}$ subtractions, since $L_{b}$ and
$\tilde{L}_{b}$ yield the same singular terms and only differ by power
corrections. A drawback of using $\tilde{L}_{b}$ here is that the Fourier
transform of $\tilde{L}_{b}^{n}$ yields complicated expressions in $q_{T}$
space, see appendix B in ref. [175], whose cumulants are not known
analytically and must be performed numerically.

##  5 Conclusions

We have studied the three-loop structure of beam and soft functions for both
$0$-jettiness $\mathcal{T}_{0}$ and transverse momentum $q_{T}$. These
functions are defined as collinear proton matrix elements and soft vacuum
matrix element, measuring the small light-cone momentum (for
$\mathcal{T}_{0}$) or total transverse momentum (for $q_{T}$) of all soft and
collinear emissions, and thus are universal objects probing the infrared
structure of QCD.

The all-order structure of the beam and soft functions is governed by
renormalization group equations, which we have employed to derive their full
three-loop structure. For the currently unknown scale-independent boundary
coefficients $I^{(3)}_{ij}(z)$ of the N3LO beam functions, we employ
consistency between different factorization limits to derive their leading
eikonal limit $I^{(3)}_{ij}(z\to 1)$, i.e. the full singular limit of the beam
functions as $z\to 1$, and estimate the size of the unknown terms beyond the
eikonal limit. All results of this paper will be made available in the C++
library SCETlib [174].

Our results provide important ingredients required for the resummation of
$\mathcal{T}_{0}$ and $q_{T}$ at N3LL′ and N4LL order. In particular, they are
important for extending the $q_{T}$ and $\mathcal{T}_{0}$ subtraction methods
to N3LO, for which we provide the complete set of differential subtraction
terms at three loops, which are for example necessary for extending the
matching of fixed-order calculations to parton showers to N3LO$+$PS. The
integrated subtraction terms are not yet fully known at three loops, but the
obtained eikonal limit allows us to provide a first approximation for a full
three-loop subtraction, and will be a useful cross check once the full $q_{T}$
and $\mathcal{T}_{0}$ beam functions become available.

#### Note added:

As discussed in the introduction, since this paper first appeared, the full
three-loop integrated subtraction terms have become available. Specifically,
results for the three-loop $\mathcal{T}_{0}$ quark beam function in the
generalized large-$N_{c}$ approximation have appeared in ref. [91] and the
complete result has been calculated in ref. [94]. The three-loop beam
functions for $q_{T}$ have been calculated in refs. [92, 93]. In all cases,
the full calculations have confirmed our predictions of the eikonal terms at
three loops.

## Acknowledgments

We thank Goutam Das for fruitful discussions, V. Ravindran for providing us
with the results of ref. [111, 14], and Hua Xing Zhu for providing us with a
numerical version of the results in ref. [82]. We also thank Marius Wiesemann
and Leandro Cieri for clarifying discussions regarding refs. [160, 18]. J.M.
and F.T. thank the MIT Center for Theoretical Physics for hospitality and M.E.
thanks the DESY theory group for hospitality. This work was supported in part
by the Office of Nuclear Physics of the U.S. Department of Energy under
Contract No. DE-SC0011090, the Alexander von Humboldt Foundation through a
Feodor Lynen Research Fellowship, the Deutsche Forschungsgemeinschaft (DFG)
under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306,
and the PIER Hamburg Seed Project PHM-2019-01. M.E. and F.T. also thank the
Mainz Institute of Theoretical Physics of the DFG Cluster of Excellence
PRISMA+ (Project ID 39083149) for hospitality while portions of this work were
performed.

##  Appendix A Plus distributions and Fourier transforms

Here, we summarize the definitions and relations for plus distributions.

###  A.1 One-dimensional plus distributions

Following ref. [107], we denote plus distributions as

$\displaystyle\mathcal{L}_{n}(x)$
$\displaystyle=\biggl{[}\frac{\theta(x)\ln^{n}x}{x}\biggr{]}_{+}=\lim_{\epsilon\to
0}\frac{\mathrm{d}}{\mathrm{d}x}\Bigl{[}\theta(x-\epsilon)\frac{\ln^{n+1}x}{n+1}\Bigr{]}\,,$
$\displaystyle\mathcal{L}^{a}(x)$
$\displaystyle=\biggl{[}\frac{\theta(x)}{x^{1-a}}\biggr{]}_{+}=\lim_{\epsilon\to
0}\frac{\mathrm{d}}{\mathrm{d}x}\Bigl{[}\theta(x-\epsilon)\frac{x^{a}-1}{a}\Bigr{]}\,,$
(A.1)

such that

$\mathcal{L}_{n}(x>0)=\frac{\ln^{n}x}{x}\,,\quad\mathcal{L}^{a}(x>0)=\frac{1}{x^{1-a}}\,,\qquad\int_{0}^{1}\\!\mathrm{d}x\,\mathcal{L}_{n}(x)=\int_{0}^{1}\\!\mathrm{d}x\,\mathcal{L}^{a}(x)=0\,.$
(A.2)

For distributions with dimensionful arguments we define

$\displaystyle\mathcal{L}_{n}(k,\mu)$
$\displaystyle=\frac{1}{\mu}\mathcal{L}_{n}\Bigl{(}\frac{k}{\mu}\Bigr{)}\,,\qquad\mathcal{L}_{n}(t,\mu^{2})=\frac{1}{\mu^{2}}\mathcal{L}_{n}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}\,.$
(A.3)

Using $\mathcal{L}^{a}(x)$ we further define the distribution

$\displaystyle\mathcal{V}_{a}(x)=\frac{e^{-\gamma_{E}a}}{\Gamma(1+a)}\bigl{[}a\mathcal{L}^{a}(x)+\delta(x)\bigr{]}\,,\qquad\mathcal{V}_{a}(k,\mu)=\frac{1}{\mu}\mathcal{V}_{a}\Bigl{(}\frac{k}{\mu}\Bigr{)}\,,$
(A.4)

which satisfies the group property

$(\mathcal{V}_{a}\mathcal{V}_{b})(k,\mu)=\int\\!\mathrm{d}k^{\prime}\,\mathcal{V}_{a}(k-k^{\prime},\mu)\,\mathcal{V}_{b}(k^{\prime},\mu)=\mathcal{V}_{a+b}(k,\mu)\,,\qquad\mathcal{V}_{0}(k,\mu)=\delta(k)\,.$
(A.5)

The $\mu$ dependence of $\mathcal{V}_{a}(k,\mu)$ is given by

$\mathcal{V}_{a}(k,\mu)=\Bigl{(}\frac{\mu^{\prime}}{\mu}\Bigr{)}^{a}\,\mathcal{V}_{a}(k,\mu^{\prime})\,,\qquad\mu\frac{\mathrm{d}}{\mathrm{d}\mu}\mathcal{V}_{a}(k,\mu)=-a\mathcal{V}_{a}(k,\mu)\,.$
(A.6)

Expanding $\mathcal{V}_{a}(k,\mu)$ in powers of $a$ we find

$\displaystyle\mathcal{V}_{a}(k,\mu)$
$\displaystyle=\delta(k)+a\,\mathcal{L}_{0}(k,\mu)+\frac{a^{2}}{2!}\bigl{[}2\mathcal{L}_{1}(k,\mu)-\zeta_{2}\delta(k)\bigr{]}$
$\displaystyle\quad+\frac{a^{3}}{3!}\bigl{[}3\mathcal{L}_{2}(k,\mu)-3\zeta_{2}\mathcal{L}_{0}(k,\mu)+2\zeta_{3}\delta(k)\bigr{]}+\mathcal{O}(a^{4})\,.$
(A.7)

The Fourier transformation of $\mathcal{V}_{a}(k,\mu)$ is given by

$\displaystyle\int\\!\mathrm{d}k\,e^{-\mathrm{i}ky}\,\mathcal{V}_{a}(k,\mu)=e^{-aL_{y}}\,,\quad\int\\!\frac{\mathrm{d}y}{2\pi}\,e^{\mathrm{i}ky}\,e^{-aL_{y}}=\mathcal{V}_{a}(k,\mu)\,,\quad
L_{y}=\ln(\mathrm{i}y\mu e^{\gamma_{E}})\,.$ (A.8)

###  A.2 Two-dimensional plus distributions for ${\vec{q}}_{T}$

$L_{b}^{n}$ | $\text{FT}^{-1}[L_{b}^{n}]$
---|---
$1$ | $\delta^{(2)}({\vec{q}}_{T})$
$L_{b}$ | $-\mathcal{L}_{0}({\vec{q}}_{T},\mu)$
$L_{b}^{2}$ | $+2\mathcal{L}_{1}({\vec{q}}_{T},\mu)$
$L_{b}^{3}$ | $-3\mathcal{L}_{2}({\vec{q}}_{T},\mu)-4\zeta_{3}\delta^{(2)}({\vec{q}}_{T})$
$L_{b}^{4}$ | $+4\mathcal{L}_{3}({\vec{q}}_{T},\mu)+16\zeta_{3}\mathcal{L}_{0}({\vec{q}}_{T},\mu)$
$L_{b}^{5}$ | $-5\mathcal{L}_{4}({\vec{q}}_{T},\mu)-80\zeta_{3}\mathcal{L}_{1}({\vec{q}}_{T},\mu)-48\zeta_{5}\delta^{(2)}({\vec{q}}_{T})$
$L_{b}^{6}$ | $+6\mathcal{L}_{5}({\vec{q}}_{T},\mu)+240\zeta_{3}\mathcal{L}_{2}({\vec{q}}_{T},\mu)+288\zeta_{5}\mathcal{L}_{0}({\vec{q}}_{T},\mu)+160\zeta_{3}^{2}\delta^{(2)}({\vec{q}}_{T})$
Table 1: Fourier transform of $L_{b}^{n}=\ln^{n}(b_{T}^{2}\mu^{2}/b_{0}^{2})$
to ${\vec{q}}_{T}$ space for $n\leq 6$, as given by eq. (A.12).

Following ref. [106], we define two-dimensional plus distributions in
${\vec{q}}_{T}$ as

$\mathcal{L}_{n}({\vec{q}}_{T},\mu)=\frac{1}{\pi\mu^{2}}\mathcal{L}_{n}\biggl{(}\frac{q_{T}^{2}}{\mu^{2}}\biggr{)}\,,$
(A.9)

where $\mathcal{L}_{n}(x)$ is defined as above in eq. (A.1), such that

$\int_{|{\vec{q}}_{T}|\leq\mu}\mathrm{d}^{2}{\vec{q}}_{T}\,\mathcal{L}_{n}({\vec{q}}_{T},\mu)=\pi\int_{0}^{\mu^{2}}\mathrm{d}q_{T}^{2}\,\frac{1}{\pi\mu^{2}}\mathcal{L}_{n}\biggl{(}\frac{q_{T}^{2}}{\mu^{2}}\biggr{)}=0\,.$
(A.10)

The cumulant for a generic cut $\lvert{\vec{q}}_{T}\rvert\leq
q_{T}^{\mathrm{cut}}$ follows to be

$\int_{|{\vec{q}}_{T}|\leq