A toolbox for qT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{T}$$\end{document} and 0-jettiness subtractions at N3LO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {N}^3\hbox {LO}$$\end{document}

We derive the leading-power singular terms at three loops for both qT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_T$$\end{document} and 0-jettiness, T0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {T}_0$$\end{document}, for generic color-singlet processes. Our results provide the complete set of differential subtraction terms for qT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_T$$\end{document} and T0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {T}_0$$\end{document} subtractions at N3LO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {N}^3\hbox {LO}$$\end{document}, which are an important ingredient for matching N3LO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {N}^3\hbox {LO}$$\end{document} 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 qT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_T$$\end{document} and T0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {T}_0$$\end{document} at N3LL′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {N}^3\hbox {LL}'$$\end{document} and N4LL\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {N}^4\hbox {LL}$$\end{document} 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 qT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_T$$\end{document} and T0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {T}_0$$\end{document} beam function coefficients at N3LO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {N}^3\hbox {LO}$$\end{document} in the z→1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z\rightarrow 1$$\end{document} threshold limit, and we also estimate the size of the unknown terms beyond threshold.


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 → 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 α s and/or the resummation of the dominant higher-order terms to all orders in α 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 N 3 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 N 3 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 τ , whose differential cross section can be factorized in the singular limit τ → 0. In this paper, we consider two such variables, 0-jettiness T 0 and the total color-singlet transverse momentum q T , and derive their singular structure to N 3 LO. Our results are necessary ingredients for carrying out the resummation for q T and T 0 at N 3 LL and N 4 LL order. For the associated q T and T 0 subtraction methods [35,38,39], we provide the complete differential N 3 LO subtraction terms in analytic form. The structure and required ingredients for the T 0 subtractions at N 3 LO were already discussed in [39]. The q T slicing method at N 3 LO was also considered in [18]. Differential 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 N 3 LO + 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 (1.1) We always use the indices a and b to label the initial states, and κ i ∈ {g, u,ū, d,d, s, . . .} denotes the parton type and flavor. When there is no ambiguity, we simply identify κ i ≡ i, e.g., we just write ab → L. The q μ a,b are lightlike Born reference (label) momenta given by (1.2) whereẑ is the beam direction and E cm is the hadronic center-of-mass energy. The precise definition of the Born momentum fractions x a,b and the associated ω a,b = x a,b E 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 = q 2 and rapidity Y of the color singlet. Other choices are possible as well, e.g., ω a,b = q ∓ ≡ Q 2 + q 2 T e ±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 τ factorizes for τ → 0 as [48,49] dσ (1. 3) The convolution structure denoted by ⊗ depends on the precise definition of τ . The key properties of Eq. (1.3) are that it captures all QCD singularities in τ and that it factorizes the dependence on the underlying process from the dependence on τ . The process dependence is carried by the hard function H ab , which describes the Born process ab → L, with the sum running over all relevant parton channels. At lowest order, H (0) ab is equivalent to the partonic Born cross section for ab → 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 → H in the m t → ∞ limit, bb → 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 τ 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 ≡ κ i , while S c only depends on the color channel of the Born process. In our case, the only possible color channels are c = {qq, gg}, which are equivalent to the color representation of the incoming partons, so we simply label it by c ≡ i = {q, g}.
The beam and soft functions do depend on the definition of τ , 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 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 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 → 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 → 1 limit for both 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 → 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 N 3 LO but fails starting at N 4 LO. Since our results capture the complete singular structure for z → 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 T 0 beam function was only known at NNLO [85][86][87][88], while only partial results were known at N 3 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 → 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 → 1 limit again provided important checks.
For q T , a similar study of the logarithmic structure at N 3 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 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 T 0 subtractions for qq and gg processes at N 3 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 T 0 in Sect. 2 and for q T in Sect. 3. The application to T 0 and q T subtractions at N 3 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.

Notation and conventions
Throughout the paper, we denote the perturbative expansion of any function F(μ) as (1.4) All anomalous dimensions γ i x (α s ) and the QCD splitting functions are expanded as (1.5) We use the following notation to abbreviate Mellin convolutions and flavor sums where i, j, k ∈ {g, u,ū, d,d, s, . . .} label parton type and flavor. We also define a corresponding identity operator as (1.7) For Fourier-type convolutions, we use the notation Here, the corresponding identity elements are simply δ(k) or δ(t).
We denote logarithmic plus distributions as For dimensionful arguments, we define (1.10) More details are given in Appendix A.

