Triple-real contribution to the quark beam function in QCD at next-to-next-to-next-to-leading order

We compute the three-loop master integrals required for the calculation of the triple-real contribution to the N$^3$LO quark beam function due to the splitting of a quark into a virtual quark and three collinear gluons, $q \to q^*+ggg$. This provides an important ingredient for the calculation of the leading-color contribution to the quark beam function at N$^3$LO.


I. INTRODUCTION
The detailed exploration of perturbative Quantum Field Theory has played an important role in collider physics during the last decade. In fact, the need to study the recently discovered Higgs boson [1,2] and the absence of any sign of physics beyond the Standard Model in LHC experiments are behind an impressive effort of particle theorists to provide predictions for important LHC observables with high precision.
Although precision physics at hadron colliders is very difficult, the LHC experiments have been performing very well, having already delivered measurements for multiple observables at the percent level and even beyond, see e.g. Refs. [3][4][5][6][7][8][9]. Comparing these experimental results with equally precise theoretical predictions, will make it possible to search for New Physics indirectly by probing energy scales far beyond the direct reach of the LHC. These considerations, augmented by an impressive experimental progress, have been continuously pushing the default standard for theoretical predictions for LHC physics from leading to next-to-leading [10][11][12][13][14][15][16][17][18] and, more recently, to next-to-next-to-leading order (NNLO) in QCD. 1 While calculations at NNLO are typically sufficient to match the foreseeable precision of present and future LHC measurements, there is a handful of interesting processes for which theoretical predictions at even higher orders of perturbative QCD (i.e. N 3 LO QCD) are warranted. This may happen for several reasons. Indeed, in some cases the convergence of the perturbative expansion in the strong coupling constant α s turns out to be so slow that even NNLO QCD predictions have a sizable uncertainty. Prominent examples of such a situation are processes where color-singlet final states are produced in gluon fusion. For the important case of Higgs boson production in gluon fusion, it was explicitly shown that N 3 LO QCD corrections are crucial to stabilize theoretical uncertainties at the few percent level [19]. In other cases, e.g. the Drell-Yan process, large statistics and clean final-state signatures led to experimental measurements with very high precision that is posed to increase further during Run III and the high-luminosity phase of the LHC. A theoretical description of the Drell-Yan process with matching or better precision remains a formidable challenge for the theory community.
The theoretical efforts aimed at extending the current computational technology to enable it to handle N 3 LO calculations have recently culminated in the computation of the N 3 LO 1 At least as far as processes with relatively simple final states are concerned.
QCD corrections to Higgs boson production in gluon fusion at the LHC [19][20][21][22]. Since these computations deal with a relatively simple final state and aim at calculating inclusive quantities, it is possible to employ the method of reverse unitarity [23] to simplify them. Although the calculation of the N 3 LO QCD corrections to the Higgs boson production cross section is a landmark in perturbative computations in collider physics, the extension of the methods used in that computation to more complicated final states and more differential observables does not appear to be straightforward and it is interesting to think about alternative options.
For definiteness, let us consider the production of a color-singlet final state V in protonproton collisions pp → V . Quite generally, the description of this process at N 3 LO in QCD requires the knowledge of the NNLO QCD corrections to the production of V together with an additional QCD jet, pp → V + j. The difference between pp → V + j at NNLO QCD and pp → V at N 3 LO QCD is that the jet in the former case can become unresolved and that the virtual corrections to pp → V have no counterpart in the pp → V + j calculation. Since the difference between the two calculations appears in the kinematic regions where the colorsinglet final state barely recoils against the QCD radiation, one can imagine partitioning the phase space into regions with and without recoil, using the NNLO QCD prediction for pp → V + j in the former region and studying the virtual corrections together with soft and collinear QCD radiation in the latter. This is the essence of a so-called slicing method.
For colorless final states, a widely used variable to slice the phase space into resolved and unresolved regions is the transverse momentum of the color-singlet V [24]. More recently, the so-called N -jettiness observable [25][26][27] has allowed to generalize this idea to cases with final-state jets. In the current paper we will focus on the latter variable and, in particular, on the case of 0-jettiness, which is required to describe the inclusive production of a color-singlet final state.
To this end, we consider the process pp → V + X, where X represents the final-state QCD radiation. We denote the momenta of the incoming and outgoing partons by p 1,2 and k 1,..,n , respectively, and write the 0-jettiness variable as In Eq. (1), Q 1,2 are the so-called hardness variables for the initial-state partons; they can be chosen in different ways and they are not relevant for the following discussion. The 0-jettiness variable T has two important properties that allow one to use it as a slicing variable. Indeed, it follows from the definition Eq. (1) that T = 0 in the absence of resolved QCD radiation, i.e. for the process pp → V . However, in the presence of any resolved QCD radiation one finds that T > 0. We can therefore introduce a cut-off T 0 and divide the phase space for V + X into two disjoint parts. We write schematically Note the NNLO subscript in the second term on the right-hand side in Eq. (2); the reason for its appearance is that by imposing the T > T 0 constraint, we exclude the situation where all final-state partons become unresolved so that the calculation for T > T 0 reduces to the computation of the NNLO QCD corrections to pp → V + j. Such calculations have already been performed for a variety of final states and we consider them to be known [27][28][29][30][31][32].
On the other hand, the first term on the right-hand side of Eq. (2) still receives contributions from those regions of phase space where the final-state radiation is fully unresolved. In general, the computation of these contributions can be as difficult as the full N 3 LO calculation itself. However, for 0-jettiness, this does not happen. Indeed, it was shown in Ref. [25] that the cross section for pp → V + X simplifies substantially in the limit T → 0 and can be written as a convolution of the hard cross section for pp → V with the so-called beam and soft functions. The cross section reads where the two functions B stand for the beam functions associated with each of the initialstate partons and S represents the soft function. The factorization of soft and collinear radiation, made apparent in Eq. (3), is the key property of the 0-jettiness variable that simplifies the calculation of the differential cross section in the small-T limit.
The cross-section formula Eq. (3) implies that, in order to employ the 0-jettiness slicing to compute the N 3 LO corrections to pp → V + X, the beam and soft functions must be known at the same perturbative order. While the soft function is a purely perturbative object and can, at least in principle, be computed order-by-order in perturbation theory, the beam-function computation requires a convolution of perturbative matching coefficients I ij with the non-perturbative parton distribution functions (pdfs) f j The computation of the N 3 LO QCD corrections to the matching coefficient I qq is the main topic of this paper. At three loops, I qq receives contributions from three classes of partonic subprocesses: the emission of three collinear partons, which we will refer to as the triple-real contribution (RRR); the one-loop corrections to the emission of two collinear partons, or the double-real single-virtual contribution (RRV); and, finally, the two-loop virtual corrections to the emission of one collinear parton, or the single-real double-virtual contribution (RVV).
In a previous paper [33], we presented the master integrals required for the calculation of the RRV contribution with two emitted gluons to the matching coefficient I qq . In this paper, we focus on the master integrals required for the computation of the RRR contribution to the matching coefficient that originate from the process where the initial-state quark emits three collinear gluons before entering the hard scattering process. We note that the same master integrals can be used to compute the N f -enhanced triple-real contribution to I qq , caused by the emission of a gluon and a quark-antiquark pair collinear to the initial-state quark.
The rest of the paper is organized as follows: in Section II we explain how to compute the RRR contribution to the matching coefficient I qq by considering collinear limits of scattering amplitudes and how reverse unitarity can be used to reduce this calculation to the computation of a large set of three-loop master integrals. We then show in Section III how these integrals can be computed using the method of differential equations. In Section IV, we explain how the calculation was validated and we present our final results in Section V.
We conclude in Section VI. The list of master integrals can be found in Appendix A. Some peculiar identities among master integrals are described in Appendix B. The results for the master integrals are provided in computer-readable format in an ancillary file, which is available at https://www.ttp.kit.edu/_media/progdata/2019/ttp19-009.tar.gz.

