N-jettiness beam functions at N3LO

We present the first complete calculation for the quark and gluon $N$-jettiness ($\Tau_N$) beam 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 class of functions appearing in the matching coefficents for all channels includes iterated integrals with non-rational kernels, thus going beyond the one of harmonic polylogarithms. Our results are a key step in extending the $\Tau_N$ subtraction methods to N$^3$LO, and to resum $\Tau_N$ distributions at N$^3$LL$^\prime$ accuracy both for quark as well as for gluon initiated processes.


Introduction
Experimental measurements at the LHC have provided remarkably precise measurements for a multitude of observables, most notably weak gauge boson production, an important benchmark for the Standard Model which has been measured at percent level accuracy [1][2][3][4]. Strong constraints on physics beyond the Standard Model are also provided by precision measurements of Higgs boson production and diboson processes [5][6][7][8][9]. To make full use of these results, it is crucial to confront them with equally-precise theory predictions, which in particular requires to include higher-order corrections in QCD.
So far, only inclusive Drell-Yan and Higgs production have been calculated at next-tonext-to-next-to-leading order (N 3 LO) in QCD [10][11][12][13][14][15][16][17], while significant progress is being made to reach the same precision for differential distributions [18,19]. A key challenge for such calculations is the cancellation of infrared divergences between real and virtual corrections, and hence a necessary prerequisite is a profound understanding of the infrared singular structure at three loops.
N -jettiness (T N ) is an infrared-sensitive N -jet resolution observable and thus provides a way to study the singular structure of QCD [20,21]. Its simplest manifestation T 0 , also referred to as beam thrust, is defined as where the sum rums over all momenta k i in the hadronic final state, q a,b are the momenta of the incoming partons projected onto the Born kinematics, and the measures Q a,b distinguish different definitions of T 0 [22,23]. A key feature of T N is that its singular structure as T N → 0 is fully captured by a factorization theorem, as shown in refs. [20,21] using softcollinear effective theory (SCET) [24][25][26][27][28]. In the simplest case, namely the production of a color-singlet final state h, the appropriate factorization theorem reads Here, Q 2 and Y are the invariant mass and rapidity of h, respectively, and we normalize by the Born partonic cross section σ 0 . In eq. (1.2), the full process dependence is given in terms of the hard function H ab , which encodes virtual corrections to the underlying hard process ab → h. The beam functions B a,b encode radiation collinear to the incoming partons. The soft function S c encodes soft radiation and only depends on the color channel c ∈ {gg, qq}, but is independent of quark flavors. Both beam and soft functions are universal and process independent. Since they are defined as gauge-invariant matrix elements in SCET, calculating them at higher orders also provides a well-defined means of separately studying the collinear and soft limits of QCD themselves. The beam functions B a,b not only appear in the factorization theorem for all T N , but also arise in the factorization theorem for the generalized threshold inclusive color-singlet production in hadronic collisions [29], and are thus of particular interest on their own. Since eq. (1.2) fully captures the singular limit of QCD, it can be employed as a subtraction scheme for higher-order calculations [30,31], in analogy to the q T subtraction method based on a similar factorization for the transverse-momentum distribution [32]. For both methods, extensions to N 3 LO have been recently proposed [18,33]. The O(T 0 /Q) corrections to eq. (1.2) have also been studied in the context of T N subtractions [34][35][36][37][38][39]. 1 These calculations are also interesting on their own as they provide insights into the infrared structure of QCD beyond leading power. T 0 subtractions are also the basis of combining NNLO calculations with a parton shower in GENEVA [41,42].
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 [65]. 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 [66][67][68][69][70], integration-by-part (IBP) identities [71,72] and the method of differential equations [73][74][75][76][77] 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. [65].
This paper is structured as follows. In section 2, we discuss our setup for calculating the beam functions based on the collinear expansion of ref. [65]. In section 3, we briefly present our results, before concluding in section 4. Our results are also provided in the form of ancillary files with this submission.

