Double-real contribution to the quark beam function at N3LO QCD

We compute the master integrals required for the calculation of the double-real emission contributions to the matching coefficients of 0-jettiness beam functions at nextto-next-to-next-to-leading order in perturbative QCD. As an application, we combine these integrals and derive the double-real gluon emission contribution to the matching coefficient Iqq (t, z) of the quark beam function.


I. INTRODUCTION
The absence of any evidence for physics beyond the Standard Model at the LHC implies a growing importance of indirect searches for new particles and interactions.An integral part of this complex endeavour are first-principles predictions for hard scattering processes in proton collisions with controllable perturbative accuracy.In recent years, we have seen a remarkable progress in an effort to provide such predictions.Indeed, robust methods for one-loop computations developed during the past decade, that allowed the theoretical description of a large number of processes with multi-particle final states through NLO QCD [1][2][3][4][5][6], were followed by the development of practical NNLO QCD subtraction and slicing schemes [7][8][9][10][11][12][13][14][15][16][17] and advances in computations of two-loop scattering amplitudes [18][19][20][21][22][23][24][25][26][27][28][29].This progress led to an opportunity to describe many 2 → 2 partonic processes relevant for the LHC physics with the NNLO QCD accuracy.These impressive developments were recently extended by a breakthrough computation of the N 3 LO QCD corrections to Higgs boson production in gluon fusion [30][31][32].Both, the total cross section and simple kinematic distributions were computed in these references.
The computational techniques employed there rely heavily on the use of reverse unitarity [33] that allows one to map complex phase space integrals to loop integrals and use the machinery of multi-loop computations to reduce the number of independent integrals that need to be computed.
It is clear that further applications of this technology will enable the computation of Z and W production cross sections and basic kinematic distributions through N 3 LO QCD as well.However, it should be also recognized that the theoretical methods employed in Refs.[30,31] limit the number of observables that can be studied with such high-order perturbative accuracy.Indeed, it is highly unlikely that complex fiducial cross sections defined at the level of decay products of the produced color-singlet particles, with additional restrictions on the QCD radiation, can be computed using these techniques.Yet, it is the knowledge of these fiducial cross sections that allows a direct connection of the refined theoretical predictions with the results of the experimental measurements.For this reason, the extension of the results of Refs.[30,31] towards fully-differential cross sections and distributions is highly desirable.
It is known from NNLO QCD computations that the calculation of an arbitrary fiducial cross section or kinematic distribution requires the development of a full-fledged subtraction scheme to identify and remove infra-red and collinear divergences.Developing such a subtraction scheme at N 3 LO QCD is a daunting task.An alternative, more practical, approach is to use slicing methods known from NLO QCD [34] and extended recently to NNLO computations [9,15,16].
To explain the general idea of the slicing method, we consider a process where a color-singlet final state V is produced in proton-proton collisions together with some accompanying QCD radiation X.We do not require the presence of any jets in the final state.It is possible to choose an infra-red safe kinematic variable, that we will refer to as ω, with the property that ω[V ] = 0 and ω[V, X] > 0, provided that the momenta of gluons and quarks in the final state X are neither soft nor collinear to the incoming protons.It follows that, if we consider the process pp → V + X and require that ω[V, X] > ω 0 > 0, we prevent the final state QCD partons from becoming fully unresolved.From the viewpoint of fixed order perturbative computations, this implies that, for ω > ω 0 , we consider a process pp → V + j, so that the final color-singlet state V is outside of the Born phase space.As the result, the desired N 3 LO QCD computation becomes a NNLO QCD computation for pp → V + j for ω > ω 0 .
Given the recent progress, such NNLO QCD computations are definitely possible.
To enable the description of the inclusive process pp → V through N 3 LO QCD, we need to supplement the NNLO QCD computation for pp → V + j described above with the computation of the contribution to the inclusive cross section from phase space region where 0 ≤ ω[V, X] < ω 0 .Note that this contribution includes ω[V, X] = 0 and, therefore, is sensitive to the fully unresolved kinematics.For finite ω 0 , the computation of the ω < ω 0 contribution is as difficult as the full computation.However, significant simplifications occur if we take ω 0 to be very small, ω 0 1.For such small values of ω < ω 0 , the allowed radiation X is either soft or collinear.This already leads to certain simplifications of the required computations since, in soft and collinear limits, the matrix elements for the partonic processes factorize into hard process-dependent matrix elements and universal splitting functions or eikonal factors.However, if this is the end of the story, the required computations are still highly non-trivial because soft and collinear divergences do overlap.
Luckily, this problem can be ameliorated by choosing a particular slicing variable ω.Indeed, for certain choices of the slicing variables, there are factorization theorems that express cross sections at small values of ω through simpler objects that separate soft and collinear dynamics and simplify the relevant computations significantly.
This is exactly what happens for the so-called jettiness variable [35].For the inclusive production of a color-singlet final state in hadron collisions, without requiring any additional jets in the final state, the relevant (zero-) jettiness variable reads Here, p 1,2 are the momenta of the incoming protons, taken to be light-like, Q 1,2 are hardness variables that can be chosen in a number of different ways, and k 1,2,...N are the momenta of the QCD partons that appear in the final state of the process pp → V + X.
The usefulness of T as the slicing variable follows from the factorization theorem [35] that states that for small values of T the differential cross section for pp → V + X can be written as the convolution of the so-called beam and soft functions and the hard cross section for The soft function can be computed order-by-order in perturbation theory.On the contrary, the beam function is determined through a convolution of perturbatively calculable matching coefficient and the non-perturbative parton distribution functions The quantity t in this formula is the so-called transverse virtuality of the quark that participates in a hard process; it is related to the jettiness variable by a simple rescaling.The index i runs over all possible partons including quarks, anti-quarks and gluons.
To arrive at the final result for the cross section in Eq. ( 2) at a particular order in perturbation theory, we need to compute the matching coefficients I ij (t, x 1 , µ), i, j ∈ {q, q, g}, and the soft function to the same order in the perturbative expansion, subtract the divergences by performing PDF renormalization and compute the relevant convolutions.Since we are interested in a fully-exclusive N 3 LO computations, we need to know the soft and the beam functions through N 3 LO QCD in perturbation theory.
We will focus on the computation of the beam function at N 3 LO QCD.The easiest way to think about the different ingredients of the computation is to realize that the matching coefficient I ij can be computed as a phase space integral of the collinear limit of the relevant scattering amplitudes squared.Integrations over multi-particle phase spaces are performed with constraints on a transverse virtuality t and the light-cone component of the fourmomentum of a parton that participates in the hard scattering process [36].
We work in third-order perturbative QCD and focus on the quark-to-quark branching process for definiteness.We then need to consider various transitions of the type q → q * + {g i }, where the number of additional gluons emitted to the final state changes from one to three.
Additional powers of the strong coupling constant, required for a particular final state to contribute at N 3 LO, are then provided by virtual diagrams that renormalize the collinear emissions of one and two real gluons.
We now define the phase space constraint more precisely.To this end, we denote the fourmomentum of the incoming quark as p, the four-momentum of a virtual quark q * that goes into the hard scattering as p * and the total momentum of the emitted gluons as We fix the component of p * along the direction of the incoming momentum p to be z.We then write where p2 = 0 and p • k ⊥ = p • k ⊥ = 0. We define the transverse virtuality t as The value of y can be computed by taking a scalar product of k tot with p, y = −(p•k tot )/(p• p).
We obtain It is also useful to express the constraint on the light-cone component of the virtual quark momentum p µ * through the momenta of the emitted gluons.We do this by considering the scalar product of p * with p and using momentum conservation p * = p − k tot .Defining Eqs. ( 6) and ( 7) provide the phase space constraints that need to be accounted for in the integration over the gluon phase space.Thus, the computation of the contribution of the N -gluon final state to the beam function is proportional to where Ĉp denotes the collinear projection of the square of the full matrix element |M| 2 .
The result is normalized to the square of tree-level matrix element |M 0 | 2 .Needless to say that, after collinear projection, all quantum effects reside on a single incoming quark line and emissions from the incoming anti-quark with momentum p decouple.
We will now discuss the three different contributions to the matching coefficient of the beam function through order O(α 3 s ).The first one is the two-loop renormalization of a single collinear gluon emission.In this contribution the kinematics of a single real emission is fully constrained and the corresponding two-loop virtual correction includes at most the two-loop three-point function with two light-like legs and one off-shell leg.In principle, such contributions are known analytically if not for the fact that they need to be computed in the light-cone gauge to achieve collinear factorization.The light-cone gauge introduces additional propagators to Feynman integrals making them unconventional and difficult to compute.
The second contribution is the single-virtual double-real emission process, i.e. a one-loop correction to the emission of two real gluons in the collinear kinematics.In this case the situation is more complex.Indeed, the one-loop virtual corrections include box diagrams that are sufficiently complicated and the kinematics of two real gluons is sufficiently unconstrained to make the computation of this double-real contribution quite a challenging task.The earlier comment about light-cone gauge also applies here.
The last contribution is the triple real emission process, which involves the integration over the three-gluon phase space subject to the jettiness constraint.Such a computation is also quite demanding.Finally, to arrive at the matching coefficient one has to perform collinear renormalization of the beam function which, at N 3 LO, is also non-trivial.
As the result, we decided to report on this computation in a few separate instalments.We will start with the discussion of the single-virtual double-real contributions to the matching coefficient of the beam function.It is schematically described by Eq. ( 8) with N = 2.