II. MATCHING COEFFICIENT
In this section we discuss how to compute the N 3 LO contributions to the matching coefficient I qq for the 0-jettiness beam function. Since the matching coefficients describe the physics of collinear emissions off the incoming partons, they can be calculated by integrating the collinear limits of the corresponding scattering amplitudes squared, over the phase space restricted by the fixed value of the 0-jettiness variable. Figure 1: The process q i → q * i + ggg and the process q i → q * i + q jqj g for i = j.
More specifically, the phase-space integration must be performed by imposing constraints on the transverse virtuality of the collinear partons and on the light-cone momentum of the parton that enters the hard-scattering process [25]. Since singular collinear emissions factorize on the external lines, the hard-scattering process decouples. The collinear emissions are described by splitting functions; for this reason, the relevant contributions to the matching coefficients can be computed by integrating these functions over a restricted phase space [34].
This observation is particularly useful since the prescription for computing the splitting functions to any order in the strong coupling constants has been laid out in Ref. [35].
Focusing on the triple-real (RRR) contribution to the matching coefficient I qq , we need to consider the tree-level splitting of a quark into a virtual quark of the same flavor and three collinear partons. These three partons can be either three gluons or a quark-antiquark pair and a gluon, so that there are two generic possibilities: q i → q * i + ggg and q i → q * i + q jqj + g. In this paper we consider the process involving three collinear gluons as well as the process involving a collinear gluon and a collinear quark-antiquark pair of a different flavor with respect to the incoming quark, i.e. i = j, see Fig. 1. The case i = j requires additional contributions that are not considered in this paper. However, it is easy to see that the neglected contributions are sub-leading in the N c → ∞ limit, where N c is the number of colors, and in the N f → ∞ limit, where N f is the number of massless quark flavors. Hence, even neglecting the i = j contributions, we can obtain the result for I qq that is valid in the large-N c or large-N f limits. In the remainder of this section, we focus our discussion on the process in Fig. 1(a) for definiteness.
We can now describe the details of the calculation. We follow the discussion in Ref. [33], where the master integrals for the double-real single-virtual contribution to I qq were computed. We consider a massless quark with momentum p which emits three collinear gluons with momenta k i , i = 1, 2, 3, and enters the hard process with momentum p * q(p) → q * (p * ) + g(k 1 ) + g(k 2 ) + g(k 3 ) , As we already explained, the relevant contribution to the matching coefficient is obtained by integrating the q → q * + ggg splitting function over the phase space of the emitted gluons with appropriate constraints. In order to write these constraints in a convenient form, we fix the component of the momentum p * along the momentum of the incoming quark p and write In Eq. (6), we used k 123 = 3 i=1 k µ i . We also introduced a light-cone momentump, which is complementary to p so thatp 2 = 0 and p · k ⊥ =p · k ⊥ = 0. The emitted gluons are on the mass shell, i.e. k 2 i = 0 for i = 1, 2, 3. With these definitions we have y = −(p · k 123 )/(p ·p). We now introduce the transverse virtuality t = −(q 2 − k 2 ⊥ ) and, using the above results, write it as Note that, in the case of collinear emissions, t ∼ T . We also impose a constraint on the light-cone component of the momentum of the quark that enters the hard process. We write it as Using Eqs. (7) and (8), we write the generic contribution of the three-gluon final state to the matching coefficient I qq in the following way In Eq. (9) the integrand P qq (p,p, {k i }) describes the q → q * + ggg splitting function. Below we explain how to compute it.
As described in Ref. [35], the q → q * + ggg splitting function can be obtained as the collinear projection of the squared scattering amplitude for the corresponding process Fig. 1(a). To this end, we generate the scattering amplitude as a sum of all diagrams that contribute to the q → q * + ggg process. The diagrams are turned into mathematical expressions with the standard QCD Feynman rules, albeit with a symbolic placeholder for the arbitrary hard-scattering process. The axial gauge is chosen for the gluons, both internal and external ones, and the light-cone vectorp from Eq. (6) is selected as the corresponding gauge-fixing vector. Squaring the amplitude, we produce a Dirac trace of the form Tr[· · ·p * Ĥp * · · · ], wherep * = γ µ p * µ and p * is the momentum that enters the hard scattering process. The Dirac matrixĤ is a symbolic representation for the (product of) gamma matrices in the hard interaction. The collinear projection of the squared scattering amplitude, schematically depicted in Fig. 2, is achieved by making the replacement which has the effect of removing all non-singular contributions in the limit where all three gluons become collinear to the incoming quark.
In practice, we generate the diagrams that contribute to the process q(p) → q * (p * ) + [36]. We perform the relevant Dirac and Lorentz algebra in FORM [37] and Mathematica in two independent implementations. Since we work in the axial gauge with the gauge-fixing vectorp, the sum over polarizations for a gluon with After applying the collinear projection in Eq. (10), the squared amplitude can be written as a linear combination of a large number of scalar phase-space integrals of the following form Here, N is a generic combination of scalar products of the parton momenta, and D j are propagators, including linear propagators that originate e.g. from the denominators in Eq. (11).
These integrals can be computed efficiently using the method of reverse unitarity [23], which allows one to turn the delta function constraints in Eq. (12) into cut propagators, mapping the problem of computing phase-space integrals onto the calculation of a large number of three-loop Feynman integrals.
We need to organize these integrals into integral families to enable the reduction to master integrals through the integration-by-parts identities (IBPs) [38][39][40]. As is often the case when dealing with phase-space integrals in the framework of reverse unitarity, this Figure 2: The collinear projection of the squared scattering amplitude for the process q → q * + ggg and the process q → q * + q q g for q = q.
step is not entirely straightforward. Indeed, a well-defined integral family requires as many propagators as the number of independent scalar products in the problem at hand. In our case there are two independent external momenta p andp and three gluon momenta k i . This implies that any integral family must contain exactly 12 independent propagators. By directly inspecting the Feynman diagrams, it is easy to see that, after accounting for the delta function from the 0-jettiness constraint, many diagrams do generate scalar integrals of the form shown in Eq. (12), but with more than 12 different propagators.
To remedy this problem, we need to use partial fractioning. For example, it may happen that an integral contains all three linear propagators 1/k i ·p with i = 1, 2, 3. However, the 0-jettiness constraint in Eq. (12) implies that the three propagators 1/k i ·p are not linearly independent. Indeed, we can write which allows us to reduce the number of propagators by one.
Unfortunately, this procedure is ambiguous, since different ways of partial fractioning can lead to different integral families and different integrals. While it is usually sufficient to use the IBP identities to remove most of this redundancy, some of the integrals that appear to be independent under IBPs can still be related by special partial fractioning identities and we need to separately account for that possibility.
Due to the ambiguity mentioned above, we find it convenient to introduce an overcomplete set of integral families in order to simplify the mapping of diagrams to topologies.
Nevertheless, performing the IBP reduction and accounting for additional identities that originate from the partial fractioning, we find that all diagrams can be expressed in terms of 91 master integrals which are drawn from 19 different topologies, see Table I. We performed the reduction to master integrals using Reduze [41] and KIRA [42], both of which support the generation and solution of IBPs for Feynman integrals with cut propagators, and we verified that the results of the two reduction codes are equivalent.
We use the following notation for the master integrals where d = 4 − 2 and the subscript 'top' indicates one of the topologies in Table I where the inverse propagators D i for each topology are defined. The integration measure for each final-state particle reads We use these notations to present the list of master integrals in Appendix A.
While the set of master integrals shown in Eq. (A1) is indeed minimal with respect to the IBPs, we were able to find two additional relations between them, that do not follow from IBPs and partial fractioning. These identities read They allow us to reduce the number of independent master integrals from 91 to 89. Nevertheless, we prefer to compute the full set of 91 master integrals and verify these identities a posteriori. We note that these identities can be proven by studying the differential equations satisfied by the four master integrals that appear in Eqs. (16) and (17) together with the direct inspection of their integral representations. We describe the proof in Appendix B.

