The double-soft integral for an arbitrary angle between hard radiators

We consider the double-soft limit of a generic QCD process involving massless partons and integrate analytically the double-soft eikonal functions over the phase-space of soft partons (gluons or quarks) allowing for an arbitrary relative angle between the three-momenta of two hard massless radiators. This result provides one of the missing ingredients for a fully analytic formulation of the nested soft-collinear subtraction scheme described in Caola et al. (Eur Phys J C 77(4):248, 2017).


Introduction
A precise description of hard processes offers an exciting opportunity to discover or constrain physics beyond the Standard Model at the LHC using indirect methods.Such a description is based on the collinear factorization framework that emphasizes the importance of understanding partonic cross sections in higher orders of perturbative QCD.Currently, it is possible to compute most processes of phenomenological interest at a fully-differential level at leading and next-to-leading orders in perturbative QCD, and 2 → 1 and 2 → 2 processes at next-to-next-to-leading order (NNLO).
Nevertheless, in spite of these successes, it is fair to say that none of the suggested schemes are fully optimal.This is unfortunate as it can limit our ability to make precise predictions for higher multiplicity processes at the LHC in the future.Hence, further developments of subtraction methods are welcome.Motivated by these considerations, two of us in collaboration with R. Röntsch have recently proposed [1] a modification of the subtraction scheme described in Refs.[22][23][24].A key element in our proposal is the double-soft limit defined as follows.
We consider the double-real emission contribution to NNLO QCD corrections to the production of an arbitrary final state X in hadron collisions.Specifically, we are interested in the X + f final state, where f are either two gluons or a quark-antiquark pair.We assign four-momenta k 4,5 to the two additional partons, and consider the double soft configuration k 4,5 → 0, with no particular hierarchy between k 4 and k 5 .It is well known that soft emissions factorize.Indeed, in the soft approximation parton emission does not change the kinematics of the final state X and does not affect infra-red safe observables.Moreover, the matrix element squared of the process ij → X +f factorizes into a color-correlated emissionless matrix element squared for the process ij → X and a universal eikonal function that depends on momenta of hard radiators that are present in either the initial or the final state, and the momenta k 4 and k 5 of the soft partons.
An important ingredient of any NNLO subtraction scheme is the integral of the double-soft eikonal function over the phase-space of the two extra partons f , subject to kinematic constraints.In the framework of Ref. [1], the following constraints on energies of the soft partons are imposed In Ref. [1], double-soft integrals with constraints as in Eq. ( 1) were computed numerically for the case when hard emittors are back-to-back.Although this is adequate for the color-singlet production processes considered in [1], in more complicated cases the numerical approach becomes cumbersome, since it requires a non-trivial continuation of phase-space integrals beyond four space-time dimensions (see [22][23][24]30] for details).Moreover, beyond the back-to-back limit, the double-soft integrals become functions of an angle between the three-momenta of the hard radiators.Since these angles change from event to event, the required numerical computations become quite expensive.
In what follows, we show how to overcome these issues and present an analytic computation of the double-soft integrals required for the description of NNLO real emission contributions to an arbitrary process.The rest of this paper is organized as follows.In Section 2 we introduce our notation, present relevant formulas for the double-soft limit and define the integrals that need to be computed.In Section 3 we discuss how to use differential equations to find the phase-space integrals.In Sec. 4 we explain how to fix the boundary conditions needed to fully reconstruct the required integrals from the differential equations.In Section 5 we present our final results for the integrals of the double-soft eikonal functions.We conclude in Section 6.