II. MATCHING COEFFICIENT
In this section, we describe the calculation of the double-real contribution to the matching coefficients for the quark beam function at N 3 LO, see Eq. (3).We follow the observation made in Ref. [36] and compute the matching coefficient as phase space integrals of the corresponding splitting functions, using the appropriate kinematic constraints, see Eq. ( 8).
The splitting functions at the relevant perturbative order can be constructed following the general method described in Ref. [37].The large number of phase-space integrals that appear in the course of the computation are calculated using the method of reverse unitarity [33].We consider a massless quark with momentum p that emits two (collinear) gluons of momenta k 1 and k 2 , before entering the hard scattering process as depicted schematically in Fig. 1.We focus on the one-loop virtual corrections to the double gluon emission, i.e. a particular contribution to N 3 LO beam function as we explained in the Introduction.Collinear dynamics on the quark line factorizes, which means that when computing the relevant amplitude squared in the collinear limit, we become insensitive to the hard scattering process.In Ref. [37], it was shown that one can perform this projection by using a physical gauge for the real collinear gluons g(k 1 ) and g(k 2 ) and by inserting p = pµ γ µ in place of the hard scattering matrix element.The vector p is a light-cone vector complementary to the light-cone vector p, see Eq. ( 4).
We find it practical to construct the collinear limit of the relevant amplitude squared Fig. 1 by considering the Feynman diagrams for the process q(p) → q(q) + g(k 1 ) + g(k 2 ) .We employ QGRAF [38] to generate all tree-level and one-loop diagrams, and use them to build up the amplitude in two independent implementations in FORM [39] and Mathematica.This requires care, since the quark q(p) is on-shell, while the quark q(q) is off-shell.This means that when one-loop corrections are generated, we must neglect all self-energy insertions on the on-shell external quark q(p), but we need to include them on the off-shell quark q(q).
As mentioned above, we need to use light-cone gauge for both real and virtual gluons.For example, when summing over real gluon polarizations, we obtain Finally, we also need to use the axial gauge for virtual gluons, with the same reference vector p.We note that the choice of p as the reference vector is convenient but not necessary.
We find 3 tree-level diagrams and 30 one-loop diagrams.Computing the interference, we obtain 90 "three-loop" phase-space diagrams.After performing the relevant Dirac algebra, each of these diagrams can be written as a linear combination of "three-loop" phase space integrals.We need to organize these integrals into integral families, in order to perform a reduction to master integrals using the well-known integration-by-parts identities [18].When dealing with phase space integrals subject to constraints, this step involves one more subtlety compared to what is done for standard multi-loop Feynman integrals.In order to understand this, we recall that a complete integral family should contain exactly as many propagators, as is the total number of independent scalar products that can be formed from the loop momenta and the external momenta.In our case, there are effectively three loop momenta and two external momenta (p and p); it follows that we can construct 12 independent scalar products.Therefore, we need to map our phase-space integrals to integral families with 12 independent propagators.This requires some extra work.Indeed, consider again Fig. 1.It is easy to convince oneself that, after Dirac algebra, the phase-space integrals stemming from these diagrams will have the following general form Here k 1 , k 2 are the momenta of two on-shell gluons and k 3 is the momentum of the virtual gluon.The D j are different, unspecified, denominators and N is a polynomial of scalar products of the loop momenta and the external momenta p, p.It should be easy to see that the diagrams obtained from Fig. 1 can generate at most 11 different propagators, including the combinations k 1 • p and k 2 • p which come from the sum over polarizations in Eq. ( 9), along with similar denominator factors coming from the gluon propagators in axial gauge.Counting four delta-functions as "propagators" (we have in mind the reverse unitarity relation δ(x) ∼ → 1/x), we get 15 propagators in total.Since this number is larger than 12, the number of independent scalar products involving the momenta k i , the propagators are linearly dependent and Eq. ( 10) therefore does not constitute an integral family.
To remedy this problem, for every contributing diagram we partial fraction some of the linearly dependent propagators.For instance, we write The partial fractioning effectively splits the diagrams into several terms, each of which contains at most 12 linearly independent propagators.At that point we can introduce well-defined integral families.
With this procedure we find that all diagrams can be expressed in terms of 19 independent integral families.The corresponding integrals can be written as where top ∈ {A1, A2, . . ., A19} and the inverse propagators D 1 , . . ., D 8 for each integral family are shown in Table I.