III. MASTER INTEGRALS
The master integrals defined in Eq. (14) depend on t, z and s = 2p ·p. However, the dependence on s and t is trivial. This becomes manifest after the simultaneous re-scaling The re-scaling has the effect of extracting powers of s and t from each integral, leaving only a non-trivial dependence on z. Explicitly, we find I top n 1 ,n 2 ,n 3 ,n 4 ,n 5 ,n 6 , Here we use the shorthand notation k ij = k i + k j and k ij = k i + k j + k .
where M = 1 + n 6 + n 7 and N = 5 − 3 − differential equations as The entries of the matrixÂ(z, ) are rational functions of z and .
The complexity of these differential equations depends strongly on the explicit form of the matrixÂ(z, ), which, in turn, depends on the choice of the master integrals. Our goal is to choose the master integrals in such a way that the matrix becomes canonical and Fuchsian [46][47][48],Â(z, ) = z 0Â z 0 z−z 0 . Note that the matricesÂ z 0 should be both zandindependent. If such a form is found, the process of solving differential equations simplifies greatly.
It turns out, however, that the system in Eq. (19) cannot be brought to a canonical Fuchsian form without replacing z with a more suitable variable. Indeed, it is easy to see that upon integration, the homogeneous terms of some of the differential equations produce the square root z(4 − z). The presence of square roots complicates substantially the problem of finding a canonical Fuchsian form. To rationalize it, we change variables from z to x according to the following equation Having removed all square roots, we can construct the appropriate transformation I(x, ) = T (x, ) I can (x, ) with the program Fuchsia [49]. As a result, we find x 0 x − x 0 I can (x, ) .
It is convenient to solve the system of differential equations Eq. (21) expanding around = 0. We write I can (x, ) as where B( ) are the integration constants. The x-dependence resides solely in the matrix M (x, ), whose elements have the form We calculate the sum over k up to and including k = 6, corresponding to O( 6 ), which is the highest order that will contribute to the finite part of the matching coefficient in the They can be evaluated numerically with the help of the program Ginac [54]. Apart from the technical difficulty in handling large expressions, the construction of the matrixM (x, ) can be done in a relatively straightforward way.
On the contrary, the determination of the integration constants B( ) in Eq. (22) is much less straightforward. We obtain them by analyzing the master integrals in the limit z → 1.
To this end, it is important to recognize that the master integrals significantly simplify in that limit. In particular, to leading order in (1 − z) we can replace the propagators 1/(p − k ij ) 2 with 1/(−2k ij · p). Note that this replacement renders the integrals uniform functions of the momenta k i so that, in the soft limit, the integral factorizes into a constant and a (1 − z)-dependent factor.
The possibility to neglect k 2 ij relative to k ij · p follows from the following argument. Let us select a frame in which the external momenta are p = 1 √ 2 (1, 0, 0, 1) andp = 1 √ 2 (1, 0, 0, −1) and introduce a Sudakov decomposition of the gluon momentum k i Since β i = k i · p = k 0 i p 0 (1 − cos θ i ) and α i = k i ·p = k 0 ip 0 (1 + cos θ i ), we conclude that all α's and β's are positive definite. According to the phase-space constraints Eq. (12), the sum α 123 = α 1 + α 2 + α 3 goes to zero in the z → 1 limit and, since all α's are positive, we conclude that each α i goes to zero in that limit at least as fast as O(1 − z). In contrast, the sum of the β i 's is constrained to be equal to one, so that up to two of them could vanish at z = 1. We write where we have used that p 2 = 0. In terms of the Sudakov parameters, k 2 ij reads where we have used k 2 i = k 2 j = 0. Assuming that, in the limit z → 1, each α i = O(1 − z) and each β i = O(1) we find k 2 ij = O(1 − z) and 2k ij · p = O(1). Hence, we can neglect k 2 ij relative to 2k ij · p. The situation does not change, should any of the α i 's vanish faster than O(1 − z). Another possibility is that both β i and β j vanish as O(1 − z), such that 2k ij · p → 0. However, in that situation k 2 ij scales as O((1 − z) 2 ) or faster, and we can again neglect it relative to 2k ij · p. Therefore, a replacement is valid in the z → 1 limit, to leading power in (1 − z).
Since the replacement in Eq. (28) implies that all propagators become uniform functions of the gluon momenta in the soft limit, the extraction of the (1 − z)-dependence of any integral becomes straightforward. We note that, in that limit, the phase-space constraints from Eq. (18) become δ k 123 · p − 1 2 δ k 123 ·p − (1−z) 2 and, upon re-scaling the momenta as of the master integrals.
It follows that in the soft limit, each integral scales as (1 − z) n−3 with an integer n that is integral-dependent. Hence, all canonical master integrals should be free of logarithmic singularities as z → 1, or equivalently as x → R ± 1 , beyond those that correspond to the expansion of (1 − z) −3 in powers of . This observation allows us to impose a regularity condition, which fixes 81 integration constants.
The remaining integration constants are obtained by an explicit computation of ten non-canonical integrals in the limit z → 1. These integrals read To present the results, it is convenient to extract the common -dependent factor, where Ω n = 2π n/2 /Γ(n/2). With this normalization, the constants C i read, up to weight six, In the following we describe the various techniques that we used for computing these constants. We discuss the integrals B 1 , B 8 , B 9 and B 10 as representative examples. All other integrals can be obtained in similar ways. We stress that all results in Eq. (31) have been checked with an independent numerical calculation, as explained in Section IV.

