A toolbox for $q_{T}$ and 0-jettiness subtractions at $\hbox {N}^3\hbox {LO}$

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 $\hbox {N}^3\hbox {LO}$, which are an important ingredient for matching $\hbox {N}^3\hbox {LO}$ 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 $\hbox {N}^3\hbox {LL}'$ and $\hbox {N}^4\hbox {LL}$ 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 $\hbox {N}^3\hbox {LO}$ in the $z\rightarrow 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\rightarrow 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 $\hbox {N}^3\hbox {LO}$. 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 $\hbox {N}^3\hbox {LO}$ 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 \rightarrow 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 $\hbox {N}^3\hbox {LO}$. Our results are necessary ingredients for carrying out the resummation for $q_T$ and $\mathcal {T}_0$ at $\hbox {N}^3\hbox {LL}'$ and $\hbox {N}^4\hbox {LL}$ order. For the associated $q_T$ and $\mathcal {T}_0$ subtraction methods [35, 38, 39], we provide the complete differential $\hbox {N}^3\hbox {LO}$ subtraction terms in analytic form. The structure and required ingredients for the $\mathcal {T}_0$ subtractions at $\hbox {N}^3\hbox {LO}$ were already discussed in [39]. The $q_T$ slicing method at $\hbox {N}^3\hbox {LO}$ was also considered in [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 $\hbox {N}^3\hbox {LO}+\hbox {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

$$\begin{aligned} \kappa _a(q_a)\, \kappa _b(q_b)\rightarrow L(q) \quad \text {with}\quad q_a^\mu + q_b^\mu = q^\mu \,. \end{aligned}$$

(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\rightarrow L$. The $q_{a,b}^\mu$ are lightlike Born reference (label) momenta given by

$$\begin{aligned} q_{a,b}^\mu = \omega _{a,b} \frac{n_{a,b}^\mu }{2} \,,\quad 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} \,, \end{aligned}$$

(1.2)

where $\hat{z}$ is the beam direction and $E_\mathrm{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 \rightarrow 0$ as [48, 49]

$$\begin{aligned} \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 ] \,. \end{aligned}$$

(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 \rightarrow 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\rightarrow 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 \rightarrow H$ in the $m_t\rightarrow \infty$ limit, $b\bar{b} \rightarrow 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 [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., [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 \rightarrow 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 [85, 86]. Here, we exploit this to obtain for the first time the three-loop beam function coefficients in the $z\rightarrow 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\rightarrow 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 $\hbox {N}^3\hbox {LO}$ but fails starting at N$^4$LO. Since our results capture the complete singular structure for $z\rightarrow 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 [85,86,87,88], while only partial results were known at $\hbox {N}^3\hbox {LO}$ [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\rightarrow 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\rightarrow 1$ limit again provided important checks.

For $q_T$, a similar study of the logarithmic structure at $\hbox {N}^3\hbox {LO}$ was performed in [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 [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 $\hbox {N}^3\hbox {LO}$.

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 Sect. 2 and for $q_T$ in Sect. 3. The application to $\mathcal {T}_0$ and $q_T$ subtractions at $\hbox {N}^3\hbox {LO}$ is discussed in Sect. 4. Readers primarily interested in this application may directly skip ahead to Sect. 4. We conclude in Sect. 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

$$\begin{aligned} F(\mu ) = \sum _{n = 0}^\infty F^{(n)}(\mu )\, \Bigl [\frac{\alpha _s(\mu )}{4\pi }\Bigr ]^n \,. \end{aligned}$$

(1.4)

All anomalous dimensions $\gamma ^i_{x}(\alpha _s)$ and the QCD splitting functions are expanded as

$$\begin{aligned} \gamma ^i_x(\alpha _s) = \sum _{n = 0}^\infty \gamma _{x\, n}^i\, \Bigl (\frac{\alpha _s}{4\pi }\Bigr )^{n+1} \,, \quad P_{ij}(z,\mu ) = \sum _{n=0}^\infty P_{ij}^{(n)}(z) \biggl [ \frac{\alpha _s(\mu )}{4\pi } \biggr ]^{n+1} \,. \end{aligned}$$

(1.5)

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

$$\begin{aligned} (I^{(m)} P^{(n)})_{ij}(z)&\equiv \sum _k \int \frac{\mathrm {d}z'}{z'}\, I^{(m)}_{ik}\Bigl (\frac{z}{z'}\Bigr ) P^{(n)}_{kj}(z') \,, \nonumber \\ [\mathcal {I}(t, \mu ) P^{(n)}]_{ij}(z)&\equiv \sum _k \int \frac{\mathrm {d}z'}{z'}\, \mathcal {I}_{ik}\Bigl (t, \frac{z}{z'}, \mu \Bigr ) P^{(n)}_{kj}(z') \,, \end{aligned}$$

(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

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

(1.7)

For Fourier-type convolutions, we use the notation

$$\begin{aligned} (FG)(k, \mu )&\equiv \int \mathrm {d}k'\, F(k - k', \mu )\,G(k', \mu ) \,,\nonumber \\ \bigl [F\,\mathcal {I}_{ij}(z)\bigr ](t, \mu ^2)&\equiv \int \mathrm {d}t'\, F(t - t', \mu ^2)\, \mathcal {I}_{ij}(t', z, \mu ^2) \,. \end{aligned}$$

(1.8)

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

We denote logarithmic plus distributions as

$$\begin{aligned} \mathcal {L}_{n}(x) = \biggl [ \frac{\theta (x) \ln ^n x}{x}\biggr ]_+ \quad \text {with}\quad \int _0^1 \mathrm {d}x \, \mathcal {L}_n(x) = 0 \,. \end{aligned}$$

(1.9)

For dimensionful arguments, we define

$$\begin{aligned} \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 ) \,. \end{aligned}$$

(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 [49, 69, 95]. 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 [69, 95]

$$\begin{aligned} \mathcal {T}_0 = \sum _i \mathrm{min}\left\{ \frac{2 q_a \cdot k_i}{Q_a},\, \frac{2 q_b \cdot k_i}{Q_b}\right\} \,, \end{aligned}$$

(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 [49, 88], are

$$\begin{aligned} \begin{aligned}&\text {leptonic:}\quad&Q_a&= Q_b = \sqrt{\omega _a \omega _b} = Q \,,\quad&\mathcal {T}_0^\mathrm{lep}&= \sum _i \min \Bigl \{ e^Y n_a \cdot k_i \,,\, e^{-Y} n_b \cdot k_i \Bigr \} \\&\text {hadronic:}\quad&Q_{a,b}&= \omega _{a,b} = Q \, e^{\pm Y} \,,&\mathcal {T}_0^\mathrm{cm}&= \sum _i \min \Bigl \{ n_a \cdot k_i \,,\, n_b \cdot k_i \Bigr \} \,. \end{aligned} \end{aligned}$$

(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]

$$\begin{aligned} \frac{\mathrm {d}\sigma }{\mathrm {d}Q^2 \mathrm {d}Y \mathrm {d}\mathcal {T}_0}&= \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 ) \nonumber \\&\quad \times \Bigl [1 + \mathcal {O}\Bigl (\frac{\mathcal {T}_0}{Q}\Bigr )\Bigr ] \,. \end{aligned}$$

(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 [49, 87].

The beam function appearing in Eq. (2.3) is the inclusive virtuality-dependent ($\hbox {SCET}_\mathrm {I}$) 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 [85,86,87,88], and they are being calculated at $\hbox {N}^3\hbox {LO}$ [83, 84].

The soft function in Eq. (2.3) is the hemisphere soft function for incoming Wilson lines. It is closely related to the hemisphere soft function for $e^+e^-\rightarrow$ 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.3) receives power corrections suppressed by $\mathcal {T}_0/Q$, as indicated. In addition, starting at N$^4$LO it also receives contributions from perturbative Glauber-gluon exchanges, which are not captured by Eq. (2.3) [104, 105].

2.2 $\mathcal {T}_0$ soft function

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

$$\begin{aligned} \mu \frac{\mathrm {d}}{\mathrm {d}\mu } S_i(k,\mu )&= \int \mathrm {d}k' \gamma _S^i(k-k',\mu )\, S_i(k',\mu ) \equiv (\gamma _S^i\, S_i)(k, \mu ) \,, \nonumber \\ \gamma _S^i(k,\mu )&= 4 \Gamma _{\mathrm{cusp}}^i[\alpha _s(\mu )]\, \mathcal {L}_0(k, \mu ) + \gamma _S^i[\alpha _s(\mu )] \delta (k) \,, \end{aligned}$$

(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.4) 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.4) 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,

$$\begin{aligned} \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 (2 m \beta _{n-m} + \gamma _{S\,n-m}^i \bigr ) S_i^{(m)}(k, \mu ) \Bigr ] , \end{aligned}$$

(2.5)

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

$$\begin{aligned} S_i^{(n+1)}(k,\mu )&= \int _{k|_+}^\mu \frac{\mathrm {d}\mu '}{\mu '} \sum _{m=0}^n \Bigl [4 \Gamma ^i_{n-m} \bigl (\mathcal {L}_0 S_i^{(m)}\bigr )(k, \mu ) + \bigl (2 m \beta _{n-m} + \gamma _{S\,n-m}^i \bigr )S_i^{(m)}(k, \mu ) \Bigr ] \nonumber \\&\quad + \delta (k)\, s_i^{(n+1)} \,, \end{aligned}$$

(2.6)

where the soft function boundary coefficients are defined by

$$\begin{aligned} S_i^{(n)}(k,\mu = k|_+) = \delta (k)\, s_i^{(n)} \quad \text {with}\quad s_i^{(0)}= 1 \,. \end{aligned}$$

(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]

$$\begin{aligned} \mathcal {L}_n(k, \mu _0 = k|_+)&= 0 \,, \nonumber \\ \delta (k) \ln ^{n+1}\frac{\mu _0}{\mu } \bigg \vert _{\mu _0 = k|_+}&= (n+1)\mathcal {L}_n(k, \mu ) \,, \nonumber \\ \int _{\mu _0 = k|_+}^\mu \frac{\mathrm {d}\mu '}{\mu '}\,\mathcal {L}_n(k, \mu ')&= -\frac{1}{n+1} \mathcal {L}_{n+1}(k, \mu ) \,. \end{aligned}$$

(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 '$ integral in Eq. (2.6), 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\ge 1$ in Eq. (2.6), we need the convolutions [107]

$$\begin{aligned} (\mathcal {L}_0 \mathcal {L}_0)(k, \mu )&= 2\mathcal {L}_1(k, \mu ) - \zeta _2\, \delta (k) \,, \nonumber \\ (\mathcal {L}_0 \mathcal {L}_1)(k, \mu )&= \frac{3}{2}\mathcal {L}_2(k, \mu ) - \zeta _2\mathcal {L}_0(k, \mu ) + \zeta _3\, \delta (k) \,, \nonumber \\ (\mathcal {L}_0 \mathcal {L}_2)(k, \mu )&= \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) \,, \nonumber \\ (\mathcal {L}_0 \mathcal {L}_3)(k, \mu )&= \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) \,. \end{aligned}$$

(2.9)

Starting from the LO result, $s_i^{(0)}= 1$, Eq. (2.6) yields up to two loops

$$\begin{aligned} S_i^{(0)}(k,\mu )&= \delta (k) \,, \nonumber \\ S_i^{(1)}(k,\mu )&= - \mathcal {L}_1(k,\mu )\,4 \Gamma _0^i - \mathcal {L}_0(k,\mu )\, \gamma ^i_{S\,0} + \delta (k)\, s_i^{(1)}\,,\nonumber \\ S_i^{(2)}(k,\mu )&= \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) \nonumber \\&\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 ] \nonumber \\&\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)}\,, \end{aligned}$$

(2.10)

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

$$\begin{aligned} S^{(3)}(k, \mu ) = \delta (k)\, s_i^{(3)}+ \sum _{\ell = 0}^5 S_{i,\ell }^{(3)}\, \mathcal {L}_\ell (k, \mu ) \end{aligned}$$

(2.11)

with the coefficients of the logarithmic distributions given by

$$\begin{aligned} S_{i,5}^{(3)}&= -8 (\Gamma ^i_0)^3 \nonumber \\ S_{i,4}^{(3)}&= -\frac{10}{3}(\Gamma ^i_0)^2\, (4\beta _0 + 3\gamma _{S\,0}^i) \nonumber \\ S_{i,3}^{(3)}&= 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 ] \nonumber \\ S_{i,2}^{(3)}&= 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} \nonumber \\&\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 ] \nonumber \\ S_{i,1}^{(3)}&= 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 ] \nonumber \\&\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)}\nonumber \\&\quad - 4 \Gamma ^i_2 - 4 \Gamma ^i_0\, s_i^{(2)}\nonumber \\ S_{i,0}^{(3)}&= 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 ] \nonumber \\&\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)}\nonumber \\&\quad - \gamma ^i_{S\,2} - (4 \beta _0 + \gamma ^i_{S\,0}) s_i^{(2)}\,. \end{aligned}$$

(2.12)

This agrees with a corresponding numerical expression in [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.3). 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,

