The fully-differential gluon beam function at NNLO

The fully-differential beam function (dBF) is a universal ingredient in resummed predictions of hadron collider observables that probe the full kinematics of the incoming parton from each colliding proton — the virtuality and transverse momentum as well as the light-cone momentum fraction x. In this paper we compute the matching coefficients between the unpolarized gluon dBF and the usual parton distribution functions (PDFs) at the two-loop order. For observables probing both the virtuality and transverse momentum of incoming gluons, our results provide the part of the NNLO singular cross section related to collinear initial-state radiation, and are required for the resummation of large logarithms through N3LL. Further to this, the dBF is closely linked to the beam function appearing in a generalized version of threshold factorization, via a simple integration. By performing this integration for the two-loop gluon matching coefficients, we also obtain the corresponding quantities for the generalized threshold beam function.


Introduction
Beam functions encode the effect of collinear initial state radiation (ISR) on the measurement of some observable T at a hadron collider [1,2]. If T is a perturbative scale (T Λ QCD ), these functions can be expressed in terms of the convolution of the standard PDFs f j and perturbatively-calculable matching coefficients I ij . The indices i, j = {q i ,q i , g} denote the different types of partons in QCD. In this paper we focus on the beam functions B i (t, x, k ⊥ ) [3], 1 which can be thought of as a parton density completely differential in all of the kinematic variables of the 'active' parton that enters the hard process. It is a main ingredient in factorized cross sections for observables sensitive to k ⊥ and t in the limit where Q represents the scale of the hard interaction. The lightcone coordinates 2 of the active parton are given in terms of the arguments of this beam function by (xP − , −t/[xP − ], k ⊥ ), with P − the large component of the proton momentum. The B i (t, x, k ⊥ ) are often referred to as fully-differential beam functions (dBFs) [3], or fully unintegrated PDFs [5,6].
The full set of dBF matching coefficients I ij was calculated at one loop in ref. [3]. In ref. [7], we computed the quark matching coefficients I qj at the two-loop level. The 1 See also ref. [4] for related beam functions defined in impact parameter space. 2 Here we use the SCET conventions for light-cone coordinates. That is, for a Lorentz vector v the components (v + , v − , v ⊥ ) are defined by v µ = v + n µ /2 + v −nµ /2 + v µ ⊥ , where n andn are oppositelypointing lightcone vectors and v µ ⊥ is a vector transverse to these directions: n ·n = 2, n 2 =n 2 = 0, n · v ⊥ =n · v ⊥ = 0. In our notation v ⊥ is a Euclidean transverse vector: v 2 ⊥ = −v 2 ⊥ .