A. Boundary integral B 1
The boundary integral B 1 is equal to the phase-space volume in the limit z → 1. The phase-space volume is simple enough to be computed directly, keeping the exact dependence on s, t, z and . The relevant integral is given by It is convenient to introduce the Sudakov decomposition as in Eq. (25) for all gluon momenta.
The change of variables from gluon momenta components to Sudakov parameters leads to We can easily integrate over k i⊥ thanks to the on-shell delta function. We obtain We re-scale the Sudakov parameters α i = (1 − z) α i and β i = t/(sz) β i , removing the dependencies on z and t. The six remaining integrations factorize into a product of parametric integrals, each of them of the form As a result, we obtain The boundary integral B 1 is extracted from this expression via Extracting the z-dependence and the common -dependent pre-factor, as in Eq. (30), we find where is the integration constant quoted in Eq. (31).

B. Boundary integral B 8
Another relatively simple example is the boundary integral B 8 , which contains two additional propagators compared to B 1 . Its integral representation reads A Sudakov decomposition of the gluon momenta would lead in this case to a non-trivial dependence on the angle between k 1 and k 3 through the propagator 1/k 2 13 . This situation can be avoided, at least for some of the boundary integrals, by introducing an auxiliary momentum Q that has the effect of factoring out an ordinary phase-space integral.
In the case of B 8 , it is convenient to choose Q = k 13 and write where B 8 (Q 2 , Q ·p, k 2 ·p) is the following integral The result in Eq. (42)