$$\begin{aligned} S_i(k, \mu ) = \int \mathrm {d}k' \, S_i(k', \mu _S) \, U_S^i(k - k', \mu _S, \mu ) \,, \end{aligned}$$

(2.13)

where the evolution factor $U_S^i(k, \mu _S, \mu )$ encodes the solution of Eq. (2.4), with $U_S^i(k, \mu _S, \mu _S) = \delta (k)$. It can be found, e.g., in [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

$$\begin{aligned} (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' \, S_i(k', \mu _S) \, U_S^i(k - k', \mu _S, \mu ) \,, \end{aligned}$$

(2.14)

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

Fig. 1

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 Fig. 1, 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 $\hbox {N}^3\hbox {LL}'$ (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 (N$^3$LL) 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 $\hbox {N}^3\hbox {LL}'$. 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 Fig. 1 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]

$$\begin{aligned} \mu \frac{\mathrm {d}}{\mathrm {d}\mu } B_i(t, x, \mu )&= \int \mathrm {d}t'\, \gamma _B^i(t-t',\mu )\, B_i(t', x, \mu ) \,, \nonumber \\ \gamma _B^i(t, \mu )&= -2 \Gamma _{\mathrm{cusp}}^i[\alpha _s(\mu )]\,\mathcal {L}_0(t, \mu ^2) + \gamma _B^i[\alpha _s(\mu )]\,\delta (t) \,, \end{aligned}$$

(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{QC}D}$, the beam function satisfies an OPE in terms of standard PDFs [49, 87]

$$\begin{aligned} 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{QC}D}^2}{t}\biggr )\biggr ] \,, \end{aligned}$$

(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]

$$\begin{aligned} \mu \frac{\mathrm {d}}{\mathrm {d}\mu }\mathcal {I}_{ij}(t,z,\mu ) = \sum _k \int \mathrm {d}t'\, \frac{\mathrm {d}z'}{z'}\, \mathcal {I}_{ik}\Bigl (t - t', \frac{z}{z'},\mu \Bigr ) \Bigl [ \gamma _B^i(t', \mu )\, \mathbf {1}_{kj}(z') -\delta (t')\, 2P_{kj}(z',\mu ) \Bigr ] , \end{aligned}$$

(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 [85, 88] to NNLO. Following the same procedure as in Sect. 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

$$\begin{aligned} \mathcal {I}_{ij}^{(n+1)}(t, z, \mu ^2)&= \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) \nonumber \\&\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 \} \,, \end{aligned}$$

(2.18)

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

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

(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

$$\begin{aligned} \mathcal {I}_{ij}^{(0)}(t,z,\mu ^2)&= \delta (t)\, \mathbf {1}_{ij}(z) \,, \nonumber \\ \mathcal {I}_{ij}^{(1)}(t,z,\mu ^2)&= \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) \,, \nonumber \\ \mathcal {I}_{ij}^{(2)}(t,z,\mu ^2)&= \mathcal {L}_3(t, \mu ^2)\, \frac{(\Gamma _0^i)^2}{2}\, \mathbf {1}_{ij}(z) \nonumber \\&\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) + 3 P^{(0)}_{ij}(z) \Bigr ] \nonumber \\&\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) \nonumber \\&\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 \} \nonumber \\&\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) \nonumber \\&\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 \} \nonumber \\&\quad + \delta (t)\, I^{(2)}_{ij}(z) \,, \end{aligned}$$

(2.20)

which agrees with [39, 85, 88]. 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 [85, 86].¹

Plugging Eq. (2.20) back into Eq. (2.18), we obtain the $\hbox {N}^3\hbox {LO}$ result

$$\begin{aligned} \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) \end{aligned}$$

(2.21)

with the coefficients

$$\begin{aligned} \mathcal {I}_{ij,5}^{(3)}(z)&= \frac{(\Gamma ^i_0)^3}{8} \mathbf {1}_{ij}(z) \nonumber \\ \mathcal {I}_{ij,4}^{(3)}(z)&= \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 ] \nonumber \\ \mathcal {I}_{ij,3}^{(3)}(z)&= \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) \nonumber \\&\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 \} \nonumber \\ \mathcal {I}_{ij,2}^{(3)}(z)&= \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} \nonumber \\&\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) \nonumber \\&\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) \nonumber \\&\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 ] \nonumber \\ \mathcal {I}_{ij,1}^{(3)}(z)&= \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 ] \nonumber \\&\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 ] \nonumber \\&\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) \nonumber \\&\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) \nonumber \\&\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) \nonumber \\ \mathcal {I}_{ij,0}^{(3)}(z)&= \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 ] \nonumber \\&\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 ] \nonumber \\&\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) \nonumber \\&\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) \nonumber \\&\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) \,. \end{aligned}$$

(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.22) using the MT package [108]. To calculate $(I^{(2)}P^{(0)})_{ij}(z)$, this required employing the identity

$$\begin{aligned} \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) \,. \end{aligned}$$

(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

$$\begin{aligned} (B_i \otimes U_B^i)_\mathrm {cut}(t_\mathrm {cut}, x, \mu ) = \int ^{t_\mathrm {cut}} \mathrm {d}t \, \int \mathrm {d}t' \, B_i(t, x, \mu _B) \, U_B^i(t - t', \mu _B, \mu ) \,. \end{aligned}$$

(2.24)

The explicit expression for the beam function evolution kernel $U_B^i(t, \mu _B, \mu )$ can be found in [87, 88].

Fig. 2

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 Fig. 2, 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 $\hbox {N}^3\hbox {LL}'$ 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 $\hbox {N}^3\hbox {LL}'$, while at lower orders, they amount to a higher-order effect. Numerical results to N$^3$LL with PDFs and $\alpha _s$ evolution at corresponding lower orders can be found in [86]. The residual dependence on $\mu _B$ is noticeably reduced by about a factor of two at $\hbox {N}^3\hbox {LL}'$ 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\rightarrow 1$ limit. As was already pointed out and exploited in the NNLO calculation in [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\rightarrow 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$,

$$\begin{aligned} q^- \equiv \bar{n}\cdot q = \sqrt{Q^2 + q_T^2}\, e^{+Y} \,, \quad q^+ \equiv n\cdot q = \sqrt{Q^2 + q_T^2}\, e^{-Y} \,, \quad x_{\mp } = \frac{q^{\mp }}{E_{\mathrm{cm}}} \,. \end{aligned}$$

(2.25)

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

$$\begin{aligned} \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 ] \,. \end{aligned}$$

(2.26)

Here, $H_{ab}$ is the same hard function as in Eq. (2.3), and $B_b$ is the same inclusive beam function as in Eq. (2.3). The threshold PDF $f^\mathrm {thr}_a(x)$ encodes the extraction of parton a from the proton for $x \rightarrow 1$.

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

$$\begin{aligned} \frac{\mathrm {d}\sigma }{\mathrm {d}q^- \mathrm {d}q^+}&= \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 ] \nonumber \\&\quad \times S^\mathrm {thr}_i(k^-, k^+, \mu ) \, \bigl [ 1 + \mathcal {O}(1-x_-, 1-x_+) \bigr ] \,. \end{aligned}$$

(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

$$\begin{aligned} S^\mathrm {thr}_i(k^-, k^+, \mu ) = S^\mathrm {thr}_i(e^{+y} k^-, e^{-y} k^+, \mu ) \,. \end{aligned}$$

(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 [14, 19, 82, 111, 117] 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.27) 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_+ \rightarrow 1$, it must reproduce Eq. (2.27). As a result, the eikonal $x\rightarrow 1$ limit of the beam function must coincide with the threshold soft function [97],

$$\begin{aligned} 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 ] \,. \end{aligned}$$

(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],

$$\begin{aligned} \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 ] \,. \end{aligned}$$

(2.30)

This relation captures all terms in $\mathcal {I}_{ij}(t, z, \mu )$ that are singular for $z \rightarrow 1$, while power corrections have at most an integrable singularity for $z \rightarrow 1$. Notably, the beam function becomes flavor diagonal as $z \rightarrow 1$, while off-diagonal 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 \rightarrow 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 [97], Eq. (2.30) is explicitly confirmed at two loops by comparison with [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

$$\begin{aligned} I_{ij}^{(1)}(z)&= \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 ] \,, \nonumber \\ I_{ij}^{(2)}(z)&= \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 ] \nonumber \\&\quad \quad + \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 ] \,, \nonumber \\ I_{ij}^{(3)}(z)&= \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 \nonumber \\&\quad \quad + \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 ] \nonumber \\&\quad \quad + \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 ] \nonumber \\&\quad \quad + \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 \nonumber \\&\quad \quad + \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 ] \nonumber \\&\quad \quad + \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} \nonumber \\&\quad \quad + \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 ] \nonumber \\&\quad \quad + \delta (1-z) \, s^{\mathrm {thr}{(3)}}_i \biggr \} + \mathcal {O}\bigl [(1-z)^0\bigr ] \,. \end{aligned}$$

(2.31)

The boundary coefficients $s^{\mathrm {thr}(n)}_i$ of the threshold soft function are given in Eq. (B.8). 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.5).

The three-loop result in Eq. (2.31) 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.22). The fact that the leading $z\rightarrow 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 [85, 118]. 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 [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 )$.² Note that in contrast to Eqs. (1.3) and (2.3), 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.15), and hence, Eq. (2.22) also holds for $\tilde{B}_i$ just with different boundary coefficients $\tilde{I}_{ij}^{(n)}(z)$. In the limit $z \rightarrow 1$, the difference between $B_i$ and $\tilde{B}_i$ becomes power suppressed in $1-z$. As a result, the $z\rightarrow 1$ limit of the modified $\tilde{I}_i^{(n)}$ is also given by Eq. (2.31).

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 Fig. 3, 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\rightarrow 0$ divergent behavior in the gluon coefficient, which is not reproduced by its $z\rightarrow 1$ expansion. If desired, it can be captured by including the leading $z\rightarrow 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.

Fig. 3

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,

$$\begin{aligned} I_{ij,\mathrm {approx}}^{(n)}(z) = 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) \,, \end{aligned}$$

(2.32)

where the NLP coefficient itself is approximated as

$$\begin{aligned} I_{ij,\mathrm {approx}}^{\mathrm {NLP}(n)}(z) = - (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 ] \,, \end{aligned}$$

(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

$$\begin{aligned} I_{ij}^{\mathrm {NLP}(n)}(z) = \sum _{k=0}^{2n-1} c_{n,k}^\mathrm {NLP} \ln ^k(1-z) \,. \end{aligned}$$

(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} = - 4 C_i$ and $c_{2,3}^\mathrm {NLP} = -8 C_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} = 4 C_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

$$\begin{aligned} \Delta = \Delta _{1} \oplus \Delta _{2} = \sqrt{\Delta _{1}^2 + \Delta _{2}^2} \,. \end{aligned}$$

(2.35)

Fig. 4

Approximate ansatzes for the NNLO (top) and $\hbox {N}^3\hbox {LO}$ (bottom) kernels, in the u-quark (left) and gluon (right) channels

In Fig. 4, we show the approximate kernel at NNLO (top) and $\hbox {N}^3\hbox {LO}$ (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 $\hbox {N}^3\hbox {LO}$, 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 $\hbox {N}^3\hbox {LO}$ 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 |\vec {q}_T| \ll Q$ was first established by Collins, Soper, and Sterman (CSS) [48, 119, 120], and was later elaborated on in [121,122,123,124]. The factorization for $q_T$ was also shown within the framework of SCET in [117, 125,126,127]. Sometimes it is also referred to as transverse-momentum dependent (TMD) factorization. We write the factorized cross section as

$$\begin{aligned} \frac{\mathrm {d}\sigma }{\mathrm {d}Q^2 \mathrm {d}Y \mathrm {d}^2 \vec {q}_T}&= \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 ] \,, \nonumber \\ \frac{\mathrm {d}\tilde{\sigma }^\mathrm {sing}(\vec {b}_T)}{\mathrm {d}Q^2 \mathrm {d}Y}&= \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 ) \,. \end{aligned}$$

(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

$$\begin{aligned} L_b = \ln \frac{b_T^2 \mu ^2}{b_0^2} \,,\quad b_0 = 2 e^{-\gamma _E} \,. \end{aligned}$$

(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 [117], which up to two loops gives results equivalent to the $\eta$ regulator of [127, 128].

The beam function appearing in Eq. (3.1) is the inclusive transverse-momentum dependent ($\hbox {SCET}_\mathrm{II}$) 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

$$\begin{aligned} \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 )} \,,\quad \zeta = \omega ^2 = (x E_{\mathrm{cm}})^2 \,, \end{aligned}$$

(3.3)

as is done, e.g., in [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 = x P^-$ 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],

$$\begin{aligned} \nu \frac{\mathrm {d}}{\mathrm {d}\nu } \tilde{B}_i(x, \vec {b}_T, \mu , \nu )&= \tilde{\gamma }_{\nu ,B}^i(b_T,\mu )\, \tilde{B}_i(x, \vec {b}_T, \mu , \nu ) \,,\nonumber \\ \nu \frac{\mathrm {d}}{\mathrm {d}\nu } \tilde{S}_i(b_T, \mu , \nu )&= \tilde{\gamma }_{\nu ,S}^i(b_T,\mu )\, \tilde{S}_i(b_T, \mu , \nu ) \,, \end{aligned}$$

(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

$$\begin{aligned} \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 ) \,, \end{aligned}$$

