Quark beam function at next-to-next-to-next-to-leading order in perturbative QCD in the generalized large-$N_c$ approximation

We present the matching coefficient for the quark beam function at next-to-next-to-next-to-leading order in perturbative QCD in the generalized large $N_c$-approximation, $N_c \sim N_f \gg 1$. Although several refinements are still needed to make this result interesting for phenomenological applications, our computation shows that a fully-differential description of simple color singlet production processes at a hadron collider at N$^3$LO in perturbative QCD is within reach.


I. INTRODUCTION
Good understanding of infra-red and collinear limits in perturbative QCD and the ability to use this understanding for an increasingly accurate description of hadron collisions is one of the key elements for the success of the future LHC physics program.Because of that, much of the current effort in theoretical collider physics focuses on achieving and advancing such understanding in a number of complementary ways, ranging from fixed-order computations, to resummations and, finally, to parton showers.Although for each of these approaches there exists a set of observables and theoretical quantities to which it is traditionally applied, there are a few cases which lie at their intersections and where progress achieved in the context of one approach has implications for the other ones.
One such theoretical quantity is the so-called beam function [1,2].Beam functions describe the dynamics of incoming partons that slightly deviate from their original direction by emitting hard quasi-collinear radiation before going into the hard process.For this reason, beam functions are important ingredients for resummation studies that aim to understand differential cross sections in the quasi-collinear region [3][4][5][6][7].
Two collinearity measures have been discussed in the literature -the total transverse momentum of the radiated partons and the 0-jettiness In Eqs. ( 1)-( 2), Q 1,2 are so-called "hardness" variables for the initial state partons (see, e.g., [1,8]), p 1,2 are the momenta of the incoming partons and k 1,...,N are the momenta of on-shell final state partons.
It is quite interesting to extend the computation of the matching coefficients to one order higher in the strong coupling constant α s .Not only will such a computation stress-test many aspects of our understanding of softcollinear dynamics in QCD, as well as many techniques of perturbative quantum field theory, but it will also provide an alternative path to next-to-next-to-next-to-leading order (N 3 LO) QCD description of color-singlet production at the exclusive level.Currently, N 3 LO QCD corrections to the inclusive cross section [22][23][24][25][26], as well as to the Higgs rapidity distribution in Higgs boson production in gluon fusion are available [27].An extension of N 3 LO computations to Drell-Yan-like processes, accounting for decays of Z and W bosons to leptons, is very desirable.
Recently, we have computed loop and phase-space integrals relevant for the so-called triple-real and doublereal single-virtual contributions to the quark-to-quark 0-jettiness matching coefficient, focusing on gluonic final states [28,29].When combined with the computation of the single-real double-virtual splitting function q * → qg described in Ref. [30], all ingredients required to obtain the N 3 LO QCD contribution to the quarkto-quark matching coefficient I qq through leading color become available.In addition, the results reported in [28,29] allow us to compute all N 3 LO contributions that scale as N 2 c N f and N c N 2 f , where N f is the number of massless quarks in the theory.
The goal of this paper is to present the N 3 LO contribution to the quark matching coefficient in the approximation N c ∼ N f 1, keeping only leading We will refer to it as the generalized large-N c or leading-color approximation.
We note that our computation of the matching coefficient I qq is restricted to generalized leading-color approximation since, so far, we have not computed all the required contributions of final states with additional quark pairs that are relevant beyond the generalized large-N c limit.In principle, the required computations are similar to what has already been done in Refs.[28,29] but, due to proliferation of integrals required for multi-quark final states, the calculations have not been finalized.
Nevertheless, we believe that the generalized large-N c N 3 LO contribution to the quark-to-quark matching coefficient is an interesting intermediate result since, at variance with our previous publications [28,29], it explicitly demonstrates how different pieces combine to produce a well-defined physical quantity at next-to-next-to-nextto-leading order in perturbative QCD.It also shows that such high-order computations, in spite of their significant complexity, appear to be doable with current computational technologies.
The rest of the paper is organized as follows.In Section II we describe how the computation of the perturbative matching coefficient is set up.In Section III we discuss how the various required ingredients are obtained.We present the result for the matching coefficient in the generalized large-N c approximation in Section IV and conclude in Section V. A number of useful formulas can be found in the Appendix.