C. Boundary integral B 9
It is not always possible to avoid non-trivial angular integrations as in the previous example; this happens in the integrals with multiple propagators of the type 1/k 2 ij . As an example, we consider the following boundary integral To calculate it, we use the Sudakov decomposition for each of the gluon momenta k i , c.f.
Eq. (25). We then remove the on-shell delta functions δ(k 2 i ) by integrating over |k i,⊥ |. Upon re-scaling α i → (1 − z)α i , we obtain an overall factor (1 − z) −1−3 while at the same time the parameters α i become constrained by δ(α 123 − 1) and are thus placed on an equal footing with the β-parameters.
Although the on-shell delta function δ(k 2 i ) fixes the length of the vector k i⊥ , its direction remains arbitrary and has to be integrated over. The required angular integrations are non-trivial. For example, the propagator 1/k 2 12 leads to an angular integral where the function F (x, y) reads Note that this function is symmetric, i.e. F (x, y) = F (y, x). The propagator 1/k 2 13 produces a similar function upon integration over the directions of k 3⊥ . As a result, we obtain where the parametric integral X 9 is given by We can use the transformation [55] 2 F 1 a, b, 2b, that simplifies the argument of the hypergeometric function in Eq. (46). We find , for x < y , Since the transformation Eq. (49) is only valid if the argument of the hypergeometric function is smaller than one, we must split the integration region into four pieces, according to the cases α 2 β 1 ≶ α 1 β 2 and α 3 β 1 ≶ α 1 β 3 . Due to the symmetry of the integrand under the simultaneous interchange of subscripts 2 ↔ 3 and α ↔ β, two of these contributions happen to be identical. The calculation of the remaining two contributions is quite similar, so that it is sufficient to describe the calculation of one of them.
Consider the contribution to X 9 that originates from the integration region defined by the conditions α 2 β 1 > α 1 β 2 and α 3 β 1 < α 1 β 3 ; we will call it X (a) 9 . After applying the transformations in Eq. (50), we find Upon changing variables β 2 → r = α 1 β 2 /(α 2 β 1 ) and β 3 → µ = α 3 β 1 /(α 1 β 3 ) and integrating over β 1 to remove the delta function, we obtain In Eq. (52) we have re-written the hypergeometric functions to make them regular in the r → 1 and µ → 1 limits. We proceed by integrating out α 2 and change the integration The integral in Eq. (53) is singular; the overlapping logarithmic singularities appear at r = 1, µ = 1, ξ = 0 and f ∈ {0, ∞}. These singularities are disentangled by performing suitable (iterated) subtractions, after which the resulting integrals are carried out using the program HyperInt [56]. The other independent contributions are obtained in a similar fashion. Upon adding all the contributions, we obtain the result for X 9 , The boundary constant C 9 is easily obtained from this result.