(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 _\mathrm{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$,

$$\begin{aligned} \mu \frac{\mathrm {d}}{\mathrm {d}\mu } \tilde{\gamma }^i_\nu (b_T,\mu ) = - 4 \Gamma _{\mathrm{cusp}}^i[\alpha _s(\mu )] \,, \end{aligned}$$

(3.6)

which predicts its all-order structure in $b_T$ and $\mu$. Similar to the $\mathcal {T}_0$ soft function in Sect. 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

$$\begin{aligned} \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 '}{\mu '}\, \tilde{\gamma }_\nu ^{i\,(m)}(b_T,\mu ') + \tilde{\gamma }_{\nu \,n+1}^i \,, \end{aligned}$$

(3.7)

where the nonlogarithmic boundary coefficients are defined as

$$\begin{aligned} \tilde{\gamma }_{\nu \,n}^i = \tilde{\gamma }_\nu ^{i\,(n)}(b_T, \mu =b_0/b_T) \,. \end{aligned}$$

(3.8)

The result up to three loops is

$$\begin{aligned} \tilde{\gamma }_\nu ^{i\,(0)}(b_T,\mu )&= -L_b\, 2\Gamma _0^i + \tilde{\gamma }_{\nu \,0}^i \,,\nonumber \\ \tilde{\gamma }_\nu ^{i\,(1)}(b_T,\mu )&= - 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 \,,\nonumber \\ \tilde{\gamma }_\nu ^{i\,(2)}(b_T,\mu )&= -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 ) \nonumber \\&\quad + \tilde{\gamma }_{\nu \,2}^i \,. \end{aligned}$$

(3.9)

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

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,

$$\begin{aligned} \mu \frac{\mathrm {d}}{\mathrm {d}\mu } \tilde{S}_i(b_T,\mu ,\nu )&= \tilde{\gamma }_S^i(\mu ,\nu )\, \tilde{S}_i(b_T,\mu ,\nu ) \,,\\ \nu \frac{\mathrm {d}}{\mathrm {d}\nu } \tilde{S}_i(b_T,\mu ,\nu )&= \tilde{\gamma }_\nu ^i(b_T, \mu )\, \tilde{S}_i(b_T,\mu ,\nu ) \,.\nonumber \end{aligned}$$

(3.10)

The perturbative structure of $\gamma _\nu ^i$ is discussed in Sect. 3.2. The $\mu$ anomalous dimension has the all-order structure

$$\begin{aligned} \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 )] \,, \end{aligned}$$

(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.10) order by order in $\alpha _s$, we obtain the coupled RGEs

$$\begin{aligned} \mu \frac{\mathrm {d}}{\mathrm {d}\mu }{\tilde{S}_i^{(n+1)}(b_T,\mu ,\nu )}&= \sum _{m=0}^{n} \Bigl ( 4 \Gamma ^i_{n-m} \ln \frac{\mu }{\nu } + 2 m \beta _{n-m} + \tilde{\gamma }_{S\,n-m}^i \Bigr ) \tilde{S}_i^{(m)}(b_T,\mu ,\nu ) \,, \nonumber \\ \nu \frac{\mathrm {d}}{\mathrm {d}\nu }{\tilde{S}_i^{(n+1)}(b_T,\mu ,\nu )}&= \sum _{m=0}^{n} \tilde{\gamma }_\nu ^{i\,(n-m)}(b_T,\mu )\, \tilde{S}_i^{(m)}(b_T,\mu ,\nu ) \,. \end{aligned}$$

(3.12)

These are easily integrated to give

$$\begin{aligned} \tilde{S}_i^{(n+1)}(b_T,\mu ,\nu )&= \sum _{m=0}^{n} \biggl [ \int _{b_0/b_T}^\mu \frac{\mathrm {d}\mu '}{\mu '}\, \Bigl ( 4 \Gamma ^i_{n-m} \ln \frac{\mu '}{\nu } + 2 m \beta _{n-m} + \tilde{\gamma }_{S\,n-m}^i\Bigr ) \tilde{S}_i^{(m)}(b_T,\mu ',\nu ) \nonumber \\&\quad + \int _{b_0/b_T}^\nu \frac{\mathrm {d}\nu '}{\nu '}\, \tilde{\gamma }_{\nu \,n-m}^i\, \tilde{S}_i^{(m)}(b_T,b_0/b_T,\nu ') \biggr ] + \tilde{s}_{i}^{(n+1)} \,, \end{aligned}$$

(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.13) are defined as

$$\begin{aligned} \tilde{s}_{i}^{(n)} = \tilde{S}_i^{(n)}(b_T, \mu = b_0/b_T, \nu = b_0/b_T) \,. \end{aligned}$$

(3.14)

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

$$\begin{aligned} L_b = \ln \frac{b_T^2 \mu ^2}{b_0^2} \,,\quad b_0 = 2 e^{-\gamma _E} \,,\quad L_\nu = \ln \frac{\mu }{\nu } \,, \end{aligned}$$

(3.15)

Eq. (3.13) yields up to two loops

$$\begin{aligned} \tilde{S}_i^{(0)}(b_T,\mu ,\nu )&= 1 \,, \nonumber \\ \tilde{S}_i^{(1)}(b_T,\mu ,\nu )&= - 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)}\,, \nonumber \\ \tilde{S}_i^{(2)}(b_T,\mu ,\nu )&= 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 ) \nonumber \\&\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 ) \nonumber \\&\quad + \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 ] \nonumber \\&\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 ] \nonumber \\&\quad + \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 \} \nonumber \\&\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)}\,. \end{aligned}$$

(3.16)

At three loops, we write the result as

$$\begin{aligned} \tilde{S}_i^{(3)}(b_T,\mu ,\nu ) = \sum _{\ell =0}^6 \tilde{S}_{i,\ell }^{(3)}(L_\nu ) \, L_b^\ell \,, \end{aligned}$$

(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

$$\begin{aligned} \tilde{S}_{i,6}^{(3)}(L_\nu )&= -\frac{(\Gamma ^i_0)^3}{48} \,,\nonumber \\ \tilde{S}_{i,5}^{(3)}(L_\nu )&= L_\nu \, \frac{(\Gamma ^i_0)^3}{4} + \frac{(\Gamma ^i_0)^2}{6} \beta _0 \,,\nonumber \\ \tilde{S}_{i,4}^{(3)}(L_\nu )&= -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 ) \,,\nonumber \\ \tilde{S}_{i,3}^{(3)}(L_\nu )&= 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 ) \nonumber \\&\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 ) \,,\nonumber \\ \tilde{S}_{i,2}^{(3)}(L_\nu )&= 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 ] \nonumber \\&\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)}\,,\nonumber \\ \tilde{S}_{i,1}^{(3)}(L_\nu )&= - 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 ) \nonumber \\&\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)}\,,\nonumber \\ \tilde{S}_{i,0}^{(3)}(L_\nu )&= -L_\nu \bigl (\tilde{\gamma }_{\nu \,2}^i + \tilde{\gamma }_{\nu \,1}^i\,\tilde{s}_i^{(1)}\bigr ) + \tilde{s}_i^{(3)}\,. \end{aligned}$$

(3.18)

Equations (3.16) and (3.18) agree with [82, 136]. 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,

$$\begin{aligned} \tilde{S}_i(b_T, \mu , \nu )&= \tilde{S}_i(b_T, \mu _S, \nu _S) \, \tilde{U}_S^i(b_T, \mu _S, \mu , \nu _S, \nu ) \,,\nonumber \\ \tilde{U}_S^i(b_T, \mu _S, \mu , \nu _S, \nu )&= \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 '}{\mu '}\, \tilde{\gamma }_S^i(\mu ',\nu )\biggr ] \,, \end{aligned}$$

(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 Fig. 5, 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 $\hbox {N}^3\hbox {LL}'$ (solid orange), normalized to the NNLL$'$ result at the central scale. The three-loop finite term is included in Fig. 5, so the NNLL$'$ and $\hbox {N}^3\hbox {LL}'$ 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 $\hbox {N}^3\hbox {LL}'$, it further stabilizes over a wider range of scales. As in Sect. 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.

Fig. 5

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

$$\begin{aligned} \mu \frac{\mathrm {d}}{\mathrm {d}\mu }{\tilde{B}_i(x,\vec {b}_T,\mu ,\nu )}&= \tilde{\gamma }_B^i(\mu ,\nu /\omega )\, \tilde{B}_i(x,\vec {b}_T,\mu ,\nu ) \,, \nonumber \\ \nu \frac{\mathrm {d}}{\mathrm {d}\nu }{\tilde{B}_i(x,\vec {b}_T,\mu ,\nu )}&= -\frac{1}{2}\tilde{\gamma }_\nu ^i(b_T,\mu )\, \tilde{B}_i(x,\vec {b}_T,\mu ,\nu ) \,, \end{aligned}$$

(3.20)

where the $\nu$ anomalous dimension is discussed in Sect. 3.2, and the $\mu$ anomalous dimension has the all-order form

$$\begin{aligned} \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 )] \,, \end{aligned}$$

(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{QC}D}$, the TMD beam function satisfies an OPE in terms of standard PDFs [48],

$$\begin{aligned} \tilde{B}_q(x,b_T,\mu ,\nu )&= \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{QC}D})\bigr ] \,,\nonumber \\ \tilde{B}^{\rho \lambda }_g(x,\vec {b}_T,\mu ,\nu )&= \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 ] \nonumber \\&\quad \times f_j\Bigl (\frac{x}{z},\mu \Bigr ) \bigl [1 + \mathcal {O}(b_T \Lambda _{\mathrm{QC}D})\bigr ] \,. \end{aligned}$$

(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.20), 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.20) by taking the evolution of the PDFs into account,

$$\begin{aligned} \mu \frac{\mathrm {d}}{\mathrm {d}\mu }{\tilde{\mathcal {I}}_{ij}(z,b_T,\mu ,\nu /\omega )}&= \sum _k \int \frac{\mathrm {d}z'}{z'}\, \tilde{\mathcal {I}}_{ik}\Bigl (\frac{z}{z'},b_T,\mu ,\nu /\omega \Bigr ) \bigl [ \tilde{\gamma }_B^i(\mu ,\nu /\omega ) \mathbf {1}_{kj}(z') - 2 P_{kj}(z',\mu ) \bigr ] \,,\nonumber \\ \nu \frac{\mathrm {d}}{\mathrm {d}\nu }{\tilde{\mathcal {I}}_{ij}(z,b_T,\mu ,\nu /\omega )}&= -\frac{1}{2}\tilde{\gamma }_\nu ^i(b_T,\mu )\, \tilde{\mathcal {I}}_{ij}(z,b_T,\mu ,\nu /\omega ) \,.\end{aligned}$$

(3.23)

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

$$\begin{aligned} \tilde{\mathcal {I}}_{ij}^{(n+1)}(z,b_T,\mu ,\nu /\omega )&= \sum _{m=0}^n \biggl [ \int _{b_0/b_T}^\mu \frac{\mathrm {d}\mu '}{\mu '} \Bigl ( 2 \Gamma ^i_{n-m} \ln \frac{\nu }{\omega } + \tilde{\gamma }_{B\,n-m}^i + 2 m \beta _{n-m} \Bigr ) \tilde{\mathcal {I}}_{ij}^{(m)}(z,b_T,\mu ',\nu /\omega ) \nonumber \\&\quad - 2 \int _{b_0/b_T}^\mu \frac{\mathrm {d}\mu '}{\mu '}\, \bigl [\tilde{\mathcal {I}}^{(m)}(b_T,\mu ',\nu /\omega ) P^{(n-m)}\bigr ]_{ij}(z) \nonumber \\&\quad - \int _{\omega }^\nu \frac{\mathrm {d}\nu '}{\nu '}\, \frac{\tilde{\gamma }_{\nu \,n-m}^i}{2}\, \tilde{\mathcal {I}}_{ij}^{(m)}(z,b_T,b_0/b_T,\nu '/\omega ) \biggr ] + \tilde{I}_{ij}^{(n+1)}(z) \,, \end{aligned}$$

(3.24)

where the nonlogarithmic boundary coefficients are defined as

$$\begin{aligned} \tilde{I}_{ij}^{(n+1)}(z) = \tilde{\mathcal {I}}_{ij}^{(n+1)}(z,b_T,\mu = b_0/b_T,\nu /\omega = 1) \,. \end{aligned}$$

(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

$$\begin{aligned} \tilde{\mathcal {I}}_{ij}^{(0)}(z,b_T,\mu ,\nu /\omega )&= \mathbf {1}_{ij}(z) \,,\nonumber \\ \tilde{\mathcal {I}}_{ij}^{(1)}(z,b_T,\mu ,\nu /\omega )&= 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) \,,\nonumber \\ \tilde{\mathcal {I}}_{ij}^{(2)}(z,b_T,\mu ,\nu /\omega )&= 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) \nonumber \\&\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 \} \nonumber \\&\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) \nonumber \\&\quad + L_\omega \, \frac{\tilde{\gamma }_{\nu \,0}^i}{2}\, P^{(0)}_{ij}(z) - P^{(1)}_{ij}(z) \nonumber \\&\quad + \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 \} \nonumber \\&\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) \,, \end{aligned}$$

(3.26)

where we abbreviated

$$\begin{aligned} L_b = \ln \frac{b_T^2\mu ^2}{b_0^2} \,, \quad b_0 = 2 e^{-\gamma _E} \,,\quad L_\omega = \ln \frac{\nu }{\omega }\,. \end{aligned}$$

(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 [136] for quark and gluon beam functions in terms of the results of [134], and were directly calculated at NNLO using the exponential regulator for the quark case in [138].

At three loops, we write

$$\begin{aligned} \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 \,, \end{aligned}$$

(3.28)

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

$$\begin{aligned} \tilde{\mathcal {I}}_{ij,3}^{(3)}(z,L_\omega )&= \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 ] \nonumber \\&\quad + \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) \nonumber \\&\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) \nonumber \\&\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) \,,\nonumber \\ \tilde{\mathcal {I}}_{ij,2}^{(3)}(z,L_\omega )&= \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 ] \nonumber \\&\quad + \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) \nonumber \\&\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) \nonumber \\&\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) \nonumber \\&\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) \,,\nonumber \\ \tilde{\mathcal {I}}_{ij,1}^{(3)}(z,L_\omega )&= \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) \nonumber \\&\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) \nonumber \\&\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) \,,\nonumber \\ \tilde{\mathcal {I}}_{ij,0}^{(3)}(z,L_\omega )&= - 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) \,, \end{aligned}$$