Beam functions from the collinear limit of cross sections
Since the T N beam function is independent of N , we calculate it from the simplest case T 0 by considering 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 n µ = (1, 0, 0, 1) ,n µ = (1, 0, 0, −1) (2.2) and carry the momenta P 1 and P 2 with the center of mass energy S = (P 1 + P 2 ) 2 . The hard probe h carries the momentum p h , and the total momentum of the hadronic final state is denoted as k. We parameterize these momenta in terms of where Q 2 and Y are the invariant mass and rapidity of the hard probe h, respectively. Eq. (2.1) receives contributions from the partonic process where i and j are the flavors of the incoming partons which carry the momenta p 1 and p 2 , and X n is a hadronic final state consisting of n partons with the momenta {−p 3 , . . . , −p n+2 }, and n = 0 at tree level. The cross section for the partonic process in eq. (2.4), differential in the variables defined in eq. (2.3), is then 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. (2.4), 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. [65]. The partonic cross section in eq. (2.5) is closely related to the beam function we are interested in. For perturbative values of T N , one can match the beam functions onto the PDFs as [20,43] (2.6) Here, I ij is a perturbative matching kernel, and t = T 0 Q a , see eq. (1.2). As shown by us in ref. [65], I ij is precisely given by the strict n-collinear limit of eq. (2.5), where all loop and real momenta are treated as being collinear to n-direction, and we refer to ref. [65] for details on how to calculate this limit: Here, we have regulated both UV and IR divergences by working in d = 4 − 2 dimensions. The renormalized matching kernel is then given by [43,44,65] Here,Ẑ αs implements the standard UV renormalization by renormalizing the bare coupling constant α b s in the MS scheme, and the convolution with the PDF counterterm Γ jk cancels infrared divergences. Explicit expressions for these ingredients are collected in appendix A.3. The remaining poles in are canceled by the convolution with the beam function counter term Z B , which in the formulation of the beam function within SCET arises as an additional UV counter term in effective theory.