D. Boundary integral B 10
The most challenging boundary integrals involve the propagator 1/k 2 123 . Their computation requires a different approach because the Sudakov decomposition of the gluon momenta does not sufficiently simplify them. To compute these integrals, we set up additional differential equations for suitable parts of their integrands, and determine the boundary constants from these differential equations. As an example, we consider the last boundary integral The z-dependence is again extracted by the re-scaling k 123 ) and integrate out the momentum k 3 to obtain In Eq. (56) we introduced the integral F 10 , that we will explicitly compute. As we indicated in Eq. (57), F 10 depends only on Q 2 since all other kinematic invariants are fixed, c.f. Eq. (56). The variable Q 2 satisfies the constraint 0 ≤ Q 2 ≤ 1. Indeed, the lower boundary appears because Q 2 = k 2 123 ≥ 0, while a Sudakov decomposition of the momentum Q gives Q 2 = 1 − Q 2 ⊥ ≤ 1. The computation of F 10 proceeds through the method of differential equations. We take the derivative of F 10 with respect to Q 2 , at fixed Q · p and Q ·p, and, after promoting the delta functions to cut propagators and performing an integration-by-parts reduction, we write the result in terms of a set of masters integrals. Performing the same steps for the other integrals that contribute to dF 10 (Q 2 )/dQ 2 , we arrive at a closed system of differential equations that contains five master integrals. They are The differential equations for these master integrals can be easily solved, but five integration constants need to be determined. We obtain these integration constants by various means. One constant follows from the calculation of J 1 (Q 2 ) at Q 2 = 1, which is closely related to the phase-space integral B 1 . Constraints on the remaining integration constants are obtained from the analysis of the solutions to the differential equations in the limits Q 2 → 0 and Q 2 → 1. For example, we require that J 5 (Q 2 ) = F 10 (Q 2 ) does not have a hard region (Q 2 ) 0 , because the integral in Eq. (56) would otherwise be ill-defined. We also find that the (Q 2 ) 0 branch of J 3 (Q 2 ) vanishes, that the (Q 2 ) − branch of J 4 (Q 2 ) vanishes and that the (1 − Q 2 ) 0 branch of J 2 (Q 2 ) is given by Putting all this information together gives us the result for F 10 (Q 2 ). Using it in Eq. (56) and integrating over Q, we obtain the boundary integral B 10 . It reads where X 10 is given by the following expression The constant C 10 is then easily extracted.

