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

TL;DR

This paper calculates complex three-loop integrals essential for advancing the precision of quark beam function predictions in QCD at N$^3$LO, focusing on the triple-real emission process.

Contribution

It provides the first computation of three-loop master integrals for the triple-real contribution to the quark beam function at N$^3$LO in QCD.

Findings

Computed three-loop master integrals for triple-real emission

Enabled next-to-next-to-next-to-leading order (N$^3$LO) calculations of the quark beam function

Contributed to precision QCD predictions for collider physics

Abstract

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.

Tables

Table 1. Table 1 : The inverse propagators D i subscript 𝐷 𝑖 D_{i} for each of the 19 topologies A 1 ​ … ​ A 19 subscript 𝐴 1 … subscript 𝐴 19 A_{1}\dotsc A_{19} . Here we use the shorthand notation k i ​ j = k i + k j subscript 𝑘 𝑖 𝑗 subscript 𝑘 𝑖 subscript 𝑘 𝑗 k_{ij}=k_{i}+k_{j} and k i ​ j ​ ℓ = k i + k j + k ℓ subscript 𝑘 𝑖 𝑗 ℓ subscript 𝑘 𝑖 subscript 𝑘 𝑗 subscript 𝑘 ℓ k_{ij\ell}=k_{i}+k_{j}+k_{\ell} .

