Transverse momentum dependent PDFs at N$^3$LO

We compute the quark and gluon transverse momentum dependent parton distribution functions at next-to-next-to-next-to-leading order (N$^3$LO) in perturbative QCD. Our calculation is based on an expansion of the differential Higgs boson and Drell-Yan production cross sections about their collinear limit. This method allows us to employ cutting edge techniques for the computation of cross sections to extract the universal building blocks in question. The corresponding perturbative matching kernels for all channels are expressed in terms of simple harmonic polylogarithms up to weight five. As a byproduct, we confirm a previous computation of the soft function for transverse momentum factorization at N$^3$LO. Our results are the last missing ingredient to extend the $q_T$ subtraction methods to N$^3$LO and to obtain resummed $q_T$ spectra at N$^3$LL$^\prime$ accuracy both for gluon as well as for quark initiated processes.

TMDPDFs measure both the longitudinal momentum fraction z and the transverse momentum q T carried by the struck parton. They are intrinsically nonperturbative objects that need to be extracted from measurements, but for perturbative | q T | they can be perturbatively related to collinear PDFs. Schematically, this matching takes the form where B i is the so-called TMD beam function for a parton of flavor i, the sum runs over all parton flavors j, I ij is the perturbative matching kernel, and f j is the collinear PDF. The precise knowledge of I ij is important for measurements dominated by transverse momenta that are small compared to the hard scale Q of the process but still perturbative, i.e. Λ QCD | q T | Q, such as Higgs and Drell-Yan production at the LHC. Precise perturbative predictions for the beam function are also essential to extract the intrinsically nonperturbative corrections, which in global fits is typically achieved through a nonperturbative model on top of eq. (1.1), see e.g. refs. [32][33][34][35][36][37].
Since TMDPDFs describe processes at small transverse momentum q T , they are intrinsically sensitive to the infrared structure of QCD. Hence, their perturbative structure is intimately related to the singular structure of QCD amplitudes. This property is employed in the q T subtraction scheme proposed in ref. [38] to achieve the cancellation of infrared divergences in next-to-next-leading order (NNLO) calculations of color-singlet cross sections. Recently, extensions of this method to next-to-next-to-next-to-leading order (N 3 LO) were discussed in refs. [39,40], which however did not include the required three-loop ingredient.
TMDPDFs are composed of the TMD beam function B i (z, q T ) and the TMD soft function S( q T ). The soft function has been known at three loops since quite some time [41][42][43], and the quark beam function has been calculated at this accuracy recently [44][45][46][47][48][49], while the gluon beam function so far is only known at two loops [46,47,50,51]. In this paper, we fill this gap by calculating the full matching coefficient for the gluon beam function at N 3 LO. We also calculate the full quark beam function at N 3 LO, where we find disagreement with the recent calculation of the corresponding result in ref. [49] in the d abc d abc color structure. These results make it possible to fully apply the q T subtraction at N 3 LO accuracy, paving the way for fully-differential cross sections of color-singlet processes at this order. Our results are also the last missing ingredient for TMD resummation at N 3 LL accuracy. They also arise in q T -dependent event shapes at hadron colliders such as the Transverse Energy-Energy Correlator (TEEC) [52], and for the azimuthal angle in vector boson +j production in the back-to-back limit [53,54].
Our computation is based on a method of expanding cross sections around the kinematic limit in which all final state radiation becomes collinear to one of the scattering hadrons [55]. This method allows one to efficiently connect technology for the computation of scattering cross sections to universal building blocks of perturbative QFT. In particular, we perform a collinear expansion of the Drell-Yan and gluon fusion Higgs boson production cross section at N 3 LO. Subsequently, we employ the framework of reverse unitarity [56][57][58][59][60], integration-by-part (IBP) identities [61,62] and the method of differential equations [63][64][65][66][67] to obtain the collinear limit of the cross sections differential in the rapidity and transverse momentum of the colorless final states. Using the connection of this limit to the desired beam functions we extract the desired perturbative matching kernels as discussed in ref. [55].
This paper is structured as follows. In section 2, we briefly review TMD factorization. In section 3, we discuss how the beam function can be calculated from the collinear limit of a color-singlet cross section using the method collinear expansions. In section 4, we briefly present our results, before concluding in section 5. Our results are also provided in the -2 -form of ancillary files with this submission.

