The Fully-Differential Quark Beam Function at NNLO

We present the first calculation of a fully-unintegrated parton distribution (beam function) at next-to-next-to-leading order (NNLO). We obtain the fully-differential beam function for quark-initiated processes by matching it onto standard parton distribution functions (PDFs) at two loops. The fully-differential beam function is a universal ingredient in resummed predictions of observables probing both the virtuality as well as the transverse momentum of the incoming quark in addition to its usual longitudinal momentum fraction. For such double-differential observables our result provides the part of the NNLO singular cross section related to collinear initial-state radiation (ISR), and is important for the resummation of large logarithms through N3LL.


Introduction
Fully-differential beam functions (dBFs) are generalized (unintegrated) PDFs that in addition to the Bjorken-variable x depend on the transverse momentum ( k ⊥ ) relative to the beam axis and the transverse virtuality (k + k − = −t < 0) of the colliding parton.
We assume both of these scales to be perturbative. More precisely, we consider the kinematic situation, where Q 2 t ∼ k 2 ⊥ Λ 2 QCD and Q ∼ k − is the scale associated with the hard partonic process. That is, the ISR forms a jet roughly along the direction of the incoming beam with the jet axis deviating from the beam axis by a small angle | k ⊥ |/k − . The dBFs are independent of the hard process and depend only on the properties of the colliding parton. They describe the effects of the collinear ISR on (double-)differential cross section measurements probing the full four-momentum k µ of the parton that enters the hard interaction. Because the total invariant mass of the ISR jet must be non-negative, k 2 ⊥ is constrained for fixed t as [1] 1 The beam functions we are concerned with in this work can be formally defined as proton matrix elements of operators in soft-collinear effective field theory (SCET) [2][3][4][5][6][7]. The first type of beam function to be defined in this way was the virtuality-dependent beam function [8,9]. This beam function was generalized to include transverse momentum dependence in ref. [10], although in that paper the beam functions are functions of the Fourier-conjugate variable to transverse momentum, i.e. the impact parameter, (and named iBFs). We will however use the momentum-space definition of the quark dBF given in ref. [1]: For t ∼ k 2 ⊥ Λ 2 QCD we then obtain the dBFs by computing the matching functions I ij (t, z, k 2 ⊥ , µ) perturbatively and convolving them with the standard PDFs f j (z, µ). Integrating this perturbative result for the dBF over the full range of k 2 ⊥ given by eq. (1.1) yields the virtuality-dependent beam function B i (t, x, µ): It is however impossible to deduce the perturbative TMD PDF B i (x, k 2 ⊥ , µ) from a simple integration of eq. (1.3) over t [1]. This is because unlike for k 2 ⊥ the t-integral is not constrained by the kinematics, eq. (1.1), and diverges indicating that the t-integration and the regularization of UV (and rapidity) divergences does not commute. The reason is that for unconstrained t we miss the, in this case, dominant contributions from soft modes with momenta ∼ Q(λ, λ, λ), which give rise to the rapidity divergences and are absent in our SCET I setup.
In ref. [1], the full set of matching coefficients I ij was calculated at one loop (correcting an earlier result in ref. [19]). Moreover it was shown, that the renormalization group (RG) evolution of the dBFs is the same as for the virtuality-dependent beam functions, which in turn equals the one of the (virtuality-dependent) jet function [9]. Since the noncusp anomalous dimension of the jet function [9,20] as well as the cusp anomalous dimension [21,22] are known to three loops, the RG running of the dBFs is known through N 3 LL. The only missing piece in the N 3 LL RG resummation kernel is the four-loop correction to the cusp anomalous dimension, which however can be expected to have an almost negligible numerical impact for processes at present colliders, see ref. [23].
On top of that, N 3 LL precision for the full differential cross section also requires the NNLO fixed-order expressions for the relevant beam, soft, hard and jet funtions. It is the aim of the present paper to determine the coefficients I ij for the (anti)quark dBF (i = q,q) at two loops. Besides the two-loop results for the virtuality-dependent beam functions [11,23] and TMD PDFs [24,25] our calculation extends the set of quark beam functions available at NNLO.
Higher order results for the dBFs may eventually help to systematically improve the initial state parton shower of Monte Carlo event generators beyond LL [26][27][28][29][30][31]. Other possible applications are precise predictions of transverse momentum distributions in Drell-Yan-like processes with a veto on hard central jets, where the jet veto is achieved by a cut on a virtuality-sensitive observable [32]. A related factorization formula is discussed in section 4.1 of ref. [1]. 1 Last but not least the quark dBF plays a prominent role for exclusive N -jet production in DIS. In ref. [34], three different (thrust-like) 1-jettiness event shape variables τ a,b,c 1 were defined, see also refs. [35][36][37], and the corresponding factorization formulae were derived. 2 The factorization formulae for τ b,c 1 involve the quark dBF. Our NNLO result for the quark dBF represents the last important ingredient to improve the corresponding resummed two-jet predictions in DIS from NNLL to N 3 LL precision.
The outline of this paper is as follows. In section 2 we sketch our two-loop matching calculation for the quark dBF and point out the main differences to our calculation of the 1 We note however that this particular factorization formula, in which the jet veto is achieved by using a cut on the global beam thrust, is incomplete and should receive leading power corrections from Glauber modes [33]. This can be avoided by a local veto on an exclusive jet-algorithm-based observable [32], where the effects of Glaubers is O(R 2 ) suppressed, with R the jet radius..  NNLO virtuality-dependent beam functions [11,23]. Section 3 contains the novel results for the dBF matching coefficients and section 4 our conclusions.