(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.26) and (3.29) with the help of the MT package [108]. In contrast to [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 $\hbox {N}^3\hbox {LO}$ cross section. They are given by

$$\begin{aligned} \tilde{\mathcal {J}}_{gj}^{(0)}(z,b_T,\mu ,\nu /\omega )&= 0 \,,\nonumber \\ \tilde{\mathcal {J}}_{gj}^{(1)}(z,b_T,\mu ,\nu /\omega )&= \tilde{J}_{gj}^{(1)}(z) = 4 C_j\, \frac{1-z}{z} \,,\nonumber \\ \tilde{\mathcal {J}}_{gj}^{(2)}(z,b_T,\mu ,\nu /\omega )&= 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 ] \nonumber \\&\quad - L_\omega \, \frac{\tilde{\gamma }_{\nu \,0}^g}{2}\, \tilde{J}_{gj}^{(1)}(z) + \tilde{J}_{gj}^{(2)}(z) \,. \end{aligned}$$

(3.30)

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

$$\begin{aligned} \tilde{J}_{gj}^{(2)}(z)&= I^{\prime (2)}_{gi}(z) \nonumber \\&= -\delta ^L C_{g \leftarrow j}^{(2;0,0)}(z) - \frac{1}{2} \tilde{s}_g^{{(1)}} \tilde{J}_{gj}^{(1)}(z) \,. \end{aligned}$$

(3.31)

In the first line of Eq. (3.31), $I^{\prime (2)}_{gi}(z)$ is the two-loop boundary term as given in [142]. In the second line of Eq. (3.31), $\tilde{s}_g^{{(1)}} = - 2 C_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 [143].

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

$$\begin{aligned} \tilde{B}_i(x, \vec {b}_T, \mu , \nu )&= \tilde{B}_i(x, \vec {b}_T, \mu _B, \nu _B) \, \tilde{U}_B^i(\omega , b_T, \mu _B, \mu , \nu _B, \nu ) \,,\nonumber \\ \tilde{U}_B^i(\omega , b_T, \mu _B, \mu , \nu _B, \nu )&= \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 '}{\mu '} \tilde{\gamma }_B^i(\mu ',\nu /\omega )\biggr ] \,. \end{aligned}$$

(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$.

Fig. 6

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 Fig. 6, we show the residual dependence on $(\mu _B, \nu _B)$ at NLL$'$ (dotted green), NNLL$'$ (dashed blue), and $\hbox {N}^3\hbox {LL}'$ 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\rightarrow 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\rightarrow 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 \rightarrow 1$ and $x_b \rightarrow 1$. As $x_{a,b} \rightarrow 1$, initial state radiation is constrained to have energy $\lesssim \lambda _- \lambda _+ Q$, where

$$\begin{aligned} \lambda _-^2 \sim 1-x_a \quad \text {and} \quad \lambda _+^2 \sim 1-x_b \end{aligned}$$

(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 [117, 148]. Some consequences of the resulting consistency relations have already been explored in [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 [74, 149] 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]

$$\begin{aligned} \frac{\mathrm {d}\tilde{\sigma }(\vec {b}_T)}{\mathrm {d}Q^2\, \mathrm {d}Y}&= \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 ] \nonumber \\&\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 ) \nonumber \\&\quad \times \Bigl [ 1 + \mathcal {O}\Bigl (\frac{1}{b_T^2 \lambda _-^2 \lambda _+^2 Q^2}, \lambda _-^2, \lambda _+^2\Bigr ) \Bigr ] \,. \end{aligned}$$

(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

$$\begin{aligned} \frac{\mathrm {d}\tilde{\sigma }(\vec {b}_T)}{\mathrm {d}Q^2\, \mathrm {d}Y}&= \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 ] \nonumber \\&\quad \times \tilde{S}^\mathrm {thr}_i(k^-, k^+, b_T, \mu )\, \bigl [ 1 + \mathcal {O}(\lambda _-^2, \lambda _+^2) \bigr ] \,. \end{aligned}$$

(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 [150], albeit in a different context, and to $\mathcal {O}(\alpha _s^3)$ in [82]. By construction, it satisfies

$$\begin{aligned} \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 ) \,, \end{aligned}$$

(3.36)

where $S^\mathrm {thr}_i(k^-, k^+, \mu )$ is the double-differential threshold soft function appearing in Eq. (2.27).

Consistency relations Consistency between Eqs. (3.1) and (3.34) implies that the $x \rightarrow 1$ limit of the $q_T$ beam function is captured by the collinear-soft function [117, 148],

$$\begin{aligned} \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 ] \,. \end{aligned}$$

(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]

$$\begin{aligned} \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 ] \,. \end{aligned}$$

(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.34) and (3.35) implies [117, 148]

$$\begin{aligned} \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 ] \,, \end{aligned}$$

(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 \rightarrow 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,

$$\begin{aligned} \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 ] \,, \end{aligned}$$

(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 [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.9) and (A.7). Up to two loops, we have

$$\begin{aligned} \tilde{I}_{ij}^{(1)}(z)&= \mathcal {O}\bigl [(1-z)^0\bigr ] \,, \nonumber \\ \tilde{I}_{ij}^{(2)}(z)&= \delta _{ij} \, \frac{\tilde{\gamma }_{\nu \,1}^i}{2} \, \mathcal {L}_0(1-z) + \mathcal {O}\bigl [(1-z)^0\bigr ] \,, \end{aligned}$$

(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

$$\begin{aligned} \tilde{I}_{ij}^{(3)}(z)&= \delta _{ij} \, \frac{\tilde{\gamma }_{\nu \,2}^i}{2} \, \mathcal {L}_0(1-z) + \mathcal {O}\bigl [(1-z)^0\bigr ] \,, \nonumber \\ \tilde{I}_{ij}^{(4)}(z)&= \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 ] \,, \nonumber \\ \tilde{I}_{ij}^{(5)}(z)&= \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 ] \,, \nonumber \\ \tilde{I}_{ij}^{(6)}(z)&= \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 ] \nonumber \\&\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 ] \,. \end{aligned}$$

(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 Sect. 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 Fig. 7, 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\rightarrow 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 toward 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\rightarrow 1$. If desired, one could also include the leading $z\rightarrow 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.

Fig. 7

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 Sect. 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 N$^{n}$LO beam coefficient,

$$\begin{aligned} \tilde{I}_{ij,\mathrm {approx}}^{(n)}(z) =&\tilde{I}_{ij}^{(n)\,\mathrm {LP}}(z) + \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) \nonumber \\&- X_3 (1-z) \tilde{I}_{ij}^{(n)\,\mathrm {LP}}(z) \,. \end{aligned}$$

(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

$$\begin{aligned} \tilde{I}_{ij,\mathrm {reg}}^{(1)}(z)&= -\delta _{ij} C_i \,, \quad \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) \,. \end{aligned}$$

(3.44)

The ansatz in Eq. (3.43) 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

$$\begin{aligned} \gamma _X^g = 3 C_A - \beta _0 \,,\quad \gamma _X^q = 10 (C_F - C_A) \,. \end{aligned}$$

(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

$$\begin{aligned} \Delta = \Delta _{1} \oplus \Delta _{2} \oplus \Delta _{3} = \sqrt{\Delta _{1}^2 + \Delta _{2}^2 + \Delta _{3}^2} \,. \end{aligned}$$

(3.46)

Fig. 8

Approximate ansatzes for the NNLO (top) and $\hbox {N}^3\hbox {LO}$ (bottom) kernels, in the u-quark (left) and gluon (right) channels

In Fig. 8, we show the approximate kernel at NNLO (top) and $\hbox {N}^3\hbox {LO}$ (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 $\hbox {N}^3\hbox {LO}$, 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 $\hbox {N}^3\hbox {LO}$ 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 [18], where the $q_T$ subtraction was first applied at $\hbox {N}^3\hbox {LO}$ 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 $\hbox {N}^3\hbox {LO}$ subtractions

The factorization in Eq. (1.3) fully describes the limit $\tau \rightarrow 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, 38, 39, 151,152,153,154,155,156]. 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 [39].

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

$$\begin{aligned} \sigma (X) = \int \mathrm {d}\tau \, \frac{\mathrm {d}\sigma (X)}{\mathrm {d}\tau } \,, \quad \sigma (X, \tau _\mathrm {cut}) = \int ^{\tau _\mathrm {cut}} \mathrm {d}\tau \, \frac{\mathrm {d}\sigma (X)}{\mathrm {d}\tau } \,, \end{aligned}$$

(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 \rightarrow 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,

$$\begin{aligned} \sigma (X) = \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 ] \,, \end{aligned}$$

(4.2)

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

$$\begin{aligned} \frac{\mathrm {d}\sigma (X)}{\mathrm {d}\tau } = \frac{\mathrm {d}\sigma ^\mathrm {sing}(X)}{\mathrm {d}\tau }\,[1 + \mathcal {O}(\tau )] \,,\quad \frac{\mathrm {d}\sigma ^\mathrm {sub}(X)}{\mathrm {d}\tau } = \frac{\mathrm {d}\sigma ^\mathrm {sing}(X)}{\mathrm {d}\tau }\,[1 + \mathcal {O}(\tau )] \,, \end{aligned}$$

(4.3)

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

$$\begin{aligned} \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}) \,. \end{aligned}$$

(4.4)

By construction, the integrand in square brackets in Eq. (4.2) contains at most integrable singularities for $\tau \rightarrow 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\rightarrow L+1$-parton process at one lower order. In practice, one always has a small IR cutoff $\delta$ on the $\tau$ integral,

$$\begin{aligned} \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 ) \,, \end{aligned}$$

(4.5)

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

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

(4.6)

which is neglected for $\delta \rightarrow 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 \rightarrow 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, 47, 157]. The differential subtractions at $\hbox {N}^3\hbox {LO}$ are a key ingredient for using this method to combine $\hbox {N}^3\hbox {LO}$ 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 [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,

$$\begin{aligned} \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}) \,. \end{aligned}$$

(4.7)

The practical advantage of the slicing method is that it allows one to readily turn an existing $L+1$-jet N$^{n-1}$LO calculation into a N$^n$LO calculation for L, and so most implementations use this approach [18, 38, 158,159,160,161,162,163,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}\rightarrow 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.2) exhibits power corrections that grow like $e^{|Y|}$ 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

$$\begin{aligned} \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) \,, \end{aligned}$$

(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

$$\begin{aligned} \frac{\mathrm {d}\sigma ^\mathrm {sing}(\Phi _0)}{\mathrm {d}\tau }&= \mathcal {C}_{-1}(\Phi _0)\,\delta (\tau ) + \sum _{n\ge 0} \mathcal {C}_n(\Phi _0)\, \mathcal {L}_n(\tau ) \nonumber \\&= \sum _{m \ge 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 \,. \end{aligned}$$

(4.9)

Their integral over $\tau \le \tau _\mathrm {cut}$ immediately follows as

$$\begin{aligned} \sigma ^\mathrm {sing}(\Phi _0, \tau _\mathrm {cut})&= \mathcal {C}_{-1}(\Phi _0) + \sum _{n\ge 0} \mathcal {C}_n(\Phi _0)\, \frac{\ln ^n \tau _\mathrm {cut}}{n+1} \nonumber \\&= \sum _{m \ge 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 \,. \end{aligned}$$

(4.10)

The differential subtractions are given by using Eq. (4.9) 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.10).

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 \rightarrow \lambda \tau$ moves contributions from $\mathcal {C}_n(\Phi _0)$ to $\mathcal {C}_{m<n}(\Phi _0)$. This freedom was used in [39] to absorb all terms with $n \ge 0$ in Eq. (4.10) 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.10) and take

$$\begin{aligned} \tau \equiv \frac{\mathcal {T}_0}{Q} \quad (\text {for } \mathcal {T}_0) \,,\quad \tau \equiv \frac{q_T^2}{Q^2} \quad (\text {for }q_T) \,. \end{aligned}$$

(4.11)

The m-loop subtraction coefficients $\mathcal {C}_n^{(m)}(\Phi _0)$ directly follow from expanding Eq. (2.3) for $\mathcal {T}_0$ or Eq. (3.1) for $q_T$ to mth order in $\alpha _s$. For the three-loop coefficients, this yields

$$\begin{aligned} \mathcal {C}_{-1}^{(3)}(\Phi _0)&= 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)} \,, \nonumber \\ \mathcal {C}_{n\ge 0}^{(3)}(\Phi _0)&= \sum _{k = 1}^3 H^{(3-k)}(\Phi _0)\, \bigl [B_a(x_a) B_b(x_b) S\bigr ]_n^{(k)} \,, \end{aligned}$$

(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 mth 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,

$$\begin{aligned} \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 ) \,. \end{aligned}$$

(4.13)