Review of q T factorization
We consider the production of a colorless hard probe h and an additional hadronic state X in a proton-proton collision, where the incoming protons are aligned along the directions and carry the momenta P 1 and P 2 with the center of mass energy S = (P 1 + P 2 ) 2 . We are interested in measuring the cross section differential in p µ h , expressed through the invariant mass Q 2 = p 2 h , rapidity Y , and transverse momentum q T . The factorization of the cross section in the limit q T Q was first established by Collins, Soper, and Sterman (CSS) [68][69][70], and was further elaborated on in refs. [71][72][73][74]. The factorization was also shown using Soft-Collinear Effective Theory (SCET) [75][76][77][78] by several groups [79][80][81][82][83][84][85][86]. The factorized cross section is typically expressed in Fourier space in two equivalent forms, For processes inclusive in h, eq. (2.3) holds up to corrections in O(q 2 T /Q 2 ). Power corrections of O(q 2 T /Q 2 ) have been firstly calculated at fixed order in perturbation theory in ref. [87], while the study of their all order structure has been initiated using SCET operator formalism [88][89][90][91][92] and their nonperturbative structure has been explored in refs. [93,94]. Eq. (2.3) receives enhanced O(q T /Q) corrections when applying fiducial cuts to h [95] that can be uniquely included in the factorization for Higgs and Drell-Yan production [96], and also receives linear corrections when radiation from massive final states is considered [97].
In eq. (2.3), we sum over all parton flavors a and b mediating the underlying partonic process ab → h. The corresponding partonic Born cross section is denoted by σ 0 , and the hard function H ab = 1 + O(α s ) encodes virtual corrections to the Born process. Thẽ B i (x, b T , µ, ν/ω) encode the probability to find a parton of flavor i with longitudinal momentum fraction x and impact parameter b T , which is Fourier-conjugate to q T . The soft functionS encodes the effect of soft radiation from either proton. In the second line of eq. (2.3), these functions are combined into the TMDPDF which is indepdent of the rapidity scale ν discussed below. Computationally, it is natural to separately consider the calculation of the beam and soft functions appearing in eq. (2.3), which can easily be combined into the TMDPDF if desired. A characteristic feature of q T factorization is the appearance of so-called rapidity divergences [68,79,82,85,[98][99][100], which require an explicit rapidity regulator. Similar to the emergence of the renormalization scale µ from UV regularization, this induces a rapidity scale, which we generically as ν. The rapidity divergences track the energy of the struck partons, encoded in the parameters There is a variety of rapidity regulators, giving rise to several formulations of the individual ingredients in eq. (2.3). However, all approaches yield the same fixed-order results for the physical cross section in eq. (2.3). Thus, we are free to choose the regulator most convenient for our calculation, and we will employ the exponential regulator of ref. [86]. Explicit definitions of the beam and soft functions in terms of gauge-invariant matrix elements formulated in SCET can be found in ref. [86] (see also refs. [48,51]), but are not required in our approach.
In this work, we focus on TMD factorization in the perturbative regime b −1 T ∼ q T Λ QCD , in which the TMD beam function and TMDPDF can be matched onto PDFs as [70] B i z, b T , µ, ν ω = 3 Beam functions from the collinear limit of cross sections We consider the contribution to eq. (2.1) from the partonic process where i and j are the flavors of the incoming partons which carry the momenta p 1 and p 2 . X n is a hadronic final state consisting of n partons with the momenta {−p 3 , . . . , −p n+2 } and total momentum k, with n = 0 at tree level. We parameterize the final-state momenta in terms of where Q 2 and Y are the invariant mass and rapidity of the hard probe h, respectively. The partonic cross section differential in these variables is defined as Here, we normalize by the partonic Born cross section σ 0 , dΦ h+n is the phase space measure of the h + X n state, and |M ij→h+Xn | 2 is the squared matrix element for the process in eq. (3.1), summed over the colors and helicities of all particles, with N ij accounting for the color and helicity average of the incoming particles. Explicit expressions for N ij and dΦ h+n can be found in ref. [55]. In ref. [55], we showed that the matching coefficient in eq. (2.6) is precisely given by the strict n-collinear limit of eq. (3.3), i.e. the limit where all loop and real momenta are treated as being collinear to n-direction, and we refer to [55] for details on its calculation: Solving the δ functions fixes w 1 and w 2 as . (3.5) The superscript naive in eq. (3.4) indicates that further steps are required to obtain the desired matching kernel. First, we note that the integral in eq. (3.4) contains the aforementioned rapidity divergences, namely divergences as x → 1 or z → 1 that are not regulated by dimensional regularization. We regulate these using the exponential regulator of ref. [86], which in fact is the only regulator in the literature compatible with our approach. Inserting the regulator factor exp(2τ e −γ E k 0 ) expressed in the above variables, we obtain the regulated kernel as [55] . (3.6) Here, ω = Qe Y is the so-called label momentum. The exponential factor in eq. (3.6) regulates divergences as x, z → 1, with τ being the rapidity regulator. UV and IR divergences are regulated by working in d = 4 − 2 dimensions, with the limit → 0 being taken after the limit τ → 0, as indicated. It is convenient to expand the exponential factor in eq. (3.6) in terms of distributions before carrying out the integral [48], and we provide more details on this in appendix A.1. It is common to express TMD beam functions in Fourier space, where convolutions in q T are replaced by simple products, which in particular greatly simplifies the resummation of large logarithms [101]. Denoting the Fourier-transform matching kernel asĨ ij , we obtain the renormalized kernel as Here, following ref. [42] we identify ν ≡ 1/τ as the rapidity renormalization scale [85]. The so-called zero-bin subtraction [102] to subtract overlap of the beam function with the soft function is implemented by dividing by the soft functionS [48]. In eq. (3.7),Ẑ αs implements the standard UV renormalization by renormalizing the bare coupling constant α b s in the MS scheme as stated in eq. (A.18). Infrared divergences are canceled through the convolution with the PDF counterterm Γ jk , which is given in eq. (A.19). The remaining poles in are canceled by the beam function counter termZ B , which in the formulation of the beam function within SCET arises as an additional UV counter term in effective theory.
Since the bare soft function is not given in the literature, we have directly calculated it from the soft limit of eq. (3.3) similar to eq. (3.6). (3.8) In the strict soft limit, both w 1 and w 2 are treated as small quantities, such the measurement δ function and the exponential regulator in eq. (3.8) are simpler than in eq. (3.4). The Fourier transform of eq. (3.8) yields the bare soft functionS(b T , , τ ) required in eq. (3.7).
We have also verified that the renormalized soft function reproduces the N 3 LO result of ref. [42]. Since the bare soft function in the exponential regulator is only given at NLO in the literature [86], we provide the bare soft function in the form of ancillary files with this submission.
The strict soft limit of the general differential partonic coefficient function is obtained by expanding the strict collinear limit in w 1 and maintaining only the first term of the generalised power series. At n th order in perturbation theory the strict soft limit of the partonic coefficient function takes the form where η (n) strict soft (x, ) is independent of w 1 and w 2 . Note that eq. (3.9) is related to the bare fully-differential soft function which measures the total soft radiation in a process. This limit is also easily related to the bare threshold soft function [103] via where z is the threshold parameter. We have used this relation as an additional check on our soft limit.

Results
In this section we report on our results for the TMD matching kernels through N 3 LO. Our computation is based on the collinear expansion of the cross sections for the production of a Higgs boson via gluon fusion and for the production of an off-shell photon (Drell-Yan) in hadron collisions. We compute the Higgs boson production cross section in the heavy top quark effective theory where the degrees of freedom of the top quark were integrated out and the Higgs boson couples directly to gluons [104][105][106][107][108][109][110][111]. We begin by computing all required matrix elements with at least one final state parton to obtain N 3 LO cross sections. All partonic cross sections corresponding to matrix elements with exactly one parton in the final state were obtained in full kinematics for the purpose of refs. [112][113][114][115] and are in part based on refs. [116][117][118]. In order to obtain the strict collinear limit we simply expand the existing results and select the required components.
To compute partonic cross sections with more than one final state parton we generate the necessary Feynman diagrams using QGRAF [119] and perform spinor and color algebra in a private code. Subsequently, we perform the strict collinear expansion of this matrix elements as outlined in ref. [55]. We make use of the framework of reverse unitarity [56][57][58][59][60] in order to integrate over loop and phase space momenta. We apply integration-by-part (IBP) identities [61,62] in order to re-express our expanded cross section in terms of collinear master integrals depending on the variables introduced in eq. (3.2). We then compute the required master integrals using the method of differential equations [63][64][65][66][67].
In order to fix all boundary conditions for the differential equations we expand the collinear master integrals further around the soft limit and integrate over phase space. The result of this procedure is then easily matched to the soft integrals that were obtained for the purpose of refs. [120][121][122][123][124].
With this we obtained all required ingredients for the bare partonic cross section expanded in the strict collinear limit of eq. (3.3). This part of the computation is the same as for the results of ref. [125]. Next, we obtain the bare matching kernel via eq. (3.6) and subsequently perform the Fourier transform over q T using eq. (A.6). We will elaborate on the details of the computation of the matching kernels in ref. [126]. Finally, the renormalized matching kernel is obtained using eq. (3.7), where the beam function counter termZ i B was predicted from the renormalization group equations (RGEs) of the beam function as shown in appendix A.4. The TMDPDF can then be straightforwardly obtained by combining the beam function with the soft function as in eq. (2.4).
We expand the matching kernelsĨ ij of the beam function and the matching kernels I TMD ij of the TMDPDF obtained in this way in powers of α s /π, can be written as a polynomial in logarithms of the appearing scales with z-dependent coefficient functions, The logarithms in eq. (4.2) are defined as where γ E is the Euler-Mascheroni constant, and we remind the reader that ζ = ω 2 is the energy of the struck parton. The logarithmic terms with m > 0 or n > 0 in eq. (4.2) encode the scale dependence of the beam function and TMDPDF, and thus their structure is determined its RGEs (see appendix A.3) in terms of its anomalous dimensions and lower-order ingredients. The genuinely new three-loop result calculated by us is the nonlogarithmic beam function boundary termĨ Remarkably, it can be expressed entirely in terms of standard plus distributions and harmonic polylogarithms [127] of argument z and transcendental weight less or equal to five.
To allow for an easy numeric implementation, we also provide a generalized power series expansion of our results around z = 0 and z = 1 with up to 50 terms in each expansion. Both power series are formally convergent within the entire unit interval but converge of course faster if the respective expansion parameter is smaller. We provide both power series for all matching kernels as well as the analytic solution in ancillary files together with the arXiv submission of this article. We performed several checks on our results. Firstly, we verified that the UV and IR subtraction as given in eq. (3.7) correctly removes all poles in . Our NLO and NNLO results for the renormalized beam function are validated against ref. [43], and the bare results through O( 4 ) at NLO and through O( 2 ) are verified against refs. [48,51]. To verify that our results obey the beam function RGE, we verified all logarithmic terms in eq. (4.2) against those predicted in ref. [40] by solving the beam function RGE. We also verified the eikonal limit which was derived in ref. [40] from consistency with joint q T and soft threshold resummation relations [86,128], and also conjectured in ref. [47]. In eq. (4.5), γ r 2 is the three-loop coefficient of the so-called rapidity anomalous dimension [85], as given in Eq. (16) of ref. [42] (see also ref. [129]), where the appropriate color structure is implicit. The rapidity anomalous dimensions is also closely related to the Collins-Soper kernel of refs. [68,69]. For the quark beam function, we also compared with the results recently obtained in ref. [49]. We find discrepancies for terms proportional to the color structure d abc d abc entering in all quark-to-quark kernels. After private communication, the authors of ref. [49] identified and resolved a minor mistake in their calculation, after which they find agreement with our result. Furthermore, we checked that the first four terms in the soft expansion of the Higgs boson production cross section correctly reproduce the collinear limit of the threshold expansion of the partonic cross section obtained for the purpose of refs. [115,130]. The inclusive cross section at N 3 LO for Drell-Yan and Higgs boson production was obtained in refs. [112,121,122,131,132]. We also reproduced the first term of the threshold expansion of all partonic initial states contributing to the collinear limit of the partonic cross sections using the collinear partonic coefficient functions obtained here after integration over phase space.
To illustrate our results, figure 1 shows the beam function boundary termsĨ ij (z) relevant for the quark beam function (left) and gluon beam function (right) as a function of z. For the purpose of this plot, we replace the occurring distributions as δ(1 − z) → 0 and L 0 (1 − z) → (1 − z) −1 . Since the different channels give rise to very different shapes and magnitudes, they are rescaled as indicated for illustration purposes only.
Next, we study the impact of our calculation on the beam function and TMDPDF themselves. We use the MMHT2014nnlo68cl PDF set from ref. [133] with α s (m Z ) = 0.118, and evaluate eq. (2.6) through an implementation of our results in SCETlib [134].
In figure 2, we show the u-quark TMDPDF (left) and the gluon TMDPDF (right) at LO (gray, dot-dashed), NLO (green, dotted), NNLO (blue, dashed) and N 3 LO (red, solid) as a function of z. We fix the impact parameter b T = 10 GeV −1 , parton energy ω = 100 GeV and renormalization scale µ = 100 GeV. Note that the LO result corresponds to the PDF itself, and thus illustrates the different shape of the beam function compared to the PDF. With the inclusion of the N 3 LO result obtain in this work, both the quark and the gluon TMDPDFs nicely show convergence over a large range of values for z.
To judge the impact of the new three-loop boundary termĨ (3) ij on resummed predictions, it is more useful to show the beam function for the canonical scale µb T = 2e −γ E and ν = ω, where all logarithms in eq. (4.2) vanish and only the boundary termĨ   In figure 3, we compare the u-quark beam function (left) and gluon beam function (right) at NLO (green, dotted), NNLO (blue, dashed) and N 3 LO (red, solid) to the corresponding PDFs, choosing canonical scales for b T = (10 GeV) −1 and ω = 100 GeV. We see that the beam function has a very different shape compared to the PDF, and that the beam function converges very well at N 3 LO. Finally, in figure 4 we show the K-factor of the N 3 LO beam function, which we define as the ratio of the N 3 LO beam function to the NNLO beam function. As before, we choose the canonical scales for b T = (10 GeV) −1 and ω = 100 GeV. We show the K factor for u quarks (red, solid), d quarks (blue, dashed) and gluons (green, dotted). In all cases, we find rather small correction of ∼ 0.2 − 0.5%, but with a notable dependence on z.
For completeness, we also present the high-energy limit z → 0 of the kernelsĨ gq (z) contributing to the gluon beam function in appendix B. The corresponding limit for the quark kernels were already presented in ref. [49], for which we find perfect agreement. These results are useful to study the small-x behavior of TMDPDFs, see e.g. refs. [135][136][137][138].

Conclusions
We have calculated the perturbative matching kernel relating transverse-momentum dependent beam functions with lightcone PDFs at N 3 LO in QCD. This provides the first results of these kernels for the gluon TMD beam function, and corrects the result in the d abc d abc color structure in the recent calculation of the quark TMD beam function in ref. [49]. After private communication, the authors of ref. [49] identified and resolved a minor mistake in their calculation, after which perfect agreement is found. This emphasizes that having two independent computations that are in accordance with each other is the most stringent cross check. Our calculation is based on a method recently developed by us to expand hadronic collinear cross sections [55], demonstrating its usefulness for the calculation of universal ingredients arising in the collinear limit of QCD. As a byproduct, we also confirmed a previous computation of the soft function for transverse momentum factorization [86]. We provide our results in the form of ancillary file with this submission, where we include the renormalized TMD beam function, as well as its expansions around z = 0 and z = 1 through 50 orders in the expansion, the renormalized TMDPDF, and the bare and renormalized soft function.
Our results have various phenomenological applications. Firstly, we provide the last missing ingredient for the fully-differential calculation of color-singlet processes at N 3 LO using the q T subtraction method [38][39][40], which can be used to obtain the first exact fullydifferential cross sections at this order. They also enable the resummation of transverse momentum distributions at hadron colliders at N 3 LL accuracy, both for gluon and quark induced hard scatterings such as Higgs-boson production and the Drell-Yan process, two key observables at the LHC.
A natural future direction of this research is the calculation of the closely related TMD fragmentation functions at N 3 LO, which are required to describe the small-q T limit semi-inclusive deep inelastic scattering and are currently only known at NNLO [47,48,51].
With the same techniques presented in this work, it would also be interesting to consider the calculation of other q T -dependent time-like collinear functions, such as the Energy-Energy Correlator (EEC) jet functions. The quark and gluon EEC jet functions enter the factorization of the EEC in the back-to-back limit for e + e − annihilation and gluon initiated Higgs decay, respectively [139], as well as that for the TEEC [52]. Their knowledge at N 3 LO is the only missing ingredient to achieve resummation of the EEC in the back-to-back limit at N 3 LL accuracy. Note that at fixed order the full angle EEC has been calculated analytically through O(α 2 s ) [140,141]. It will also be interesting to study the collinear expansions beyond leading power to shed light on the structure of TMD factorization at subleading power, in particular on the structure of rapidity divergences at subleading power [87,142,143].

Acknowledgments
We thank the Ming-xing Luo, Tong-Zhi Yang, Hua Xing Zhu and Yu Jiao Zhu for private communication concerning the results of ref. [49] prior to the publication of this article that allowed both groups to agree on our results. We also thank Johannes Michel, Iain Stewart and Frank Tackmann for useful discussions.

A Ingredients for the calculation of the beam function
In this appendix, we provide more details on the regularization and renormalization of the beam function kernels. Details of the calculation of all required integrals will be presented in ref. [126].

A.1 Rapidity regularization
In practice, it is useful to expand the regulator in eq. (3.6) in terms of plus distributions, which allows one to take the limit τ → 0 before evaluating the integral in eq. (3.6), see also [48] for a discussion at two loops. In general we want to expand a power divergence using the relation (A.1) Here, f (x) is a suitable test function that is holomorphic at x = 0 and a represents a generalized power, typically proportional to the dimensional regulator . We find that We furthermore require the two-dimensional generalization of the above integral, The regularization of an integral of the type of the above equation including a test function follows along the lines of eq. (A.1).
In order to compute the soft function we use the following relations

A.2 Fourier transform
The Fourier transform required when going from eq. (3.6) to eq. (3.7) can be conveniently evaluated using where κ = τ q 2 T /ω, is the loop-order, and we express all resulting logarithms in terms of In eq. (A.6) we divide by the angular factor 1 2 Ω 1−2 = π 1− /Γ(1 − ) and q −2 T to account for the fact that q T is defined as the magnitude of the (2 − 2 )-dimensional vector, and that the associated 2 − 2 -dimensional solid angle has already been integrated over in I ij .

A.3 Renormalization group equations
The beam function depends on the renormalization scale µ and the rapidity scale ν, and thus obeys two coupled RGEs [85,86] whereγ i ν is the so-called rapidity anomalous dimension [85]. Its prefactor or −1/2 arises becauseγ i ν is defined as the ν-anomalous definition of the soft function. The beam anomalous dimension has the all-order form where Γ i cusp (α s ) andγ i B (α s ) are the cusp and the beam noncusp anomalous dimensions in the appropriate color representation i = q or i = g, but are independent of the quark flavor.
The RGE for the matching kernel follows from eqs. (2.6) and (A.8) together with the DGLAP equation It is given by The rapidity anomalous dimension itself is governed by an RGE in µ, which can be solved as (A.14) This anomalous dimension appears in the eikonal limit in eq. (4.5), and is related to the notation of ref. [42] by γ r 2 = 2γ i ν 2 , where the color representation i is implicit in γ r 2 .

A.4 Structure of the beam function counterterm
The beam function counterterm can be predicted from eq. (A.8) using The all-order form of the counterterm is given by (see also ref. [144]) where β(α s , ) = −2 α s + β(α s ) is the QCD beta function in d = 4 − 2 dimensions. Expanding eq. (A.16) systematically in α, we obtain the result through three loops as Here, L ω = ln(ν/ω), and the γ n are the coefficients of the corresponding anomalous dimensions at O[(α s /4π) n ]. Explicit expressions for all anomalous dimensions in the convention of eq. (A.17) are collected in ref. [40]. The required three-loop results for Γ cusp and β were calculated in refs. [145][146][147] and [148,149], respectively. The coefficients ofγ B follow from consistency with the anomalous dimensions of the hard and soft functions in eq. (2.3), which can be obtained from the quark and gluon anomalous dimensions of the corresponding form factors calculated in refs. [150][151][152][153][154][155][156]. Our calculation explicitly confirms the beam anomalous dimension obtained from these relations.

A.5 α s renormalization and IR counterterms
The bare strong coupling constant α b s is renormalised as The mass factorisation counter term can be expressed in terms of the splitting functions [146,147] as Here, we suppress the argument z of the splitting functions on the right hand side and keep the summation over repeated flavor indices implicit. The convolution in eq. (A.19) is defined as (A.20) B High-energy limit of the beam function kernels The high-energy limit z → 0 of the kernelsĨ gq (z) contributing to the gluon beam function is given by Here, the color factors C A and C F are only used for compactness of the result and should be replaced with their expressions in terms of n c . The corresponding limit for the quark kernels were already presented in ref. [49], for which we find perfect agreement. Note that the expressions for the high energy limit z → 0 up to O(z 50 ), as well as that for the threshold limit z → 1 up to O((1 − z) 50 ), can be found for all channels in electronic form in the ancillary files of this work.