The double-soft current and its integration
In this section, we consider the double-soft limit of a generic scattering process.It is well known that soft emissions factorize.We now recall basic features of this factorization, following closely Ref. [31].Interested readers should consult Ref. [31] for further details.
In QCD, soft emissions involve non-trivial color correlations.It is then convenient to introduce a color basis |c 1 , ..., c n , and write a generic scattering amplitude as where c i are the color indices.It is also useful to associate a color charge T i with the emission of soft gluons off a parton i.Its action is defined as where a is the gluon color index (a = 1, ..., N 2 c −1) and Here f abc and t c ab are the generators of the SU(N c ) Lie algebra in the adjoint and fundamental representations, respectively.The color charge operators satisfy with ) if i is a quark or an antiquark.Note also that each vector |M(p 1 , ..., p n ) is a color-singlet state, which implies Using this notation, the matrix element squared for the process ij → X + f (k 4 , k 5 ) in the double-soft limit k 4 , k 5 → 0 can be written as [31] |M(g 4 , g if f (k 4 , k 5 ) are two gluons, and if f (k 4 , k 5 ) are a quark and an antiquark.In Eqs.(6,7), "≈" means that we only consider the most singular contribution in the double-soft limit, and We also used g s,b to denote the bare QCD coupling constant and T R = 1/2.The sums in Eqs.(6,7) run over all pairs of hard color-charged radiators.The functions S ij and I ij read All other eikonal functions are defined as follows [31] where ǫ = (4 − d)/2 with d being the space-time dimensionality, and According to the computational scheme described in Ref. [1], the doublesoft matrix elements in Eqs.(6,7) should be integrated over the three-momenta of soft partons, subject to the constraints in Eq. (1).It follows that a doublesoft contribution to any differential cross section can be constructed if the following integrals are known where the factor -2 in S S (q q) ij is introduced for convenience.In Eq. ( 13) we introduced the short-hand notation As explained in Ref. [1], the energy ordering E 4 > E 5 accounts for the 1/2! symmetry factor relevant for gg emission.We find it convenient to use the same phase-space parametrization for q q emission as well.Since I ij (k 4 , k 5 ) = I ij (k 5 , k 4 ), the full result in the q q case is twice the result that is obtained by imposing the E 4 > E 5 ordering.
We satisfy constraints in Eqs.(1,13) by choosing the following parametrization of the energies of the two soft partons We note that integrals in Eq. ( 13) depend on the relative angles between the three-momenta of the hard partons but not on their energies.The first integral S S ij,kl in Eq. ( 13) is easy to compute since it is a product of two single-gluon eikonal functions.For further convenience, we introduce the following notation where n 2 i = 1.Moreover, we will use in what follows.
To compute S S ij,kl , we integrate over ξ and z and obtain The angular integrals over the emission angles of the gluons g 4 and g 5 completely factorize and can be easily performed.One obtains [32] where Finally, we obtain It is straightforward to obtain the expansion of the hypergeometric function through O(ǫ 4 ) using existing computer algebra packages [33].
The non-trivial part of the computation requires the calculation of the correlated emission terms S S (gg/q q) ij as a function of the scattering angle between the two hard partons.We describe such a calculation in the next section.

Double soft integrals
In this section, we describe the calculation of the correlated contributions to the eikonal integrals S S (gg/q q) ij defined in Eq. (13).For definiteness, we focus on the S S (gg) ij computation.The calculation of S S (q q) ij proceeds in a similar fashion.
We use the fact that S ij (k 4 , k 5 ) is a homogeneous function of the soft momenta k 4,5 and of the hard momenta p i,j .This implies that if we use the parametrization of the four-momenta as in Eq. ( 16), we can integrate over the variable ξ.After the ξ integration, we separate the integration over z and write The function G ij (z) is defined as and the four-momentum P is a time-like vector P = (1, 0).Although the final result for double soft integrals does not depend on the normalization of the four-momenta of hard radiators, when computing individual contributions we will use p i,j = 1/2(1, n i,j ).The two δ-functions in Eq. ( 23) provide constraints on the energies of the two gluons k 4,5 ; their arguments are chosen to make them "propagator-like" for reasons that will become clear later.
To calculate G ij (z) we need to integrate S ij (k 4 , k 5 ) over the phase-space of the two gluons with energy constraints shown in Eq. (23).We do this by mapping these phase-space integrals onto loop integrals following Ref.[34].After defining integral families, we apply the integration-by-part identities to reduce the number of independent integrals that need to be computed and to derive differential equations that these integrals satisfy.
We identify 19 master integrals to be calculated.To display them, we introduce seven propagator-like structures and define With this notation, we require the following integrals1 These master integrals are functions of the energy fraction z and of the relative angle θ between the two hard radiators i and j.To compute them, we use differential equations.
In principle, we can write differential equations for master integrals in both z and θ.As it is easy to see from their definition, the two integrals I 18 and I 19 are homogeneous in z; this implies that the z-differential equation does not give any non-trivial information in this case.Therefore, we computed these two integrals by solving the differential equation with respect to the scattering angle.The boundary conditions for these differential equations were determined from the values of I 18,19 computed in a situation when the three-momenta of the radiators are back-to-back, i.e. θ = π.We find the following results − 16xG 1,0,0 (x) + 8xG 1,1,0 (x) + 48G 0,0,0,0 (x) − 32G 0,0,1,0 (x) − 32G 0,1,0,0 (x) + 16G 0,1,1,0 (x) − 32G 1,0,0,0 (x) + 24G 1,0,1,0 (x) In writing the expressions for I 18,19 , we used the following expression for the normalization factor Also, x is the sine squared of half the relative angle between the threemomenta of hard radiators x = sin 2 δ, δ = θ/2, and G a 1 ,a 2 ,...,am (x) are the standard Goncharov polylogarithms.
The situation with the remaining seventeen integrals is rather different.Indeed, many of them couple to each other and the majority of them are not homogeneous functions of z.Although it is possible to use differential equations w.r.t. the relative energy and angle to determine the integrals also in this case, we found it more convenient to consider the differential equation in the energy fraction z and to determine the full dependence on the angle between the hard radiators by computing boundary conditions as functions of θ.We did not use a canonical form [35] for the z differential equation.In fact, it is relatively straightforward to achieve a canonical form for the first eleven integrals but after that it becomes much more difficult to do so.However, we managed to re-write the system of differential equations in such a way that integrating it by expanding master integrals order-by-order in ǫ becomes possible.In principle, this is absolutely sufficient for solving the system of differential equations.A possible drawback of this approach is that intermediate results tend to be quite cumbersome.This is, however, easy to deal with once all the expressions for the integrals are substituted to obtain the physical result.
The differential equations in z are of the following form where Â(0, z, δ) is a triangular matrix with vanishing diagonal elements.To integrate these differential equations, it is important to expose the dependence of the matrix Â on inverse powers of the monomials of z.This dependence is characterized by elements of the list shown below z, (1 + z), (sin The integration of the system of differential equations is greatly simplified if its coefficients are rational functions of the integration variables.To achieve this, we rationalize the square root in Eq. ( 30) using the following change of variables It leads to In addition to making the square root rational, the variable transformation Eq. ( 31) also maps all other z-dependent monomials in Eq. ( 30) onto rational functions of t.We obtain where As the result of the z → t mapping, we obtain a system of linear differential equations for the seventeen integrals with rational coefficients Since the matrix B is a rational function of t, integration over t can be performed in terms of Goncharov polylogarithms in a straightforward manner.This gives the result up to an integration constant that must be determined by matching to appropriate boundary conditions.We discuss the computation of the latter in the next section.

Boundary conditions
As explained in the previous section, we only integrate the t differential equation, without considering a differential equation in the scattering angle θ.We then need to compute the master integrals at a given value of t (or z) as a function of θ.It is natural to consider the boundary condition at z = 0, which corresponds to the situation where one of the two soft particles is much softer than the other.Not only is this the simplest kinematic point where such a computation can be performed, but it is also very useful for the subsequent integration over z in Eq. ( 22) since that integration is, in fact, singular at z = 0.The computation of boundary conditions is relatively straightforward for the majority of the master integrals but there are a few of them that require some effort.We will illustrate the relevant techniques by considering two representative examples.
The simplest master integral is the phase-space itself.It can be computed in a straightforward way by first integrating over energies and then over emission angles.The result reads where the normalization factor N ǫ is defined in Eq. (28).A significantly more complex integral is I 13 , which reads (cf.Eqs.(26,24)) Upon integrating over gluon energies, we obtain where we introduced the normalized solid angle integration measure as By inspecting Eq. ( 38), it is easy to see that, at small z, the master integral I 13 scales as z −1 .Therefore, we need to compute it to first subleading power to determine the integration constant.To accomplish this, we first combine denominators using Feynman parameters where ρ 4,5 η = 1 − n 4,5 • η and η = n i x + (1 − x) n j .We then use this representation in Eq. ( 38) and integrate over directions of the gluon g 4 .We obtain We still need to integrate the right hand side of Eq. (41) over x and direction of the vector n 5 to obtain I 13 , i.e.
It is quite obvious that such computations simplify dramatically if the expansion in small z is possible at early stages of the computation.Unfortunately, the hypergeometric function in Eq. ( 42) can not be expanded in powers of z because the maximal value of its argument in the z → 0 limit is one and the hypergeometric function in Eq. ( 42) is non-analytic there.To transform the integrand in Eq. ( 42) to a suitable form, we use the standard transformation for hypergeometric functions that connects F 21 (..., y) with F 21 (..., 1 − y).
We also note that since η is invariant under the replacement x ↔ (1 − x), one can replace 2x with 1 in the integrand in Eq. ( 42) without affecting the value of the integral I 13 .Splitting the integral into two contributions as the consequence of the hypergeometric transformation, we write with where κ = zρ 5η .In principle, the hypergeometric functions in Eq. ( 44) can be directly expanded in powers of z, since the goal of the transformations described above has already been achieved.However, the remaining integrations over x and the directions of the vector n 5 would have been quite difficult in this case.Fortunately, there exists another transformation of the hypergeometric function that reduces the complexity of the remaining integrations dramatically.It reads As we will see, this transformation completely removes square roots from the computation of the boundary conditions.
We begin by applying this relation to the computation of T (a) 13 .We note that T (a) 13 ∼ O(1) in the z → 0 limit, so that it contributes directly to the subleading term in the z-expansion of I 13 .Therefore, we are allowed to set z = 0 in the computation of T To compute this integral, we write the function F 21 as hypergeometric series and integrate over x using 1 − η 2 = 4x(1 − x) sin 2 δ.We find a very simple result T The computation of T (b) 13 is somewhat more complex, primarily because the z → 0 limit can not be taken directly.Using the transformation Eq. ( 45), we obtain To proceed further we note that if we write the hypergeometric function in Eq. ( 48) as the standard hypergeometric series, we can take the z → 0 limit in all but the first two terms of the expansion.Hence, we consider the contribution of these two terms separately.We write We begin with the computation of T (b),Σ 13 . We use the series representation of the hypergeometric function in Eq. (48), set z to zero, integrate over x and n 5 and obtain Next, we compute the contribution T (b), 2 13 that arises if we take the second term in the series representation of the hypergeometric function in Eq. ( 48).This term can be written as , it is impossible to expand the integrand in Eq. ( 42) in Taylor series in z.Nevertheless, to obtain an approximation to W (b),2 at small values of z by expanding the integrand, we can follow ideas about asymptotic expansions of Feynman diagrams known as the "strategy of regions" [36].
The integral in Eq. (51) has, obviously, three regions: The first two (soft) regions give identical contributions; we consider one of them and multiply the result by two.We refer to the third region as the "hard region".We therefore write The contribution of the hard region is obtained upon expanding the integrand in Eq. ( 51) in Taylor series in powers of z.Since we are interested in the O(z 0 ) term only, we obtain To compute the soft contribution, we focus on the region x ∼ z.We expand the integrand assuming x ∼ z ≪ 1 and extend the upper x integration boundary to infinity [36].We obtain The calculation of T (b), 1 13 proceeds along the same lines.The only complication is that, since T (b), 1 13 ∼ O(z −1 ), we need to expand the soft contribution in Taylor series in x ∼ z ≪ 1 to first subleading term.Apart from additional algebraic complexity, this does not lead to any conceptual complications.We provide the results of the calculation for completeness.We write where and Finally, we assemble the full result for the boundary condition of the integral I 13 using Eqs.(43,47,49-57).
We note that the computation of the other boundary integrals is performed following similar steps; in fact, the manipulations of the hypergeometric functions and the "expansion by regions" are similar for all complicated master integrals that one has to compute.We believe that the discussion of the boundary condition of the integral I 13 provides enough insight into how to deal with them and we do not discuss other integrals for that reason.We note that boundary conditions for all the seventeen integrals I 1,..,17 are given in the Appendix.

Final results
With the boundary conditions known and the system of linear equations rationalized, it is straightforward to integrate the equations in terms of Goncharov polylogarithms and to match the result of the integration to the boundary conditions.This allows us to write the function G ij (z) in Eq. ( 22) as a linear combination of master integrals where R (k) ij (z) are the reduction coefficients to master integrals.The integration over z in Eq. ( 22) is straightforward -we change variables z → t using Eq. ( 31) and integrate from t = (1 − sin δ)/ cos δ, that corresponds to z = 1, to t = cos δ that corresponds to z = 0.The only subtlety is that the integration over z diverges at z = 0. To overcome this problem, we write where Gij (z) ∼ z −1−2ǫ describes the non-integrable behavior of the function G ij (z) at small z.This function can be extracted from the computed boundary conditions for the master integrals and the small-z expansion of the reduction coefficients R (k) ij (z).Finally, since Gij (z) ∼ z −1−2ǫ , the last term in Eq. ( 59) can be trivially integrated over z and, since the first term is not singular at z = 0, the integrand can be expanded in ǫ and, after changing variables from z to t, the integration over t can be performed in a relatively straightforward way.
We note that after performing this final z (or, rather, t) integration, we obtain the result given by a linear combination of Goncharov polylogarithms up to weight four with indices drawn from the following set a ± , b ± , e ±iδ , cos δ, where a ± , b ± are given in Eq. (34).The arguments of these Goncharov polylogarithms are either (1 −sin δ)/ cos δ or cos δ.The result of the z integration (or rather t integration) appears to be very large and complex.However, it can be simplified using the (by now standard ) symbol techniques [37][38][39].
Computing the symbol of the result, simplifying it and integrating the result back, we arrive at the following expressions for the double-eikonal integrals Here, s = sin δ, c = cos δ and δ = θ/2 is half the angle between the threemomenta of the hard radiators i and j.The Clausen functions are defined as and G a 1 ,a 2 ,...,am (x) are the standard Goncharov polylogarithms.We have checked these analytic results by computing the functions S S (gg,q q) ij numerically for a few values of δ using the Mellin-Barnes representation for the original eikonal integrals and the MB.m routine for numerical integration [40].We found agreement within the expected numerical precision of the latter. 2he results shown in Eqs.(61,62) describe the integrals S S (gg,q q) ij as a function of the relative angle between the hard emittors.A useful special case corresponds to back-to-back radiators; this kinematic situation is relevant for the description of color singlet production and decay.In the back-to-back limit  (65) We have checked the back-to-back results Eqs.(64,65) against the numerical values used in Ref. [1], and found perfect agreement.

Conclusions
We computed integrals of the double-soft eikonal functions over phase-spaces of two soft gluons or a soft q q pair in the case when the three-momenta of the hard massless radiators are at an arbitrary angle to each other.Within the framework of a nested soft-collinear subtraction scheme [1], our results will allow for an analytic treatment of the double-soft contribution to NNLO QCD corrections to a generic process with arbitrary number of hard massless color-charged particles.Our results for the integrated double-soft functions are compact; they are expressed in terms of ordinary and harmonic polylogarithms which ensures that they can be evaluated numerically fast and efficiently.We look forward to the applications of these results in NNLO QCD computations.