top	$D_1$	$D_2$	$D_3$	$D_4$	$D_5$	$D_6$	$D_7$
$A_1$	${\left(p-k_1\right)}^2$	${\left(p-k_2\right)}^2$	$\left(p-k_{12}\right)$	${\left(p-k_{13}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_2$
$A_2$	${\left(p-k_1\right)}^2$	${\left(p-k_2\right)}^2$	$\left(p-k_{12}\right)$	${\left(p-k_{13}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_3$
$A_3$	${\left(p-k_1\right)}^2$	${\left(p-k_2\right)}^2$	$\left(p-k_{12}\right)$	${\left(p-k_{13}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_2$	$\overline{p}\cdot k_3$
$A_4$	$k_{12}^2$	${\left(p-k_1\right)}^2$	$\left(p-k_{12}\right)$	${\left(p-k_{13}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_2$
$A_5$	$k_{12}^2$	${\left(p-k_1\right)}^2$	$\left(p-k_{12}\right)$	${\left(p-k_{13}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_3$
$A_6$	$k_{12}^2$	${\left(p-k_1\right)}^2$	$\left(p-k_{12}\right)$	${\left(p-k_{13}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_2$	$\overline{p}\cdot k_3$
$A_7$	${\left(p-k_1\right)}^2$	${\left(p-k_2\right)}^2$	$\left(p-k_{13}\right)$	${\left(p-k_{23}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_2$
$A_8$	$k_{12}^2$	$k_{13}^2$	${\left(p-k_2\right)}^2$	${\left(p-k_{23}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_2$
$A_9$	$k_{12}^2$	$k_{13}^2$	${\left(p-k_2\right)}^2$	${\left(p-k_{23}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_3$
$A_{10}$	$k_{12}^2$	${\left(p-k_1\right)}^2$	${\left(p-k_3\right)}^2$	${\left(p-k_{13}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_2$
$A_{11}$	$k_{12}^2$	${\left(p-k_1\right)}^2$	${\left(p-k_3\right)}^2$	${\left(p-k_{13}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_2$	$\overline{p}\cdot k_3$
$A_{12}$	$k_{12}^2$	$k_{13}^2$	$k_{123}^2$	${\left(p-k_2\right)}^2$	${\left(p-k_{12}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_2$
$A_{13}$	$k_{12}^2$	$k_{13}^2$	$k_{123}^2$	${\left(p-k_2\right)}^2$	${\left(p-k_{12}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_3$
$A_{14}$	${\left(p-k_1\right)}^2$	${\left(p-k_2\right)}^2$	${\left(p-k_{12}\right)}^2$	${\left(p-k_{13}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_{13}$
$A_{15}$	$k_{12}^2$	${\left(p-k_1\right)}^2$	${\left(p-k_{12}\right)}^2$	${\left(p-k_{13}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_{12}$
$A_{16}$	$k_{12}^2$	$k_{13}^2$	${\left(p-k_2\right)}^2$	${\left(p-k_{23}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_{12}$
$A_{17}$	$k_{12}^2$	$k_{13}^2$	${\left(p-k_2\right)}^2$	${\left(p-k_{23}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_{13}$
$A_{18}$	$k_{12}^2$	$k_{13}^2$	${\left(p-k_{12}\right)}^2$	${\left(p-k_{13}\right)}^2$	${\left(p-k_{123}\right)}^2$	$\overline{p}\cdot k_1$	$\overline{p}\cdot k_{12}$
$A_{19}$	$k_{12}^2$	$k_{13}^2$	$k_{123}^2$	${\left(p-k_2\right)}^2$	${\left(p-k_{12}\right)}^2$	$\overline{p}\cdot k_{13}$	$\overline{p}\cdot k_3$

Equations

$$\mathcal{T}=\sum _{j=1}^n{\min }_{i\in \{1,2\}}\left[\frac{2p_i\cdot k_j}{Q_i}\right].$$

$${\sigma }_{pp\to V+X}^{{\mathrm{N}}^3\mathrm{LO}}={\sigma }_{pp\to V+X}^{{\mathrm{N}}^3\mathrm{LO}}\left(\mathcal{T}\le {\mathcal{T}}_0\right)+{\sigma }_{pp\to V+X}^{\mathrm{NNLO}}\left(\mathcal{T}>{\mathcal{T}}_0\right).$$

$$\underset{{\mathcal{T}}_0\to 0}{\lim }d{\sigma }_{pp\to V+X}^{{\mathrm{N}}^3\mathrm{LO}}\left(\mathcal{T}\le {\mathcal{T}}_0\right)\sim B\otimes B\otimes S\otimes d{\sigma }_{pp\to V}^{{\mathrm{N}}^3\mathrm{LO}},$$

$$B_i=\sum _{\mathrm{partons}\text{leavevmode}\mathrm{j}}{I}_{ij}\otimes f_j,\text{mbox}wherei,j=\{q,\bar{q},g\}.$$

$$q(p)\to q^{\ast }(p^{\ast })+g(k_1)+g(k_2)+g(k_3),p^{\ast }=p-k_1-k_2-k_3.$$

$$p^{\ast \mu }=z{p}^{\mu }+y{\bar{p}}^{\mu }+k_{\perp }^{\mu },k_{123}^{\mu }=(1-z)p^{\mu }-y{\bar{p}}^{\mu }-k_{\perp }^{\mu }.$$

$$t=-zy2p\cdot \bar{p}=z2p\cdot k_{123}.$$

$$s(1-z)=2\bar{p}\cdot k_{123},\text{mbox}withs=2p\cdot \bar{p}.$$

$$I_{qq}(t,s,z)\sim \int$$

$$\times P_{qq}(p,\bar{p},\{k_i\}).$$

$$\mathrm{Tr}[\cdots {\hat{p}}^{\ast }\hat{\mathcal{H}}{\hat{p}}^{\ast }\cdots ]\to \mathrm{Tr}[\cdots {\hat{p}}^{\ast }\hat{\bar{p}}{\hat{p}}^{\ast }\cdots ],$$

$$\sum _{\mathrm{pol}}{ϵ}_i^{\mu }(k_i){\left({ϵ}_i^{\nu }(k_i)\right)}^{\ast }=-g^{\mu \nu }+\frac{k_i^{\mu }{\bar{p}}^{\nu }+k_i^{\nu }{\bar{p}}^{\mu }}{k_i\cdot \bar{p}},\text{mbox}fori=1,2,3.$$

$$\mathcal{I}=\int \prod _{i=1}^3\left[\frac{d^d{k}_i}{(2\pi {)}^{d-1}}{\delta }^{+}(k_i^2)\right]\frac{\delta \left(p\cdot k_{123}-\frac{t}{2z}\right)\delta \left(\bar{p}\cdot k_{123}-\frac{s(1-z)}{2}\right)\mathcal{N}}{D_1^{n_1}\cdots D_t^{n_t}}.$$

$$\frac{1}{k_1\cdot \bar{p}{k}_2\cdot \bar{p}{k}_3\cdot \bar{p}}=\frac{2}{1-z}\left[\frac{1}{k_1\cdot \bar{p}{k}_2\cdot \bar{p}}+\frac{1}{k_1\cdot \bar{p}{k}_3\cdot \bar{p}}+\frac{1}{k_2\cdot \bar{p}{k}_3\cdot \bar{p}}\right],$$

$$\displaystyle\mathcal{I}^{\rm top}_{n_{1},n_{2},n_{3},n_{4},n_{5},n_{6},n_{7}}=\int\mathfrak{D}^{d}k_{1}\mathfrak{D}^{d}k_{2}\mathfrak{D}^{d}k_{3}\,\frac{\delta\big{(}p\cdot k_{123}-\frac{t}{2z}\big{)}\delta\big{(}\bar{p}\cdot k_{123}-\frac{s(1-z)}{2}\big{)}}{D_{1}^{n_{1}}D_{2}^{n_{2}}D_{3}^{n_{3}}D_{4}^{n_{4}}D_{5}^{n_{5}}D_{6}^{n_{6}}D_{7}^{n_{7}}}\,,$$

$${\mathfrak{D}}^d{k}_i=\frac{d^d{k}_i}{(2\pi {)}^{d-1}}{\delta }^{+}(k_i^2).$$

$${\mathcal{I}}_{1,1,1,0,1,0,1}^{A_{10}}+\frac{1}{2}\frac{(1-z)}{z}{\mathcal{I}}_{1,1,0,0,1,1,1}^{A_6}=0,$$

$${\mathcal{I}}_{0,0,1,0,0,1,0}^{A_1}+\frac{2(1-2ϵ)}{ϵ}\frac{1}{z}{\mathcal{I}}_{0,0,1,0,1,0,0}^{A_{12}}=0.$$

$$\displaystyle\begin{aligned} \mathcal{I}^{\rm top}_{n_{1},n_{2},n_{3},n_{4},n_{5},n_{6},n_{7}}&=\int\mathfrak{D}^{d}k_{1}\mathfrak{D}^{d}k_{2}\mathfrak{D}^{d}k_{3}\,\frac{\delta\big{(}p\cdot k_{123}-\frac{t}{2z}\big{)}\delta\big{(}\bar{p}\cdot k_{123}-\frac{s(1-z)}{2}\big{)}}{D_{1}^{n_{1}}D_{2}^{n_{2}}D_{3}^{n_{3}}D_{4}^{n_{4}}D_{5}^{n_{5}}D_{6}^{n_{6}}D_{7}^{n_{7}}}\\ &=s^{-\mathcal{M}}\,t^{\mathcal{N}}\int\mathfrak{D}^{d}k_{1}\mathfrak{D}^{d}k_{2}\mathfrak{D}^{d}k_{3}\,\frac{\delta\big{(}p\cdot k_{123}-\frac{1}{2}\big{)}\delta\big{(}\bar{p}\cdot k_{123}-\frac{(1-z)}{2}\big{)}}{D_{1}^{n_{1}}D_{2}^{n_{2}}D_{3}^{n_{3}}D_{4}^{n_{4}}D_{5}^{n_{5}}D_{6}^{n_{6}}D_{7}^{n_{7}}}\,\,,\end{aligned}$$

$$\frac{d}{dz}\vec{\mathcal{I}}(z,ϵ)=\hat{A}(z,ϵ)\vec{\mathcal{I}}(z,ϵ).$$

$$z=\frac{(1+x{)}^2}{x}.$$

$$\frac{d}{dx}{\vec{\mathcal{I}}}_{\text{can}}(x,ϵ)=ϵ\left(\sum _{x_0}\frac{{\hat{A}}_{x_0}}{x-x_0}\right){\vec{\mathcal{I}}}_{\text{can}}(x,ϵ).$$

$${\vec{\mathcal{I}}}_{\text{can}}(x,ϵ)=\hat{M}(x,ϵ)\vec{B}(ϵ),$$

$$M_{ij}(x,ϵ)=\sum _{k\ge 0}\sum _{\vec{w}\in W(k)}{c}_{i,j,k,\vec{w}}{ϵ}^{k}G(\vec{w};x).$$

$$G(w_1,w_2,\dots ,w_n;x)=\underset{0}{\overset{x}{\int }}dt\frac{G(w_2,w_3,\dots ,w_n;t)}{t-w_1}.$$

$$k_i^{\mu }={\alpha }_i{p}^{\mu }+{\beta }_i{\bar{p}}^{\mu }+k_{i\perp }^{\mu }.$$

$$\frac{1}{(p-k_{ij}{)}^2}=\frac{1}{k_{ij}^2-2k_{ij}\cdot p},$$

$$k_{ij}^2={\alpha }_i{\beta }_j+{\alpha }_j{\beta }_i-2\sqrt{{\alpha }_i{\beta }_j{\alpha }_j{\beta }_i}\cos {\theta }_{ij}\text{and}2k_{ij}\cdot p={\beta }_i+{\beta }_j={\beta }_{ij},$$

$$\frac{1}{(p-k_{ij}{)}^2}\to \frac{1}{-2k_{ij}\cdot p},$$

$$B_1$$

$$B_2$$

$$B_3$$

Full text

††arxiv: 1904.02433

1 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 Aad:2012tfa ; Chatrchyan:2012xdj 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. Aad:2016naf ; Chatrchyan:2014mua ; Aaboud:2018hip ; Sirunyan:2018goh ; Aaboud:2019gxl ; Sirunyan:2019bez ; Aaboud:2019pkc . 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 Kubar:1980zv ; Nason:1987xz ; Ohnemus:1991gb ; Frixione:1992pj ; Ossola:2006us ; Berger:2008sj ; Cascioli:2011va ; Alwall:2014hca ; Actis:2016mpe and, more recently, to next-to-next-to-leading order (NNLO) in QCD.111 At least as far as processes with relatively simple final states are concerned. 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. N3LO QCD) are warranted. This may happen for several reasons. Indeed, in some cases the convergence of the perturbative expansion in the strong coupling constant $\alpha_{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 N3LO QCD corrections are crucial to stabilize theoretical uncertainties at the few percent level Anastasiou:2015ema . 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 N3LO calculations have recently culminated in the computation of the N3LO QCD corrections to Higgs boson production in gluon fusion at the LHC Anastasiou:2015ema ; Anastasiou:2016cez ; Mistlberger:2018etf ; Dulat:2018bfe . 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 Anastasiou:2002yz to simplify them.222Recently, an approximated N3LO differential calculation for Higgs production has been completed using the $q_{T}$-subtraction formalism Cieri:2018oms . Although the calculation of the N3LO 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 proton-proton collisions $pp\to V$. Quite generally, the description of this process at N3LO in QCD requires the knowledge of the NNLO QCD corrections to the production of $V$ together with an additional QCD jet, $pp\to V+j$. The difference between $pp\to V+j$ at NNLO QCD and $pp\to V$ at N3LO QCD is that the jet in the former case can become unresolved and that the virtual corrections to $pp\to V$ have no counterpart in the $pp\to V+j$ calculation. Since the difference between the two calculations appears in the kinematic regions where the color-singlet 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\to 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$ Catani:2007vq . More recently, the so-called $N$-jettiness observable Stewart:2010tn ; Gaunt:2015pea ; Boughezal:2015dva 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 [math]-jettiness, which is required to describe the inclusive production of a color-singlet final state.

To this end, we consider the process $pp\to 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 [math]-jettiness variable as

[TABLE]

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 [math]-jettiness variable $\mathcal{T}$ has two important properties that allow one to use it as a slicing variable. Indeed, it follows from the definition eq. 1 that $\mathcal{T}=0$ in the absence of resolved QCD radiation, i.e. for the process $pp\to V$. However, in the presence of any resolved QCD radiation one finds that $\mathcal{T}>0$. We can therefore introduce a cut-off $\mathcal{T}_{0}$ and divide the phase space for $V+X$ into two disjoint parts. We write schematically

[TABLE]

Note the $\rm NNLO$ subscript in the second term on the right-hand side in eq. 2; the reason for its appearance is that by imposing the $\mathcal{T}>\mathcal{T}_{0}$ constraint, we exclude the situation where all final-state partons become unresolved so that the calculation for $\mathcal{T}>\mathcal{T}_{0}$ reduces to the computation of the NNLO QCD corrections to $pp\to V+j$. Such calculations have already been performed for a variety of final states and we consider them to be known Boughezal:2015dra ; Ridder:2015dxa ; Boughezal:2015aha ; Boughezal:2015dva ; Chen:2016zka ; Gehrmann-DeRidder:2017mvr .

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 N3LO calculation itself. However, for 0-jettiness, this does not happen. Indeed, it was shown in Ref. Stewart:2010tn that the cross section for $pp\to V+X$ simplifies substantially in the limit $\mathcal{T}\to 0$ and can be written as a convolution of the hard cross section for $pp\to V$ with the so-called beam and soft functions Berger:2010xi ; Gaunt:2014xga ; Gaunt:2014cfa . The cross section reads

[TABLE]

where the two functions $B$ stand for the beam functions associated with each of the initial-state partons and $S$ represents the soft function. The general factorization formula for $N$-jettiness was originally derived in SCET Bauer:2000ew ; Bauer:2000yr ; Bauer:2001ct ; Bauer:2001yt ; Bauer:2002nz . 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-$\mathcal{T}$ limit.

The cross-section formula eq. 3 implies that, in order to employ the [math]-jettiness slicing to compute the N3LO corrections to $pp\to 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}$

[TABLE]

The computation of the N3LO 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 Melnikov:2018jxb , 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 2 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 3 how these integrals can be computed using the method of differential equations. In section 4, we explain how the calculation was validated and we present our final results in section 5. We conclude in section 6. 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.

2 Matching coefficient

In this section we discuss how to compute the N3LO 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.

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 Stewart:2010tn . 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 Ritzmann:2014mka . 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. Catani:1999ss .

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}\to q_{i}^{*}+ggg$ and $q_{i}\to q_{i}^{*}+q_{j}\bar{q}_{j}+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\neq 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}\to\infty$ limit, where $N_{c}$ is the number of colors, and in the $N_{f}\to\infty$ 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. Melnikov:2018jxb , 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^{*}$

[TABLE]

As we already explained, the relevant contribution to the matching coefficient is obtained by integrating the $q\to 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

[TABLE]

In eq. 6, we used $k_{123}=\sum\limits_{i=1}^{3}k_{i}^{\mu}$. We also introduced a light-cone momentum $\bar{p}$, which is complementary to $p$ so that $\bar{p}^{2}=0$ and $p\cdot k_{\perp}=\bar{p}\cdot k_{\perp}=0$. The emitted gluons are on the mass shell, i.e. $k_{i}^{2}=0$ for $i=1,2,3$. With these definitions we have $y=-(p\cdot k_{123})/(p\cdot\bar{p})$.

We now introduce the transverse virtuality $t=-((p^{*})^{2}-k_{\perp}^{2})$ and, using the above results, write it as

[TABLE]

Note that, in the case of collinear emissions, $t\sim\mathcal{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

[TABLE]

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

[TABLE]

In eq. 9 the integrand $P_{qq}(p,\bar{p},\{k_{i}\})$ describes the $q\to q^{*}+ggg$ splitting function. Below we explain how to compute it.

As described in Ref. Catani:1999ss , the $q\to 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\to 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 vector $\bar{p}$ from eq. 6 is selected as the corresponding gauge-fixing vector. Squaring the amplitude, we produce a Dirac trace of the form $\mathrm{Tr}[\dotsm\,\hat{p}^{*}\,\hat{\mathcal{H}}\,\hat{p}^{*}\,\dotsm]$, where $\hat{p}^{*}=\gamma^{\mu}p^{*}_{\mu}$ and $p^{*}$ is the momentum that enters the hard scattering process. The Dirac matrix $\hat{\mathcal{H}}$ 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

[TABLE]

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)\to q^{*}(p^{*})+g(k_{1})+g(k_{2})+g(k_{3})$ with QGRAF Nogueira:1991ex . We perform the relevant Dirac and Lorentz algebra in FORM Vermaseren:2000nd and Mathematica in two independent implementations. Since we work in the axial gauge with the gauge-fixing vector $\bar{p}$, the sum over polarizations for a gluon with momentum $k_{i}$ reads

[TABLE]

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

[TABLE]

Here, $\mathcal{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 Anastasiou:2002yz , 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) Tkachov:1981wb ; Chetyrkin:1981qh ; Laporta:2001dd . As is often the case when dealing with phase-space integrals in the framework of reverse unitarity, this 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$ and $\bar{p}$ 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 [math]-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}\cdot\bar{p}$ with $i=1,2,3$. However, the [math]-jettiness constraint in eq. 12 implies that the three propagators $1/k_{i}\cdot\bar{p}$ are not linearly independent. Indeed, we can write

[TABLE]

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 1. We performed the reduction to master integrals using Reduze vonManteuffel:2012np and KIRA Maierhoefer:2017hyi , 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

[TABLE]

where $d=4-2\epsilon$ and the subscript ‘top’ indicates one of the topologies in table 1 where the inverse propagators $D_{i}$ for each topology are defined. The integration measure for each final-state particle reads

[TABLE]

We use these notations to present the list of master integrals in appendix A.

While the set of master integrals shown in eq. 72 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

[TABLE]

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.

3 Master integrals

The master integrals defined in eq. 14 depend on $t$, $z$ and $s=2p\cdot\bar{p}$. However, the dependence on $s$ and $t$ is trivial. This becomes manifest after the simultaneous re-scaling $k_{i}\to k_{i}\,\sqrt{t}$, $p\to p\,\sqrt{t}$ and $\bar{p}\to\bar{p}\,s/\sqrt{t}$. 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

[TABLE]

where $\mathcal{M}=1+n_{6}+n_{7}$ and $\mathcal{N}=5-3\epsilon-\sum\limits_{i=1}^{5}n_{i}$. As a consequence, we can set $s=t=1$ everywhere and focus only on the $z$-dependence of the master integrals.

We determine the $z$-dependence of the master integrals with the method of differential equations Kotikov:1990kg ; Remiddi:1997ny ; Gehrmann:1999as . To this end, we differentiate each of the master integrals with respect to $z$ and express the result in terms of master integrals using integration-by-parts identities. We collect the master integrals into a vector $\vec{\mathcal{I}}(z,\epsilon)$ and write the resulting closed system of differential equations as

[TABLE]

The entries of the matrix $\hat{A}(z,\epsilon)$ are rational functions of $z$ and $\epsilon$.

The complexity of these differential equations depends strongly on the explicit form of the matrix $\hat{A}(z,\epsilon)$, 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 Kotikov:2012ac ; Henn:2013pwa ; Lee:2014ioa , $\hat{A}(z,\epsilon)=\epsilon\sum\limits_{z_{0}}\frac{\hat{A}_{z_{0}}}{z-z_{0}}$. Note that the matrices $\hat{A}_{z_{0}}$ should be both $z$- and $\epsilon$-independent. 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 $\sqrt{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

[TABLE]

Having removed all square roots, we can construct the appropriate transformation $\vec{\mathcal{I}}(x,\epsilon)=\hat{T}(x,\epsilon)\,\vec{\mathcal{I}}_{\text{can}}(x,\epsilon)$ with the program Fuchsia Gituliar:2017vzm . As a result, we find

[TABLE]

The differential equations have singularities drawn from the list $x_{0}\in\{-1,0,1,R_{1}^{\pm},R_{2}^{\pm},R_{3}^{\pm}\}$, which in turn correspond to singularities in $z$ given by $z_{0}\in\{0,\infty,4,1,-1,2\}$. The symbols $R_{1}^{\pm}$, $R_{2}^{\pm}$ and $R_{3}^{\pm}$ represent the two roots of each of the quadratic polynomials $P_{1}=1+x+x^{2}$, $P_{2}=1+3x+x^{2}$ and $P_{3}=1+x^{2}$, respectively.

It is convenient to solve the system of differential equations eq. 21 expanding around $\epsilon=0$. We write $\vec{\mathcal{I}}_{\text{can}}(x,\epsilon)$ as

[TABLE]

where $\vec{B}(\epsilon)$ are the integration constants. The $x$-dependence resides solely in the matrix $\hat{M}(x,\epsilon)$, whose elements have the form

[TABLE]

We calculate the sum over $k$ up to and including $k=6$, corresponding to ${\mathcal{O}}(\epsilon^{6})$, which is the highest order that will contribute to the finite part of the matching coefficient in the $\epsilon\to 0$ limit. For a given $k$, the inner sum in eq. 23 runs over $W(k)$, containing all vectors $\vec{w}$ of the length $k$ with components drawn from the set of roots $\{-1,0,1,R_{1}^{\pm},R_{2}^{\pm},R_{3}^{\pm}\}$. The functions $G(\vec{w};x)$ are the Goncharov polylogarithms Kummer ; Goncharov:1998kja ; Remiddi:1999ew ; Goncharov:2001iea

[TABLE]

They can be evaluated numerically with the help of the program Ginac Vollinga:2004sn . Apart from the technical difficulty in handling large expressions, the construction of the matrix $\hat{M}(x,\epsilon)$ can be done in a relatively straightforward way.

On the contrary, the determination of the integration constants $\vec{B}(\epsilon)$ in eq. 22 is much less straightforward. We obtain them by analyzing the master integrals in the limit $z\to 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}\cdot 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_{ij}^{2}$ relative to $k_{ij}\cdot p$ follows from the following argument. Let us select a frame in which the external momenta are $p=\frac{1}{\sqrt{2}}(1,0,0,1)$ and $\bar{p}=\frac{1}{\sqrt{2}}(1,0,0,-1)$ and introduce a Sudakov decomposition of the gluon momentum $k_{i}$

[TABLE]

Since $\beta_{i}=k_{i}\cdot p=k_{i}^{0}p^{0}(1-\cos\theta_{i})$ and $\alpha_{i}=k_{i}\cdot\bar{p}=k_{i}^{0}\bar{p}^{0}(1+\cos\theta_{i})$, we conclude that all $\alpha$’s and $\beta$’s are positive definite. According to the phase-space constraints eq. 12, the sum $\alpha_{123}=\alpha_{1}+\alpha_{2}+\alpha_{3}$ goes to zero in the $z\to 1$ limit and, since all $\alpha$’s are positive, we conclude that each $\alpha_{i}$ goes to zero in that limit at least as fast as $\mathcal{O}(1-z)$. In contrast, the sum of the $\beta_{i}$’s is constrained to be equal to one, so that up to two of them could vanish at $z=1$. We write

[TABLE]

where we have used that $p^{2}=0$. In terms of the Sudakov parameters, $k_{ij}^{2}$ reads

[TABLE]

where we have used $k_{i}^{2}=k_{j}^{2}=0$. Assuming that, in the limit $z\to 1$, each $\alpha_{i}=\mathcal{O}(1-z)$ and each $\beta_{i}=\mathcal{O}(1)$ we find $k_{ij}^{2}=\mathcal{O}(1-z)$ and $2k_{ij}\cdot p=\mathcal{O}(1)$. Hence, we can neglect $k_{ij}^{2}$ relative to $2k_{ij}\cdot p$. The situation does not change, should any of the $\alpha_{i}$’s vanish faster than $\mathcal{O}(1-z)$. Another possibility is that both $\beta_{i}$ and $\beta_{j}$ vanish as $\mathcal{O}(1-z)$, such that $2k_{ij}\cdot p\to 0$. However, in that situation $k_{ij}^{2}$ scales as $\mathcal{O}((1-z)^{2})$ or faster, and we can again neglect it relative to $2k_{ij}\cdot p$. Therefore, a replacement

[TABLE]

is valid in the $z\to 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 \delta\big{(}k_{123}\cdot p-\tfrac{1}{2}\big{)}\delta\big{(}k_{123}\cdot\bar{p}-\tfrac{(1-z)}{2}\big{)} and, upon re-scaling the momenta as $k_{i}\to k_{i}\sqrt{1-z}$, $\bar{p}\to\bar{p}\sqrt{1-z}$ and $p\to p/\sqrt{1-z}$, we extract the overall $(1-z)$-dependence of the master integrals.

It follows that in the soft limit, each integral scales as $(1-z)^{n-3\epsilon}$ with an integer $n$ that is integral-dependent. Hence, all canonical master integrals should be free of logarithmic singularities as $z\to 1$, or equivalently as $x\to R_{1}^{\pm}$, beyond those that correspond to the expansion of $(1-z)^{-3\epsilon}$ in powers of $\epsilon$. 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\to 1$. These integrals read

[TABLE]

To present the results, it is convenient to extract the common $\epsilon$-dependent factor,

[TABLE]

where $\Omega_{n}=2\pi^{n/2}/\Gamma(n/2)$. With this normalization, the constants $\widetilde{C}_{i}$ read, up to weight six,

[TABLE]

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 4.

3.1 Boundary integral $B_{1}$

The boundary integral $B_{1}$ is equal to the phase-space volume in the limit $z\to 1$. The phase-space volume is simple enough to be computed directly, keeping the exact dependence on $s,t,z$ and $\epsilon$. The relevant integral is given by

[TABLE]

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

[TABLE]

We can easily integrate over $k_{i\perp}$ thanks to the on-shell delta function. We obtain

[TABLE]

We re-scale the Sudakov parameters $\alpha_{i}=(1-z)\widetilde{\alpha}_{i}$ and $\beta_{i}=t/(sz)\widetilde{\beta}_{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

[TABLE]

As a result, we obtain

[TABLE]

The boundary integral $B_{1}$ is extracted from this expression via

[TABLE]

Extracting the $z$-dependence and the common $\epsilon$-dependent pre-factor, as in eq. 30, we find

[TABLE]

where

[TABLE]

is the integration constant quoted in eq. 31.

3.2 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

[TABLE]

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_{13}^{2}$. 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

[TABLE]

where $\widetilde{B}_{8}(Q^{2},Q\cdot\bar{p},k_{2}\cdot\bar{p})$ is the following integral

[TABLE]

The result in eq. 42 is most easily obtained by computing the integral in the rest frame of the vector Q=\big{(}Q_{0},\vec{0}\,\big{)} and expressing the result of the integration in the Lorentz-invariant way by replacing $Q_{0}\bar{p}_{0}$ with $Q\cdot\bar{p}$ and $Q_{0}^{2}$ with $Q^{2}$. Upon inserting the result for the integral into eq. 41, one can proceed by introducing the Sudakov decomposition for the remaining momenta $k_{2}$ and $Q$. Carrying out the resulting parametric integrations yields the desired result

[TABLE]

3.3 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_{ij}^{2}$. As an example, we consider the following boundary integral

[TABLE]

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 $\delta(k_{i}^{2})$ by integrating over $|k_{i,\perp}|$. Upon re-scaling $\alpha_{i}\to(1-z)\alpha_{i}$, we obtain an overall factor $(1-z)^{-1-3\epsilon}$ while at the same time the parameters $\alpha_{i}$ become constrained by $\delta(\alpha_{123}-1)$ and are thus placed on an equal footing with the $\beta$-parameters.

Although the on-shell delta function $\delta(k_{i}^{2})$ fixes the length of the vector $\vec{k}_{i\perp}$, its direction remains arbitrary and has to be integrated over. The required angular integrations are non-trivial. For example, the propagator $1/k_{12}^{2}$ leads to an angular integral

[TABLE]

where the function $F(x,y)$ reads

[TABLE]

Note that this function is symmetric, i.e. $F(x,y)=F(y,x)$. The propagator $1/k_{13}^{2}$ produces a similar function upon integration over the directions of $\vec{k}_{3\perp}$. As a result, we obtain

[TABLE]

where the parametric integral $X_{9}$ is given by

[TABLE]

We can use the transformation abramowitzstegun

[TABLE]

that simplifies the argument of the hypergeometric function in eq. 46. We find

[TABLE]

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 $\alpha_{2}\beta_{1}\lessgtr\alpha_{1}\beta_{2}$ and $\alpha_{3}\beta_{1}\lessgtr\alpha_{1}\beta_{3}$. Due to the symmetry of the integrand under the simultaneous interchange of subscripts $2\leftrightarrow 3$ and $\alpha\leftrightarrow\beta$, 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 $\alpha_{2}\beta_{1}>\alpha_{1}\beta_{2}$ and $\alpha_{3}\beta_{1}<\alpha_{1}\beta_{3}$; we will call it $X_{9}^{(a)}$. After applying the transformations in eq. 50, we find

[TABLE]

Upon changing variables $\beta_{2}\to r=\alpha_{1}\beta_{2}/(\alpha_{2}\beta_{1})$ and $\beta_{3}\to\mu=\alpha_{3}\beta_{1}/(\alpha_{1}\beta_{3})$ and integrating over $\beta_{1}$ to remove the delta function, we obtain

[TABLE]

In eq. 52 we have re-written the hypergeometric functions to make them regular in the $r\to 1$ and $\mu\to 1$ limits. We proceed by integrating out $\alpha_{2}$ and change the integration variables $\alpha_{1}\to\xi=\alpha_{1}/(1-\alpha_{3})$ and $\alpha_{3}\to f=\alpha_{3}/(\mu(1-\alpha_{3}))$. We obtain

[TABLE]

The integral in eq. 53 is singular; the overlapping logarithmic singularities appear at $r=1,\mu=1,\xi=0$ and $f\in\{0,\infty\}$. These singularities are disentangled by performing suitable (iterated) subtractions, after which the resulting integrals are carried out using the program HyperInt Panzer:2014caa . The other independent contributions are obtained in a similar fashion. Upon adding all the contributions, we obtain the result for $X_{9}$,

[TABLE]

The boundary constant $\widetilde{C}_{9}$ is easily obtained from this result.

3.4 Boundary integral $B_{10}$

The most challenging boundary integrals involve the propagator $1/k_{123}^{2}$. 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

[TABLE]

The $z$-dependence is again extracted by the re-scaling $k_{i}\to k_{i}\sqrt{1-z}$, $\bar{p}\to\bar{p}\sqrt{1-z}$ and $p\to p/\sqrt{1-z}$. We introduce $1=\int d^{d}Q\,\delta^{d}(Q-k_{123})$ and integrate out the momentum $k_{3}$ to obtain

[TABLE]

In eq. 56 we introduced the integral $F_{10}$,

[TABLE]

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\leq Q^{2}\leq 1$. Indeed, the lower boundary appears because $Q^{2}=k_{123}^{2}\geq 0$, while a Sudakov decomposition of the momentum $Q$ gives $Q^{2}=1-Q_{\perp}^{2}\leq 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\cdot p$ and $Q\cdot\bar{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

[TABLE]

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}\to 0$ and $Q^{2}\to 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})^{-\epsilon}$ branch of $J_{4}(Q^{2})$ vanishes and that the $(1-Q^{2})^{0}$ branch of $J_{2}(Q^{2})$ is given by

[TABLE]

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

[TABLE]

where $X_{10}$ is given by the following expression

[TABLE]

The constant $\widetilde{C}_{10}$ is then easily extracted.

4 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\rightarrow 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. Melnikov:2018jxb 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. Anastasiou:2013srw . 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. Melnikov:2018jxb , we only give a short overview of the steps required for the numerical checks.

For reasons explained in Ref. Melnikov:2018jxb , in order to perform the numerical evaluation of the phase-space integrals, it is preferable to consider the decay process $q^{*}\to q+ggg$ instead of the production process $q\to q^{*}+ggg$. We accomplish this by formally changing the four-momenta $p\rightarrow-p,\bar{p}\rightarrow-\bar{p}$ in the definition of the master integrals. We obtain

[TABLE]

where we introduced $\kappa=t/z$. It follows from the above equation that we need to take $z\geq 1$ and $t\leq 0$, or otherwise the integrals would identically vanish. The analytic expression for the integral in the decay kinematics, that we refer to as $I_{\text{decay}}(\kappa,z)$, can be determined from the solutions in the production channel $I_{\text{production}}(\kappa,z)$ by an analytic continuation to the region $z\geq 1$ and $\kappa\leq 0$. Note that since the propagators are positive definite both in the production and in the decay kinematics, both integrals $I_{\text{production}}(\kappa,z)$ and $I_{\text{decay}}(\kappa,z)$ are real-valued. This consideration provides a useful constraint on the results of the analytic continuation.

As the next step, we set $\kappa=z-2$ and write

[TABLE]

Note that in the second step in eq. 63 we used the fact that $I_{\text{decay}}(\kappa,z)$ vanishes outside the region $\kappa\leq 0,z\geq 1$.

Equation 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_{\text{production}}(\kappa,z)$, continued to the decay region $z\geq 1$, $\kappa\leq 0$. On the other hand, $W$ can be written as a MB integral, following the discussion in Ref. Anastasiou:2013srw . Indeed, we consider the right-hand side of eq. 63 and write the integral as

[TABLE]

where $D_{j}$ are the propagators of the particular integral, c.f. table 1. To proceed further, we may use the Mellin-Barnes representation

[TABLE]

to re-write propagators of the form

[TABLE]

into integrals of products of $k_{i}\cdot k_{j}$, $p\cdot k_{i}$ and $p\cdot k_{j}$. Upon doing so, we obtain integrals that are identical to the ones studied in Ref. Anastasiou:2013srw and we can follow that reference to construct the Mellin-Barnes representation for those integrals. The resulting Mellin-Barnes integrals are finally computed numerically with the package MBtools MBTools . The two results for the quantity $W$ in eq. 63 must agree and we, therefore, get an indirect check of the results for the master integrals. We have performed this comparison for the master integrals and found agreement within the numerical errors. Furthermore, we note that we can use the same procedure to compute the soft limits of all integrals, checking the boundary values for all of them through weight six.

5 Results

The analytic expressions for the 91 master integrals $\mathcal{I}_{\vec{n}}^{\mathrm{top}}$ listed in eq. 72 are the main result of this paper. To present them we choose the normalization such that

[TABLE]

where the powers $\mathcal{M}$ and $\mathcal{N}$ depend on the index vector $\vec{n}$, as explained in eq. 18. With this normalization, the integral related to the phase-space volume becomes

[TABLE]

In general, the integrals ${\rm INT}(\mathrm{top},\vec{n})$ depend on the variable $x$, which is related to the longitudinal momentum fraction $z$ via eq. 20. We did not express all the master integrals in terms of the variable $z$ since, if one does this, square roots of $z$ appear. Explicit expressions for the integrals ${\rm INT}(\mathrm{top},\vec{n})$ are provided in an ancillary file, which may be downloaded from https://www.ttp.kit.edu/_media/progdata/2019/ttp19-009.tar.gz.

To illustrate the usefulness of the integrals presented in this paper, we construct the RRR contribution to the $I_{qq}$ matching coefficient at N3LO in QCD in the large-$N_{c}$ and the large-$N_{f}$ limits. Interestingly, upon inserting our results for the master integrals, we find that all $x$-dependent multiple polylogarithms as well as the rational functions of $x$ combine in such a way, that the final result is expressible in terms of rational functions of $z$ and harmonic polylogarithms of $z$ only. The required mappings from $G(\vec{w};x)$ to $H(\vec{w},z)$ were obtained by expressing all harmonic polylogarithms up to weight $6$ in terms of $G(\vec{w};x)$ with the program HyperInt Panzer:2014caa and subsequently inverting the (underdetermined) system of linear equations. The resulting $s$-independent contributions can be written as

[TABLE]

where $g_{s}$ is the strong coupling constant. The subscript $i$ is either $N_{c}^{3}$ to indicate the leading-color 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,

[TABLE]

in order to properly extract the collinear singularities. As a consequence, $\mathcal{A}_{i}(z,\epsilon)$ is needed up to first order in $\epsilon$. In turn, $\mathcal{A}_{i}(z,\epsilon)$ contains soft singularities, which are extracted by writing $(1-z)^{-1-3\epsilon}$ in terms of plus distributions. The results are rather lengthy, so we choose to only display their soft limits. They read

[TABLE]

6 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 N3LO in QCD due to the splitting of an incoming quark $q$ into a virtual quark of the same flavor and three collinear gluons, $q\to 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\to q^{*}+q\bar{q}g$. 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 N3LO 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. Melnikov:2018jxb 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. Melnikov:2018jxb .

Appendix A List of master integrals

In this appendix we list the $91$ master integrals.

[TABLE]

The definition of the topologies $A_{1}$ through $A_{19}$ may be found in table 1.

Appendix B Additional relations among the master integrals

In this appendix we prove two simple relations among some of the master integrals, eqs. 16 and 17. We define the two quantities

[TABLE]

Below we show that $F_{1}(z,\epsilon)=F_{2}(z,\epsilon)=0$.

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

[TABLE]

The solutions to these equations are

[TABLE]

In the limit $z\to 1$ these solutions for $F_{1}$ and $F_{2}$ behave as $(1-z)^{0}$ and $(1-z)^{-\epsilon}$. However, we have argued in the main body of the paper that all master integrals in the soft $z\to 1$ limit should be proportional to $(1-z)^{n-3\epsilon}$ for some integer $n$. The only way to make this scaling compatible with eq. 77 is to choose $c_{1}(\epsilon)=c_{2}(\epsilon)=0$ which implies that $F_{1,2}$ vanish identically. This proves the identities among master integrals shown in eqs. 16 and 17.

Acknowledgements.

We thank A. von Manteuffel for a useful advice concerning the calculation of boundary integrals. L.T. would like to acknowledge the Mainz Institute for Theoretical Physics (MITP), which is part of the DFG Cluster of Excellence PRISMA+ (Project ID 39083149), for its partial support during the completion of this work. The research of K.M. and R.R. was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762 - TRR 257. The research of L.T. was supported by the ERC starting grant 637019 “MathAm”. The research of Ch.W. was supported in part by the BMBF project No. 05H18WOCA1.