III. MASTER INTEGRALS
Having expressed the amplitude in terms of master integrals, we proceed with the discussion of their evaluation.The master integrals depend on the three quantities s = 2p • p, t and z, cf.Eq. ( 14).However, some of these dependencies are quite simple.Indeed, we re-define the gluon momenta as k i = ki t/z, p = p t/z and p = p s z/t.Since p • p = s/2, we find p • p = 1/2.Applying these transformations to integrands of master integrals, we find Note that the second step in the above equation implies that one can get the expression for the re-scaled integral upon taking the definition of the original integral with its full s, t and z dependence and then evaluating it for s = 1 and t = z.
In order to calculate the master integrals I top n 1 ,...,n 8 (z) on the right-hand side of Eq. ( 17), we compute their derivative with respect to z and re-express the result in terms of master integrals by means of integration-by-parts reduction, thus producing a closed system of differential equations.We write it as where I(z, ) contains the 128 master integrals in Eq. ( 15) with s = 1 and t = z, and Â(z, ) is a 128 × 128 matrix.We use the program Fuchsia [41] to transform the system of equation to the -form using the algorithm described in Ref. [42].We find where the Âk are constant matrices.The vector I can (z, ) contains the canonical master integrals; it is related to the original master integrals I(z, ) = T (z, ) I can (z, ) by a transformation matrix T (z, ).
The system of differential equations in Eq. ( 19) is solved iteratively for the coefficients of I can (z, ) in an expansion around = 0.The result of this procedure can be written as where the 128 × 128 matrix M (z, ) contains elements of the form The inner sum runs over vectors w of length k whose components are drawn from a set {−1, 0, 1, 2}, which are the singular points in the differential equations Eq. ( 19).The G( w; z) are multiple polylogarithms [22,25,43,44], which are defined iteratively as with the special cases The constant vector B( ) in Eq. ( 20) is fixed by the boundary conditions required for the solution of the system of first-order differential equations.We extract the boundary condition from the behaviour of the master integrals in the limit z → 1.
In principle, it is possible to extract the boundary conditions by writing and taking the limit z → 1.The multiple polylogarithms G( w; z) inside M (z, ), appearing on the right-hand side of Eq. ( 24), are, in general, logarithmically divergent in this limit.
This divergence should be extracted by writing the multiple polylogarithms in the form Computing a sufficient number of integrals I(z, ) that appear on the left-hand side of Eq. ( 24) in the z → 1 limit allows to determine the constants B( ).
In performing the computation of the master integrals we should remember that the beam function matching coefficients should be treated as distributions in (1 − z).Indeed, the matching coefficients have to be convoluted with the PDFs and all divergences in the limit z → 1 have to be regularized, extracted and renormalized in order to obtain the final finite result in terms of quantities like δ(1−z), [1/(1−z)] + , etc.To extract these distributions, one needs to have the master integrals written in the "resummed" form, i.e. the dependence in the limit z → 1 has to be made explicit in the form of (1 − z) a+b powers.It is convenient, therefore, to extract the boundary constants in a similar way, namely by matching equal powers of (1 − z) on both sides of Eq. (24).
In order to do that, the right-hand side of Eq. ( 24) should first be written in the "resummed" form ,m c ,m ( )(1 − z) n+ log m (1 − z), which can be easily achieved by solving the differential equations in Eq. ( 19) in the limit z → 1 whose solution is given by a matrix exponential The new constants in H( ) can be expressed as linear combinations of the constants in B( ) by equating the right-hand side of Eq. ( 20) taken in the limit z → 1 and the right-hand side of Eq. ( 26) expanded in .Multiplying both sides in Eq. ( 26) by the transformation matrix T (z, ) in the limit z → 1 and performing the matrix exponentiation gives the leading-power behaviour of the master integrals The coefficients C i, ,m ( ) are in general linear combinations of the 128 constants in B( ).
Therefore, we need to compute sufficiently many linearly independent C i, ,m ( ) to determine all boundary constants.Fortunately, we can fix many boundary constants by identifying the constants C i, ,m ( ) that must vanish in order to produce the correct behaviour of the master integrals in the limit z → 1.
As an explicit example, consider the second master integral I (2) (z → 1) = I A1 1,0,0,0,0,1,0,0 | z→1 .On the one hand, the differential equations predict that in the limit z → 1 the integral has the form On the other hand, its integral representation suggests that its leading behaviour scales as (1 − z) 1−2 , which implies that C 2,0,0 ( ) = C 2,−2,0 ( ) = 0. Let us verify this by inspecting the integral representation The integral over k 3 is straightforward to evaluate.The result is proportional to ) and the on-shellness p 2 = 0 constraints.The last delta-function in the numerator of Eq. ( 29) fixes the projection of the total emitted momentum k 12 ≡ k 1 + k 2 to be small in the limit z → 1.To capture that, it is convenient to make a Sudakov decomposition of k 1 and k 2 We use this decomposition to compute 2p . Moreover, the fact that the momenta k 1 and k 2 are both on-shell and have positive energy, implies that and scale as O ((1 − z) 2− ).Three of the delta functions in Eq. ( 29) scale as O ((1 − z) −1 ).Putting everything together, we find that . Therefore, there are no contributions to the integral that ).We conclude that the coefficients of these regions vanish: After finding all the coefficients C i, ,m ( ) that must vanish because of similar arguments, we acquire enough relations to express 104 of the boundary constants in B( ) in terms of a remaining set of 24 constants.To determine the latter constants we performed explicit computations of non-vanishing regions of selected master integrals in the limit z → 1.
The boundary integrals that we have calculated are B 2 = I A1 1,0,0,1,0,0,0,0 s=1,t=z,z≈1 = C 3,−3,0 ( ) the computations that led to the expressions listed below.For convenience we extract a common -dependent pre-factor, The results for the constants, up to weight six, are C 10,−2,0 ( ) = In the following subsection we provide some examples of the calculation of some of the constants above.All other constants can be obtained by suitable extensions of the computations presented below.