Calculation
Our calculation of the two-loop quark dBF matching coefficients closely follows our calculation for the virtuality-dependent quark beam function in ref. [11], see also ref. [23] for further details. Due to the large overlap of the calculations, we will restrict ourselves to discussing the important differences rather than going through the whole calculation in detail.
Since QCD is charge conjugation invariant the antiquark matching coefficients can easily be obtained from the quark ones according to Iq q = I qq and Iq g = I qg (q = u, d, s, . . .), so we will only consider the quark coefficients in the remainder of this section. We begin by computing the bare two-loop quark beam function in a partonic state j (j = q,q, g), which is defined by the matrix element of the same bare dBF operator in eq. (1.2), but with the incoming proton replaced by the parton j. We denote this by B bare q/j (t, z, k 2 ⊥ ). In accordance with refs. [9,11,23] we denote the light-cone minusmomentum fraction by z rather than x when using a partonic state j. The kinematic constraint eq. (1.1) obviously also holds on the partonic level and hence for both x and z.
To obtain the B bare q/j (t, z, k 2 ⊥ ) we calculate the discontinuity of two-loop diagrams like the ones shown in figure 1. The complete set of diagrams relevant in axial-gauge and dimensional regularization is given in ref. [11], where now the bilocal operator represented by the two ⊗ symbols also measures k 2 ⊥ according to eq. (1.2). The bare partonic dBF is related to the renormalized one with Z i B the same renormalization factor as for the virtuality-dependent beam function [1]. This relation holds on the operator level, i.e. independently from the state j and hence also for the physical proton dBF B i (t, z, k 2 ⊥ , µ). Finally the matching coefficients I ij may be extracted using the partonic analog of eq. (1.3). The PDFs in the partonic state j are given up to the two-loop order we need in eqs. (2.22) and (2.23) of ref. [11]. Upon integration over k ⊥ , both the bare and the renormalized partonic dBFs as well as our final results for the I ij must yield the respective results for the virtuality-dependent beam function. At each step this serves us as a strong cross check of the independently obtained parts of the two calculations.
We evaluate the Feynman diagrams together with taking their discontinuity (i.e. performing the unitarity cut) using two methods -the "On-Shell Diagram Method" and the "Dispersive Method", which are described in detail in refs. [11,23]. We also use 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. As in the previous calculations, the calculation of the bare beam function is done in two stages -first we compute the beam function away from z = 1, B bare q/j (t, z, k 2 ⊥ ), and then we add the endpoint z → 1 contribution, δ(1 − z)D q/j (t, k 2 ⊥ ). However, in this calculation we do not need to calculate the endpoint contribution explicitly -we can extract it from our previously-calculated bare virtuality-dependent beam function as follows. Since the endpoint contribution is proportional to δ(1 − z), it must also be proportional to δ( k 2 ⊥ ) according to the constraint eq. (1.1) (z playing the role of The integral of the fully differential beam function over k 2 ⊥ must give the virtuality-dependent beam function, also at the bare level, leading to All terms in this equation apart from D q/j (t) are known, so we can use this equation to extract D q/j (t) and therefore the endpoint. Note that in our previous calculation of the bare virtuality-dependent quark beam function, the integrations over the transverse components of the loop momenta and the expansion in the dimensional were performed before the integrations over the loop minus components. This means that we cannot trivially obtain the dBF by taking our previous calculation and undoing the last integration. However, many of the integral results obtained in that calculation could be re-used in the present context and no additional tool was needed to carry out the integrations for the dBF.
One result required to simplify the piece of our bare partonic dBF proportional to δ(t) is the following distributional identity (which holds if the function f is integrable): Note that the θ-function on the left hand side of eq. (2.3) is present in all terms of the bare partonic dBF (that are regular in the argument of this θ). It technically originates from the unitarity cut through the two-loop diagrams, similarly to the θ(1 − z) in the virtuality-dependent beam function calculation, and enforces the constraint analogous to eq. (1.1). Together with the δ(t) the θ-function restricts the integration range for k 2 ⊥ and k 2 ⊥ itself to zero on the left hand side of eq. (2.3). The t → 0 limit of the last three factors on the left hand side therefore gives a δ( k 2 ⊥ ) normalized by the integral on the right hand side of eq. (2.3). The correctness of eq. (2.3) can be verified by integrating both sides of the equation over k 2 ⊥ . The δ(t) piece of the bare partonic dBF as obtained directly from the two-loop calculation outlined above has the form of the left hand side of eq. (2.3). We therefore conclude that it is the same as the δ(t) piece of the bare partonic virtuality-dependent beam function, up to a factor of δ( k 2 ⊥ )/π. It is actually not surprising that one can use the virtuality-dependent beam function to predict both the δ(t) and the δ(1 − z) pieces of the dBF, given the similar role of t and (1 − z) in the constraint eq. (1.1) (for x = z). In fact eq. (2.3) also holds for δ(t) replaced by δ(1 − z) on both sides of the equation.
The renormalization scale (µ) dependent terms in the matching functions I ij (t, z, k 2 ⊥ , µ) are fixed by solving the corresponding renormalization group equation (RGE), where the standard Mellin convolution in the minus-momentum fraction is denoted by ⊗ z and defined in eq. (B.4). The function is the full (virtuality-dependent) beam function anomalous dimension, and P kj (z, µ) is the QCD splitting function. The anomalous dimension γ i B (t, µ) equals the jet function anomalous dimension γ i B (t, µ) = γ i J (t, µ) [9]. The various cusp (Γ i n ) and non-cusp contributions (γ i B n ) are collected up to three loops (n = 2) for the quark case (i = q) in appendix A.1 of ref. [11]. The terms in the perturbative expansion of the splitting function P (n) ij can be found up to NLO (n = 1) in appendix A.3 of ref. [11]. (We also list the LO expression (n = 0) in appendix B.1.) Let us expand the matching coefficient as follows (note the overall 1/π factor compared to our corresponding expansion for the integrated version of I ij in refs. [11,23]): The tree level and one-loop terms, I Solving the RGE, eq. (2.4), iteratively to NNLO yields the master formula for the two-loop matching coefficient: defines the usual plus distributions. All ingredients in eq. (2.7) were explained and given for the quark case (i = q) in ref. [11], except for the functions J ij (t, z, k 2 ⊥ ) and the zconvolutions in the L 1 (t/µ 2 ) term. The results for the latter convolutions are presented for i = q in appendix B.2. The I (2) ij (z) functions here are the same ones appearing in the matching coefficient for the virtuality-dependent beam function. This is because the δ(t) pieces of B bare q/j (t, z, k 2 ⊥ ) and B bare q/j (t, z) are equal (up to a factor of δ( k 2 ⊥ )/π) as shown above and the respective renormalization factor Z q B is the same. The functions J (2) ij (t, z, k 2 ⊥ ) for i = q are the novel outcome of our dBF calculation. Despite being the coefficient of L 0 (t/µ 2 ) we cannot predict the J (2) ij from the RGE, only their integral over k 2 ⊥ . The reason for this is that when we differentiate the L 0 (t/µ 2 ) term in eq. (2.7) with respect to µ, we obtain [using eq. (C.2)]: where in the second line we use eq. (2.3). We see that by comparing this term to the corresponding δ(t) term on the right hand side of the RGE, eq. (2.4), we can only extract the integral of J (2) ij (t, z, k 2 ⊥ ) over k 2 ⊥ . We present results for the J ij (t, z, k 2 ⊥ ) in the case i = q (and i =q) in the next section.
The diagonal components (i = j = q) of our master formula, eq. (2.7), and also the results for the J (2) ij contain terms that might appear ill-defined as z → 1 or t → 0 at first sight. To obtain a meaningful result in these limits one is however forced to perform the integration over r ≡ k 2 ⊥ /t analogous to eq. (2.3) first. This will generate terms that exactly cancel the ones ill-defined in the z → 1 or t → 0 limits leaving only regular contributions and well-defined distributions. As an example consider the ominous terms in eq. (2.7), ], see appendix B.1. In the limit z → 1 (or t → 0) we have and the term in square brackets in eq. (2.10) vanishes. Hence we are free to replace qq (z) multiplying this term without spoiling the behaviour of eq. (2.10) for z → 1.
Using the notation for the plus-distribution with the boundary condition at r max = (1 − z)/z, as defined in appendix B of ref. [38], we could also compactly express the terms in square brackets in eq. (2.10) as where r ≡ k 2 ⊥ /t. For the sake of simplicity we however refrain from introducing another type of plus distribution and only use the one defined in eq. (2.8), which has the boundary condition at 1, for the presentation of our results in the next section.
For simplicity we have suppressed the overall θ t(1−z)/z − k 2 ⊥ factor already mentioned in section 2 (in all terms regular in the argument of this θ-function). The combination of this function, the θ(z) in eq. (3.1) and the support of the PDFs in eq. (1.3) enforces the kinematic constraints (1 − x)/x ≥ r ≡ k 2 ⊥ /t ≥ 0 and 1 ≥ x ≥ 0. Also, we emphasize again, that due to the overall θ-function the proper limits z → 1 and t → 0 of the above results for the J (2) ij require the integration over k 2 ⊥ (or equivalently r).
The expressions for the J (2) ij in eqs. (3.2)-(3.5) as well as for the I (2) ij in eqs. (4.2)-(4.5) of ref. [11] are also available in electronic form upon request to the authors.

Conclusions
In this paper we have presented the first NNLO calculation of a fully-unintegrated parton distribution, namely the (anti)quark dBF. We have computed at two-loop order the matching coefficients I (2) qj (t, z, k 2 ⊥ , µ) between the dBF B q (t, x, k 2 ⊥ , µ) and the PDFs f j (x, µ). We have checked our computation by using two different gauges, Feynman and axial lightcone gauge, and two different methods for taking the discontinuities of the operator diagrams that are required to obtain the partonic dBF matrix elements. Integration of B q (t, x, k 2 ⊥ , µ) over the transverse momentum k ⊥ yields the virtuality-dependent beam function B q (t, 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 quarks.