IV. NUMERICAL CHECKS OF MASTER INTEGRALS
We have performed several checks to ensure the correctness of the master integrals computed in the previous section. First, we inserted the master integrals into the system of differential equations from which they were derived and checked that the differential equations are indeed satisfied. Second, some of the boundary constants for z → 1 have been computed in several different ways. Nevertheless, a completely independent check of the integrals is desirable. Unfortunately, contrary to standard Feynman integrals, there exists no automated code to evaluate phase-space integrals numerically and therefore we have to proceed differently.
In Ref. [33] we have considered similar phase-space integrals, albeit in a situation where one of the gluons was virtual and two were real. The double-virtual single-real master integrals in that paper were checked numerically using the Mellin-Barnes (MB) integration, following the discussion in Ref. [57]. We employ the same approach to check the triple-real integrals computed in the current paper; since there are significant similarities with the calculation described in Ref. [33], we only give a short overview of the steps required for the numerical checks.
For reasons explained in Ref. [33], in order to perform the numerical evaluation of the phase-space integrals, it is preferable to consider the decay process q * → q + ggg instead of the production process q → q * + ggg. We accomplish this by formally changing the four-momenta p → −p,p → −p in the definition of the master integrals. We obtain where we introduced κ = t/z. It follows from the above equation that we need to take z ≥ 1 and t ≤ 0, or otherwise the integrals would identically vanish. The analytic expression for the integral in the decay kinematics, that we refer to as I decay (κ, z), can be determined from the solutions in the production channel I production (κ, z) by an analytic continuation to the region z ≥ 1 and κ ≤ 0. Note that since the propagators are positive definite both in the production and in the decay kinematics, both integrals I production (κ, z) and I decay (κ, z) are real-valued. This consideration provides a useful constraint on the results of the analytic continuation.
As the next step, we set κ = z − 2 and write Note that in the second step in Eq. (63) we used the fact that I decay (κ, z) vanishes outside the region κ ≤ 0, z ≥ 1.
Eq. (63) can be used to check our integrals numerically. Indeed, on the one hand, the first integral in Eq. (63) can be calculated using the analytic solution I production (κ, z), continued to the decay region z ≥ 1, κ ≤ 0. On the other hand, W can be written as a MB integral, following the discussion in Ref. [57]. Indeed, we consider the right-hand side of Eq. (63) and write the integral as where D j are the propagators of the particular integral, c.f. Table I. To proceed further, we may use the Mellin-Barnes representation to re-write propagators of the form into integrals of products of k i · k j , p · k i and p · k j . Upon doing so, we obtain integrals that are identical to the ones studied in Ref. [  To illustrate the usefulness of the integrals presented in this paper, we construct the RRR contribution to the I qq matching coefficient at N 3 LO in QCD in the large-N c and the large- The resulting s-independent contributions can be written as where g s is the strong coupling constant. The subscript i is either N 3 c to indicate the leadingcolor contribution, or N f to indicate the contribution proportional to N f . The t-dependence factorizes by construction, since we computed the leading contribution in the collinear limit.
In fact, this factor will eventually be expanded in terms of plus distributions,

VI. CONCLUSIONS
In this paper, we computed the master integrals required to describe the real-emission contribution to the matching coefficient of a quark beam function at N 3 LO in QCD due to the splitting of an incoming quark q into a virtual quark of the same flavor and three collinear gluons, q → q * + ggg. We used reverse unitarity and integration-by-parts identities to derive differential equations satisfied by the master integrals. We solved the differential equations and fixed the boundary conditions for the master integrals using both regularity requirements and the explicit computation of a small subset of integrals in the soft limit.
Our final results for the master integrals are expressed in terms of Goncharov polylogarithms up to weight six.
The master integrals computed in this paper allow us to obtain the triple-real contribution to the matching coefficient I qq in the large-N c and large-N f limits. To extend this calculation to include terms that are sub-leading in N c , we have to account for processes where an incoming quark q splits into a quark-antiquark pair of the same flavor and a gluon, q → q * + qqg. The contribution of this process to I qq requires additional master integrals. We expect their computation to be feasible using the techniques described in this paper.
As we pointed out in the Introduction, there are three N 3 LO QCD contributions to I qq , the triple-real, the double-real single-virtual and the single-real double-virtual, that need to be calculated. We studied the double-real single-virtual contribution in Ref. [33] and the triple-real contribution in this paper. The so far unattended contribution is the single-real double-virtual one; its computation will require us to understand how to compute a massless two-loop three-point function in an axial gauge. Although such a computation appears to be quite challenging, we believe that it can be dealt with using calculational methods developed both in this paper and in Ref. [33].
Using the result for the differential equations for the master integrals, we find that F 1 and F 2 satisfy the following homogeneous differential equations The solutions to these equations are In the limit z → 1 these solutions for F 1 and F 2 behave as (1 − z) 0 and (1 − z) − . However, we have argued in the main body of the paper that all master integrals in the soft z → 1 limit should be proportional to (1 − z) n−3 for some integer n. The only way to make this scaling compatible with Eq. (B5) is to choose c 1 ( ) = c 2 ( ) = 0 which implies that F 1,2 vanish identically. This proves the identities among master integrals shown in Eqs. (16) and (17).