Results
In this section we report on our results for the 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 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 [78][79][80][81][82][83][84][85].
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. [17,[86][87][88] and are in part based on refs. [89][90][91]. 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 [92] 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. [65]. We make use of the framework of reverse unitarity [66][67][68][69][70] in order to integrate over loop and phase space momenta. We apply integration-by-part (IBP) identities [71,72] in order to re-express our expanded cross section in terms of collinear master integrals depending on the variables introduced in eq. (2.3). We then compute the required master integrals using the method of differential equations [73][74][75][76][77]. 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. [10,15,[93][94][95].
This yields all required ingredients for the bare partonic cross section expanded in the strict collinear limit of eq. (2.5). This part of the computation is the same as for the results of ref. [96]. Next, we perform the Fourier transform over t and make use of eq. (2.7) to obtain the matching kernel through N 3 LO in QCD perturbation theory. We will elaborate on the details of the computation of the matching kernels in ref. [97]. Finally, we subtract poles in as given in eq. (2.8) to obtain the renormalized matching kernel I ij (t, z, µ) through N 3 LO in QCD perturbation theory. This is carried out in Fourier (y) space, where the convolution in t becomes a simple product, and the Fourier-transformed counter term Z B can be easily predicted from the known renormalization group equation (RGE) of the beam function. We collect the required formulas in appendix A.2. It straightforward to Fourier transform back to t space after the UV renormalization, and we will provide results in both spaces.
We express the perturbative matching kernels in terms of harmonic polylogarithms [98] and Goncharov polylogarithms [99] as well as a set of iterated integrals. We define the iterated integrals recursively via J a 1 (z), a 2 (z), . . . , a n (z) = z 0 dx a 1 (x) J a 2 (x), . . . , a n (x) , (3.1) with the prescription to regularize logarithmic singularities as We refer to the arguments of the iterated integrals as letters. The explicit end point of the iterated integration used for our iterated integrals is always the variablez = 1 − z. In order to express our matching kernels we require the following set of letters (or alphabet): It is possible to rationalise the square root in A by introducing the variable transformation z → (y + 1) 2 /y as noted in ref. [49] and to rewrite the iterated integrals in terms of Goncharov polylogarithms using well known techniques, see for example refs. [100][101][102][103].
Studying the letters of our alphabet and the singularities appearing in our matching kernels we see that they contain logarithmic singularities at the boundaries of the physical interval z ∈ [0, 1]. In order to provide a representation of our perturbative matching kernels that is suitable for numeric evaluation we perform a generalised power series expansion around two different points z = 0 and z = 1 up to 50 terms in the 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 have calculated the matching kernel in Fourier (y) space, where its renormalization becomes simpler. As it is more commonly used in momentum (t) space, we provide results in both spaces. The corresponding kernels are expanded in powers of α s /π, The coefficientsĨ ij can be further expanded as where the logarithm L y and the distribution L m are defined as where the [· · · ] + denotes the standard plus distribution. Note that there is no one-to-one correspondence between theĨ ( ,m) ij (z) and I ( ,m) ij (z), as the Fourier transform induces a nontrivial mixing. For explicit relations for the Fourier transform, see e.g. ref. [104].
The logarithmic terms in eq. (3.5) encode the scale dependence of the beam function, and thus their structure is fully determined by its renormalization group equation (see appendix A.1) in terms of its anomalous dimensions and lower-order ingredients. The genuinely new three-loop results calculated by us are the nonlogarithmic boundary terms I  ij (z). We performed several checks on our results. Firstly, we verified that all poles in cancel after applying UV renormalization and IR subtraction as given in eq. (2.8), where the beam function counterterm was predicted from its RGE as shown in appendix A.2. To check that our results obey the beam function RGE, we verified all logarithmic terms in eq. (3.5) against those predicted in ref. [33] by solving the beam function RGE. We also checked that our results for I ij (z) agree with the eikonal limit lim z→1 I (3) ij (z) that was predicted in ref. [33] using a consistency relation with the threshold soft function [29], and that our results agree with the generalized large-n c approximation n c ∼ n f 1 obtained in ref. [49]. Furthermore, we checked that the first four terms in the soft expansion of the Higgs boson production cross section reproduce correctly the collinear limit of the threshold expansion of the partonic cross section obtained for the purpose of refs. [19,88]. The inclusive cross section at N 3 LO for Drell-Yan and Higgs boson production was obtained in refs. [10,11,14,17,94]. We confirm that we can reproduce 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 I 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 n (1 − z) → ln n (1 − z)/(1 − z). Since the different channels give rise to very different shapes and magnitudes, they are rescaled as indicated for illustration purposes only.
To study the impact of our calculation on the beam function itself, we consider the cumulative beam function where we distinguish both quantities only by their arguments. As indicated, this always involves the sum over all flavors j contribution the desired beam function of flavor i. We use the MMHT2014nnlo68cl PDF set from ref. [105] with α s (m Z ) = 0.118, and evaluate eq. (3.7) through an implementation of our results in SCETlib [106].
In figure 2, we compare the u-quark beam function (left) and gluon beam function (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 t cut = (10 GeV) 2 and µ = 100 GeV and rescale the beam functions by z. Note that the LO result corresponds to the PDF itself, and thus  illustrates the different shape of the beam function compared to the PDF. While we observe a notable effect of the N 3 LO corrections, the beam functions show good convergence overall.
To judge the impact of the new three-loop boundary term I (3) ij on resummed predictions, it is more useful to show the beam function B i (t cut , z, µ) at its canonical scale µ = √ t cut , where all distributions L m in eq. (3.5) vanish and only the boundary term I (3) ij contributes. In figure 3, we show the cumulative beam functions at the canonical scale with µ = √ t cut = 30 GeV, showing the relative difference of the u-quark beam function (left) and the gluon beam function (right) at NLO (green, dotted), NNLO (blue, dashed) and N 3 LO (red, solid) to the corresponding PDF itself. We observe that the shape of the beam functions differ significantly from the shape of the PDF for large z, but tend to converge towards the PDF for small z 10 −1 . As before, we see good convergence at N 3 LO, but still a notable effect of the N 3 LO corrections itself.
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 to the NNLO beam function. As before, we choose the canonical scales µ = √ t cut = 30 GeV as relevant for a resummed calculation, We show the K factor for u quarks (red, solid), d quarks (blue, dashed) and gluons (green, dotted). In all cases, we see corrections of ∼ 1 − 2% with a sizable dependence on z.
For completeness, we also show the high-energy limit z → 0 of the kernels I ij (z) in appendix B. This limit is for example interesting because the small-T 1 region is known to grow at small z in deep inelastic scattering [107,108].

Conclusions
We have calculated the perturbative matching kernel relating N -jettiness beam functions with lightcone PDFs in all partonic channels for the first time at N 3 LO in QCD. Our calculation is based on a method recently developed by us to expand hadronic collinear cross sections [65], demonstrating its usefulness for the calculation of universal ingredients arising in the collinear limit of QCD.
We provide our results in the form of ancillary file with this submission, where we include the renormalized N -jettiness beam function in both momentum (t) and Fourier (y) space. For the t space result, we also provide its expansions around z = 0 and z = 1 through 50 orders in the expansion.
In contrast to the TMD beam functions, which are based on the same collinear limit and at N 3 LO can be entirely expressed in terms of harmonic polylogarithms up to weight 5 [96,109], the T N beam functions have a much richer structure of the appearing functions and are expressed in terms of Goncharov polylogarithm, as well as iterated integrals with letters that involve square roots. It will be interesting to better understand the source of this difference.
Our results have various phenomenological applications. Firstly, we provide a key ingredient to extend the N -jettiness subtraction method [30,31] to N 3 LO, which can be used to obtain exact fully-differential cross sections at this order. They are also crucial to extend the resummation of T N to N 3 LL and N 4 LL accuracy, and for matching N 3 LO calculations to parton showers based on T 0 resummation [41,42].

Acknowledgments
We 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. [97].

A.1 Renormalization group equations
In t space, the beam function B i (t, z, µ) obeys the RGE [20,43] where the anomalous dimension γ i B has the all-order form Here, Γ i cusp (α s ) and γ i B (α s ) are the cusp and beam noncusp anomalous dimensions, which both depend on the color representation i = q or i = g only, but are independent of the quark flavor. The RGE for the matching kernel follows from eqs. (2.6) and (A.1) and the DGLAP equation It is given by [43] µ d dµ (A.4)

A.2 Structure of the beam function counterterm
We define the Fourier transformation of a function f as f (y, · · · ) = dt e −ity f (t, · · · ) , f (t, · · · ) = dy 2π e ityf (y, · · · ) . (A.5) The Fourier transform of the bare kernel I ij (t, z, ) can be conveniently evaluated using Here, L y is the canonical logarithm in Fourier space, and γ E is the Euler-Mascheroni constant. In Fourier space, the renormalization of the bare matching kernel in eq. (2.8) becomes multiplicative in y, and the countertermZ i B follows from the RG eq. (A.2) in y space, Solving eq. (A.8), we can predict the all-order pole structure ofZ i B as (see also ref. [120]) where β(α s , ) = −2 α s + β(α s ) is the QCD beta function in d = 4 − 2 dimensions. Expanding eq. (A.9) systematically in α, we obtain the result through three loops as Here, 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.10) are collected in ref. [33]. The required three-loop results for Γ cusp and β were calculated in refs. [121][122][123] and refs. [124,125], respectively. The beam noncusp anomalous dimension were originally determined in refs. [43,44], see also refs. [63,64].

A.3 α s renormalization and IR counterterms
The bare strong coupling constant is renormalised as The mass factorisation counter term can be expressed in terms of the splitting functions P ij [122,123] 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.12) is defined as (A.13)

B High-energy limit of the beam function kernels
Here, we present the high-energy limit z → 0 of the beam function I C 3 A ln 5 (z) + ln 4 (z) C 2 A C F ln 5 (z) + ln 4 (z)  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 . 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 in electronic form in the ancillary files.