II. PERTURBATIVE MATCHING COEFFICIENT
In this section we explain how the perturbative matching coefficient is computed.The starting point is the relation between beam functions and parton distribution functions where the sign ⊗ z stands for the convolution1 The proportionality coefficients between the beam functions and the parton distribution functions, I ik (t, z, µ) in Eq. (3), are the matching coefficients.The sum in Eq. ( 3) runs over all species of partons that are found in the proton for a particular value of the factorization scale µ.The parameter t is the so-called transverse virtuality, which is related to the 0-jettiness variable T in Eq. (2) and will be defined below in Eq. ( 14).For √ t Λ QCD , the matching coefficient I ik can be calculated in perturbative QCD.To this end, we replace the non-perturbative parton distributions with their perturbative counter-parts, calculate the partonic beam function and extract the matching coefficient by comparing the two sides of Eq. (3).Similar to parton distribution functions, this can be done for any combination of an incoming parton j and the parton i that eventually goes into the hard scattering.We therefore write In contrast to Eq. ( 3), all quantities in Eq. ( 5) admit an expansion in the strong coupling constant α s .Writing and defining the leading-order quantities through , we solve Eq. ( 5) to express the matching coefficients through the partonic beam function.We find Perturbative parton distribution functions in various orders in α s are obtained as iterative solutions of the Altarelli-Parisi equation with the boundary condition given above.We note that since in Eq. ( 3) the parton distribution functions are the MS ones, the perturbative parton distribution functions that we need can only contain poles in the dimensional regularization parameter .Explicit results for f (1,2,3) ij in terms of the splitting functions P ik are given in the Appendix.
Eq. ( 7) allows us to iteratively compute the matching coefficients once the perturbative beam functions become available.However, a beam function computed directly from the quasi-collinear limits of the relevant scattering amplitudes is what one refers to as a bare beam function, because it contains both soft and collinear divergences.Soft divergences must be removed by a dedicated MSsubtraction that, schematically, is given by the following formula In Eq. ( 9) the convolution with respect to t is defined by the equation To compute the quark-to-quark matching coefficient, we require the renormalization constant Z q .Similar to other renormalization constants, Z q satisfies a renormalization group equation [1] where the anomalous dimension reads The anomalous dimensions γ q B and Γ q cusp are known through O(α 3 s ) [2,[32][33][34].Here, L 0 (t/µ 2 ) is the modified plus-distribution L 0 (t/µ 2 ) = µ −2 [µ 2 /t] + with the (regularized) singularity at t/µ 2 = 0 rather than at t/µ 2 = 1.In practice, we construct the renormalization constant Z q in the MS-scheme from Eq. ( 11) by expanding the various quantities in the strong coupling constant, see e.g.Eq. ( 29), and inserting an ansatz for Z q in terms of t-distributions.The ansatz is constructed following an observation that Z q must have the same t-dependence as the bare beam function in order to cancel the soft divergences.We then use Eq. ( 9) to obtain the renormalized partonic beam function from the bare one.Finally, we employ Eq. ( 7) to derive the desired matching coefficient.Explicit formulas for various steps described above are given in the Appendix.
We note that since the partonic PDFs are singular in the → 0 limit, f (n) ij ∼ −n , it follows from Eq. ( 7) that the matching coefficients I (1,2) ij need to be known to higher powers in the dimensional regularization parameter .The relevant computation was performed in Ref. [35] and we borrow the results from there.
It remains to discuss the computation of the bare beam function.We do that in the next section.

III. COMPUTATION OF THE BARE 0-JETTINESS QUARK BEAM FUNCTION
It is clear that the major challenge for computing matching coefficients through third order in perturbative QCD is the calculation of the bare beam functions.We can obtain the bare quark beam function from any physical process that features a quark in the initial state, by extracting the leading collinear-enhanced contributions.Since leading collinear singularities factorize into products of universal splitting functions and hard matrix elements, one can organize the calculation in a processindependent way.Indeed, in physical gauges, collinear splitting functions can be obtained by considering QCD radiation off a single external line [36], for example the incoming quark line in our case.It is important that the emissions, both real and virtual, that originate from any other incoming lines, do not contribute to leading collinear singularities and, for this reason, can be ignored.The splitting functions so obtained must be integrated over the particular phase space for real emission(s) that is constrained in such a way as to keep the momentum fraction z and the transverse virtuality t of the incoming quark that goes into the hard scattering process fixed [37].
The bare quark beam function at N 3 LO is then computed by adding such collinear-enhanced contributions with up to three real partons in the final state, with the number of virtual loops required to provide the O(α 3 s ) correction to the leading-order transition q → q.Hence, we need to consider a tree-level contribution where a quark splits into a virtual quark that goes into a hard process and three real partons, a one-loop correction to a process where a quark splits into a virtual quark and two real partons and a two-loop correction to the q → q * + g splitting.
Since in this paper we focus on the generalized large-N c contribution to the quark beam function, where the number of colors and the number of flavors are taken to be large N c ∼ N f 1, it is sufficient to consider gluons in the final state as well as quarks that exclusively originate from a final-state gluon splitting.Other final states are sub-leading in the generalized large-N c approximation.Fig. 1 illustrates which types of quark-antiquark final states have been included and which types have been excluded from our calculation.
We schematically write the O(α 3 s ) contribution to the bare beam function of a quark in the following way where the label Rn R V n V refers to processes with n R real partons and n V virtual loops.The quantities where p is the four-momentum of the incoming parton, p is the complementary collinear direction, denotes the n V -loop contribution to the collinear splitting functions that describes the q → q * + g 1 + ... + g n R process or, if n R ≥ 2, the q → q * + q + q + g 3 + ... + g n R process.We note that the functions B b,Rn R V n V qq (t, z) scale uniformly with the transverse virtuality, i.e.
This observation will be important for the discussion below where we describe the computation of the doublevirtual single-real contribution B b,R1V2 qq .The calculation of the triple-real and double-real single-virtual contributions B b,R3V 0 qq and B b,R2V 1 qq was discussed in Refs.[28,29], respectively.We will briefly summarize these discussions here.
Although, as we already said, the collinear splitting functions in Eq. ( 14) are universal objects, they are not available in closed form beyond NNLO.Since, as shown in Eq. ( 14), our goal is not only to construct the splitting functions, but also to integrate them over the realemission phase space, it is important to have an algorithm that allows us to perform both of these tasks in a concerted way.We achieve this by following the procedure outlined in Ref. [36] that describes how to extract splitting functions by considering emissions off a single external line and by employing relevant projection operators.An important ingredient in this construction is the use of physical gauges for both virtual and real glu-ons that, unfortunately, complicates the computations significantly.In Ref. [36] this procedure was used to explicitly construct all tree-level splitting functions at NNLO in QCD.Here, we just use this procedure to find a suitable expression for the collinear splitting functions P (Rn R V n V ) qq (p, p, {k i }) that may involve unintegrated momenta of both real and virtual gluons.Once such a representation for P (Rn R V n V ) qq (p, p, {k i }) is available, we apply reverse unitarity [38] to map phase-space integrals onto loop integrals.We then use integration-by-parts technology [39,40] to express each particular contribution to B bare qq in terms of master integrals and to derive the differential equations that these integrals satisfy [41][42][43][44].
A detailed discussion of how the master integrals are computed from the relevant differential equations can be found in Refs.[28,29].Here, we just note that the use of physical gauges makes their computation much more difficult, in that it introduces additional propagator-like structures that arise from polarization sums of real and virtual gluons.Unfortunately, this leads to a proliferation of integrals that need to be calculated.Another interesting point is that the master integrals, that describe triple-real emissions, are initially written as linear combinations of generalized polylogarithms of a complexvalued variable which arises during the rationalization of the differential equations, see Ref. [29].Curiously, as we will see from the final result, the dependence on x disappears once the complete triple-real emission contribution to the beam function is constructed.
In principle, one can compute the B b,R1V 2 qq contribution to the beam function using a similar approach.This would require the calculation of the two-loop correction to the process q → q * + g in a physical gauge; such computation is, currently, not available.Fortunately, there is a way out.The contribution we are interested in can be extracted from the two-loop amplitude of the process q(p)q(p) → V + g(k 1 ) in the limit when the gluon is emitted along the direction of the incoming quark q.To see this, consider the Mandelstam variables T = (p − k 1 ) 2 , U = (p − k 1 ) 2 and S = 2p • p that are needed to describe this process.Then, from the phase-space constraints in Eq. ( 14), we find T = −t/z, U = −s(1 − z).Therefore, we can obtain the required splitting function by studying the T → 0 limit of the NNLO QCD contribution to the amplitude squared for the process q(p)q(p) → V + g(k 1 ), and by extracting the contribution with the appropriate T −1−2 scaling. 2 The calculation of the 0 → q qV g scat-tering amplitude in the T → 0 limit is available [30], so that the splitting function P (R1V 2) qq can be extracted from that reference.An analytic continuation is required to obtain the initial-state splitting function from the finalstate one; this can be done following the discussion in Ref. [30].For the correct regularisation of the soft limit z → 1 it is important to keep also factors of (1 − z) −a unexpanded in , which fortunately is the case in that reference.Finally, we note that the remaining integration over the single-gluon phase space is straightforward since the phase-space constraints restrict the gluon kinematics to a point that, in fact, no non-trivial integration is needed.The integration over the singular limits of the single-real emission phase space introduces up to two additional powers of −1 so that, in order to correctly obtain the 0 term of the bare beam function, the first six orders of the expansion in of the splitting function have to be known.Ref. [30] contains the first five orders of the splitting function, but the sixth order is only necessary for the soft limit z → 1, so that it can be reconstructed from the soft current calculated in Ref. [45], see Ref. [30] for more details.
In addition to the two-loop virtual corrections to the q → q * + g process, the square of the one-loop correction to the single-gluon emission process has to be included into the calculation of B b,R1V 2 qq .We obtained this contribution by adapting the computation of the NNLO QCD bare beam function to higher orders in dimensional regularization parameter , as reported in Ref. [35].

IV. RESULT FOR THE MATCHING COEFFICIENT
We are now in a to present the N 3 LO contribution to the quark matching coefficient in the generalized large-N c approximation.To this end, we write the O(α n s ) contribution to the matching coefficient, as defined in Eq. ( 6), in the following way where and the plus-distributions δ (z) are referred to as "hard".We therefore write As we already mentioned, the NLO and NNLO contributions to the matching coefficient I (1),( 2) qq are fully known [19,21].Recently, in Ref. [8], it was shown how to extract the soft contributions to N 3 LO matching coefficient described by the constants C (3) k , k = −1, . . ., 5 from known results in the literature [27,[46][47][48][49][50].Also, by using the renormalization group equations for the matching coefficient, all functions F (3,k) + (z) were calculated in that reference.These results, especially the ones for the soft constants, provide an important check on the correctness of our computation.Indeed, we have verified that our results reproduce the constants C (z) reported in Ref. [8] in the limit The new result of this paper is the contribution of hard collinear gluons to the function F (3) δ (z) in the generalized large-N c limit.The result turns out to be remarkably simple.It is expressed in terms of harmonic polylogarithms of the variable z of up to weight five.To present the result in a compact form, we use a notation for harmonic polylogarithms (HPLs) introduced in Ref. [51] and extended in Ref. [52].To this end, we explicitly list the right-most zeros of an HPL index but, starting from the first non-vanishing entry, we do not display trailing zeros in an index anymore.Instead, we add one to the absolute value of the index entry per trailing zero and continue doing so until the next non-zero entry is reached.For example, in the formulas below, H 1,2,1,0 means H(1, 0, 1, 1, 0, z) whereas H 4,1 is H(0, 0, 0, 1, 1, z) etc. Armed with this understanding, we present the result for the hard contribution to I (3) qq (t, z) in the generalized large-N c approximation.To this end, we write where T R = 1/2.We note that all other contributions are subleading either in N f or in N c and are thus neglected.The three functions read We note that the NLO, NNLO and N 3 LO contributions to the matching coefficient I qq can be found in an ancillary file attached to this submission.In addition to the functions F (1,2,3) (z), also the functions F (3) k can be found there, in a computerreadable form.

V. CONCLUSIONS
In this paper, we presented the N 3 LO matching coefficient for the 0-jettiness quark beam function in the large-N c large-N f approximation.We have compared our results for the matching coefficient I qq with the results in the literature [8] and found perfect agreement for all terms that are available.The new result of this paper is the hard contribution to the matching coefficient I qq given in Eqs. ( 20)- (22).The full matching coefficient with soft terms and t-dependent plus-distributions can be found in an ancillary file provided with this article.
Although our large-N c large-N f result is, perhaps, not quite suitable for phenomenology per se, we believe it is an important milestone in the computation of beam functions through N 3 LO QCD.Indeed, it clearly shows that computations of complete matching coefficients for quark and gluon beam functions at N 3 LO are within reach.In fact, although only planar Feynman diagrams are needed for computations in the large-N c limit, we already have all the ingredients for gluonic final states to go beyond this approximation.We are in the process of computing all relevant integrals to describe q → q * + q q (+g) transitions; once these integrals are obtained, going beyond the generalized large-N c approximation will be quite straightforward.The NLO and NNLO coefficients agree with Ref. [37].

FIG. 1 :
FIG.1: Top: example of a triple-real emission amplitude with a quark-antiquark pair in the final state which contributes to the bare beam function in the leading-color approximation and therefore has been included in our computation.Bottom: example of a similar amplitude which is sub-leading in Nc and therefore is not included in our computation.The box labeled H denotes the hard scattering process.