The Quark Beam Function at Two Loops

In differential measurements at a hadron collider, collinear initial-state radiation is described by process-independent beam functions. They are the field-theoretic analog of initial-state parton showers. Depending on the measured observable they are differential in the virtuality and/or transverse momentum of the colliding partons in addition to the usual longitudinal momentum fraction. Perturbatively, the beam functions can be calculated by matching them onto standard quark and gluon parton distribution functions. We calculate the inclusive virtuality-dependent quark beam function at NNLO, which is relevant for any observables probing the virtuality of the incoming partons, including N-jettiness and beam thrust. For such observables, our results are an important ingredient in the resummation of large logarithms at N3LL order, and provide all contributions enhanced by collinear t-channel singularities at NNLO for quark-initiated processes in analytic form. We perform the calculation in both Feynman and axial gauge and use two different methods to evaluate the discontinuity of the two-loop Feynman diagrams, providing nontrivial checks of the calculation. As part of our results we reproduce the known two-loop QCD splitting functions and confirm at two loops that the virtuality-dependent beam and final-state jet functions have the same anomalous dimension.


Introduction
At the LHC, we are mainly interested in exploring the high-energy, short-distance processes that produce Higgs particles and possible beyond Standard Model particles. The complication we face in studying these processes is that they occur inside a QCD environment of initial-state radiation (ISR), soft interactions, final-state radiation (FSR), and multiparton interactions due to the the fact that we are colliding bags of colored particles, commonly known as protons. This is illustrated in figure 1.
If we are interested in sufficiently inclusive measurements, for example the total cross section for pp → L + X, where L is the colorless final state of interest (e.g. W , Z, or H) and we make no restrictions on the remaining hadronic final state X, then a factorization theorem exists to express the cross section in terms of the partonic short-distance cross section and the parton distribution functions (PDFs) up to corrections that are suppressed by Λ 2 QCD /Q 2 , with Q 2 the scale of the hard interaction. In this case, the PDFs describe all of the initial-state radiation effects, and have only two arguments -the fraction x of f Figure 1. Schematic depiction of a generic LHC collision. An energetic parton from each proton emits initial-state radiation (I) before the two collide in a hard process (H), producing here a lepton pair and colored particles that emit final state radiation in jets (J i ). In addition, there are soft interactions connecting the initial-state and final-state colored particles (S). Also further pairs of initial-state partons may interact (not shown). Figure taken from ref. [1].
light-cone momentum of the proton carried by the parton entering the hard process, and the factorization scale µ [2][3][4]: where i, j = {g, u,ū, d, ...}. The PDFs absorb nonperturbative as well as perturbative effects and so need to be fitted to data, while their evolution with respect to the scale µ can be predicted according to the DGLAP equation.
If one wishes to make more differential measurements, and if the scale of the extra measurement is substantially lower than the scale of the hard process, then eq. (1.1) no longer applies and we either have to derive a new factorization formula, or it is possible that no factorization formula exists. If a new factorization formula can be derived, it will have a more complex structure than eq. (1.1), involving separate functions to describe the effects of the soft radiation, energetic FSR and ISR. In particular, the most appropriate description of energetic ISR may involve parton density objects that depend on more variables than the standard PDFs, referred to as beam functions [5] (they are also referred to as unintegrated PDFs in the literature). The factorization formula differential in some infrared (IR) sensitive observable T has the generic form where ⊗ denotes a convolution of some sort, and there is a correspondence between the objects on the right hand side of eq. (1.2) and the different parts of figure 1: H contains the hard process, while the soft function S, the jet functions J i , and the beam functions B a,b describe the contributions to the measurement of T from soft radiation, energetic FSR, and energetic ISR respectively. A variety of beam functions have appeared in the literature and have been applied in factorization formulae of various kinds. An extensively studied case [6][7][8][9][10][11][12][13][14] is that of the inclusive p T -dependent beam function B i (p 2 T , x, µ), or transverse momentum dependent PDF (TMD PDF), which measures the total transverse momentum p T of the parton i entering the hard process in addition to its light-cone momentum fraction x. It appears for example in the cross section for pp → L + X with L uncharged under color in which the transverse momentum of L is measured. We are interested in the case of the inclusive virtuality (t)-dependent beam function B i (t, x, µ) [1,5]. Here, the argument t measures the component of momentum of the parton i in the light-cone direction opposite to the respective proton momentum, or equivalently the (space-like, transverse) virtuality of the parton i entering the hard interaction. The virtuality-dependent beam functions appear in any cross section in which the incoming parton virtualities are probed by the measurement performed on the final state. They have appeared in a variety of processes -e.g. in J/ψ photoproduction [15], hadron-hadron N -jettiness [16,17], including the Higgs and Drell-Yan 0-jettiness (or beam thrust) [5,18,19] and the Higgs + 1-jet cross section [20], the 1-jettiness DIS cross section [21,22], and 1-jettiness in nuclear dynamics [23,24]. There are also more complex beam functions, that depend on both p T and t [21,[25][26][27][28][29][30], or that involve more complicated jet-algorithm dependent measurements [31][32][33][34][35][36].
The standard PDFs can be defined as proton matrix elements of renormalized operators Q i composed of partonic fields, where |p n (P − ) denotes the external proton state with lightlike momentum P µ = P − n µ /2. The matrix elements are always averaged over proton spins, which we suppress in our notation. Similarly, a generic beam function B(x, T , µ) that depends on the extra dimensionful quantity T (where T = p 2 T , t, . . .) can be expressed in soft-collinear effective theory (SCET) [37][38][39][40][41][42] as the proton matrix elements of renormalized operators composed of partonic fields with an additional dependence on T , (1.4) Like the PDFs, the B i contain nonperturbative information and cannot be predicted entirely from theory. We could fit them from data just like the PDFs -however, this would be a very difficult undertaking due to the paucity of appropriate data and the fact that they depend on a further argument T . However, as long as we are in a region in which T is large enough, i.e. T Λ 2 QCD , we can exploit the fact that in this region T is predominantly generated perturbatively (up to corrections suppressed by powers of 1/T ). In this region, the beam functions can be calculated as the convolution of the usual collinear PDFs, which provide the nonperturbative information, and a perturbatively calculable matching coefficient I ij (T , x, µ) that describes the T dependence due to perturbative ISR: This equation can be regarded as an operator product expansion (OPE) in SCET [5,15]. In almost all of the above-mentioned examples, the matching coefficients are known at next-to-leading order (NLO) in perturbation theory. For the inclusive p T -dependent beam function or TMD PDF, the coefficients at next-to-next-to-leading order (NNLO) have been obtained in ref. [43] using the NNLO calculations in refs. [44,45] (except for gluon-spincorrelated contributions). The quark-to-quark I qq (p 2 T , x, µ) matching coefficient has been directly calculated at NNLO using SCET in ref. [46], and the remaining partonic channels are being calculated [47].
We are interested in the inclusive virtuality-dependent beam function, which for brevity in the following we refer to simply as the beam function. It has been calculated at NLO in refs. [1,19]. In this paper, we perform the two-loop calculation for the quark beam function in SCET, i.e. we calculate the matching coefficients for the quark and antiquark beam functions, I qj (t, x, µ) and Iq j (t, x, µ), at NNLO. They are relevant for processes where the hard interaction is initiated by quarks. The gluon beam function at NNLO will be presented in ref. [48]. The NNLO beam functions are a necessary ingredient for the N 3 LL prediction of the N -jettiness cross section and all other observables involving the virtuality-dependent beam functions. From a fixed-order perspective, the NNLO beam function describes the full two-loop singular contributions differential in t and x from collinear initial-state t-channel singularities.
This paper is organized as follows. In section 2 we discuss the general setup of the matching calculation, including the relevant operator definitions and renormalization equations. For the calculation itself, we utilize two different methods, which are described in section 3. Our results for the NNLO quark matching coefficients I qj are presented in section 4. We conclude in section 5.

General setup and matching at NNLO
Our calculation of the (virtuality-dependent) quark beam function follows the general setup in ref. [1]. For a detailed discussion we refer the reader there. Here, we summarize the relevant definitions of the bare and renormalized beam function and PDF operators in section 2.1, and in section 2.2 we explain the necessary steps to extract the matching coefficients I ij and give all relevant expressions at NNLO.

Operator definitions of beam functions and PDFs
We use the usual light-cone (Sudakov) decomposition for four-vectors, writing an arbitrary four-vector A µ as with n 2 =n 2 = 0 , n ·n = 2 , Quarks and gluons with momentum p are n-collinear if their momentum components scale as (p + , p − , p ⊥ ) ∼ p − (λ 2 , 1, λ), where λ 1 is the power expansion parameter of SCET. The effective size of λ is set by the measurement of interest, e.g., λ 2 T N /Q for N -jettiness [16].
In SCET, n-collinear quarks and gluons are described by composite quark and gluon field operators χ n and B µ n⊥ , where ξ n is the n-collinear quark field and iD µ n⊥ = P µ n⊥ + gA µ n⊥ is the covariant derivative involving the n-collinear gluon field A µ n . The Wilson line accounts for arbitrary emissions of n-collinear gluons, which are all of O(λ 0 ) in the power counting. The P n and P n⊥ are the SCET label momentum operators [39]. The χ n and B µ n⊥ fields are gauge invariant with respect to collinear gauge transformations [38,39]. They are defined after the field redefinition [40] decoupling soft gluons from collinear particles and do not interact with soft gluons through their Lagrangian (at leading order in the power counting). This means they do not transform under soft gauge transformation, such that operators built from them are gauge invariant under both soft and collinear gauge transformations. The soft interactions with collinear particles are factorized into a product of soft Wilson lines, whose matrix element gives the soft function in eq. (1.2).
The bare quark, antiquark, and gluon beam function operators are defined in terms of the fields in eq. (2.3) as [ Here, P n acts within the square brackets and returns the sum of the large minus momentum components of all fields in χ n and B n⊥ . Hence, the effect of the δ(ω − P n ) operator is to set the total minus momentum of the composite quark/gluon field to ω. Thep + is the momentum operator of the small plus momentum, and the effect of δ(t − ωp + ) is to set the total plus momentum of all initial-state radiation (i.e. of any intermediate state inserted between the fields) to t/ω. The bare quark, antiquark, and gluon PDF operators are defined in terms of collinear SCET fields as In contrast to the operators in eq. (2.5), they only measure the large minus components of the fields. The beam functions and PDFs are defined as in Eqs. (1.4) and (1.3) as the proton matrix elements of the corresponding renormalized operators O i and Q i , which are related to the bare operators by Here and in the following we never sum over repeated parton indices unless explicitly stated otherwise. The renormalization constants Z i B (t − t , µ) and Z f ij (ω/ω , µ) are defined to remove the ultraviolet (UV) divergences in the bare beam function and PDF matrix elements, respectively. In this paper, we always renormalize using the MS scheme with dimensional regularization (which yields the standard renormalized PDFs). The renormalization of the PDF operators is well known. In contrast to the PDFs, the renormalization of the beam function operator in SCET depends on t but not on ω. This structure was proven to all orders in perturbation theory in ref. [1], with the renormalization constant (or equivalently the anomalous dimension) being identical to that of the SCET jet function.
As already mentioned, the collinear fields in the operators include a field redefinition decoupling them from soft interactions. Hence, their matrix elements (at leading order in the power counting) are computed using the purely n-collinear sector of the SCET Lagrangian, which has the same form as the full QCD Lagrangian in a boosted frame. Therefore, as is well known, the matrix elements of boost-invariant collinear operators, such as those of the beam functions and PDFs, can be computed using QCD Feynman rules. This renders the calculation more compact because it avoids having to use the more complex SCET vertices.
As discussed in detail in ref. [1], the definitions of the PDFs in terms of SCET fields are equivalent to the standard operator definitions in QCD. For example, the quark PDF in terms of the full QCD quark field ψ(x) in position space is defined as [49] where we used square brackets to denote the renormalized operator and ϕ(x) comprises a quark field attached to a lightlike Wilson line: The inclusion of the Wilson line renders the product of fields separated along then direction in eq. (2.9) gauge invariant. The SCET fields in eq. (2.6) involve a Fourier transform of ψ in y + , and the Wilson line in eq. (2.10) corresponds to the W n contained in the definitions of χ n and B µ n⊥ . In general, the definition of the n-collinear SCET fields explicitly excludes the soft region where all their momentum components become small. In practice, this is implemented by carrying out all integrations in the calculation of collinear SCET matrix elements over the full momentum range (including the soft region) and then subtracting the zero-bin contributions where the collinear momenta become soft [50]. As is well known, the soft region does not contribute to the PDFs, which means that no zero-bin subtractions are needed for the collinear SCET fields appearing in the PDF operators. Therefore, as long as one uses the same renormalization scheme, the QCD and SCET definitions of the PDF are precisely equivalent at both the calculational and formal level.
The situation for the beam function operators is more complicated. First, the explicit dependence on the smallp + momentum corresponds to a large separation in y − in position space (as can be seen from the corresponding definitions in ref. [1]). Thus the quark beam function operator expressed in terms of full QCD fields should contain a quark and an antiquark field operator separated in both the n andn directions, Fourier-transformed with respect to both of these separations. Since the fields have to be separated in both light-cone directions, it is a priori not clear how to obtain an unambiguous gauge-invariant definition for the beam function in terms of full QCD fields, because different paths for the Wilson lines connecting the fields are not equivalent. Second, the zero-bin subtractions on the collinear SCET fields are required in the case of the beam function operators. Practically, they involve scaleless integrals and vanish. Formally, the scaleless zero-bin subtraction is still crucial as it converts 1/ soft IR divergences into 1/ UV divergences that can be renormalized. Nevertheless it is interesting to note that performing the calculation of the beam function matrix elements in SCET using QCD Feynman rules (with a vanishing zerobin subtraction) is operationally equivalent to evaluating the matrix elements of the QCD operator

Calculation of the quark matching coefficients
The explicit form of the matching equation for the beam functions is [1] Our aim is to calculate the quark matching coefficients I qj in eq. (2.12) at O(α 2 s ). Since the I ij only depend on short-distance physics they can be calculated in perturbation theory at the scale µ 2 ∼ t. Furthermore, the matching holds at the operator level in the sense that it is independent of the external state, so to perform the matching calculation we can choose whatever external state is most convenient, provided that state has some overlap with the parton state j [51]. We choose to use on-shell massless quarks and gluons that are n-collinear with momentum p µ = p − n µ /2 as the external states. The corresponding partonic beam functions and PDFs with parton j in the external state are denoted as and analogously for the bare versions. Here, z = ω/p − is now defined as the fraction of parton j's minus momentum carried by parton i. The matrix elements are always understood to be averaged over the color and spin of parton j.
In terms of the partonic matrix elements, the matching equation becomes where in the second line we defined the shorthand notation ⊗ z to denote the convolution in z and the integration limits are implicit in the support of the functions, z < 1 and z/z < 1. The partonic matrix elements B i/j and f i/j in eq. (2.14) are renormalized quantities. However, unlike the proton matrix elements B i and f i , they are unphysical and only serve as a tool to perform the matching. Using on-shell external states makes the calculation simpler algorithmically because there are fewer scales involved. This comes at the cost of introducing explicit IR divergences in both B i/j and f i/j due to which they are ill-defined in four dimensions even after UV renormalization without further IR regularization. We use dimensional regularization in d = 4 − 2 dimensions to also regulate the IR divergences alongside the UV divergences. Therefore, we need to perform the entire calculation in d = 4 − 2 dimensions, right up until the point at which we obtain the matching coefficient I ik [see eq. (2.20) below]. We then set d = 4 to obtain the final result for I ik .
Our two-loop calculation (see section 3) yields the O(α 2 s ) piece of the bare partonic quark beam function B bare q/j . The bare and renormalized partonic matrix elements are related through the renormalization constants Z j B in eq. (2.7), Order by order in the perturbative expansion, Z i B (t, µ, g) can be derived from the known beam function anomalous dimension γ i B (t, µ) [1]. Here, we have explicitly denoted the dependence on the bare and renormalized coupling, g 0 and g, for which we need the oneloop relation We are interested in the part of eq. (2.15) proportional to α 2 s . Let us define the perturbative expansions of B i/j (either bare or renormalized) and Z i B as Then the term in B i/j we are going to compute is B (2) i/j , and eq. (2.15) gives Having obtained B q/j from eq. (2.18), we can extract the two loop I qj by expanding eq. (2.14) to second order in α s . Let us define the expansions of f i/j and I ij as Then, the part of eq. (2.14) proportional to α 2 s yields the matching relation All individual terms on the right-hand side are IR divergent, and as mentioned earlier, have to be consistently evaluated in d dimensions. In particular, in the last term we need the one-loop contribution to I ij in d = 4 − 2 dimensions, which can be straightforwardly obtained from Appendix C of ref. [1]. The IR divergences cancel between the terms such that we obtain an IR-finite result for I (2) ij (z, t, µ), and we can take the limit → 0. The renormalized PDF matrix elements f i/j (z, µ) are related to the bare ones by [see eq. (2.8)] In pure dimensional regularization, all loop corrections to the bare partonic PDF matrix elements are scaleless and vanish. Hence, the MS renormalized f (n) i/j (for n ≥ 1) are given by a pure counterterm contribution. In particular f (1) i/j and f (2) i/j , which are needed in eq. (2.20), are expressed in terms of the well-known one-and two-loop splitting functions P (0) ij and P (1) ij as Note that since f i/j only contains 1/ n -poles and I ij contains none, knowing B (2) q/j and the one-loop quantities in eq. (2.20) allows us to extract both I qj and f (2) q/j , the latter of which allows us to calculate the two-loop splitting functions P (1) qj via eq. (2.23). Hence, from our calculation we get an independent determination of P (1) qj "for free", which should of course agree with the known results [52][53][54], serving as a very useful cross check of our calculation. Formally, the fact that the beam function calculation reproduces the complete set of IR divergences in Eqs. (2.22) and (2.23), such that it yields an IR-finite result for I ij , shows that SCET correctly reproduces the IR structure of QCD at two loops as it should. Since we use dimensional regularization for both UV and IR, this statement relies on knowing the UV renormalization of the beam function. Alternatively, we can take it for granted that SCET reproduces the IR structure of QCD, in which case our calculation provides an explicit confirmation at two loops that the beam function renormalization is identical to that of the jet function. (Having both checks at the same time would require using different regulators in the UV and IR in order to separate UV and IR divergences as was done at one loop in ref. [1].) The µ-dependent logarithmic terms in I ij (t, z, µ) can be obtained from the known two-loop renormalization group equation (RGE), as has been done for the j = g case in ref. [19]. The procedure is as follows. We start from the all-order RGEs for the beam functions, 24) and the PDFs,

(2.27)
Solving this equation recursively to two-loop order we obtain where we denote the plus distributions as All terms on the right hand side of eq. (2.28) are known apart from the two-loop matching functions I ij (z), which are the novel output of our calculation performed in this paper. We present our results for the quark I

Method of calculation
As explained in ref. [1], the bare partonic beam function matrix elements may be calculated by computing the time-ordered partonic matrix elements We evaluate the diagrams together with taking the discontinuity in eq. (3.2) using two different methods, which we shall refer to as "On-Shell Diagram Method" and "Dispersive Method". They are described in sections 3.1 and 3.2 below. When performing the calculation using the On-Shell Diagram Method, we use axial gaugen · A (n) = 0, while we use Feynman gauge in the calculation with the Dispersive Method. Both calculations yield the Figure 2. Diagrams contributing to the calculation of the NNLO matching coefficients I qiqj (a-g), I qiqj (h) and I qig (i-o) using dimensional regularization. Left-right mirror graphs (up to the fermion flow) are not displayed. We also do not display the diagram analogous to (g) but with the fermion flow of the external quarks reversed, which contributes to the flavor-singlet part of I qiqj . The blob in diagram (f) represents the full one-loop gluon self-energy. The graphs can either be interpreted as standard QCD diagrams or as SCET diagrams with collinear quark and gluon lines. Using axial gauge, this set of nontrivial diagrams is complete when using QCD Feynman rules, while it has to be supplemented by diagrams involving vertices of four collinear particles when using SCET Feynman rules. In Feynman gauge, additional diagrams with Wilson line connections as shown in figure 3 contribute.
same result providing us with a strong cross check. 2 We emphasize that the gauge choice is independent of the method and any gauge could have been used with either method. In both methods the calculation of B bare,(2) q/q , which is singular in the limit z → 1 (just like P qq ), is divided into two stages. First, B bare,(2) q/q (z) is calculated for z < 1, and then at the endpoint for z → 1. The advantage of this is that in the z < 1 calculation we can avoid one light-cone divergence, thereby avoiding having to expand master integrals to an extra order in . In the endpoint calculation, we can take the limit z → 1 in an appropriate fashion, which simplifies master integrals and allows us to more easily extend them to the extra order in required in this limit. Another possible way to calculate B bare,(2) q/q (z) near the endpoint is to replace the quark lines by eikonal (Wilson) lines in direction n (following the direction of the quark), which then allows one to use web techniques [55,56] to calculate the graphs. In this way we checked the Abelian C 2 F endpoint piece of the qq beam function, which can be calculated entirely from the one-loop result given in ref. [1] using the web technique.
The fact that the endpoint contribution of B bare,(2) q/q can be calculated by replacing the quark line by an eikonal line can be used to make some definite statements about the form of the endpoint contributions, even without doing any explicit calculations. The eikonal diagrams for the endpoint involve a Wilson line in the n direction connected to a Wilson line in then direction on both sides of the cut, with all possible gluon connections between these. Four-momentum conservation requires that the summed momenta of the emitted gluons crossing the cut must have an n component equal to (1 − z)P − . Then component is fixed by the measurement, cf. eq. (2.5), to be t/(zP − ) → t/P − for z → 1. Due to the symmetry between n andn in the endpoint diagrams the result must be symmetric in (1 − z)P − and t/P − . Because t is the only dimensionful Lorentz-invariant quantity involved, we can argue on dimensional grounds that B bare,(2) q/q must be proportional to t −1−2 , since t has mass dimension 2 and at two loops we have a factor µ 4 . Therefore, the overall endpoint contribution must be equal to ((1 − z)t) −1−2 multiplied only by a function of . This result is useful to check our brute-force calculation of the endpoint contributions, which in individual diagrams gives terms like (1−z) −1− t −1−2 F ( ) that must cancel between diagrams in the final result. Alternatively, it could be used to predict/check the δ(t)L n (1 − z) for n = 0 terms in the beam function from the δ(1 − z)L n (t) terms for n = 0, whose form is in turn dictated by the anomalous dimension γ B (see also ref. [57]).
Just as in the one-loop case [1], the zero-bin subtractions in this two-loop calculation are scaleless and vanish. Their effect is to convert some 1/ n divergences in B bare,(2) q/j from IR divergences to UV divergences, which are then renormalized by the counterterm Z q B and contribute to the anomalous dimension γ q B .

On-Shell Diagram Method
This method employs the Cutkosky rules [58,59] for computing the discontinuities, and closely mirrors the method detailed in the paper by Ellis and Vogelsang [54]. That is, each possible cut of the Feynman diagrams in figure 2 is evaluated separately with the particles crossing the cut being put on shell from the very beginning. The integral over the phase space of the cut (real) partons as well as any integral over a virtual uncut loop momentum is performed for each diagram and each allowed cut.
The nonzero diagrams to be evaluated using this method can be divided into two classes: real-real diagrams in which two lines are cut and there are no loops either side of the cut, and real-virtual diagrams in which only one line is cut and there is a loop on one side of the cut. The real-virtual diagrams are computed with the aid of the integrals given in Appendix A of ref. [54], except converted to use dimensional regularization to regulate the IR divergences, and extended to one higher order in . For the virtual integral with a light-cone singularity in the center of the minus momentum integration region, corresponding to eq. (A.15) in ref. [54], it is necessary to use a further regulator during the calculation, e.g. principal value, or modifying the power of the light-cone divergence. The dependence on the regulator drops out in the final result.
The calculation of the real-real diagrams involves first using the on-shell constraints to perform the integrals over plus momenta. Then the integrals over transverse momenta are done, with the delta function for t being used to fix the magnitude of one transverse momentum. Finally, the integrals over minus momenta are done, with the other measurement delta function being used to fix one minus momentum. Just as in ref. [54], a suitable change of transverse variables is performed to facilitate the integration over transverse momenta.

Dispersive Method
This method is a direct extension of the approach followed in ref. [1] from one to two loops. Here, the matrix element in eq. (3.1) is obtained using standard perturbation theory in terms of loop diagrams like the ones shown in figures 2 and 3 without cuts. The discontinuity is then taken after performing the integrations over most of the components of the loop momenta.
More precisely, we proceed as follows. Letting k and l be the two loop momenta, two light-cone components, say l + and l − , can be directly fixed by the (measurement) delta functions in the beam function operator, eq. (2.5). The integration over the unconstrained k + momentum component we perform using residues. Nonzero contributions require poles in both halves of the complex k + plane, which consequently restricts the integration range for the k − momentum component to a finite interval. After partial fraction decomposition of the resulting amplitude expressions we identify a set of master integrals. For the master integrals we introduce Feynman parameters to combine the propagator denominators and carry out the Euclidean transverse momentum integrations (k ⊥ , l ⊥ ) in d − 2 dimensions. We now take the discontinuity according to eq. (B.5). After that we perform the integration over the Feynman parameter(s) taking into account possible constraints on the integration range due to the θ-function in eq. (B.5). Finally we integrate the amplitudes over k − and express the results as a function of z = ω/p − = l − /p − . Depending on the complexity of the master integral, it is convenient to expand the integrand in (to high enough order) before or after the Feynman parameter integration. In either case, we use the distributional identity in eq. (B.2) to consistently treat the 1/ n poles induced by the Feynman parameter and/or k − integrations.
After the expansion the k − integration is finite for most of the two-loop diagrams. Feynman gauge diagrams with a triple gluon vertex and connections to Wilson lines, like the one in figure 3 a, however, exhibit light-cone singularities within the k − integration region. As in the On-Shell Diagram method these divergences require further regularization. Since in the Dispersive Method the analogy to the two-loop PDF calculation in ref. [54] is somewhat obscured and the distinction between real-real and real-virtual contributions is often not clear we find it most straightforward to modify the power of the light-cone divergence in order to regulate it. Namely, we chose the η-regulator proposed in ref. [8], which assigns a noninteger 1 + η power to each Wilson line propagator. After the k − integration, the divergent 1/η terms cancel for each two-loop diagram individually and we can safely take the η → 0 limit to obtain a regulator-independent result.
In is similar for both methods. For more details on technical issues in the two-loop beam function calculation we refer to an upcoming publication [48].

Results
Here we present our new results for the two-loop quark matching functions I (2) ij (z) entering in eq. (2.28) for i = q and i =q. Since QCD is charge conjugation invariant we can easily obtain the antiquark coefficients I (2) qj (z) from our results for I (2) qj (z). We decompose the two-loop coefficients as follows where q i (q i ) denotes the (anti)quark of flavor i and find the following results: For simplicity we have suppressed the overall θ(1 − z) multiplying the regular terms. To write the above results in a compact form, we have defined the auxiliary functions, which all vanish for z → 1 at least like 1 − z.
Following the discussion in section 2.2 we can extract the two-loop quark splitting functions P (1) qi from our calculation, and we find agreement with the results of refs. [52][53][54] using both axial and Feynman gauge.
Our results are also the last important ingredient to obtain the beam function at N 3 LL order, which accounts for all collinear ISR effects at N 3 LL. The all-order structure of the beam function RGE has been discussed in ref. [1]. In addition to our two-loop matching results, the other required ingredients at this order are the three-loop noncusp anomalous dimension, which is known from the analysis in ref. [1], and the four-loop cusp anomalous dimension. The latter is not yet known, but can be expected to have an almost negligible numerical impact (similar to what has been observed in thrust [60,61]).

Conclusions
We have calculated at two-loop order the perturbative matching coefficients I qj (t, z, µ) between the virtuality-dependent quark beam function B q (t, z, µ) and the PDFs f j (z, µ). We have performed the calculation using two different methods and in two different gauges -covariant Feynman and axial gauge -with both methods and gauges yielding the same result. The two methods differ in their procedure for taking the discontinuities of the operator diagrams that are required to obtain the partonic beam function matrix elements: in the first method the discontinuity is taken immediately using the Cutkosky rules following refs. [52,54] whilst in the second the discontinuity is taken after most of the loop integrals have been performed, following ref. [1]. The calculational effort to determine the NNLO matching coefficients is similar for both methods.
Our calculation provides an explicit verification at two loops of the all-orders result [1] that the beam and jet function anomalous dimensions are equal. Conversely, relying on this fact, we are able to extract the two-loop quark splitting functions, P qi , and find agreement with the well-known results [53,54]. Our results are an important ingredient to obtain the NNLO singular contributions as well as the N 3 LL resummation for observables that probe the virtuality of the colliding partons, such as N -jettiness.

Acknowledgments
We like to thank Wouter Waalewijn, Jonathan Walsh, and Iain Stewart for helpful discussions. The Feynman diagrams in this paper have been drawn using JaxoDraw [62]. This work was supported by the DFG Emmy-Noether Grant No. TA 867/1-1.

A.1 Anomalous dimensions
We define the expansion of the beam function noncusp anomalous dimension as For the quark beam function in MS we have [1] γ q B 0 = 6C F , The coefficients of the MS cusp anomalous dimension to three loops are [63,64] Γ q i = C F Γ i , Γ g i = C A Γ i , (for i = 0, 1, 2) , Γ 0 = 4 , T F n f = 4 3 (4 − π 2 )C A + 5β 0 ,