The Gluon Beam Function at Two Loops

The virtuality-dependent beam function is a universal ingredient in the resummation for observables probing the virtuality of incoming partons, including N-jettiness and beam thrust. We compute the gluon beam function at two-loop order. Together with our previous results for the two-loop quark beam function, this completes the full set of virtuality-dependent beam functions at next-to-next-to-leading order (NNLO). Our results are required to account for all collinear ISR effects to the N-jettiness event shape through N^3LL order. We present numerical results for both the quark and gluon beam functions up to NNLO and N^3LL order. Numerically, the NNLO matching corrections are important. They reduce the residual matching scale dependence in the resummed beam function by about a factor of two.


Introduction
In differential measurements at hadron colliders, collinear initial-state radiation is described and can be resummed by process-independent beam functions [1]. In this paper, we are concerned with the two-loop virtuality-dependent beam functions B i (t, x, µ), which are defined formally as the proton matrix element of operators in soft collinear effective field theory (SCET) [2][3][4][5][6][7] in refs. [1,8].
The beam functions B i (t, x, µ) are an integral component for predictions of observables in hadronic processes that probe the virtuality of the incoming partons via a measurement performed on the hadronic final state, such as N -jettiness [9] and beam thrust [1]. The B i (t, x, µ) are a type of unintegrated parton density that, loosely speaking, gives the probability to find a parton i in the initial state carrying a fraction x of the proton's lightcone momentum, accompanied by initial-state radiation that causes the parton i to be off shell with virtuality −t. 1 The SCET definition of the virtuality-dependent gluon beam function (i = g) relevant for this work reads (x ≡ ω/P − ) B g (t, x, µ) = p n (P − ) (−ω)θ(ω) B c n⊥µ (0) δ(t − ωp + ) δ(ω − P n )B µc n⊥ (0) p n (P − ) , (1.1) where p n (P − ) denotes the incoming (spin-averaged) proton state with lightlike momentum P µ = P − n µ /2, P n is the SCET minus-momentum label operator [4],p + is the plusmomentum operator, and B µ n⊥ denotes the gauge-invariant gluon field strength operator in SCET. For more details on the SCET notations and conventions, we refer to refs. [1,8,13].
For t Λ QCD , the beam function can be obtained as the convolution of the collinear parton density functions (PDFs) f i (x, µ) with perturbatively calculable matching functions I ij (t, z, µ) [1,8,14] B i (t, x, µ) = (1. 2) The matching coefficients were calculated to one-loop order for the quark case (i = q) in ref. [8] and for the gluon case (i = g) in ref. [15]. In a previous publication [13] we calculated the (anti)quark matching coefficients at two-loop order (NNLO). In this paper, we complete the full set of matching coefficients through NNLO.
In section 2, we give a brief summary of the calculational approach, referring to ref. [13] for more details on the setup of the matching calculation. In the appendices, we give several technical details on the calculation. In section 3, we present our results for the gluon matching coefficients at two loops. In section 4, we present numerical results for both the quark and gluon beam functions up to NNLO and N 3 LL order. We conclude in section 5.