B. Boundary integral B 1
Boundary integral B 1 is one of the simplest integrals and can be computed exactly in .Its integral representation is given by The integral over k 3 is performed first.In this case it is a simple one-loop bubble integral given in Eq. (A1).For convenience, we also introduce the abbreviation After these steps the boundary integral is written as where we have used the on-shell conditions p 2 = k 2 2 = 0.The prefactor C bub ( ) is defined below Eq. (A1).
We proceed by writing the remaining integration measure in the form which is convenient for extracting the leading behaviour in the limit z → 1.The expression in Eq. ( 83) is obtained by inserting the Sudakov decomposition Eq. ( 30) into Eq.( 81) and integrating out the length of the vector k i⊥ .Accordingly, one should set 2k i • p = α i , ) in the integrand when using the measure in Eq. ( 83).
where the bounds on the integral over Q 2 ⊥ are dictated by the conditions As a result, we find Other boundary integrals, that can be computed via the same method, lead to a two-particle massless phase-space integral of the form One way to calculate the phase-space integral in Eq. ( 93) is to work in the rest frame Q = (Q 0 , 0), to carry out the resulting angular integrations (see eq. ( 49) in Ref. [46]), and to write the expression back in a Lorentz invariant form.The result is D. Boundary integral B 8 As the next example, we consider the boundary integral B 8 .It is given by The result of the one-loop one-mass box integral is given as Box1 in Eq. (A3).In order to find the behaviour of Box1 in the limit z → 1, it is convenient to rewrite the last two hypergeometric functions in Eq. (A3) using the identity As a result, the last two terms in Eq. (A3) combine to produce a contribution that is subleading with respect to the first term in the limit z → 1.We are left with the calculation of the following integral Here the integrand depends on k 2 12 but, unlike in the previous example, inserting 1 = d d Q δ d (k 12 − Q) will not be helpful, because that would lead to a two-particle phasespace integral whose integrand contains the hypergeometric function in Eq. ( 97).In such a situation there is no choice but to perform a non-trivial angular integration directly.In this example we demonstrate how to carry out such an integral.
Following the above discussion, we write the integral B 8 as the sum of two terms The contribution from case (i) to the integral is where α 2 = 1 − α 1 and β 2 = 1 − β 1 .After a change of variables β 1 → r with r = α 2 β 1 α 1 β 2 and, subsequently, α 1 → t with α 1 = t r+t−rt , it becomes After applying identities for hypergeometric functions to simplify their argument and extract their singularities at the endpoints of the integration, we obtain The integrand has singularities at points r = 1 and t = 0, 1.The integral may be carried out by performing suitable subtractions at this points that enables expansion of complicated integrals in .This procedure is tedious but relatively standard and its explanation is thus omitted here.The calculation of X (ii) 8 can be performed along the same lines.The final result for this boundary integral reads