JHEP07(2020)234
goal of the present paper is to extend this computation to the case of the gluon matching coefficients I gj . In an unpolarized proton, there are two varieties of gluon dBFs: the unpolarized and linearly polarized dBFs (similar to what is observed for the transverse momentum dependent beam functions, or TMD PDFs [8]). We limit our attention here to the computation of the former quantity at two loops, which is sufficient for NNLO or N 3 LL resummed predictions of the most important gluon-initiated processes like Higgs production, see section 2.
Our motivation for performing this calculation is the large and growing list of potential applications for the dBF: • The gluon dBF is required to make precise resummed predictions of the Higgs transverse momentum distribution in the presence of a jet veto, where the jet veto is imposed via a cut on a virtuality-sensitive observable [3,[9][10][11]. In Higgs measurements, veto constraints of this kind are often imposed to enhance the Higgs signal over the backgrounds [12,13], or categorise the cross section in terms of different initial states [14].
• Beam functions are a core component of the Geneva methodology for combining parton showers with higher order resummation and fixed order computations [15][16][17][18]. The two-loop gluon dBF would be a necessary ingredient in constructing a shower in the context of (for example) Higgs production that was simultaneously NNLL accurate in both the Higgs p T and beam thrust T 0 , and matched to NNLO.
• The two-loop gluon dBF could be utilised in multi-differential extensions of the Njettiness subtraction scheme [19,20] in computations of NNLO gluon-initiated processes. This is discussed in section 3.4 of ref. [19].
• Finally, in ref. [21] a generalized threshold factorization theorem was derived for colour singlet production that holds when the momentum fraction of one of the incoming partons approaches unity. This corresponds to the limit of large rapidity but generic invariant mass of the produced colour singlet. A beam functionB j (t, x) appears in this factorization theorem that can be obtained from the dBF as follows: Thus one can straightforwardly obtain the two-loop matching coefficients for the generalized threshold beam function from our two-loop results for the gluon dBF. In this paper we actually perform this exercise, and include the results below. Note that only the unpolarized part of the dBF contributes in eq. (1.1).
The paper is organised as follows. In section 2 we give the operator definition of the dBFs and recall their key properties as well as their distributional structure. For most of the calculational details we refer to previous work. Section 3 contains the novel results for the two-loop gluon dBF matching coefficients. In section 4 (and appendix C) we present the corresponding results for the generalized threshold beam function. Finally we conclude in section 5.

JHEP07(2020)234
2 Definitions and properties of the dBF Our calculation of the two-loop gluon dBF is performed in close analogy to our calculation of the two-loop quark dBF [7] and largely based on the two-loop calculation of the virtualitydependent beam functions [22,23]. In particular the integrand of the NNLO gluon beam function is taken from ref. [23]. The (loop) integrations are then performed with the additional constraint on the transverse momentum crossing the unitarity cut as in ref. [7]. For details on the calculational methods we therefore refer the reader to the above three references and shall only review the generic structure and some important properties of the dBFs here.
The bare gluon dBF is defined by the operator matrix element [3] B µν in soft-collinear effective theory (SCET) [24][25][26][27][28][29]. Here p n (p − ) denotes the incoming spinaveraged proton state with lightlike momentum p µ = p − n µ /2, x ≡ k − /p − and B µ n⊥ is the gauge-invariant gluon field strength operator in SCET. The measurement delta functions involve the SCET label momentum operators P n⊥ and P n ≡n · P n [26] as well as the usual (derivative) plus-momentum operatorp + . The label momentum operators act to the right on the SCET fields in B µ n⊥ , whilep + also acts on the external proton state. For further details and explanations on the relevant SCET operators, see e.g. refs. [1,2,22].
The gluon dBF is a Lorentz tensor transverse to the lightcone directions n µ andn µ , i.e. n µ B µν g =n µ B µν g = 0, as a consequence of the transverse polarization of the physical gluon field operators generating it. As it depends on the external two-dimensional transverse momentum vector k ⊥ the renormalized gluon dBF (with renormalization scale µ) naturally decomposes into two orthogonal tensor structures 3 in four dimensions: The two scalar coefficients B g and B g depend on k ⊥ only through k 2 ⊥ and can be regarded as independent (unpolarized and linearly polarized) beam functions. In particular, the corresponding projections of the dBF operator do not mix under renormalization. We can thus choose different regularization schemes for the computation of the respective bare functions. While a momentum space calculation of B g requires a regularization scheme that treats k ⊥ as strictly two-dimensional vector, there is in fact no such restriction for B g . That is because we may replace in the dBF matrix element B g ≡ B µ gµ , eq. (2.1), which is not possible for B g since the second term in eq. (2.2) depends on k ⊥ rather than k 2 ⊥ . We can thus employ conventional JHEP07(2020)234 dimensional regularization (CDR) for the computation of B g by interpreting k ⊥ on the r.h.s. of eq. (2.3) as a (d − 2) dimensional vector with d = 4 − 2 . This in turn is equivalent to replacing in our calculation of the bare B g and in analogy to what we did in our calculation of the two-loop quark dBF [7]. While intermediate (bare) results may differ depending on the dimensional regularization scheme, the renormalized results do not. The argument closely follows the one given in ref. [30] in the context of the transverse momentum dependent soft function. It is based on the observation that all 1/ n poles associated with ultraviolet, infrared, or rapidity divergences are, after the subtraction of subdivergences in the renormalization or matching procedure, proportional to δ( k 2 ⊥ ) = πδ (2) ( k ⊥ ) and therefore cannot promote the O( ) differences in the transverse momentum measurement to a finite contribution. Throughout this work we use the MS renormalization scheme.
For all gluon-initiated processes that are insensitive to the polarization of the incoming gluons, like Higgs or (azimuthally averaged) top pair production (cf. refs. [31,32]) only the NLO result for B g [3] contributes to NNLO (leading power) or N 3 LL resummed predictions. On the other hand this precision does require the NNLO result of B g . The reason is that due to the orthogonality of the two tensor structures in eq. (2.2) only diagonal terms ∝ B g ⊗ B g or ∝ B g ⊗ B g and no mixed terms occur in the factorized cross section for these processes. At LO (tree-level) B g vanishes (unlike B g ) and thus only the NLO ⊗ NLO term of B g ⊗ B g type survives at NNLO. 4 It is the aim of this paper to compute the so far unknown NNLO correction to B g , i.e. more precisely the matching coefficient In practice we determine the I gj coefficients from a perturbative calculation along the lines of ref. [7], where the incoming proton is replaced by a single parton on both sides of eq. (2.5). The relevant two-loop diagrams (in lightcone axial gauge) are the same as the ones shown in ref. [23]. Note that k 2 ⊥ is constrained by [ where z is the large lightcone momentum fraction argument of I gj as in eq. (2.5). That is because the total invariant mass of the collinear radiation emitted in the transition of parton j to the active gluon must be non-negative. As z ≥ x ≥ 0, eq. (2.6) also holds for z → x. Integrating over the full range of k 2 ⊥ yields the virtuality (t) dependent gluon beam JHEP07(2020)234 function. For the corresponding matching coefficients this means with the I gj (t, x, µ) known to NNLO [23,36]. 5 As explained in ref. [7] we can use eq. (2.7) to fix the δ(1 − z) δ( k 2 ⊥ ) term in I gg (t, z, k 2 ⊥ , µ) provided we know the result for z < 1. It equals the corresponding term in I qq (t, z, k 2 ⊥ , µ) upon replacing C F → C A in the latter, see ref. [23] for an explanation. There is no δ(1 − z) (endpoint) term in I gq .
As a consequence of eq. (2.6) the matching functions I ij (t, z, k 2 ⊥ , µ) obey the same renormalization group equation (RGE) as the I ij (t, z, µ) [3]: 9) and the collinear QCD splitting function P ij (z, µ). The latter are completely known to three loops [37,38]. We give the LO expressions in appendix B.1. The (virtuality-dependent) beam function anomalous dimension γ i B (t, µ) equals the jet function anomalous dimension, , which was derived for i = g to three-loop order in refs. [36,39]. Recently, the computation of the universal cusp anomalous dimension Γ i cusp to four loops in QCD has been completed [40][41][42][43][44][45][46]. Thus all pieces in the dBF anomalous dimension necessary for N 3 LL resummation are now available. The relevant cusp and non-cusp contributions to γ g B (t, µ) are collected up to three loops in the appendices of refs. [22,23] using our notation.
The first two terms (n = 0, 1) [3] in the expansion of the dBF matching coefficient, JHEP07(2020)234 with the usual plus distributions The anomalous dimension coefficients in the expression for I with n f the number of (massless) quark flavours. The coefficients of the δ(t)δ( k 2 ⊥ ) terms at one and two loops, I ij (z) and I (2) ij (z), agree by virtue of eq. (2.7) with the δ(t) coefficients in the virtualitydependent beam function and were computed in refs. [23,36]. The I (1) gj (z) can also be found in appendix A. We explicitly carry out the convolutions involving the one-loop splitting functions in the next-to-last line of eq. (2.11) for i = g in appendix B.2.
The novel results of this paper are the functions J gj (t, z, k 2 ⊥ ). We present the explicit expressions in the next section. Although they multiply the µ-dependent distribution L 0 (t/µ 2 ) in eq. (2.11) they cannot be obtained by solving the RGE, because ij (t, z, t r) .

(2.13)
Hence, only the k ⊥ integral of the J ij (t, z, k 2 ⊥ ), but not the function itself, is fixed by the dBF anomalous dimension. In eq. (2.13) we used the distributional identity [7] which holds for any integrable function f ( k 2 ⊥ ) as can be verified by integration over k 2 ⊥ . The θ-function in eq. (2.14) implements the kinematic constraint in eq. (2.6) and is implicit in the J

JHEP07(2020)234
3 Results for the dBF Our calculation of the two-loop dBF perfectly reproduces the known terms in eq. (2.11), which serves as a strong cross check, and yields the following new results for the NNLO coefficient functions J (2) gj (t, z, k 2 ⊥ ). Our notation is in close analogy to the one for the NNLO coefficients of the virtuality-dependent gluon beam function in ref. [23]. We define 6 where q i (q i ) denotes the (anti)quark of flavor i. Together with the finite support of the PDFs in eq. (2.5) the two θ-functions in eqs. (3.1) and (3.2) enforce the kinematic constraints (1 − x)/x ≥ r ≡ k 2 ⊥ /t ≥ 0 and 1 ≥ x ≥ 0. They also imply that taking the limits z → 1 or t → 0 of the J (2) gi coefficients requires to integrate over the full range of k 2 ⊥ (or equivalently r) first, see ref. [7] for details. The explicit expressions for J (2) ggA , J ggF , and J Throughout this work θ(x)δ(x) = δ(x) is understood. 7 All expressions for the J (2) ij as well as for the I (2) ij , which appear in eq. (2.11), are available in electronic form upon request to the authors.

Results for the generalized threshold beam function
The beam function for generalized threshold factorization,B i (t, x), is introduced and defined in ref. [21]. Fort Λ 2 QCD it can be written as the convolution of perturbative matching coefficientsĨ ij (t, z) and the PDFs. The gluon matching coefficientsĨ gj (t, z) can be computed at two loops from the dBF results above via eq. (1.1). The two-loop structure ofĨ ij (t, z) is analogous to the one of I ij (t, z) [22,23] and given in appendix E of ref. [21]. It is completely fixed by the RGE and (known) lower-order ingredients, apart from the term proportional to δ(t): 9 ij (t, z) + terms fixed by RGE and lower-order pieces .
(4.1) 8 Not to be confused with the proper (regulated) splitting function Pgg(z) in eq. (B.3). 9 Note that our definition ofĨ (2) ij (t, z) is in analogy to I ij (t, z) of refs. [22,23] and differs by a conventional factor of 4 compared to the corresponding quantity in ref. [21].

JHEP07(2020)234
All terms in eq. (2.11), when integrated according to eq. (1.1), contribute toĨ (2) ij (t, z), except those terms proportional to L n (t/µ 2 )δ( k 2 ⊥ ). Summing all contributions, one obtains The quantity I (2) ij (z) is the same one that appears in eq. (2.11). The term ∆Ĩ (2) ij (z) corresponds to the contribution from integrating the terms ∝ L 0 (t/µ 2 )J ij and ∝ L 1 (t/µ 2 )P ik ⊗ P kj in eq. (2.11). We computed the ∆Ĩ (2) gj (z) using our results from section 3 and appendix B.2. The ∆Ĩ (2) gj (z) are regular functions in z, i.e. they do not involve distributions. We give the explicit expressions, expressed in terms of ordinary polylogarithms up to weight three, in appendix C. We emphasize that the linearly polarized term ∝ B g in eq. (2.2) vanishes upon the k ⊥ integration in eq. (1.1). Our two-loop results for the gluon generalized threshold beam function are therefore complete.

Conclusions
The fully-differential beam function B j (t, x, k 2 ⊥ , µ) has many potential applications in highenergy collider physics, from resummations of multi-differential observables, to IR subtraction schemes and improved parton showers. In this paper we utilised the methods we developed in refs. [7,22,23] to compute the two-loop matching coefficients I (2) gj (t, z, k 2 ⊥ , µ) between the unpolarised gluon dBF B g (t, x, k 2 ⊥ , µ) and the usual PDFs f j (x, µ). Our results are an important ingredient to obtain the full NNLO singular contributions as well as the NNLL and N 3 LL resummation for observables that probe both the virtuality and the transverse momentum of the colliding gluons. From our two-loop dBF results, we also obtained, for the gluon case, the complete two-loop matching coefficients between the beam function appearing in generalized threshold factorisation [21] and the PDFs.

Acknowledgments
JRG acknowledges the University of Freiburg for hospitality, where a part of this work was completed.

A Tree-level and one-loop matching coefficients
We expand the beam function matching coefficient as At tree level we have

Splitting functions
We expand the PDF anomalous dimensions (γ f ij = 2P ij ) in the MS scheme as The one-loop terms read with the usual one-loop (LO) quark and gluon splitting functions 10 Here Iij(z) ≡ I (1,δ) ij (z) in the notation of refs. [2,36].
Here we employed the definitions in eqs. (3.6) and (B.3). Note that the two terms in eq. (B.7) combine to an expression regular in the limit r → 0.