Factorization
The factorization for N -jettiness, T N , has been derived using SCET in [49,69,95]. Here, we focus on 0-jettiness, 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] 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 leptonic: For our present purposes, the precise choice of the Q i is not important, so we will simply use the symbol T 0 . The factorization for T 0 is given by [49] dσ Explicit definitions of the beam and soft functions for 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 (SCET I ) beam function. It appears in the N -jettiness factorization for any N [95], including deepinelastic 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 N 3 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 − → 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 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].

T 0 soft function
The beam thrust soft function satisfies the all-order RGE [49,87] where i cusp (α s ) and γ i S (α s ) are the cusp and soft noncusp anomalous dimensions. The RGE in Eq. (2.4) fully predicts the structure of S i (k, μ) in k and μ to all orders in perturbation theory. By solving it recursively order by order in α s , we can derive this structure at any given fixed order. Expanding both sides of Eq. (2.4) to fixed order in α s (μ) and accounting for the running of α s (μ), we obtain a relation for the (n + 1)-loop term in terms of the terms up to n loops, where we used the short-hand notation in Eq. (1.8) for the convolution in k. This can be integrated to give where the soft function boundary coefficients are defined by Here, we have used distributional scale setting μ 0 = k| + [106], which is defined such that it effectively allows us to treat the μ dependence of the logarithmic distributions like ordinary logarithms. In particular, it satisfies [106] L n (k, μ 0 = k| + ) = 0 , The first relation is used in Eq. (2.7) to define the boundary coefficients as the coefficients of the δ(k) by setting all logarithmic distributions in S (n) i (k, μ) to zero. The other two relations allow us to easily perform the μ integral in Eq. (2.6), essentially turning a δ(k) into a L 0 (k) and a L n (k) into a L n+1 (k). In addition, to evaluate the cross terms for m ≥ 1 in Eq. (2.6), we need the convolutions [107] (L 0 L 0 )(k, μ) = 2L 1 (k, μ) − ζ 2 δ(k) , (2.9) Starting from the LO result, s (0) i = 1, Eq. (2.6) yields up to two loops which agrees with [39]. Evaluating Eq. (2.6) at the next order, we obtain the three-loop result, with the coefficients of the logarithmic distributions given by 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, μ) has an explicit dependence on μ, which cancels against that of the hard and beam functions in Eq. (2.3). Therefore, simply varying the scale where the evolution factor U i S (k, μ S , μ) encodes the solution of Eq. (2.4), with U i S (k, μ S , μ S ) = δ(k). It can be found, e.g., in [87,88]. Formally, the μ S dependence on the right-hand side cancels, but when the starting condition S i (k, μ 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 (2.14) for which the distributions turn into ordinary logarithms of T cut . In Fig. 1, we take as an example T cut = μ = 20 GeV and show the residual μ S dependence of the resummed soft function when varying μ S around the canonical central value μ S = T cut at NLL (dotted green), NNLL (dashed blue), and N 3 LL (solid orange). For the latter, we set the currently unknown three-loop constant term s 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 α s , which formally amounts to a higher-order effect at (N)NLL . For the quark soft function (left panel), the μ S dependence is more than halved going from NLL to NNLL , and roughly halved again at N 3 LL . In the gluon case (right panel), the μ 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 μ S dependence due to its α 3 s (μ S ) prefactor, which, however, should not change the general picture. We stress that the residual μ 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.

T 0 beam function
The beam function B i (t, x, μ) obeys the all-order RGE [49,87] where i cusp (α s ) and γ i B (α s ) are the cusp and beam noncusp anomalous dimensions. For t QCD , the beam function satisfies an OPE in terms of standard PDFs [49,87] (2. 16) where the I i j (t, z, μ) are perturbatively calculable matching coefficients. Taking into account the evolution of the PDFs, they obey the RGE [87] μ d dμ (2.17) where 1 i j (z) ≡ δ i j δ(1 − z) and 2P i j (z, μ) are the PDF anomalous dimensions.
By solving the RGE in Eq. (2.17) recursively order by order, we can derive the complete structure of I i j (t, z, μ) 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 where the μ-independent boundary coefficients are defined as Starting from the LO result, , we obtain up to two loops i j (z) i j (z) 20) which agrees with [39,85,88]. The NLO and NNLO boundary coefficients I with the required Mellin convolutions (P (0) P (0) ) i j (z) and (I (1) P (0) ) i j (z) can be found in [85,86] with the coefficients We caution that the functions P i j (z) and I i j (z) in [39,85,86] are expanded in powers of α s /(2π), while here we expand them in powers of α s /(4π). (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) ) i j (z), this required employing the identity (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 The explicit expression for the beam function evolution kernel U i B (t, μ B , μ) can be found in [87,88].
In Fig. 2, we show the residual μ B dependence of the resummed integrated beam function at fixed representative values of x = 10 −2 and √ t cut = μ = 30 GeV. We again show the relative difference to the NNLL central result at μ B = √ t cut at NLL (dotted green), NNLL (dashed blue), and N 3 LL with the unknown three-loop I (3) i j (z) = 0 (solid orange). We use the MMHT2014nnlo68cl [109] NNLO PDFs and four-loop running of α s everywhere. These evolution orders are sufficient to ensure the formal cancellation of the μ B dependence at N 3 LL , while at lower orders, they amount to a higher-order effect. Numerical results to N 3 LL with PDFs and α s evolution at corresponding lower orders can be found in [86]. The residual dependence on μ B is noticeably reduced by about a factor of two at N 3 LL compared to NNLL . The missing three-loop constant terms will again add an additional source of μ B dependence due to its α 3 s (μ B ) prefactor and also the scale dependence of the PDFs, which, however, should not change the qualitative picture.

Beam function coefficients in the eikonal limit
We now obtain the beam function coefficients I (n) i j (z) in the z → 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 (n) i j (z → 1) to the threshold soft function to all orders in α s . Consistency relations of this kind generically arise from different factorization theorems that apply in different limits of the same multidifferential 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 ∓ and corresponding momentum fractions x ∓ , As recently shown in [97], in the generalized threshold limit x − → 1 but generic x + , the inclusive color-singlet cross section differential in q ± factorizes as (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 thr a (x) encodes the extraction of parton a from the proton for x → 1.
On the other hand, in the well-known and stronger soft threshold limit, where both x − → 1 and x + → 1, the cross section factorizes as [110][111][112][113][114] dσ The new ingredient is the threshold soft function S thr i (k − , k + , μ). It describes the processindependent contribution of soft emissions with total lightcone momenta k + = n · k and k − =n · k. It also only depends on the color representation c ≡ 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 (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 dσ/dQ 2 , which is known to O(α 3 s ) [11,12]. In Appendix B, we review this procedure and give explicit results for S thr i (k − , k + , μ) 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 + → 1, it must reproduce Eq. (2.27). As a result, the eikonal x → 1 limit of the beam function must coincide with the threshold soft function [97], which is justified at leading power in 1 − z, yields the corresponding relation for the matching coefficients [97], This relation captures all terms in I i j (t, z, μ) that are singular for z → 1, while power corrections have at most an integrable singularity for z → 1. Notably, the beam function becomes flavor diagonal as z → 1, while off-diagonal channels are 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 → 1, they become invariant under a simultaneous rescaling t → e +y t and 1 − z → e −y (1 − z). In other words, they are symmetric in t/ω and ω(1 − z) such that the dependence on ω 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 δ(k − ) in the threshold soft function upon identifying Including the one-loop and two-loop results for reference, we find (2.31) The boundary coefficients s . 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 → 1 terms must be symmetric in t/ω and ω(1 − z) allows one to directly determine (or check) the δ(t)L n (1 − z) terms from the RGE-predicted L n (t)δ(1 − z) terms, which was already noted in [85,118]. However, the δ(t)δ(1−z) coefficient cannot be predicted in this way, and Eq. (2.30) explicitly identifies it with the threshold soft function coefficients s 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 functionB i (t, x, μ). 2 Note that in contrast to Eqs. (1.3) and (2.3), here the difference between q ± and (Q, Y ) matters. The RGE forB i (t, x, μ) is the same as for B i (t, x, μ) in Eq. (2.15), and hence, Eq. (2.22) also holds forB i just with different boundary coefficients I (n) i j (z). In the limit z → 1, the difference between B i andB i becomes power suppressed in 1 − z. As a result, the z → 1 limit of the modifiedĨ (n) i is also given by Eq. (2.31).

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 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 ( 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 α s /(4π) at each order, so the overall normalization shows the percent impact relative to the LO result. For definiteness, we choose μ = 30 GeV for the scale entering the PDFs and α s . The eikonal approximation reproduces the correct divergent behavior of the full flavordiagonal 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 → 0 divergent behavior in the gluon coefficient, which is not reproduced by its z → 1 expansion. If desired, it can be captured by including the leading z → 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 flavordiagonal channel and also much flatter in x, since they only start at NLP.
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, where the NLP coefficient itself is approximated as 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 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 NLP 1,1 = −4C i and c NLP 2,3 = −8C 2 i for both quarks and gluons. Here, the additional double logarithm is determined by the same power of ( i 0 ) n as at LP, and this pattern can be expected to hold at higher orders. The second term in Eq. (2.33) generates a next-to-leading logarithmic NLP series. We fix the central value for X 1 = 1 to reproduce the NLL constant term at NLO c 1,0 = 4C i . Interestingly, we find that this choice also reproduces all NNLO coefficients c 2,k very well, typically to within 10%, for both quarks and gluons and also independently of the choice of n f . This provides a very nontrivial check and so we expect that Eq. (2.33) provides a very good model of the true NLP structure also at higher orders. To estimate the uncertainties, we vary X 1 by ±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 ±1.
Since X i probe independent structures, we can consider them as uncorrelated. Hence, we add the impacts i on the final result of their variation in quadrature (2. 35) In Fig. 4, we show the approximate kernel at NNLO (top) and N 3 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 N 3 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 N 3 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.

Factorization theorem
The factorization of the q T distribution in the limit q T ≡ | q T | 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 It receives power corrections suppressed by q 2 T /Q 2 as indicated. As is common, we consider the factorized singular cross section in Fourier-conjugate b T space, where convolutions in q T space turn into simple products. In particular, Fourier transforming the q T -dependent plus distributions L n ( q T , μ) turns them into powers of the canonical b T -space logarithms, which we denote as 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 ν in Eq. (3.1). We use the exponential regulator of [117], which up to two loops gives results equivalent to the η regulator of [127,128].
The beam function appearing in Eq. (3.1) is the inclusive transverse-momentum dependent (SCET II ) beam function, which also appears in the q T factorization of Z + j and γ + 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].

Rapidity anomalous dimension
The ν dependence of the beam and soft functions is encoded in their rapidity RGEs [127], whereγ i ν,B andγ i ν,S 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 ν, they are related bỹ and we will simply refer toγ i ν (b T , μ) as the rapidity anomalous dimension. An important property ofγ i ν (b T , μ) 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 −1 T 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 μ, which predicts its all-order structure in b T and μ. Similar to the T 0 soft function in Sect. 2.2, it can be solved recursively order by order in α s . Expanding both sides of Eq. (3.6) to fixed order in α s (μ) and accounting for the running of α s (μ), the (n + 1)-loop term is related to the lower-order terms bỹ (3.7) where the nonlogarithmic boundary coefficients are defined as The result up to three loops is (3.9) The boundary coefficientsγ i ν n are known up to three loops [82,136,141] and are summarized in Eq. (D.10).

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 μ and ν RGEs, The perturbative structure of γ i ν is discussed in Sect. 3.2. The μ anomalous dimension has the all-order structureγ where i cusp (α s ) andγ i S (α s ) are the cusp and the soft noncusp anomalous dimensions. Expanding both sides of Eq. (3.10) order by order in α s , we obtain the coupled RGEs (3.12) These are easily integrated to givẽ where we first integrated the ν RGE at fixed μ = b 0 /b T and then the μ RGE at arbitrary ν.
In this way, the rapidity anomalous dimension reduces to its boundary coefficients γ i ν n . The soft boundary coefficients in Eq. (3.13) are defined as (3.14) Starting from the LO result,s i = 1, and expressing the results in terms of Eq. (3.13) yields up to two loops At three loops, we write the result as where theS (3) i,k coefficients themselves are polynomials in L ν . Insertingγ i S 0 =γ i ν 0 = 0 for brevity, they are given bỹ Numerical impact The soft functionS i (b T , μ, ν) has an explicit dependence on the scales μ and ν that cancels against that of the hard and beam functions in Eq. (3.1). Therefore, varying μ and ν is not very meaningful for illustrating the numerical impact of the scale-dependent three-loop terms. Instead, we consider the resummed soft function, where we have chosen to first evolve in ν and then in μ.
To probe the full set of terms in the fixed-order expansion ofS i (b T , μ S , ν S ), we consider simultaneous variations of (μ S , ν S ) around the canonical central scales μ S = ν S = μ = ν = 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 GeV at NLL (dotted green), NNLL (dashed blue), and N 3 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 N 3 LL results do not coincide at the central scales. We use four-loop running of α 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 N 3 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.

Beam function
The beam function obeys the coupled RGEs where the ν anomalous dimension is discussed in Sect. 3.2, and the μ anomalous dimension has the all-order formγ where i cusp (α s ) andγ i B (α s ) are the cusp and the beam noncusp anomalous dimensions. For perturbative b 0 /b T QCD , the TMD beam function satisfies an OPE in terms of standard PDFs [48], (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 piecẽ I g j and the polarization-dependent pieceJ g j , where g ρλ ⊥ = g ρλ − (n ρnλ +n ρ n λ )/2 is the transverse metric and b ρ ⊥ is the transverse four vector with b 2 ⊥ = − b 2 T . Due to the multiplicative structure of Eq. (3.20), bothĨ g j andJ g j obey the same RGE, and in the following, we will only consider the RGEs forĨ i j .
TheĨ i j are perturbatively calculable matching coefficients, whose RGEs follow from Eq. (3.20) by taking the evolution of the PDFs into account, (3.23) Similar to the soft function, these coupled RGEs can be solved recursively as where the nonlogarithmic boundary coefficients are defined as Starting from the LO result,Ĩ , we obtain up to two loops where we abbreviated (3.27) Note that L ω differs from the characteristic logarithm of the soft function in the previous section. TheĨ (n) i j (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 writẽ and usingγ i ν 0 = 0 for brevity, the coefficients arẽ where the three-loop boundary coefficientsĨ (3) i j (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 kernelsJ g j have a simpler structure than theĨ i j 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 N 3 LO cross section. They are given bỹ The two-loop termsJ (2) g j have recently been calculated in [142] using the exponential regulator and in [143] using the δ regulator. They can be converted to our convention via the relatioñ (3. 31) In the first line of Eq. (3.31), I gi (z) is the two-loop boundary term as given in [142]. In the second line of Eq. (3.31),s (1) g = −2C A ζ 2 is the soft function constant at one loop and δ L C (2;,0,0) g← j (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 functioñ (3.32) For i = g, we restrict to the polarization-independent pieceĨ g j and writeB g ≡ −g ⊥,ρλB ρλ g for short. As for the soft function, we restrict to simultaneous variations of μ B and ν B . In Fig. 6, we show the residual dependence on (μ B , ν B ) at NLL (dotted green), NNLL (dashed blue), and N 3 LL with the unknown I  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Ĩ

Beam function coefficients in the eikonal limit
We now proceed to extract the three-loop beam function coefficients in the z → 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 dσ pp→L /dQ 2 dY dq 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 → 1 and x b → 1. As x a,b → 1, initial state radiation is constrained to have energy λ − λ + Q, where 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 λ − λ + 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 λ − λ + ∼ 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 O(q 2 T /Q 2 ), but captures the exact dependence on x a,b via the beam functions. q T /Q λ − λ + 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 1. The leftover radiation in this limit is described by intermediate collinear-soft modes [74,149] in terms of n a,b -collinear-soft functionsS i (k, b T , μ, ν). They are matrix elements of collinearsoft Wilson lines and depend on the small additional momentum k = k ∓ 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] dσ (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.
In this regime, the threshold constraint dominates and all radiation is forced to be soft. The recoil against soft radiation with transverse momentum k T = − q T is encoded in the fully differential threshold soft function S thr i (k − , k + , k T ). In terms of its Fourier transform with respect to k T ,S thr i (k − , k + , b T ), the factorization theorem reads 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 O(α 2 s ) in [150], albeit in a different context, and to O(α 3 s ) in [82]. By construction, it satisfies where S thr i (k − , k + , μ) 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 → 1 limit of the q T beam function is captured by the collinear-soft function [117,148], 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] 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] S thr (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 functionS i (b T , μ, ν) (which determines the scheme). Furthermore, the scheme ambiguity amounts to moving terms from the soft function boundary coefficients into the coefficient of δ(1 − z) in the beam function coefficients. Since δ(1 − z) is a leading-power contribution as z → 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Ĩ i j in the exponential regulator scheme, where the plus distribution 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 δ-regulator in [137]. It is straightforward to expand Eq. (3.40) in α 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 havẽ in agreement with the full two-loop result [136], and where we have used thatγ i ν 0 = 0. Including terms up to six loops for illustration, we find 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 δ(1 − z) in the beam function matching coefficient due to the inverse Fourier transform to k ± back from the conjugate b ± space, where the regularization procedure is applied.

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 i j ⊗ 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 α s /(4π) 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 μ = 30 GeV.
In both flavor-diagonal contributions, denoted as qqV and gg, the eikonal limit correctly reproduces the divergent behavior as x → 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 (z), which is not captured by the expansion around z → 1. If desired, one could also include the leading z → 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.
Similar to the 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 L n (1 − z) with n > 0, but only 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, Eur. Phys. J. Plus (2021) 136:214 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 comparisoñ Here,Ĩ i j,reg refers to the full regular (non-eikonal) piece of the beam coefficient at O(α n s ). 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 usẽ 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 γ i X needed to achieve this is given by γ (3.45) The size of these additional logarithms can be probed by varying the coefficients X 1,2 by ±1 around the central choice. Furthermore, we add the eikonal limitĨ (n) LP i j suppressed by one power of (1 − z) to estimate the pure NLP constant. Its coefficient X 3 is varied by ±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 i on the final result of their variation in quadrature (3.46) In Fig. 8, we show the approximate kernel at NNLO (top) and N 3 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 N 3 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 L 0 (1 − z) with a rather small coefficient γ i ν 2 , while the NLP now contains up to ln 4 (1 − z). The fact that the uncertainty bands are of similar size at NNLO and N 3 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 N 3 LO for Higgs production. There, the employed approximation wasĨ N 3 δ(1−z), withC N 3 fixed such that the inclusive cross section is correctly reproduced. This effectively absorbs the averaged effect of the actual z dependence into an effective δ(1 − z) coefficient. From our results, we know the exact δ(1 − z) coefficient, and so our approximate results give an independent estimate of the actual rapidity dependence and total size of these unknown terms.

N 3 LO subtractions
The factorization in Eq. (1.3) fully describes the limit τ → 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 τ 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 τ or as a global τ 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 τ , where the second relation defines the cumulant in τ 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 τ → 0, the cross sections scales like ∼ 1/τ , so performing the τ integral requires knowing the full analytic distributional structure involving δ(τ ) and L n (τ ), which encodes the cancellation of real and virtual IR divergences. To separate out the singular structure in τ , we introduce a subtraction term, where dσ sub (X )/dτ captures all singularities for τ → 0, and σ sub (X, τ off ) is the integrated subtraction term, By construction, the integrand in square brackets in Eq. (4.2) contains at most integrable singularities for τ → 0 and so the integral can be performed numerically. Hence, the full cross section dσ (X )/dτ is only ever evaluated at finite τ > 0, which means it can be obtained from a calculation of the corresponding ab → L + 1-parton process at one lower order. In practice, one always has a small IR cutoff δ on the τ integral, where the last term contains the integral over τ ≤ δ, which is neglected for δ → 0. The above is a differential τ -subtraction scheme, where the parameter τ off ∼ 1 determines the range over which the subtraction acts. The key advantage of formulating the subtractions in terms of a physical resolution variable τ , 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 τ → 0, which is also the basis for the resummation in τ . 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, τ off corresponds to the point where the τ 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 T 0 subtractions are used in this way in the Geneva Monte Carlo to combine the fully differential NNLO calculation together with the NNLL T 0 resummation with a parton shower [46,47,157]. The differential subtractions at N 3 LO are a key ingredient for using this method to combine N 3 LO calculations with parton showers.
In contrast to a fully local subtraction scheme, all singularities are projected onto the resolution variable τ , so the subtractions are local in τ 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 T 0 resummation [74,75] offers the possibility to construct doubledifferential q T − T 0 subtractions.
The key point of the differential subtraction is that δ can in principle be made arbitrarily small, because the integrand of the τ integral is nonsingular, which also means that the numerically expensive small τ region does not need to be sampled with weight 1/τ . On the other hand, by letting δ = τ cut be a small but finite cutoff and setting τ off = τ cut , Eq. (4.5) turns into a global τ subtraction or slicing, (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 τ . This makes the L+1-jet calculation very demanding, both in terms of computational expense and numerical stability, because the 1/τ -divergent integral of dσ (X )/dτ must be computed with sufficient accuracy down to sufficiently small τ cut , which in practice limits how small one can take τ cut . Since the integral is divergent, one cannot let τ cut → 0 even in principle, so one always has a leftover systematic uncertainty from the neglected power corrections σ (X, τ cut ).
The numerical efficiency of the subtractions can be improved by including the power corrections in the subtractions for both 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 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 O(τ ) to O( √ τ ) [173].

Subtraction terms
The singular terms needed for the subtractions only depend on the Born phase space, so we can write them as dσ sing (X ) where 0 ≡ 0 (κ a , κ b , ω a , ω b ) denotes the full Born phase space, including the parton labels κ a,b , the total color-singlet momentum q μ parametrized in terms of ω a,b as in Eqs.
(1.1) and (1.2) as well as the internal phase space of L. The X ( 0 ) denotes the measurement function that implements the measurement X on a Born configuration. The singular terms are defined such that their τ dependence is minimal and given by (4.9) Their integral over τ ≤ τ cut immediately follows as (4.10) The differential subtractions are given by using Eq. (4.9) for τ > 0, which amounts to dropping the C −1 ( 0 )δ(τ ) term and using L n (τ > 0) = ln n−1 (τ )/τ . The integrated subtractions are directly given by Eq. (4.10). The precise definition of the C n ( 0 ) coefficients depends on the normalization of the dimensionless variable τ or equivalently on the boundary condition of the L n (τ ). Rescaling τ → λτ moves contributions from C n ( 0 ) to C m<n ( 0 ). This freedom was used in [39] to absorb all terms with n ≥ 0 in Eq. (4.10) into a C −1 ( 0 , T off ) by taking τ ≡ T 0 /T off . Here, we prefer to keep the cutoff dependence explicit as in Eq. (4.10) and take The m-loop subtraction coefficients C (m) n ( 0 ) directly follow from expanding Eq. (2.3) for T 0 or Eq. (3.1) for q T to mth order in α s . For the three-loop coefficients, this yields where for simplicity we have suppressed the dependence on μ and the distinction of the and thus are needed for the integrated subtraction terms but not the differential ones. Note also that most of the process and 0 dependence resides in the hard coefficients, while the soft/collinear contributions only depend on x a,b and the parton types, The results for the subtraction coefficients C n ( 0 ) in Eq. (4.12) up to three loops for both 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 T 0 requires rescaling and convolving the plus distributions in the beam and soft functions, as discussed in [39]. For q T , expanding the b T -space result dσ sing ( b T ) yields powers of the b T -space logarithm L n b up to n ≤ 6. Their Fourier transform, given in Table 1 in Appendix A.2, yields simple δ( q T ) and L n ( q T , μ), which are easily rescaled to δ(τ ) and L n (τ ).
Note also that the original q T subtraction method in [35] was based on the q T resummation framework of [175], where the canonical b T -space logarithms are replaced by This form is also used, e.g., in [18,160]. While usingL b has certain advantages in the context of q T resummation, it is unnecessary for the purpose of q T subtractions, since L b andL b yield the same singular terms and only differ by power corrections. A drawback of using L b here is that the Fourier transform ofL n b yields complicated expressions in q T space, see Appendix B in [175], whose cumulants are not known analytically and must be performed numerically.

Conclusions
We have studied the three-loop structure of beam and soft functions for both 0-jettiness 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 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 i j (z → 1), i.e., the full singular limit of the beam functions as z → 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 T 0 and q T at N 3 LL and N 4 LL order. In particular, they are important for extending the q T and T 0 subtraction methods to N 3 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 fixedorder calculations to parton showers to N 3 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 T 0 beam functions become available.
Note added: As discussed in the introduction, since this paper first appeared, the full threeloop integrated subtraction terms have become available. Specifically, results for the threeloop 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.

A Plus distributions and Fourier transforms
Here, we summarize the definitions and relations for plus distributions.
A.1 One-dimensional plus distributions Following [107], we denote plus distributions as For distributions with dimensionful arguments, we define Using L a (x), we further define the distribution which satisfies the group property The μ dependence of V a (k, μ) is given by Expanding V a (k, μ) in powers of a, we find The Fourier transformation of V a (k, μ) is given by dk e −iky V a (k, μ) = e −aL y , dy 2π e iky e −aL y = V a (k, μ) , L y = ln(iyμe γ E ) .
(A.8) Table 1 Fourier transform of L n b = ln n (b 2 T μ 2 /b 2 0 ) to q T space for n ≤ 6, as given by Eq. (A.12) A.2 Two-dimensional plus distributions for q T Following [106], we define two-dimensional plus distributions in q T as where L n (x) is defined as above in Eq. (A.1), such that The cumulant for a generic cut | q T | ≤ q cut T follows to be The Fourier transformation of L n ( q T , μ) and its inverse are [106] where L b is the usual logarithm in Fourier space and the coefficients R Up to N 3 LO, we require the Fourier transforms of L n b with n ≤ 6, which are summarized in Table 1.

B Threshold soft function
Here, we discuss the double-differential threshold soft function S thr i (k − , k + , μ), which appears in the soft threshold factorization for the inclusive cross section in Eq. (2.27) and determines the eikonal limit of the T 0 beam function in Eq. (2.31). We give its complete N 3 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.

B.1 Plus distribution basis
A key property of the threshold soft function is that is invariant under the simultaneous rescaling k − → k − e +y and k + → k + e −y , see Eq. (2.28). To make this property manifest, we define a basis of plus distributions in k ± that individually satisfy this property, Note that the leading δ(k − , k + ) term multiplies a double pole in a. The second line implicitly defines the L n (k − , k + , μ) 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 The threshold soft function satisfies the all-order RGE 3) Expanding the threshold soft function in α s as and suppressing all arguments for brevity, S because the hard function is the same in all cases. Here, γ i f (α s ) is the coefficient of δ(1 − z) in the PDF anomalous dimension Eq. (D.14). Solving Eq. (B.6) for γ i thr (α s ), we find where the soft anomalous dimension coefficients γ i S n are given in Eq. (D.7). The boundary coefficients s thr(n) i , which are defined as the coefficients of δ(k − , k + ) in Eq. (B.5), are given by [11,12]  (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 ).

B.3 Extraction method
The double-differential threshold soft function depends on the total lightcone momentum components k ± of the soft hadronic final state. Equivalently, its Fourier transform This implies that at any given order in perturbation theory,Ŝ thr The relevant Fourier transforms between L n thr and L n (k − , k + , μ) follow from the onedimensional 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 dσ/dQ 2 , where the corresponding inclusive threshold soft function S thr i (k 0 , μ) only depends on the total energy k 0 of soft radiation. In particular, S thr i (k 0 , μ) is the process-independent soft contribution to the inclusive partonic cross section σ ab (z) in the soft-virtual limit z → 1, where 1 − z = 2k 0 /Q. In position space, the inclusive threshold soft functionŜ thr i (b 0 , μ) is defined in terms of Wilson lines separated by b 0 (n μ +n μ )/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Ŝ thr i (b 0 , b 0 , μ) =Ŝ thr i (b 0 , μ) .

(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Ŝ thr i (b + , b − , μ) 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 positionspace boundary coefficients of the double-differential soft function at L thr = 0, i.e., at μ = μ * ≡ +i2e −γ 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, 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 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Ŝ thr i (b + , b − , b T , μ) = d 4 k e +ib·k S thr i (k − , k + , k T , μ) , Correspondingly, we define the Fourier transform of S i (k ± , b T , μ, ν) with respect to its lightcone momentum argument k ± aŝ The coefficients up to three loops in the MS scheme are [176,177] β 0 = 11 3 C A − 4 3 T F n f , The coefficients of the MS cusp anomalous dimension to three loops are [89,90,178] i where C i = C F for i = q and C i = C A for i = g.

D.2 Ingredients for T 0
The quark beam function noncusp anomalous dimension coefficients to three loops are [87] γ q B 0 = 6C F , 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] γ g B 0 = 2β 0 , Eur. Phys. J. Plus (2021) 136:214 (K K ) qqV (z) = (K qqV K qqV )(z) + (K qqV K qqV )(z) , (K K ) qq S (z) = K qq S K qqV + K qqV (z) + K qqV + K qqV K qq S (z) + 2n f K qq S K qq S (z) + K qg K gq (z) , (K K ) qq S (z) = K qq S (K qqV − K qqV ) (z) + (K qqV − K qqV )K qq S (z) + 2n f K qq S K qq S (z) , (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, (D.14) We perturbatively expand them in powers of α s /4π, 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 P to distinguish them from our P