IV. NUMERICAL CHECKS OF MASTER INTEGRALS
The calculation of the master integrals required many non-trivial steps and therefore it is good to have a completely independent check of our results for the integrals.There are various public codes that can evaluate loop integrals numerically, but none as of yet that can compute phase-space integrals, especially of the type that we consider in this paper.The complication arises from integration over angles of the emitted gluons since it is challenging to find a suitable parametrization for the angular degrees of freedom.
One possibility, pointed out in Ref. [30], is to use the Mellin-Barnes (MB) representation for this purpose.The idea is to split complex denominators that appear in integrals into integrals of products of simpler scalar products, perform ensuing angular integrals analytically using results of Ref. [46] and compute the resulting MB integrals numerically using available MB packages [47].
Taking, as an example, a propagator ((p In Eq. ( 105) the contour has to be chosen in such a way that the poles of Γ(−z) are to the right and the poles of Γ(λ + z) are to the left of the contour that runs along the imaginary axis.However, if either x or y is negative in Eq. ( 105), the numerical evaluation of the right hand side of Eq. ( 105) may become unstable because of the exponential increase of (−x−i0) z as Im[z] → ∞.This implies that if we would split the denominator (k into MB integrals, numerical integration may become unstable. 1t is possible to get around this problem by considering a decay process instead of the production process.Indeed, our phase-space integrals correspond to an incoming parton emitting two collinear particles before entering a hard process; since the incoming parton has zero invariant mass, the off-shellness of a quark line becomes negative after gluon emissions. If, on the other hand, a quark with positive virtuality leaves the hard process and decays to a zero-virtuality final state quark by emitting gluons, virtual quark lines at intermediate stages have positive virtualities for which numerical integration of the relevant MB representations is straightforward.In order to get from a production kinematics to a decay kinematics, we need to change the four-momenta p → −p, p → −p.The constraint on the longitudinal momentum of a virtual quark in the decay kinematics becomes k 12 • p = (z − 1).Since this quantity should be positive, we have to take z ≥ 1.The virtuality constraint reads k 12 • p = −t/(2z).Since k 12 • p is positive definite, we have to take t ≤ 0. For the sake of example, consider an integral in the decay kinematics where κ = t/z and 2p • p = 1.Its analytic expression can be found from our solutions for integrals in the production channel by an analytic continuation of κ and z to the region κ ≤ 0, z ≥ 1.2 Note that the variable dependence s = 2p • p is unchanged when moving from production to decay kinematics.
The propagators in the decay kinematics, Eq. ( 106), are given by sums of positive-definite quantities and, for this reason, are more suitable for the MB integration.It is therefore more convenient to numerically compute integrals in decay kinematics, Eq. ( 106), and compare them with analytically-continued integrals computed in the production channel.
To implement this in practice, we note that if specific combination δ(1 − 2k 1...n • (p + p)) n i=1 δ + (k 2 i ) of delta functions appears in the integrand, it is known how to perform phase-space integrals with the MB method [30].In our case, different delta-functions appear in integrands but we may produce such a combination of delta functions by integrating the function I decay (κ, z) over the variables κ and z, both from −∞ to +∞, with an extra delta function insertion that imposes a further constraint κ = z − 2. We obtain N = ( The second equality in Eq. (108) follow from the fact that the decay kinematics imposes that I decay (κ, z) is exactly zero outside of the region κ ≤ 0, z ≥ 1.As we already mentioned, the first line in Eq. ( 108) can be evaluated starting from the analytic solution for I production (t, z) constants C i ( ) multiplied by Taylor expansions in that contain multiple polylogarithms G a (z) up to weight 6.The constants C i ( ) are provided separately as Laurent expansions in .
Our results are ingredients to the computation of the double-real contribution to third-order matching coefficients I ij (t, z, µ) of quark and gluon beam functions in perturbative QCD.
To illustrate this, we focus on the quark-to-quark branching process, i.e. i = q and j = q, and write the perturbative expansion as The N 3 LO term I Its dependence on t, z and µ factorises as expected.The t-dependent factor is expanded in terms of plus distributions according to the formula The non-trivial z-dependent factor on the right-hand side of Eq. ( 112) may be split into a part that diverges in the soft limit z → 1 and a finite remainder With our results for the master integrals we find that, for instance, the divergent part of the matching coefficient is given by the following compact expression We remark that Eq. ( 115) does not contain the colour factor C 3 F ; this is similar to the observation in Ref. [48] that eikonal factors are not renormalized by one-loop QED corrections.

VI. CONCLUSIONS
We computed the master integrals for the double-real contribution to the matching coefficient of the quark beam function in third-order perturbative QCD.The matching coefficients are obtained as collinear limits of QCD amplitudes, with additional constraints on the phase space that fix both light-cone components of the real radiation.We calculated the resulting non-standard phase-space integrals with the methods of reverse unitarity and differential equations, and obtained the boundary conditions from explicit computations of suitable integrals in the soft limit.We provide the master integrals as Laurent series in the dimensional regulator , which contain multiple polylogarithms up to weight six.The result of this paper is an important ingredient for the computation of quark and gluon beam functions through order O(α 3 s ) in QCD.The completion of that task requires additionally the results for two-loop corrections to the single-real emission process as well as the triple-real emission process, which should be feasible using techniques employed in this paper.

Figure 1 .
Figure 1.Double-real diagrams for the calculation of the quark beam function at N 3 LO.

( 3 )
qq (t, z, µ) receives three contributions: single-real, double-real and triplereal.The master integrals calculated in this paper can be used to compute the double-real contribution to the matching coefficient.Restoring the normalization, we writeI (3),RRV qq (t, z, µ)