The results for the subtraction coefficients $\mathcal {C}_n(\Phi _0)$ in Eq. (4.12) 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.12) for $\mathcal {T}_0$ requires rescaling and convolving the plus distributions in the beam and soft functions, as discussed in [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 \le 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 [35] was based on the $q_T$ resummation framework of [175], where the canonical $\vec {b}_T$-space logarithms are replaced by

$$\begin{aligned} L_b \rightarrow \tilde{L}_b \equiv \ln \Bigl (\frac{b_T^2 \mu ^2}{b_0^2} + 1 \Bigr ) \,. \end{aligned}$$

(4.14)

This form is also used, e.g., in [18, 160]. 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 [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 $\hbox {N}^3\hbox {LO}$ beam functions, we employ consistency between different factorization limits to derive their leading eikonal limit $I^{(3)}_{ij}(z\rightarrow 1)$, i.e., the full singular limit of the beam functions as $z\rightarrow 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 $\hbox {N}^3\hbox {LL}'$ and $\hbox {N}^4\hbox {LL}$ order. In particular, they are important for extending the $q_T$ and $\mathcal {T}_0$ subtraction methods to $\hbox {N}^3\hbox {LO}$, 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 $\hbox {N}^3\hbox {LO}$ $+$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 [91], and the complete result has been calculated in [94]. The three-loop beam functions for $q_T$ have been calculated in [92, 93]. In all cases, the full calculations have confirmed our predictions of the eikonal terms at three loops.

Appendices

A Plus distributions and Fourier transforms

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

1.1 A.1 One-dimensional plus distributions

Following [107], we denote plus distributions as

$$\begin{aligned} \mathcal {L}_{n}(x)&= \biggl [ \frac{\theta (x) \ln ^n x}{x}\biggr ]_+ = \lim _{\epsilon \rightarrow 0} \frac{\mathrm {d}}{\mathrm {d}x} \Bigl [ \theta (x-\epsilon ) \frac{\ln ^{n+1} x}{n+1} \Bigr ] \,, \\ \mathcal {L}^{a}(x)&= \biggl [ \frac{\theta (x)}{x^{1-a}} \biggr ]_+ = \lim _{\epsilon \rightarrow 0} \frac{\mathrm {d}}{\mathrm {d}x} \Bigl [ \theta (x-\epsilon ) \frac{x^a - 1}{a} \Bigr ] \,,\nonumber \end{aligned}$$

(A.1)

such that

$$\begin{aligned} \mathcal {L}_n(x> 0) = \frac{\ln ^n x}{x} \,, \quad \mathcal {L}^a(x > 0) = \frac{1}{x^{1-a}} \,, \quad \int _0^1 \mathrm {d}x \, \mathcal {L}_n(x) = \int _0^1 \mathrm {d}x \, \mathcal {L}^a(x) = 0 \,. \end{aligned}$$

(A.2)

For distributions with dimensionful arguments, we define

$$\begin{aligned} \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 ) \,.\end{aligned}$$

(A.3)

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

$$\begin{aligned} \mathcal {V}_a(x) = \frac{e^{-\gamma _E a}}{\Gamma (1+a)} \bigl [ a \mathcal {L}^a(x) + \delta (x) \bigr ] \,, \quad \mathcal {V}_a(k, \mu ) = \frac{1}{\mu } \mathcal {V}_a \Bigl ( \frac{k}{\mu } \Bigr ) \,, \end{aligned}$$

(A.4)

which satisfies the group property

$$\begin{aligned} (\mathcal {V}_a \mathcal {V}_b)(k, \mu ) = \int \mathrm {d}k'\, \mathcal {V}_a(k - k',\mu ) \, \mathcal {V}_b(k',\mu ) = \mathcal {V}_{a+b}(k,\mu ) \,, \quad \mathcal {V}_0(k, \mu ) = \delta (k) \,. \end{aligned}$$

(A.5)

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

$$\begin{aligned} \mathcal {V}_a(k, \mu ) = \Bigl ( \frac{\mu '}{\mu } \Bigr )^a \, \mathcal {V}_a(k, \mu ') \,, \quad \mu \frac{\mathrm {d}}{\mathrm {d}\mu } \mathcal {V}_a(k, \mu ) = -a \mathcal {V}_a(k, \mu ) \,. \end{aligned}$$

(A.6)

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

$$\begin{aligned} \mathcal {V}_a(k, \mu )&= \delta (k) + a\,\mathcal {L}_0(k, \mu ) + \frac{a^2}{2!}\bigl [2\mathcal {L}_1(k, \mu ) - \zeta _2\delta (k) \bigr ] \nonumber \\&\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) \,. \end{aligned}$$

(A.7)

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

$$\begin{aligned} \int \mathrm {d}k \, e^{-\mathrm {i}k y} \, \mathcal {V}_a(k, \mu ) = e^{-a L_y} \,, \quad \int \frac{\mathrm {d}y}{2\pi }\, e^{\mathrm {i}k y} \, e^{-a L_y} = \mathcal {V}_a(k, \mu ) \,, \quad L_y = \ln (\mathrm {i}y \mu e^{\gamma _E}) \,. \end{aligned}$$

(A.8)

1.2 A.2 Two-dimensional plus distributions for ${\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 \le 6$, as given by Eq. (A.12)

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

$$\begin{aligned} \mathcal {L}_n(\vec {q}_T, \mu ) = \frac{1}{\pi \mu ^2} \mathcal {L}_n\biggl (\frac{q_T^2}{\mu ^2}\biggr ) \,, \end{aligned}$$

(A.9)

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

$$\begin{aligned} \int _{|\vec {q}_T| \le \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 \,. \end{aligned}$$

(A.10)

The cumulant for a generic cut $|\vec {q}_T| \le q_T^\mathrm {cut}$ follows to be

$$\begin{aligned} \int _{|\vec {q}_T| \le q_T^\mathrm {cut}} \mathrm {d}^2\vec {q}_T\, \mathcal {L}_n(\vec {q}_T,\mu ) = \frac{\theta (q_T^\mathrm {cut})}{n+1} \ln ^{n+1} \frac{(q_T^\mathrm {cut})^2}{\mu ^2} \,. \end{aligned}$$

(A.11)

The Fourier transformation of $\mathcal {L}_n(\vec {q}_T, \mu )$ and its inverse are [106]

$$\begin{aligned} \int \mathrm {d}^2\vec {q}_T\, e^{-i \vec {q}_T\cdot \vec {b}_T} \mathcal {L}_n(\vec {q}_T,\mu )&= \frac{1}{n+1} \sum _{k=0}^{n+1} (-1)^k \left( {\begin{array}{c}n+1\\ k\end{array}}\right) R_2^{(n+1-k)} L_b^k \,, \nonumber \\ \int \frac{\mathrm {d}^2\vec {b}_T}{(2\pi )^2}\, e^{i \vec {q}_T\cdot \vec {b}_T} L_b^n&= \sum _{k=0}^{n-1} (-1)^{k+1} n \left( {\begin{array}{c}n-1\\ k\end{array}}\right) R_2^{(n-k-1)} \mathcal {L}_k(\vec {q}_T,\mu ) + R_2^{(n)} \delta ^{(2)}(\vec {q}_T) \,, \end{aligned}$$

(A.12)

where $L_b$ is the usual logarithm in Fourier space

$$\begin{aligned} L_b = \ln (b_T^2 \mu ^2 / b_0^2) \,,\quad b_0 = 2 e^{-\gamma _E} \,, \end{aligned}$$

(A.13)

and the coefficients $R_2^{(n)}$ in Eq. (A.12) are given by

$$\begin{aligned} R_2^{(n)} = \frac{\mathrm {d}^n}{\mathrm {d}a^n}\, e^{2\gamma _E a}\, \frac{\Gamma (1+a)}{\Gamma (1-a)}\, \bigg |_{a=0} \,. \end{aligned}$$

(A.14)

Up to $\hbox {N}^3\hbox {LO}$, we require the Fourier transforms of $L_b^n$ with $n \le 6$, which are summarized in Table 1.

B Threshold soft function

Here, we discuss the double-differential threshold soft function $S_i^\mathrm {thr}(k^-, k^+, \mu )$, which appears in the soft threshold factorization for the inclusive cross section in Eq. (2.27) and determines the eikonal limit of the $\mathcal {T}_0$ beam function in Eq. (2.31). We give its complete $\hbox {N}^3\hbox {LO}$ result in Appendix B.2 in terms of a convenient plus distribution basis defined in Appendix B.1. In Appendix B.3, we discuss how the three-loop coefficients are extracted from the known three-loop results for the closely related inclusive threshold soft function.

1.1 B.1 Plus distribution basis

A key property of the threshold soft function is that is invariant under the simultaneous rescaling $k^- \mapsto k^- e^{+y}$ and $k^+ \mapsto k^+ e^{-y}$, see Eq. (2.28). To make this property manifest, we define a basis of plus distributions in $k^\pm$ that individually satisfy this property,

$$\begin{aligned} \frac{\theta (k^-)\theta (k^+)}{\mu ^2} \Bigl (\frac{k^- k^+}{\mu ^2}\Bigr )^{-1+a}&= \biggl [\frac{\delta (k^-)}{a} +\sum _{n=0}^{\infty } \frac{a^n}{n!} \mathcal {L}_n(k^-,\mu )\biggr ] \biggl [\frac{\delta (k^+)}{a} +\sum _{m=0}^{\infty } \frac{a^m}{m!} \mathcal {L}_m(k^+,\mu )\biggr ] \nonumber \\&\equiv \frac{\delta (k^-, k^+)}{a^2} + \sum _{n=0}^{\infty } \frac{a^{n-1}}{n!}\mathcal {L}_{n}(k^-,k^+,\mu ) \,. \end{aligned}$$

(B.1)

Note that the leading $\delta (k^-, k^+)$ term multiplies a double pole in a. The second line implicitly defines the $\mathcal {L}_n(k^-, k^+, \mu )$ by the expansion of the first line in powers of a. They are by construction invariant under rescaling, because the left-hand side is. Explicitly, they are given by

$$\begin{aligned} \delta (k^-, k^+)&= \delta (k^-) \, \delta (k^+) \,, \nonumber \\ \mathcal {L}_{n}(k^-, k^+, \mu )&= \delta (k^-)\, \mathcal {L}_{n}(k^+,\mu ) + \mathcal {L}_{n}(k^-,\mu )\, \delta (k^+) \nonumber \\&\quad + n \sum _{m = 0}^{n-1} \left( {\begin{array}{c}n-1\\ m\end{array}}\right) \mathcal {L}_{m}(k^-, \mu ) \, \mathcal {L}_{n-1-m}(k^+, \mu ) \,. \end{aligned}$$

(B.2)

1.2 B.2 Three-loop result

The threshold soft function satisfies the all-order RGE

$$\begin{aligned} \mu \frac{\mathrm {d}}{\mathrm {d}\mu } S^\mathrm {thr}_i(k^-, k^+, \mu )&= \int \mathrm {d}\ell ^- \mathrm {d}\ell ^+ \, \gamma _\mathrm {thr}^i(k^- - \ell ^-, k^+ - \ell ^+, \mu ) \, S^\mathrm {thr}_i(\ell ^-, \ell ^+, \mu ) \,, \nonumber \\ \gamma _\mathrm {thr}^i(k^-, k^+, \mu )&= -2 \Gamma _{\mathrm{cusp}}^i[\alpha _s(\mu )] \, \mathcal {L}_0(k^-, k^+, \mu ) + \gamma _\mathrm {thr}^i[\alpha _s(\mu )] \, \delta (k^-, k^+) \,. \end{aligned}$$

(B.3)

Expanding the threshold soft function in $\alpha _s$ as

$$\begin{aligned} S^\mathrm {thr}_i(k^-, k^+, \mu ) = \sum _{n=0}^\infty S^{\mathrm {thr}(n)}_i(k^-, k^+, \mu )\, \Bigl [\frac{\alpha _s(\mu )}{4\pi }\Bigr ]^n \,, \end{aligned}$$

(B.4)

and suppressing all arguments for brevity, $S_{i}^{\mathrm {thr}(n)} \equiv S_{i}^{\mathrm {thr}(n)}(k^-, k^+, \mu )$, $\mathcal {L}_n \equiv \mathcal {L}_n (k^-, k^+, \mu )$, $\delta \equiv \delta (k^-, k^+)$, the three-loop solution of Eq. (B.3) takes the form

$$\begin{aligned} S_i^{\mathrm {thr}{(0)}}&= \delta \,, \nonumber \\ S_i^{\mathrm {thr}{(1)}}&= \mathcal {L}_1 \, \Gamma _0^i - \mathcal {L}_0 \, \frac{\gamma _{\mathrm {thr}\,0}^{i}}{2} + \delta \, s^{\mathrm {thr}{(1)}}_i \,, \nonumber \\ S_i^{\mathrm {thr}{(2)}}&= \mathcal {L}_3 \, \frac{(\Gamma _0^i)^2}{2} - \mathcal {L}_2 \, \frac{\Gamma _0^i}{2} \Bigl (\beta _0 + \frac{3}{2}\gamma _{\mathrm {thr}\,0}^{i} \Bigr ) \nonumber \\&\quad + \mathcal {L}_1 \Bigl [ - 2\zeta _2 (\Gamma _0^i)^2 + \Bigl (\beta _0 + \frac{\gamma _{\mathrm {thr}\,0}^{i}}{2} \Bigr ) \frac{\gamma _{\mathrm {thr}\,0}^{i}}{2} + \Gamma _1^i + \Gamma _0^i \, s^{\mathrm {thr}{(1)}}_i \Bigr ] \nonumber \\&\quad + \mathcal {L}_0 \Bigl [ \Gamma _0^i\bigl (2\zeta _3 \Gamma _0^i + \zeta _2 \gamma _{\mathrm {thr}\,0}^{i} \bigr ) - \frac{\gamma _{\mathrm {thr}\,1}^{i}}{2} - \Bigl (\beta _0 + \frac{\gamma _{\mathrm {thr}\,0}^{i}}{2}\Bigr ) s^{\mathrm {thr}{(1)}}_i \Bigr ] + \delta \, s^{\mathrm {thr}{(2)}}_i \,, \nonumber \\ S_i^{\mathrm {thr}{(3)}}&= \mathcal {L}_5\, \frac{(\Gamma ^i_0)^3}{8} - \mathcal {L}_4\, \frac{5}{8}(\Gamma ^i_0)^2 \Bigl (\frac{2}{3}\beta _0 + \frac{\gamma _{\mathrm {thr}\,0}^i}{2} \Bigr ) \nonumber \\&\quad + \mathcal {L}_3\, \Gamma ^i_0 \Bigl [ - 2\zeta _2 (\Gamma ^i_0)^2 + \frac{\beta _0^2}{3} + \Bigl (\frac{5}{3}\beta _0 + \frac{\gamma _{\mathrm {thr}\,0}^i}{2}\Bigr ) \frac{\gamma _{\mathrm {thr}\,0}^i}{2} + \Gamma _1^i + \frac{\Gamma _0^i}{2} s_i^{\mathrm {thr}{(1)}} \Bigr ] \nonumber \\&\quad + \mathcal {L}_2\, \biggl \{ (\Gamma ^i_0)^2 \Bigl [ 5\zeta _3 \Gamma ^i_0 + 3\zeta _2 (\beta _0 + \gamma ^i_{\mathrm {thr}\,0}) \Bigr ] - \Bigl (\beta _0 + \frac{3}{4}\gamma ^i_{\mathrm {thr}\,0}\Bigr )\Bigl (\beta _0\frac{\gamma ^i_{\mathrm {thr}\,0}}{2} + \Gamma ^i_1\Bigr ) \nonumber \\&\quad - \frac{(\gamma ^i_{\mathrm {thr}\,0})^3}{16} - \frac{\Gamma ^i_0}{2} \Bigl [ \beta _1 + \frac{3}{2}\gamma ^i_{\mathrm {thr}\, 1} + \Bigl (4\beta _0 + \frac{3}{2}\gamma ^i_{\mathrm {thr}\,0} \Bigr ) s_i^{\mathrm {thr}{(1)}} \Bigr ] \biggr \} \nonumber \\&\quad + \mathcal {L}_1\, \biggl \{ (\Gamma ^i_0)^2 \bigl [4\zeta _4\Gamma ^i_0 - \zeta _3(6\beta _0 + 4\gamma ^i_{\mathrm {thr}\,0}) \bigr ] - \zeta _2 \Gamma ^i_0 \bigl [ (3\beta _0 + \gamma ^i_{\mathrm {thr}\,0})\gamma ^i_{\mathrm {thr}\,0} + 4\Gamma ^i_1 \bigr ] \nonumber \\&\quad + \beta _0 \gamma ^i_{\mathrm {thr}\,1} + \frac{\gamma ^i_{\mathrm {thr}\,0}}{2} (\beta _1 + \gamma ^i_{\mathrm {thr}\,1}) + \Gamma ^i_2 \nonumber \\&\quad + \Bigl [ -2\zeta _2(\Gamma ^i_0)^2 + 2\beta _0^2 + \Bigl (3\beta _0 + \frac{\gamma ^i_{\mathrm {thr}\,0}}{2}\Bigr ) \frac{\gamma ^i_{\mathrm {thr}\,0}}{2} + \Gamma ^i_1 \Bigr ] s_i^{\mathrm {thr}{(1)}} + \Gamma ^i_0\, s_i^{\mathrm {thr}{(2)}} \biggr \} \nonumber \\&\quad + \mathcal {L}_0\, \biggl \{ (\Gamma ^i_0)^2 \bigl [ - \Gamma ^i_0 (8 \zeta _2 \zeta _3 - 6 \zeta _5) + 2\zeta _4 (\beta _0 - \gamma ^i_{\mathrm {thr}\,0}) \bigr ] + \zeta _3\Gamma ^i_0 \Bigl [ \Bigl (\beta _0 + \frac{\gamma ^i_{\mathrm {thr}\,0}}{2} \Bigr ) \gamma ^i_{\mathrm {thr}\,0} \nonumber \\&\quad + 4\Gamma ^i_1 \Bigr ] + \zeta _2 \bigl (\gamma ^i_{\mathrm {thr}\,0} \Gamma ^i_1 + \Gamma ^i_0 \gamma ^i_{\mathrm {thr}\,1} \bigr ) - \frac{\gamma ^i_{\mathrm {thr}\,2}}{2} + \Bigl [ (\Gamma ^i_0)^2 2\zeta _3 + \Gamma ^i_0\zeta _2 (2\beta _0 + \gamma ^i_{\mathrm {thr}\,0} ) \nonumber \\&\quad - \Bigl (\beta _1 + \frac{\gamma ^i_{\mathrm {thr}\,1}}{2}\Bigr ) \Bigr ] s_i^{\mathrm {thr}{(1)}} - \Bigl (2\beta _0 + \frac{\gamma ^i_{\mathrm {thr}\,0}}{2}\Bigr ) s_i^{\mathrm {thr}{(2)}} \biggr \} + \delta \, s^{\mathrm {thr}{(3)}}_i \,. \end{aligned}$$

(B.5)

Consistency of the factorization theorems in Eqs. (2.3), (2.26), and (2.27) implies

$$\begin{aligned} 2\gamma _B^i(\alpha _s) + \gamma _S^i(\alpha _s) = 2\gamma _f^i(\alpha _s) + \gamma _\mathrm {thr}^i(\alpha _s) = \gamma _f^i(\alpha _s) + \gamma _B^i(\alpha _s) \,, \end{aligned}$$

(B.6)

because the hard function is the same in all cases. Here, $\gamma _f^i(\alpha _s)$ is the coefficient of $\delta (1-z)$ in the PDF anomalous dimension Eq. (D.14). Solving Eq. (B.6) for $\gamma _\mathrm {thr}^i(\alpha _s)$, we find

$$\begin{aligned} \gamma _\mathrm {thr}^i(\alpha _s) = - \gamma _S^i(\alpha _s) \,, \quad \gamma _{\mathrm {thr}\,n}^i = - \gamma _{S\,n}^i \,, \end{aligned}$$

(B.7)

where the soft anomalous dimension coefficients $\gamma _{S\,n}^i$ are given in Eq. (D.7).

The boundary coefficients $s_i^{\mathrm {thr}(n)}$, which are defined as the coefficients of $\delta (k^-,k^+)$ in Eq. (B.5), are given by [11, 12]³

$$\begin{aligned} s^{\mathrm {thr}{(1)}}_i&= -C_i\, 2\zeta _2 \,, \nonumber \\ s^{\mathrm {thr}{(2)}}_i&= C_i \biggl [ C_i\, 21\zeta _4 + C_A \Bigl (\frac{208}{27} - 4 \zeta _2 - 10 \zeta _4 \Bigl ) + \beta _0 \Bigl (\frac{164}{27} - 5 \zeta _2 - \frac{10 \zeta _3}{3}\Bigr ) \biggr ] \,, \nonumber \\ s^{\mathrm {thr}{(3)}}_i&= C_i \biggl [ C_i^2 \Bigl (\frac{640}{3}\zeta _3^2 - \frac{499}{6}\zeta _6 \Bigr ) + C_i C_A \Bigl ( -\frac{416}{27}\zeta _2 - \frac{512}{9}\zeta _3 + \frac{188}{3}\zeta _4 + 224 \zeta _3^2 - 77 \zeta _6 \Bigr ) \nonumber \\&\quad + C_i \beta _0 \Bigl ( -\frac{328}{27}\zeta _2 - \frac{448}{9}\zeta _3 + \frac{235}{3}\zeta _4 + \frac{308}{3} \zeta _2 \zeta _3 - 64 \zeta _5 \Bigr ) \nonumber \\&\quad + C_A^2 \Bigl ( \frac{115895}{324} - \frac{45239}{486}\zeta _2 - \frac{23396}{81}\zeta _3 - \frac{334}{3}\zeta _4 + 240 \zeta _2 \zeta _3 - 224 \zeta _5 + \frac{1072}{9}\zeta _3^2 + \frac{4348}{27}\zeta _6 \Bigr ) \nonumber \\&\quad + C_A \beta _0 \Bigl ( - \frac{363851}{2916} + \frac{1043}{486}\zeta _2 - \frac{140}{81}\zeta _3 + \frac{230}{9}\zeta _4 - \frac{164}{3}\zeta _2 \zeta _3 + \frac{632}{9}\zeta _5 \Bigr ) \nonumber \\&\quad + \beta _0^2 \Bigl ( - \frac{64}{729} - \frac{34}{3} \zeta _2 - \frac{20}{27}\zeta _3 + \frac{31}{3}\zeta _4 \Bigr ) \nonumber \\&\quad + \beta _1 \Bigl ( \frac{42727}{972} - \frac{275}{18}\zeta _2 - \frac{1636}{81}\zeta _3 - \frac{76}{9} \zeta _4 + \frac{40}{3} \zeta _2 \zeta _3 - \frac{112}{9}\zeta _5 \Bigr ) \biggr ] \,.\end{aligned}$$

(B.8)

We have also checked that inserting the above coefficients into Eq. (B.5) and expanding against the Drell–Yan hard function, we reproduce the three-loop soft-virtual partonic cross section in [14, 111] in terms of $1 - z_a = k^-/(Qe^{+Y})$ and $1 - z_b = k^+/(Qe^{-Y})$.

1.3 B.3 Extraction method

The double-differential threshold soft function depends on the total lightcone momentum components $k^\pm$ of the soft hadronic final state. Equivalently, its Fourier transform

$$\begin{aligned} \hat{S}_i^\mathrm {thr}(b^+, b^-, \mu ) = \int \mathrm {d}k^- \mathrm {d}k^+ \, e^{+\mathrm {i}(k^-b^+ + k^+b^-)/2} \, S_i^\mathrm {thr}(k^-, k^+, \mu ) \,, \end{aligned}$$

(B.9)

depends on the time-like separation $(b^-n^\mu + b^+\bar{n}^\mu )/2$ between the Wilson lines in the soft matrix element.

Importantly, $\hat{S}_i^\mathrm {thr}(b^+, b^-, \mu )$ only depends on the product $b^+ b^-$ by the rescaling relation Eq. (2.28), and thus only depends on $b^+ b^- \mu ^2$ by dimensional analysis. On the other hand, the dependence on $\mu$ is fully predicted by the RGE Eq. (B.3), which in position space reads

$$\begin{aligned} \mu \frac{\mathrm {d}}{\mathrm {d}\mu } \hat{S}_i^\mathrm {thr}(b^+, b^-, \mu ) = \Bigl \{ 2\Gamma _{\mathrm{cusp}}^i[\alpha _s(\mu )] \, L_\mathrm {thr}(b^+, b^-, \mu ) + \gamma ^i_\mathrm {thr}[\alpha _s(\mu )] \Bigr \}\, \hat{S}_i^\mathrm {thr}(b^+, b^-, \mu ) \,. \end{aligned}$$

(B.10)

This implies that at any given order in perturbation theory, $\hat{S}_i^\mathrm {thr}(b^+, b^-, \mu )$ is a polynomial in

$$\begin{aligned} L_\mathrm {thr}(b^+, b^-, \mu ) \equiv \ln \Bigl ( -\frac{b^+ b^- \mu ^2 e^{2\gamma _E}}{4} - \mathrm {i}0 \Bigr ) \,. \end{aligned}$$

(B.11)

The relevant Fourier transforms between $L_\mathrm {thr}^n$ and $\mathcal {L}_n(k^-, k^+, \mu )$ follow from the one-dimensional Fourier transforms in Appendix B of [106], accounting for the relative factors of $-1/2$ in the Fourier exponent in Eq. (B.9).

A factorization analogous to Eq. (2.27) holds for the inclusive cross section $\mathrm {d}\sigma / \mathrm {d}Q^2$, where the corresponding inclusive threshold soft function $S_i^\mathrm {thr}(k^0, \mu )$ only depends on the total energy $k^0$ of soft radiation. In particular, $S_i^\mathrm {thr}(k^0, \mu )$ is the process-independent soft contribution to the inclusive partonic cross section $\sigma _{ab}(z)$ in the soft-virtual limit $z \rightarrow 1$, where $1-z = 2k^0/Q$. In position space, the inclusive threshold soft function $\hat{S}_i^\mathrm {thr}(b^0, \mu )$ is defined in terms of Wilson lines separated by $b^0 (n^\mu + \bar{n}^\mu )/2$, i.e., strictly along the time axis. This is a special case of Eq. (B.9), so the two position-space threshold soft functions are simply related by

$$\begin{aligned} \hat{S}_i^\mathrm {thr}(b^0, b^0, \mu ) = \hat{S}_i^\mathrm {thr}(b^0, \mu ) \,. \end{aligned}$$

(B.12)

This is of course equivalent to integrating over the longitudinal momentum $k^3$ of soft radiation. We stress that Eq. (B.12) cannot be used to approximate $\hat{S}_i^\mathrm {thr}(b^+, b^-, \mu )$ by taking $b^+ = b^-$ in general. This is because in Eq. (2.27) the $k^+$ and $k^-$ dependences are separately convolved with the PDFs, and thus, the rescaling property Eq. (2.28) is lost at the level of the cross section. See also Appendix D of [97] for further discussion of this point.

Inserting Eq. (B.12) into Eq. (B.10) implies that both threshold soft functions have the same noncusp anomalous dimension given by Eq. (B.7). Moreover, the position-space boundary coefficients of the double-differential soft function at $L_\mathrm {thr}= 0$, i.e., at $\mu = \mu _*\equiv +\mathrm {i}2 e^{-\gamma _E}/b^0$, are equal to the inclusive ones at the same scale. Hence, the double-differential threshold soft function can be constructed from the inclusive one.

The inclusive threshold soft function was calculated to three loops in [11, 12]. Here, we use the results of [12], where the three-loop soft function for $i = g$ is reported in exponentiated form,

$$\begin{aligned} \hat{S}_i^\mathrm {thr}( b^0, \mu _*) = \exp \biggl \{ \frac{C_i}{C_A} \Bigl [ \frac{\alpha _s(\mu _*)}{4\pi } c_{1\,\text {[12]}}^H + \frac{\alpha _s^2(\mu _*)}{(4\pi )^2} \Delta c_{2\,\text {[12]}}^H + \frac{\alpha _s^3(\mu _*)}{(4\pi )^3} \Delta c_{3\,\text {[12]}}^H \Bigr ] + \mathcal {O}(\alpha _s^4) \biggr \} . \end{aligned}$$

(B.13)

We have also exploited Casimir scaling to three loops to restore the dependence on $C_i$. Comparing Eq. (B.13) to the position-space solution of Eq. (B.10) at $L_\mathrm {thr}= 0$, we obtain Eq. (B.8) for the momentum-space boundary coefficients after performing the inverse Fourier transform.

C Collinear-soft function for the exponential regulator

In this appendix, we derive the all-order expression for the collinear-soft function using the exponential regulator, which leads to Eq. (3.40) in the main text.

We start by defining the complete Fourier transform of the fully differential threshold soft function

$$\begin{aligned} \hat{S}^\mathrm {thr}_i(b^+, b^-, b_T, \mu ) = \int \mathrm {d}^4 k\, e^{+\mathrm {i}b\cdot k}\, S^\mathrm {thr}_i(k^-, k^+, k_T, \mu ) \,, \end{aligned}$$

(C.1)

where $b^\mu = (b^+, b^-, \vec {b}_T)$ is the four-vector Fourier conjugate of $k^\mu = (k^+, k^-, \vec {k}_T)$ with $b \cdot k = b^+ k^-/2 + b^- k^+/2 - \vec {b}_T \cdot \vec {k}_T$. Correspondingly, we define the Fourier transform of $\tilde{\mathcal {S}}_i(k^\pm , b_T, \mu , \nu )$ with respect to its lightcone momentum argument $k^\pm$ as

$$\begin{aligned} \hat{\mathcal {S}}_i(b^+, b_T, \mu , \nu ) = \int \mathrm {d}k^- \, e^{+\mathrm {i}k^- b^+/2} \, \tilde{\mathcal {S}}_i(k^-, b_T, \mu , \nu ) \,, \end{aligned}$$

(C.2)

and analogously for $b^- \leftrightarrow b^+$ and $k^+ \leftrightarrow k^-$. Fully in position space, the consistency relation Eq. (3.39) reads

$$\begin{aligned} \hat{S}_i^\mathrm {thr}(b^+, b^-, b_T, \mu ) = \hat{\mathcal {S}}_i(b^+, b_T, \mu , \nu ) \, \hat{\mathcal {S}}_i(b^-, b_T, \mu , \nu ) \, \tilde{S}_i(b_T, \mu , \nu ) \, \Bigl [ 1 + \mathcal {O}\Bigl (\frac{b^+ b^-}{b_T^2}\Bigr ) \Bigr ] \,. \end{aligned}$$

(C.3)

In the exponential regulator scheme, the regulated $q_T$ soft function is defined as [82, 117]⁴

$$\begin{aligned} \tilde{S}_i(b_T, \mu , \nu ') = \lim _{\nu ' \rightarrow \infty } \hat{S}^\mathrm {thr}_i\Bigl ( \frac{\mathrm {i}b_0}{\nu '}, \frac{\mathrm {i}b_0}{\nu '}, b_T, \mu \Bigr ) \,, \end{aligned}$$

(C.4)

where we use $\nu '$ to distinguish it from the scale at which we later wish to evaluate the collinear-soft function. The prescription for taking the limit is to keep all nonvanishing terms. In particular, a logarithmic dependence of the right-hand side on $\nu '$ is to be kept. Inserting Eq. (C.3), we have

$$\begin{aligned} \tilde{S}_i(b_T, \mu , \nu ')&= \lim _{\nu ' \rightarrow \infty } \Bigl [ \hat{\mathcal {S}}_i\Bigl (\frac{\mathrm {i}b_0}{\nu '}, b_T, \mu , \nu \Bigr ) \, \hat{\mathcal {S}}_i\Bigl ( \frac{\mathrm {i}b_0}{\nu '}, b_T, \mu , \nu \Bigr ) \, \tilde{S}_i(b_T, \mu , \nu ) + \mathcal {O}\Bigl (\frac{1}{\nu '^2 b_T^2}\Bigr ) \Bigr ] \nonumber \\&= \tilde{S}_i(b_T, \mu , \nu ) \lim _{\nu ' \rightarrow \infty } \Bigl [ \hat{\mathcal {S}}_i\Bigl (\frac{\mathrm {i}b_0}{\nu '}, b_T, \mu , \nu \Bigr ) \, \hat{\mathcal {S}}_i\Bigl ( \frac{\mathrm {i}b_0}{\nu '}, b_T, \mu , \nu \Bigr ) \Bigr ] \nonumber \\&= \tilde{S}_i(b_T, \mu , \nu ) \, \hat{\mathcal {S}}_i\Bigl (\frac{\mathrm {i}b_0}{\nu '}, b_T, \mu , \nu \Bigr ) \, \hat{\mathcal {S}}_i\Bigl ( \frac{\mathrm {i}b_0}{\nu '}, b_T, \mu , \nu \Bigr ) \,. \end{aligned}$$

(C.5)

In the second line, we moved the $q_T$ soft function out of the limit, since it does not depend on $\nu '$, and dropped the power corrections. On the third line, we used that all dependence of the $\hat{\mathcal {S}}_i$ on $\nu '$ is logarithmic, so the limit is trivial. Because the exponential regulator is symmetric under an interchange of collinear-soft directions, we find

$$\begin{aligned} \hat{\mathcal {S}}_i^2\Bigl (\frac{\mathrm {i}b_0}{\nu '}, b_T, \mu , \nu \Bigr ) = \frac{\tilde{S}_i(b_T, \mu , \nu ')}{\tilde{S}_i(b_T, \mu , \nu )} = \exp \Bigl [ \tilde{\gamma }_\nu ^i(b_T, \mu ) \ln \frac{\nu '}{\nu } \Bigr ] \,,\end{aligned}$$

(C.6)

where the second equality follows from solving the rapidity RGE of the soft function between $\nu$ and $\nu '$ at fixed $\mu$. Assuming we are dealing with the $n_a$-collinear-soft function that depends on $b^+$, we can analytically continue back to $\nu ' = \mathrm {i}b_0/b^+ = 2 \mathrm {i}/(b^+ e^{\gamma _E})$, leaving

$$\begin{aligned} \hat{\mathcal {S}}_i(b^+, b_T, \mu , \nu ) = \exp \Bigl [ -\frac{1}{2} \tilde{\gamma }_\nu ^i(b_T, \mu ) \ln (-\mathrm {i}b^+ \nu e^{\gamma _E}/2) \Bigr ] \,. \end{aligned}$$

(C.7)

Evaluating the inverse Fourier transform using Eq. (A.8), we find the following all-order relation for the momentum-space $n_a$-collinear-soft function in the exponential regulator scheme,

$$\begin{aligned} \tilde{\mathcal {S}}_i(k^-, b_T, \mu , \nu ) = \mathcal {V}_{\tilde{\gamma }_\nu ^i(b_T, \mu )/2}(k^-, \nu ) \,, \end{aligned}$$

(C.8)

and identically for the $n_b$-collinear one as a function of $k^+$. In other words, the collinear-soft function in the exponential regulator scheme is simply given by the rapidity RG evolution between its canonical rapidity scale $\nu _\mathcal {S}\sim k^-$ and $\nu$, with trivial boundary condition at $\nu _\mathcal {S}$. Inserting this result into Eq. (3.38) leads to Eq. (3.40) in the main text.

D Perturbative ingredients

1.1 D.1 Anomalous dimensions

We expand the QCD $\beta$ function as

$$\begin{aligned} \mu \frac{\mathrm {d}\alpha _s(\mu )}{\mathrm {d}\mu } = \beta [\alpha _s(\mu )] \,, \quad \beta (\alpha _s) = -2\alpha _s\sum _{n = 0}^\infty \beta _n \Bigl ( \frac{\alpha _s}{4\pi } \Bigr )^{n+1} \,. \end{aligned}$$

(D.1)

The coefficients up to three loops in the $\overline{\text {MS}}$ scheme are [176, 177]

$$\begin{aligned} \beta _0= & {} \frac{11}{3}\,C_A - \frac{4}{3}\,T_F\,n_f \,, \nonumber \\ \beta _1= & {} \frac{34}{3}\,C_A^2 - 2T_F\,n_f \Bigl (\frac{10}{3}\, C_A + 2 C_F\Bigr ) \,, \nonumber \\ \beta _2= & {} \frac{2857}{54}\,C_A^3 + 2T_F\,n_f \Bigl (- \frac{1415}{54}\,C_A^2 - \frac{205}{18}\,C_F C_A + C_F^2 \Bigr ) + 4T_F^2\,n_f^2 \Bigl (\frac{79}{54}\, C_A + \frac{11}{9}\, C_F \Bigr ).\nonumber \\ \end{aligned}$$

(D.2)

The cusp anomalous dimension and all noncusp anomalous dimensions are expanded as

$$\begin{aligned} \Gamma _{\mathrm{cusp}}^i(\alpha _s) = \sum _{n = 0}^\infty \Gamma _n^i \Bigl ( \frac{\alpha _s}{4\pi } \Bigr )^{n+1} \,,\quad \gamma (\alpha _s) = \sum _{n = 0}^\infty \gamma _n \Bigl ( \frac{\alpha _s}{4\pi } \Bigr )^{n+1} \,. \end{aligned}$$

(D.3)

The coefficients of the $\overline{\mathrm {MS}}$ cusp anomalous dimension to three loops are [89, 90, 178]

$$\begin{aligned} \Gamma _0^i&= 4 C_i \,,\nonumber \\ \Gamma _1^i&= 4 C_i \biggl [ C_A \Bigl (\frac{67}{9} - 2\zeta _2 \Bigr ) - \frac{20}{9}\, T_F\, n_f \biggr ] \,,\nonumber \\ \Gamma _2^i&= 4 C_i \biggl \{ C_A^2 \Bigl (\frac{245}{6} - \frac{268}{9}\zeta _2 + \frac{22}{3}\zeta _3 + 22\zeta _4 \Bigr ) \nonumber \\&\quad \quad + 2T_F\,n_f \biggl [ C_A \Bigl (-\frac{209}{27} + \frac{40}{9}\zeta _2 - \frac{28}{3}\zeta _3 \Bigr ) + C_F \Bigl (-\frac{55}{6} + 8\zeta _3 \Bigr ) \biggr ] - \frac{16}{27}\, T_F^2\, n_f^2 \biggr \} \,,\end{aligned}$$

(D.4)

where $C_i=C_F$ for $i=q$ and $C_i=C_A$ for $i=g$.

1.2 D.2 Ingredients for $\mathcal {T}_0$

The quark beam function noncusp anomalous dimension coefficients to three loops are [87]

$$\begin{aligned} \gamma _{B\,0}^q&= 6 C_F \,, \nonumber \\ \gamma _{B\,1}^q&= 2C_F \biggl [ C_A \Bigl (\frac{73}{9} - 40 \zeta _3\Bigr ) + C_F \Bigl (\frac{3}{2} - 12\zeta _2 + 24 \zeta _3 \Bigr ) + \beta _0 \Bigl (\frac{121}{18} + 2\zeta _2 \Bigr ) \biggr ] \,, \nonumber \\ \gamma _{B\,2}^q&= 2 C_F \biggl [ C_A^2 \Bigl ( \frac{52019}{162} - \frac{1682}{27}\zeta _2 - \frac{2056}{9}\zeta _3 - \frac{820}{3}\zeta _4 + \frac{176}{3}\zeta _2\zeta _3 + 232 \zeta _5 \Bigr ) \nonumber \\&\quad + C_A C_F \Bigl ( \frac{151}{4} - \frac{410}{3}\zeta _2 + \frac{844}{3}\zeta _3 - \frac{494}{3}\zeta _4 + 16\zeta _2\zeta _3 + 120 \zeta _5 \Bigr ) \nonumber \\&\quad + C_F^2 \Bigl ( \frac{29}{2} + 18\zeta _2 + 68 \zeta _3 + 144\zeta _4 - 32\zeta _2\zeta _3 - 240 \zeta _5 \Bigr ) \nonumber \\&\quad + C_A\beta _0 \Bigl ( -\frac{7739}{54} + \frac{650}{27}\zeta _2 - \frac{1276}{9}\zeta _3 + \frac{617}{3}\zeta _4 \Bigr ) \nonumber \\&\quad + \beta _0^2 \Bigl ( -\frac{3457}{324} + \frac{10}{3}\zeta _2 + \frac{16}{3}\zeta _3 \Bigr ) + \beta _1 \Bigl ( \frac{1166}{27} - \frac{16}{3}\zeta _2 + \frac{52}{9}\zeta _3 - \frac{82}{3}\zeta _4 \Bigr ) \biggr ] \,. \end{aligned}$$

(D.5)

They have been confirmed recently by an explicit three-loop calculation of the jet function [179], see also [180].

The gluon beam function noncusp anomalous dimension coefficients to three loops are [88]

$$\begin{aligned} \gamma _{B\,0}^g&= 2 \beta _0 \,, \nonumber \\ \gamma _{B\,1}^g&= 2C_A \biggl [ C_A \Bigl (\frac{91}{9} - 16\zeta _3 \Bigr ) + \beta _0 \Bigl (\frac{47}{9} - 2\zeta _2 \Bigr ) \biggr ] + 2\beta _1 \,, \nonumber \\ \gamma _{B\,2}^g&= 2C_A \biggl [ C_A^2 \Bigl ( \frac{49373}{162} - \frac{944}{27}\zeta _2 - \frac{2260}{9}\zeta _3 - 144\zeta _4 + \frac{128}{3}\zeta _2\zeta _3 + 112\zeta _5 \Bigr ) \nonumber \\&\quad + C_A\, \beta _0 \Bigl ( - \frac{6173}{54} - \frac{376}{27}\zeta _2 + \frac{140}{9}\zeta _3 + 117\zeta _4 \Bigr ) + \beta _0^2 \Bigl ( -\frac{493}{81} - \frac{10}{3}\zeta _2 + \frac{28}{3}\zeta _3 \Bigr ) \nonumber \\&\quad + \beta _1 \Bigl ( \frac{1765}{54} - 2\zeta _2 - \frac{152}{9}\zeta _3 - 8\zeta _4 \Bigr ) \biggr ] + 2 \beta _2 \,. \end{aligned}$$

(D.6)

The soft noncusp anomalous dimension coefficients to three loops follow from consistency by $\gamma _S^i(\alpha _s) = - 2\gamma _B^i(\alpha _s) -4\gamma _C^i(\alpha _s)$, where the $\gamma _C^i(\alpha _s)$ are taken from [27]. They are the hard noncusp anomalous dimensions and are known up to three loops from the quark and gluon form factors [50,51,52,53, 55,56,57]. We obtain,

$$\begin{aligned} \gamma _{S\,0}^i&= 0 \,, \nonumber \\ \gamma _{S\,1}^i&= 2C_i \biggl [ C_A \Bigl (-\frac{64}{9} + 28\zeta _3 \Bigr ) + \beta _0 \Bigl (-\frac{56}{9} + 2\zeta _2 \Bigr ) \biggr ] \,, \nonumber \\ \gamma _{S\,2}^i&= 2C_i \biggl [ C_A^2 \Bigl ( -\frac{37871}{162} + \frac{620}{27}\zeta _2 + \frac{2548}{9}\zeta _3 + 144\zeta _4 - \frac{176}{3}\zeta _2\zeta _3 - 192\zeta _5 \Bigr ) \nonumber \\&\quad + C_A\, \beta _0 \Bigl ( \frac{4697}{54} + \frac{484}{27}\zeta _2 + \frac{220}{9}\zeta _3 - 112\zeta _4 \Bigr ) + \beta _0^2 \Bigl ( \frac{520}{81} + \frac{10}{3}\zeta _2 - \frac{28}{3}\zeta _3 \Bigr ) \nonumber \\&\quad + \beta _1 \Bigl ( - \frac{1711}{54} + 2\zeta _2 + \frac{152}{9}\zeta _3 + 8\zeta _4 \Bigr ) \biggr ] \,. \end{aligned}$$

(D.7)

Finally, the soft function coefficients to two loops are [49, 100,101,102,103]

$$\begin{aligned} s_i^{(0)}&= 1 \,, \nonumber \\ s_i^{(1)}&= C_i\, 2\zeta _2 \,, \nonumber \\ s_i^{(2)}&= C_i \biggl [ - C_i\, 27 \zeta _4 + C_A \Bigl (-\frac{640}{27} + 8\zeta _2 + 44 \zeta _4 \Bigr ) + \beta _0 \Bigl (-\frac{20}{27} - \frac{37}{3}\zeta _2 + \frac{58}{3} \zeta _3 \Bigr ) \biggr ] \,.\end{aligned}$$

(D.8)

The three-loop coefficient is still unknown.

1.3 D.3 Ingredients for $q_T$

In the exponential regulator, the noncusp anomalous dimension $\tilde{\gamma }_S^i$ of the $q_T$ soft function is equal to that of the threshold soft function $\gamma _\mathrm {thr}^i$, which in turn is the negative of the $\mathcal {T}_0$ soft anomalous dimension $\gamma _S^i$. As a result, we have

$$\begin{aligned} \begin{aligned} \tilde{\gamma }_S^i(\alpha _s)&= \gamma _\mathrm{thr}^i(\alpha _s) = -\gamma _S^i(\alpha _s) \,, \quad&\tilde{\gamma }_{S\,n}^i&= - \gamma _{S\,n}^i \,, \\ \tilde{\gamma }_B^i(\alpha _s)&= \gamma _B^i(\alpha _s) + \gamma _S^i(\alpha _s) \,, \quad&\tilde{\gamma }_{B\,n}^i&= \gamma _{B\,n}^i + \gamma _{S\,n}^i \,. \end{aligned} \end{aligned}$$

(D.9)

The result for $\tilde{\gamma }^i_B$ follows from RG consistency and the fact that the hard anomalous dimension is the same for $q_T$ and $\mathcal {T}_0$. The $\gamma _{S\,n}^i$ and $\gamma _{B\, n}^i$ coefficients are given in Eqs. (D.5), (D.6) and (D.7) above.

The rapidity anomalous dimensions coefficients, which enter the fixed-order expansion of $\tilde{\gamma }_\nu ^i$ in Eq. (3.9), are known up to three loops [82, 136, 141]. They are given by

$$\begin{aligned} \tilde{\gamma }_{\nu \,0}^i&= 0 \,,\nonumber \\ \tilde{\gamma }_{\nu \,1}^i&= 2 C_i \biggl [ C_A \Bigl ( -\frac{64}{9} + 28\zeta _3 \Bigr ) - \frac{56}{9} \beta _0 \biggr ] \,,\nonumber \\ \tilde{\gamma }_{\nu \,2}^i&= 2 C_i \biggl [ C_A^2 \Bigl ( -\frac{37871}{162} + \frac{620}{27} \zeta _2 + \frac{2548}{9} \zeta _3 + 144 \zeta _4- \frac{176}{3} \zeta _2 \zeta _3 - 192 \zeta _5 \Bigr ) \nonumber \\&\quad \quad + C_A \beta _0 \Bigl ( \frac{3865}{54} + \frac{412}{27} \zeta _2 + \frac{220}{9} \zeta _3 - 50 \zeta _4 \Bigr ) \nonumber \\&\quad \quad + \beta _0^2 \Bigl ( - \frac{464}{81} - 8 \zeta _3 \Bigr ) + \beta _1 \Bigl ( -\frac{1711}{54}+ \frac{152}{9} \zeta _3 + 8 \zeta _4 \Bigr ) \biggr ] \,. \end{aligned}$$

(D.10)

The soft function coefficients are also known up to three loops [82, 136], and are given by

$$\begin{aligned} \tilde{s}_i^{(0)}&= 1 \,,\nonumber \\ \tilde{s}_i^{(1)}&= -C_i\, 2\zeta _2 \,,\nonumber \\ \tilde{s}_i^{(2)}&= C_i \biggl [C_i\, 5 \zeta _4 + C_A \Bigl ( \frac{208}{27} - 4 \zeta _2 + 10 \zeta _4 \Bigr ) + \beta _0 \Bigl ( \frac{164}{27} - 5 \zeta _2 - \frac{14}{3} \zeta _3 \Bigr ) \biggr ] \,,\nonumber \\ \tilde{s}_i^{(3)}&= C_i \biggl [ - C_i^2\, \frac{35}{6} \zeta _6 + C_i C_A \Bigl ( - \frac{416}{27}\zeta _2 + 20 \zeta _4 - 35 \zeta _6 \Bigr ) + C_i \beta _0 \Bigl (-\frac{328}{27}\zeta _2 + 25 \zeta _4 + \frac{28}{3} \zeta _2 \zeta _3 \Bigr ) \nonumber \\&\quad + C_A^2 \Bigl ( \frac{115895}{324} - \frac{51071}{486}\zeta _2 - \frac{23396}{81}\zeta _3 - 58 \zeta _4 + 240 \zeta _2 \zeta _3 - 224 \zeta _5 + \frac{928}{9}\zeta _3^2 - \frac{3086}{27}\zeta _6 \Bigr ) \nonumber \\&\quad + C_A \beta _0 \Bigl ( -\frac{363851}{2916} + \frac{2987}{486}\zeta _2 - \frac{428}{81}\zeta _3 + \frac{830}{9}\zeta _4 - \frac{220}{3} \zeta _2 \zeta _3 + \frac{1388}{9}\zeta _5 \Bigr ) \nonumber \\&\quad + \beta _0^2 \Bigl ( - \frac{64}{729} - \frac{34}{3}\zeta _2 - \frac{140}{27} \zeta _3 - \frac{11}{3}\zeta _4 \Bigr ) \nonumber \\&\quad + \beta _1 \Bigl ( \frac{42727}{972} - \frac{275}{18}\zeta _2 - \frac{1744}{81} \zeta _3 - \frac{76}{9}\zeta _4 + \frac{40}{3} \zeta _2 \zeta _3 - \frac{112}{9}\zeta _5 \Bigr ) \biggr ] \,. \end{aligned}$$

(D.11)

1.4 D.4 Mellin kernels and splitting functions

We decompose the flavor dependence of a generic Mellin-convolution kernel $K_{ij}(z)$ as

$$\begin{aligned} K_{q_i q_j}(z) = K_{\bar{q}_i \bar{q}_j}(z)&= \delta _{ij} K_{qqV}(z) + K_{qqS}(z) + K_{qq\Delta S}(z) \,,\nonumber \\ K_{q_i \bar{q}_j}(z) = K_{\bar{q}_i q_j}(z)&= \delta _{ij} K_{q\bar{q}V}(z) + K_{qqS}(z) - K_{qq\Delta S}(z) \,,\nonumber \\ K_{q_i g}(z) = K_{\bar{q}_i g}(z)&= K_{q g}(z) \,, \nonumber \\ K_{gg}(z)&= K_{gg}(z) \,,\nonumber \\ K_{g q_i}(z) = K_{g \bar{q}_i}(z)&= K_{g q}(z) \,. \end{aligned}$$

(D.12)

This decomposition is sufficient and unique to all orders by the flavor and charge symmetries of QCD. The $K_{qqV}$ and $K_{gg}$ contributions are already present at tree level, the $K_{qg}$ and $K_{gq}$ channels start at one loop, the $K_{qqS}$ and $K_{q\bar{q}V}$ channels open up at two loops, and the $K_{qq\Delta S}$ channel only receives contributions from topologies at three loops and beyond. This decomposition also makes it straightforward to evaluate and iterate sums over intermediate partons. For example, for the convolution of two generic kernels $(K K')_{ij}(z)$, we have

$$\begin{aligned} (K K')_{gg}(z)&= (K_{gg} K'_{gg})(z) + 2n_f (K_{gq} K'_{qg})(z) \,, \nonumber \\ (K K')_{qg}(z)&= \bigl [ \bigl ( K_{qqV} + K_{q\bar{q}V} + 2n_f K_{qqS} \bigr ) K'_{qg} \bigr ](z) + (K_{qg} \, K'_{gg})(z) \,, \nonumber \\ (K K')_{gq}(z)&= \bigl [ K_{gq} \bigl ( K'_{qqV} + K'_{q\bar{q}V} + 2n_f K'_{qqS} \bigr ) \bigr ](z) + (K_{gg} \, K'_{gq})(z) \,, \nonumber \\ (K K')_{qqV}(z)&= (K_{qqV} K'_{qqV})(z) + (K_{q\bar{q}V} K'_{q\bar{q}V})(z) \,, \nonumber \\ (K K')_{q\bar{q}V}(z)&= (K_{qqV} K'_{q\bar{q}V})(z) + (K_{q\bar{q}V} K'_{qqV})(z) \,, \nonumber \\ (K K')_{qqS}(z)&= \bigl [ K_{qqS} \bigl ( K'_{qqV} + K'_{q\bar{q}V}\bigr ) \bigr ](z) + \bigl [ \bigl ( K_{qqV} + K_{q\bar{q}V}\bigr ) K'_{qqS} \bigr ](z) \nonumber \\&\quad + 2 n_f \bigl (K_{qqS} \, K'_{qqS}\bigr )(z) + \bigl (K_{qg} \, K'_{gq}\bigr )(z) \,, \nonumber \\ (K K')_{qq\Delta S}(z)&= \bigl [ K_{qq\Delta S} (K'_{qqV} - K'_{q\bar{q}V})\bigr ](z) +\bigl [ (K_{qqV} - K_{q\bar{q}V}) K'_{qq\Delta S} \bigr ](z) \nonumber \\&\quad + 2 n_f \bigl ( K_{qq\Delta S} K'_{qq\Delta S}\bigr )(z) \,, \end{aligned}$$

(D.13)

where $n_f$ is the number of active flavors, and the outer brackets on the right-hand side indicate Mellin convolutions without flavor sums.

The DGLAP splitting functions are defined as the anomalous dimension of the PDFs,

$$\begin{aligned} \mu \frac{\mathrm {d}}{\mathrm {d}\mu } f_i(x,\mu ) = 2 \sum _j \int \frac{\mathrm {d}z}{z}\, P_{ij}(z,\mu )\, f_j\Bigl (\frac{x}{z}, \mu \Bigr ) \,. \end{aligned}$$

(D.14)

We perturbatively expand them in powers of $\alpha _s/4\pi$, see Eq. (1.5), and decompose their flavor dependence as in Eq. (D.12). The DGLAP kernels have been calculated at three loops in [89, 90]. Denoting the results of [89, 90] by a calligraphic $\mathcal {P}$ to distinguish them from our $P_{ij}^{(n)}$, we can relate the two notations by

$$\begin{aligned} \begin{aligned} P_{qqV}^{(n)}(z)&= \frac{1}{2} \Bigl [ \mathcal {P}_{ns}^{(n)+}(z) + \mathcal {P}_{ns}^{(n)-}(z) \Bigr ] \,, \quad&P_{gg}^{(n)}(z)&= \mathcal {P}_{gg}^{(n)}(z) \,,\\ P_{q\bar{q}V}^{(n)}(z)&= \frac{1}{2} \Bigl [ \mathcal {P}_{ns}^{(n)+}(z) - \mathcal {P}_{ns}^{(n)-}(z) \Bigr ] \,, \quad&P_{g q}^{(n)}(z)&= \mathcal {P}_{gq}^{(n)}(z) \,,\\ P_{qqS}^{(n)}(z)&= \frac{1}{2 n_f} \mathcal {P}_{ps}^{(n)}(z) \,, \quad&P_{q g}^{(n)}(z)&= \frac{1}{2 n_f} \mathcal {P}_{qg}^{(n)}(z) \,,\\ P_{qq \Delta S}^{(n)}(z)&= \frac{1}{2 n_f} \mathcal {P}_{ns}^{(n)s}(z) \,. \end{aligned} \end{aligned}$$

(D.15)