Calculation
We calculate the NNLO gluon beam function using the same methodology as in ref. [13]. That is, we compute the bare (unrenormalized) partonic beam function B bare g/j (t, z, g 0 ) from the discontinuity with respect to t of the Feynman diagrams shown in figure 1. From these we then extract the coefficients I gj (t, z, µ) using the matching equations given in section 2.2 of ref. [13]. We evaluate the diagrams together with taking the discontinuity using two different methods, the 'On-Shell Diagram Method' and the 'Dispersive Method', which are explained in ref. [13]. We also perform the calculation using two different gauges, namely light-cone axial (n · A n = 0) gauge and Feynman gauge. The two different methods and gauge choices gave the same final results, hence providing us with a strong cross check.
In the appendices, we provide some further calculational details: in Appendix B, we provide the change of transverse variables in the On-Shell Diagram Method that is used to compute the cuts of the diagrams. In Appendix C, we analyze the three-point integral that contains the only (light-cone) divergence in the calculation that is not regulated by dimensional regularization (related to the discussion in section 3.1 of ref. [13]). Finally in Appendix D, we complete the list of required SCET Feynman rules in Feynman gauge (see ref. [15]) by giving the expression for the triple gluon vertex associated with the gluon field strength operator B µ n⊥ , which is first needed at NNLO, and to our knowledge has not been given in the literature before.
To compute the endpoint (z → 1) contributions to B bare g/g , we deviate slightly from the procedure taken in [13]. In contrast to the quark case, we do not need to calculate . Diagrams contributing to the calculation of the NNLO matching coefficients I gg (a-l) and I gq (m-s) when using dimensional regularization. Left-right mirror graphs and graphs with reversed fermion flow in the loop are not displayed. The blob in diagrams (h,i,p) represents the full one-loop gluon self-energy. The graphs can be computed using either standard QCD Feynman rules or SCET Feynman rules with collinear quark and gluon lines. Using axial gauge and QCD Feynman rules, this is the complete set of nontrivial diagrams. Using SCET Feynman rules, it has to be supplemented by diagrams involving vertices of four collinear particles. In Feynman gauge, additional diagrams with Wilson line connections (see e.g. figure 2) or ghost loops contribute. the endpoint of B bare g/g directly. Rather, we appeal to the argument given in ref. [13], that the endpoint contributions can be obtained by replacing the incoming lines by (collinear) Wilson lines, this time in the adjoint representation rather than the fundamental representation as was appropriate for the quark case. This is illustrated in figure 2 for an example diagram with a nonzero endpoint contribution in Feynman gauge. Then the quark and gluon endpoint diagrams become identical, up to a replacement of fundamental color ma- trices by adjoint color matrices for connections along the Wilson lines. At the two-loop order, this means that the quark and gluon endpoints are actually equal, up to the replacement C F → C A in going from the quark to the gluon endpoint. Therefore, we can obtain the gluon endpoint contribution simply by taking the quark endpoint contribution already calculated and making the replacement C F → C A . 2 We regulate divergences associated with light-cone propagators (in light-cone axial gauge) or Wilson line propagators (in Feynman gauge) by dimensional regularizationexcept for one particular case concerning diagram (c) that is discussed further in Appendix C. Other ways to consistently treat these light-cone divergences, employed in the past for the two-loop (axial gauge) calculations of the QCD splitting functions, are the principal value [16][17][18], or the Mandelstam-Leibbrandt prescription [19][20][21][22]. Using dimensional regularization, we find in light-cone gauge that the swordfish diagrams (f) and (g) are zero, whilst the diagram (d) contributes only a finite part, with no poles. Therefore our prescription gives zero for the contributions of diagrams (d), (f), and (g) to the gluon-gluon splitting function, which is related to the infrared divergent part of the diagrams. In this respect, dimensional regularization behaves similarly to the principal value prescription, which in contrast to the Mandelstam-Leibbrandt prescription [23] was observed to give zero contribution of the diagrams (d), (f), and (g) to the splitting function [18]. In Feynman gauge, the diagrams (d), (f), and (g) contribute both to the poles and the finite pieces of B bare g/g .

Results
We expand the matching coefficient I ij in a perturbative series as At tree level and one loop we have Iteratively solving the renormalization group equation for I ij to O(α 2 s ) yields for the twoloop matching coefficient [13] denotes the usual plus distributions. The new results of our calculation are the two-loop gluon δ(t)-terms I (2) gj (z) in the last line of eq. (3.3). All remaining ingredients in eq. (3.3) are known and, for the case i = g, have been given in ref. [15]. They are collected in Appendix A for completeness.
We write the I gj (z) as and we find and For simplicity we have suppressed the overall θ(1 − z) multiplying the regular terms. The auxiliary functions, all vanish for z → 1 at least like 1 − z and are identical to those for the quark case [13]. As in the quark beam function calculation, one can extract from the poles of the bare beam function either the two-loop anomalous dimension for the gluon beam function γ g B1 (t, µ) , or the two-loop splitting functions P (1) gi (z), assuming the other quantity is known (alternatively one may extract both using the sum rules for the splitting functions [12]). We extracted these functions and found agreement with the known results [8,17,18], which serves as an additional check of our calculation.

Numerics
To illustrate the numerical impact of the NNLO corrections to the beam functions associated with the I (2) ij (z) computed here and in ref. [13], we consider the integrated beam function In all our numerical results, we pick a representative value of t max = (30 GeV) 2 . The qualitative features in the numerical results only depend very little on the value of t max . We use the MSTW 2008 PDFs [24] with their corresponding α s (m Z ).
In figure 3 ∝ ln n+1 (t max /µ 2 B ) from integrating the plus distributions L n (t/µ 2 B ) in eq. (3.3) vanish. Hence, with this scale choice, the n-loop correction to B i is directly given by where For each parton i = u, d,d, g, figure 3 shows the pure one-loop correction, B i / B i / B (0) i (orange) in percent relative to the tree level result, as a function of the minus-momentum fraction x. Since here we care about the size of the terms in the perturbative series of the matching coefficients, we use the same NNLO PDFs everywhere. For each order, we show three curves corresponding to the contributions from the diagonal (dashed lines), the off-diagonal (dotted lines), and the sum of all parton channels (solid lines). The diagonal contributions (q → d, u,d) to B d,u,d include the sum of all possible (anti)quark-(anti)quark channels (q i ,q i → d, u,d). Similarly, the off-diagonal contribution (q → g) to B g includes the sum over all (anti)quark-to-gluon contributions (q i ,q i → g). Numerically, the corrections from the q i → q j channels in B (2) q j with i = j however turn out to be completely negligible.
In all cases in figure 3 we observe sizable negative total O(α 2 s ) corrections compared to the positive total O(α s ) corrections. For the (anti)quark beam functions, the two-loop corrections are about half the size of the one-loop corrections. The reason is that the g → q mixing contribution is sizeable and always negative. As a result, at NLO it partially compensates the diagonal q → q contribution, reducing the absolute size of the total O(α s ) correction. In contrast, at NNLO it adds to the diagonal contribution, enhancing the absolute size of the total O(α 2 s ) correction. For the gluon beam function, the corrections are much larger than for the quark case, as expected from the larger color factor for gluons. Here, the off-diagonal mixing contributions q → g have a relatively smaller effect and add to the diagonal g → g contribution (for most of the relevant x range). For x → 1 the corrections become large due to the presence of threshold logarithms α n s ln m (1 − x) in the diagonal terms. For x 0.2, the corrections are largely independent of x, and their size complies with the typical pattern expected for perturbative QCD corrections at the considered scales.
To study the uncertainties in the perturbative series for the beam function, we cannot use a simple variation of µ B in eq. (4.1), since the fixed-order beam function has an explicit  dependence on µ B , containing Sudakov double logarithms in µ 2 B /t max . Instead, we can consider the resummed beam function, where the evolution factor U i B only depends on i = q or i = g and can be found in refs. [8,15].
The perturbative ingredients entering at a given resummation order are summarized in Table 1. The one-loop matching enters at NLL and NNLL, while the two-loop matching enters at NNLL and N 3 LL. The explicit expressions for the resummation kernels to N 3 LL are taken from ref. [25]. The noncusp anomalous dimensions are known [8,15] from the three-loop results in refs. [26][27][28][29]. The beta function is known up to four loops [30][31][32]. The cusp anomalous dimension is known at present to three loops [26,33]. For its fourloop coefficient formally needed at N 3 LL, we use the Pade approximation Γ 3 = Γ 2 2 /Γ 1 .  In all cases we show the correction in percent relative to the fixed NNLO result at the central scale The numerical effect of varying Γ 3 by a factor of ±3 is much smaller than the effect induced by the known three-loop noncusp coefficient. To study the µ B dependence in the resummed beam function, we use a consistent set of PDFs as required by the matching order, together with their corresponding α s (m Z ) value, see Table 1. The α s running deserves some comment. Strictly speaking, the PDFs require three-loop (two-loop) running at NNLO (NLO), which is the same running order as required by the NNLL (NLL) RGE running. This means, the employed α s running is formally fully consistent between PDFs and RGE at NNLL (NLL ) order, while at N 3 LL (NNLL) order, the resummation requires the α s running at four (three) loops, i.e., one order higher than the PDFs. In these cases, we use the higher α s running, to have a fully consistent resummation, together with the numerical α s (m Z ) value of the PDFs, which is the dominant effect as far as numerical consistency with the PDFs goes. 3 The resummed B i (t max , x, µ) explicitly depends on the arbitrary scale µ but is formally independent of the matching scale µ B , with the µ B dependence canceling between the fixed-order B i (t max , x, µ B ) and the evolution factor U B (µ B , µ) to the order one is working at. Hence, in eq. (4.3), we can use the residual dependence on the matching scale µ B as an indication of the uncertainties due to missing higher-order corrections in the perturbative series of the beam function as long as µ 2 B t max (so there are no large unresummed logarithms in the fixed-order series for B i (t max , x, µ B )). In figure 4 we show the residual matching scale dependence, varying µ B between √ t max /4 and 4 √ t max , at different resummation orders, and for a representative value of x = 0.01. In these plots, we choose µ 2 = t max such that the central value at µ B = √ t max is equivalent to the pure NNLO or NLO result from eq. (4.2). The purpose of the resummation here is thus not to resum large logarithms of µ/µ B , but rather as a means to have a meaningful way to estimate the perturbative uncertainties from the residual matching scale dependence. All lines in these plots are shown as the percent change relative to the NNLO result at the central scale.
The NNLL (NLL) evolution factor cancels the explicit logarithmic µ B dependence in the fixed-order B i (t max , x, µ B ) to NNLO (NLO). Therefore, at NNLL (NLL ), shown by the solid lines in figure 4, the µ B dependence comes from both the residual µ B dependence of the nonlogarithmic NNLO (NLO) matching corrections as well as the higher-order logarithmic corrections resummed by the evolution kernel. (The cancellation of the PDF scale dependence inside B i (t max , x, µ B ) by both the diagonal and off-diagonal matching corrections also plays a nontrivial role.) By going to N 3 LL (NNLL), shown by the dashed lines, one includes an additional uncanceled µ B dependence from the three-loop (two-loop) noncusp anomalous dimension (as well as higher-order β function and cusp pieces).
For the quark beam function, the effect of the two-loop noncusp is large and increases the overall µ B variation at NNLL compared to NLL . At N 3 LL the opposite happens. Here, the three-loop noncusp corrections turn out to be tiny due to an accidental but almost perfect numerical cancellation in the combination γ q B 2 − γ q B 1 β 1 /β 0 that appears in the RGE solution. The effect of the four-loop cusp is tiny, which is not unusual. This overall pattern is consistent with that seen in ref. [25] for thrust to N 3 LL, which has an equivalent resummation structure. For the gluon beam function, the same cancellation does not happen, so the effect of the three-loop noncusp at N 3 LL is visible and consistent with the size of the NNLL effect suppressed by an additional power of α s . (The numerical effect of the four-loop cusp is tiny here as well.) Overall, the scale dependence is larger in the gluon case due to the larger color factor for gluons than quarks.
Including the two-loop matching corrections reduces the matching scale dependence by a factor of two, from 2 − 3% to 1 − 1.5% for quarks and ∼ 8% to ∼ 4% for gluons. In complete resummed cross sections, the perturbative uncertainties due to the beam function component are often evaluated by separately varying the beam function scale, see e.g. refs. [15,[34][35][36][37][38][39][40][41]. This is typically implemented through profile scale variations [25,42], which in the resummation region correspond to the canonical beam scales we have used here. In this context, this can be an important source of uncertainties. For example, in ref. [15] the gluon beam function gives the largest uncertainty in the resummation regime. We emphasize that final conclusions concerning the perturbative convergence and uncertainties can of course only be drawn by looking at the cross section for physical observables. Nevertheless, the overall reduction in the matching scale dependence in the resummed beam function at NNLL and N 3 LL gives a good indication of the possible reduction in uncertainties in resummed predictions.

Conclusions
In this paper, we have completed the calculation of the NNLO virtuality-dependent beam functions B i (t, x, µ), by computing the gluon matching coefficients I gj (t, z, µ) for the gluon beam function onto the PDFs at two loops. These 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 the virtuality of the colliding partons, such as beam thrust and N -jettiness.
The methodology used here is the same as in our previous calculation of the twoloop quark matching coefficients in ref. [13]. As in the quark case, we have checked our calculation by using two different gauges -Feynman and axial light-cone gauge -and two different methods for taking the discontinuities of the operator diagrams that are required to obtain the partonic beam function matrix elements. Our calculation provides an explicit verification at two loops of the all-orders result [8] that the beam and jet function anomalous dimensions are equal also for the gluon case. Conversely, relying on this fact, we are able to extract the two-loop gluon splitting functions, P gi , and find agreement with the well-known results [17,18].
We have also presented numerical results for the (anti)quark and gluon beam functions at NNLO as well as for the resummed beam function to N 3 LL. We find that the numerical effects of the two-loop corrections are important. They are about half the size of the oneloop corrections but with opposite sign (except near x = 1). For the resummed beam function, the residual dependence on the matching scale µ B gives an indication of the perturbative uncertainties due to missing higher order corrections. It reduces by roughly a factor of two when our new two-loop matching corrections are included.

A Perturbative ingredients
In this appendix we summarize the additional perturbative ingredients for the gluon beam function. These have been given previously in ref. [15] and are repeated here for completeness.
The coefficients of the cusp, noncusp, and PDF anomalous dimensions are defined according to The MS anomalous dimension coefficients for the gluon beam function up to three loops are The coefficients of the cusp anomalous dimension and beta function are given in ref. [13]. The one-loop gluon matching coefficients appearing in eq. (3.2) are written as with the one-loop matching functions 4 The one-loop PDF anomalous dimension in the MS scheme are with the LO gluon splitting functions At two loops we write where the NLO gluon splitting functions are given by [17,18] P (1) After the angular integral, we have two terms, one of which corresponds to r 1 < r 2 and the other of which corresponds to r 2 < r 1 . The integrals over one of the magnitudes r 1 or r 2 in these terms can be performed using the delta function of t, whilst the other can be performed using a straightforward variable transform (often only as an expansion in ). Then all that remains are the integrals over the components z 1 and z 2 , which can be performed using standard techniques.
The change of transverse variables in eq. (B.4) is sufficient to evaluate all other nontrivial integrals for the ladder diagram, and indeed the integrals for all other topologies in the two-loop calculation. where y =n · l /n · p and w = −y(1 − z)/z here. The first integral in curly brackets is regulated by the w − factor. The second integral requires further regulation, because we integrate 1/y over the origin, but this is a very simple integral. We can use any one of the standard regulators in the 1/y factor: principal value, multiplying the integrand by an infinitesimal negative power of |y|, adding a small imaginary part of either sign to the denominator 5 (as is done e.g. in ref. [12], section 3.1), cutting out a small symmetric region from the integration either side of the origin, etc.. After finally setting the regulator to zero, we will get the same result for any regulator.

D Triple gluon field strength vertex
The Feynman gauge expressions for the vertices with one and two external legs associated with the gluon field strength operator B µ n⊥ , which is part of the operator definition of B g in eq. (1.1), can e.g. be found in ref. [15]. For completeness we also give the Feynman rule for the B µ n⊥ vertex with three external gluons needed in the Feynman gauge calculation of I − g 2 2n αnβnγ n·p an ·p cn ·(p a + p b )n·(p b + p c )n·(p a + p b + p c ) tr T a T b T c T d (D.1) × p µ a⊥n ·p an ·(p a + p b ) − p µ b⊥n ·(p a + p b )n·(p b + p c ) + p µ c⊥n ·p cn ·(p b + p c ) + perms.
Here 'cycl.' and 'perms.' stand for additional terms generated by cyclic (two more terms) and full permutations (five more terms) in {(p a , α, a), (p b , β, b), (p c , γ, c)}, respectively. 5 In Feynman gauge this corresponds to consistently assigning a ±i prescription to the collinear Wilson line propagators, which in SCET is a priori not fixed by causality. Prescription-dependent terms cancel once the complex conjugated (i.e. the left-right mirror) graph of figure 6 is added.
Note that in our two-loop calculation of the partonic beam function with on-shell transverse polarized incoming gluons the second term ∝n αnβnγ does not contribute.