Analytic double-soft integrated subtraction terms for two massive emitters in a back-to-back kinematics

We consider the double-soft limit of QCD amplitudes with two massive quarks in a back-to-back kinematics accompanied by two soft partons. We integrate analytically the respective double-soft eikonal functions over the phase space of the two soft partons. Within the context of the nested soft-collinear subtraction scheme, our results may serve as one of the integrated subtraction terms needed for the analytic and fully-differential description of next-to-next-to-leading order (NNLO) QCD corrections to colour-singlet decay into massive partons or to heavy-quark pair production.

Despite the large number of available subtraction and slicing schemes, it is fair to say that an optimal subtraction scheme, capable of dealing with complex processes, is yet to be designed. A set of criteria that should be considered when attempting the construction of a subtraction scheme may include physical transparency, scalability and locality as well as analyticity and efficiency. With these considerations in mind, the nested softcollinear subtraction scheme was introduced in ref. [49], building on the sector-improved residue subtraction scheme [29][30][31]. There it was shown that subtractions applied to gaugeinvariant scattering amplitudes, rather than to individual Feynman diagrams, can be done in a nested fashion, yielding a somewhat simpler description.
As the name suggests, subtraction schemes handle infrared singularities of real corrections by designing suitable subtraction terms for soft and collinear divergences. Properly constructed differences of real emission contributions and subtraction terms become integrable over the full phase-space in four dimensions. The subtraction terms, however, still need to be integrated in d = 4 − 2 dimensions where soft and collinear singularities manifest themselves as 1/ poles.
When the nested soft-collinear subtraction scheme is applied to massless partons, there are two genuinely double-unresolved limits that need to be addressed. These are the doublesoft limit, where two emitted partons become soft, and the triple-collinear limit, where momenta of three partons become collinear to each other. For both of these cases, integrals over the phase space of unresolved partons, subject to specific energy constraints dictated by the setup of the nested soft-collinear subtraction scheme, have been computed [50,51]. These results facilitate an analytic and fully-differential description of colour-singlet production [52], colour-singlet decay [53] and DIS-like processes [54]. We note that these "dipole-like" building blocks should enable a fully-differential NNLO QCD description of arbitrary processes.
The structure of IR singularities changes if massive quarks are involved in a partonic process. Indeed, since there are no collinear singularities related to massive external legs, only soft singularities need to be considered. They can be subtracted using appropriate soft eikonal functions that have to be integrated over the unresolved phase space. The goal of this article is to start exploring the subtraction terms that arise in NNLO QCD calculations for processes involving massive quarks in the context of the nested soft-collinear subtraction scheme. Specifically, we compute the integrated subtraction terms which are required to describe double-soft emissions off two radiators of the same mass in a back-to-back kinematics.
The remainder of this paper is organized as follows. In section 2, we describe the nested soft-collinear subtraction scheme and the single-and the double-soft functions for massive radiators. In section 3, we integrate the single-soft and the double-soft eikonal functions over the respective unresolved phase space. We discuss results in section 4 and conclude in section 5.

JHEP07(2020)011 2 Preliminary remarks
In this section, we specify a physical setup, describe the idea behind the nested soft-collinear subtraction scheme and establish notations by writing down the factorisation formulae for QCD amplitudes in the single-soft and the double-soft limits. We conclude the section by defining sets of single-and double-soft emission integrals that need to be computed.
Our goal is to describe infrared (IR) singularities that arise in NNLO QCD calculations for processes involving massive quarks. In this work, we focus on the IR singularities that appear in the double-real contribution to the decay of a colour-singlet particle X to massive quarks, i.e. a tree-level process In eq. (2.1), Q stands for a massive quark while f andf denote a pair of massless partons (gluons or a quark-antiquark pair). Following ref. [49], we write the contribution of the partonic process in eq. (2.1) to the decay rate as In eq. (2.2), N includes normalization and symmetry factors, dLips AB;12 denotes the Lorentz invariant phase-space measure of the massive quark system, including the energymomentum conserving δ-function, and F is the measurement function of an arbitrary infrared-safe observable. All of these quantities are then absorbed into the function F LM in eq. (2.2). Note that, following the original formulation of the nested soft-collinear subtraction scheme, we have introduced an energy ordering for the radiated partons k 1 and k 2 , i.e. we require E 1 > E 2 in eq. (2.2). We define the phase-space element of a massless parton as In eq. (2.3), we introduced an energy cut-off E max which is arbitrary but must be large enough so that it does not change the value of the integral in eq. (2.2), see ref. [49] for details. The need for such a cut-off parameter will become clear later when the double-soft limit of eq. (2.2) is discussed. For now, we only note that E max breaks Lorentz invariance but leaves rotational invariance intact.

The nested soft-collinear subtraction scheme
Infrared divergences of QCD amplitudes can be regulated by introducing appropriate subtraction terms for all relevant kinematic configurations. Within the nested soft-collinear subtraction scheme, such subtractions are constructed in an iterative manner, starting from the double-soft limit.

JHEP07(2020)011
The double-soft limit describes kinematic configurations where energies of both emissions in eq. (2.2), E 1 and E 2 , vanish at a comparable rate. To describe this limit, we introduce a double-soft projection operator S S. For a generic amplitude M involving momenta k 1 and k 2 , we consider the scaling E 1 ∼ E 2 ∼ λ and define (2.4) We then split the double-real contribution in eq. (2.2) into double-soft regulated and unresolved parts, i.e. (2.5) The first term on the right-hand side is not divergent in the double-soft limit. This term still contains singularities in the single-soft limit, where E 2 → 0, and in the collinear limit, where the two emitted partons become collinear to each other. Deriving a full subtraction would require us to remove these singularities as well. However, since integrated subtraction terms in these two limits can be obtained in a rather straightforward manner (see e.g. ref. [55]), in this paper we focus on the second term on the right-hand side of eq. (2.5) and its integration over the double-unresolved phase space. It reads We note that in the double-soft limit the momenta k 1 and k 2 completely decouple from the hard matrix element, from the energy-momentum conserving δ-function and from the measurement function F . This allows us to obtain integrals over the double-unresolved phase space in a universal manner. After the decoupling from the energy-momentum conservation, integrals over dE 1 and dE 2 in eq. (2.6) are only limited by the cut-off parameter E max introduced in eq. (2.3). As a consequence of the factorization in the double-soft limit, the reduced matrix element describes the Born-like process and the momenta p A,B are back-to-back in the rest frame of the decaying particle X. In this kinematic situation and the heavy-quark momenta are on the mass shell

JHEP07(2020)011
We note that the quark energy E is half the mass of the decaying particle, E = M X /2, the vector n describes the direction of flight of the heavy quark in the rest frame of the decaying colour singlet and (2.10) The threshold limit E = m implies β = 0.

Eikonal functions for single-and double-soft emissions
Soft factorization formulas for generic QCD tree-level amplitudes involving massless radiators and up to two soft partons were studied, for example, in ref. [56]. This result was extended to cover massive radiators in ref. [30] using the observation that eikonal currents are identical for massive and massless emitters and that emitters' masses become relevant only when eikonal currents are squared. We begin with the single-gluon emission. The limit of an amplitude that contains a gluon with a soft momentum k readŝ . . , p n } and the sum runs over all n hard emitters. The operatorŜ k extracts the leading asymptotic behaviour of the matrix element in the soft limit, . (2.12) The colour correlations in eq. (2.11) are encoded in the reduced matrix element 1 The double-soft function that describes emission of two gluons with momenta k 1 and k 2 reads where the additional colour correlated matrix element is defined as 15) and the notation {·, ·} stands for an anticommutator in colour space. The first term on the right-hand side of eq. (2.14) is the abelian contribution. It is simply a product of singleeikonal factors defined in eq. (2.12). Due to its factorized form, it is particularly easy to JHEP07(2020)011 integrate this term over the soft-gluons phase space. The second, non-abelian contribution is proportional to the colour factor C A . It is given by the function S ij (k 1 , k 2 ) which reads where we note that the term in square brackets explicitly depends on the squared masses of the emitters m 2 i and m 2 j . Both S 0 ij (k 1 , k 2 ) and S m ij (k 1 , k 2 ) implicitly depend on the masses. The first term in eq. (2.16), S 0 ij (k 1 , k 2 ), also appears in the factorization formula for massless emitters [56]; it reads .
( 2.17) The other two contributions in eq. (2.16) are only relevant for massive hard emitters. The function S m ij (k 1 , k 2 ) is given by [30] S m ij (k 1 , k 2 ) = − Note that we use an abbreviation k 12 = k 1 + k 2 in eqs. (2.17) and (2.18). When a soft quark-antiquark pair is emitted, the soft limit of the matrix element is described by where T F = 1/2. The soft function I ij (k 1 , k 2 ) is given by As we already mentioned, in the soft limit the dependence on the soft gluon momenta drops out from the matrix element as well as from the momentum conserving δ-function. For this reason, the eikonal factors in eqs. (2.14) and (2.19) can be integrated over the soft-gluons phase space, irrespective of matrix elements that describe the underlying hard process.

JHEP07(2020)011
In the following, we explain how to do that in the case of two equal mass emitters whose momenta p A and p B are back-to-back. To simplify notations, we introduce the single-emission phase-space integral where ij ∈ {AA, AB, BA, BB} and the phase-space measure and eikonal functions are defined in eqs. (2.3) and (2.12), respectively. For the double-emission phase-space integrals, we distinguish between emissions of gluons and quarks and define where, again, ij ∈ {AA, AB, BA, BB}. We note that in case of back-to-back kinematics, integrated subtraction terms BB and AA, as well as BA and AB, are equal to each other.
Therefore, in what remains, we will only consider cases ij = AA and ij = AB.
The integrals in eqs. (2.21) and (2.22) fully describe the integrated soft subtraction terms in the decay process of eq. (2.1) and are an important ingredient for more complex processes, such as heavy-quark pair production. The computation of phase-space integrals in eq. (2.22) is the main goal of this paper. We describe the details of the computation in the following section. We note that a similar calculation was performed in ref. [57], however, the unresolved phase space in that paper was subject to a slightly different constraint.

Phase-space integrals
In this section we present details of the calculation of the integrals defined in eqs. (2.21) and (2.22). We start with the single-soft emission to clarify notation and then proceed to the double-emission case. The results of the latter calculation are discussed in section 4.

Single-emission integrals
We start with a brief discussion of single-emission integrals. The first integral reads where we have parametrised the gluon four momentum as k = E(1, n k ). Further, we choose the reference frame in such a way that the z-axis points in the n direction. This yields (1 ± βn · n k ) = (1 ± β cos θ) and, after introducing η = (1 − cos θ)/2, we obtain

JHEP07(2020)011
where Ω (n) = 2π n/2 / Γ (n/2) denotes the surface of a unit sphere embedded in n dimensions. The integral in eq. (3.2) can be written as a hypergeometric function of the type , which further simplifies to [58] We find The second integral, for a self-correlated emission, reads Note that the hypergeometric function which appears in eqs. (3.4) and (3.5) features an expansion in powers of in terms of classical polylogarithms with arguments that involve square roots of β. In order to simplify the expansion, we again rewrite the hypergeometric function [58] and find Using HypExp [59], the hypergeometric function in eq. (3.6) can be expanded as The results shown in eqs. (3.4) and (3.5) were derived earlier in the literature [60,61].

Double-emission integrals
We now turn to the calculation of the double-soft subtraction terms. We need to compute the four functions GG AA , GG AB , QQ AA and QQ AB in eq. (2.22). To this end we employ reverse unitarity [62] that has been previously used for the computation of other integrated subtraction terms [50,51].
Computational setup. The integration measure for the two energy-ordered emissions in eq. (2.22) reads

JHEP07(2020)011
Note that all integrands, S ij (k 1 , k 2 ) and I ij (k 1 , k 2 ), are homogeneous under uniform rescaling of E 1 and E 2 . For this reason we parametrise the energies as and integrate over x to obtain where n i = (1, n i ), and the angular integration measure reads dΩ . It remains to carry out angular and z integrations in eqs. (3.10) and (3.11). However, the gluon emission case exhibits an additional singularity in the strongly-ordered limit, where the gluon with momentum k 2 is much softer than the gluon with momentum k 1 . Such behaviour results in a logarithmic divergence in the z integration at z = 0, which prevents us from a naive Taylor expansion of the integrand in . The problem can be ameliorated by using endpoint subtraction at z = 0. To accomplish this, we extract the divergent part using the following formula We note that it is beneficial to perform such a subtraction at the level of the full integrand since the resulting expression fully accounts for gauge properties of QCD amplitudes and, in variance to individual integrals, does not exhibit unphysical singularities. Note that an emission of a soft quark-antiquark pair does not exhibit the z → 0 singularity and, for this reason, does not require additional subtraction.
To perform angular integrals in eqs. (3.10) and (3.11) we proceed as follows. In the spirit of reverse unitarity [62], we rewrite δ-functions through cut propagators. To this end, we first rewrite the angular integration measures for both emissions as with ξ 1 = 1 and ξ 2 = z. By applying Cutkosky rules [63] backwards, we define cut loop integrals (3.14)

JHEP07(2020)011
We note that the variable z appears only in one of the cut propagators and plays the role of an internal mass. We use the definitions of eq. (3.14) in eqs. (3.10) and (3.11) and write ij (z, β, ) After mapping angular integrals onto ordinary loop integrals with cut propagators, we employ standard techniques of loop calculations to compute the integrals that appear in eqs. (3.15) and (3.16).
IBP reduction. We apply integration-by-parts (IBP) techniques [64] to the integrands of eqs. (3.15) and (3.16) to express them in terms of a few master integrals. The integrands consist of two-loop cut integrals where the propagators to be cut are given by The variables α i in eq. (3.17) refer to powers of propagators in integrals in a certain topology T a 1 ,a 2 ,a 3 . The prefactor in eq. (3.17) was chosen to render integrals dimensionless.
To express all integrals in eq. (3.14) through these topologies, we use the following list of linear relations between propagators where the last two equations follow from the cut constraints. We use Reduze2 [65] to express integrals shown in eq. (3.14) through master integrals. We write where R X X ij (z, β, ) are vectors of reduction coefficients and I(z, β, ) stands for a vector constructed out of thirteen master integrals grouped into five topologies. The first integral is the phase-space volume
Differential equations. Having obtained a set of master integrals we employ the method of differential equations [66][67][68] to compute them. To this end, we derive a closed system of first order partial differential equations for the master integrals I as functions of variables β and z. We then cast the differential equations into the -homogeneous form [69] by changing the basis of master integrals HereT can is the transformation that brings master integrals into their so-called canonical basis J . In general, finding a canonical basis or, equivalently, constructing a transformation T can is a complicated task. In our case, we accomplish this by using the algorithmic approach suitable for multi-scale problems proposed in ref. [70] and implemented in the CANONICA package [71] for Mathematica. This transformation can also found using the approach of ref. [72] implemented in a private Mathematica tool Libra. 2 In this case, a sequential application of the algorithm of ref. [72] is required.
In the canonical basis J , differential equations take the -homogeneous form with x ∈ {z, β}. The matricesM z andM β feature simple poles and can be written aŝ

JHEP07(2020)011
In eq. (3.26), the residue matricesm x i are composed of rational numbers and the poles x i are drawn from the two alphabets Thanks to the -homogeneous form of the differential equations in eq. (3.25), the expansion of the functions J (z, β) can be obtained by recursive integration of the righthand side. Since matricesM x contain only simple poles, the result can be expressed in terms of linear combinations of Goncharov Polylogarithms (GPLs) [73] that depend on z and β and constants of integration. Note that, since we are interested in a final integration over the variable z in eqs. (3.15) and (3.16), it is beneficial to write master integrals in such a way that z appears only as an argument of the GPLs. For this reason, at each order in , we first integrate the system of differential equations with respect to z. A constant of integration in this case is an unspecified function of β. To determine this function, we substitute the solution into the differential equations in β, and explicitly check that the resulting differential equations are z independent. After integration over β, all master integrals are expressed in terms of GPLs, G({ z 0 }; z) and G({ β 0 }; β), where the elements of z 0 are drawn from the alphabet A z , cf. eq. (3.27), and elements in β 0 belong to the z-independent part of the alphabet A β in eq. (3.28), i.e.Ã β = {0, −1, +1}.
This concludes the computation of master integrals up to constants of integration. These constants are determined by calculating suitable boundary conditions as we discuss in the next section.
Boundary conditions. We find it suitable to determine constants of integration by computing master integrals in the threshold limit β → 0. This limit is particularly convenient, since many of the integrals simplify. This happens because in that limit the dependencies of all scalar products on quark momenta disappear. For example (3.29) By inspecting master integrals in eq. (3.23), we observe that Moreover, we find that all entries, except for the first diagonal element of the canonical transformation matrixT −1 can are suppressed as O(β) and therefore vanish in the threshold limit. The transformation matrix in the threshold limit reads This means that to fix all integration constants we only need the phase-space master integral I 1 , which is straightforward to compute, cf. eq. (3.22).

JHEP07(2020)011
After fixing all the integration constants using boundary conditions, we transform master integrals J into the original basis I. We check the resulting expressions numerically for several values of β and z.
Integration over z. Having computed the required master integrals, we obtain the integrands in eqs. (3.15) and (3.16) and perform the integration over z. The masters integrals I of eqs. (3.22) and (3.23) allow us to express the functions E X X ij (z, β, ) in terms of rational functions of z, β and GPLs of z and β with z-independent letters. Such a representation enables the final z integration in eqs. (3.15) and (3.16) in a straightforward manner. We note that after obtaining the primitive, the z → 0 limit features spurious 1/z n poles and needs to be taken with care. We use PolyLogTools [74] to expand all GPLs around z = 0 up to the order required to cancel these 1/z n poles and facilitate z-integration over the interval 0 < z < 1. We report results for the functions GG AA , GG AB , QQ AA and QQ AB in the next section.

Results
In this section, we present some results for the integrated double-soft subtraction terms, cf. eq. (2.22). We write with Ω (n) defined after eq. where we used the abbreviations Even though expressions for the finite parts of the functions f gg,qq ij (β, ) are rather long, they simplify in certain limits. In what follows, we present the expansions in the threshold limit, β → 0, and the high-energy limit, β → 1.
We begin with the threshold limit, where the energies of the emitting quarks are close to their masses, i.e. E ≈ m, which implies β 1. We perform a Taylor expansion in small JHEP07(2020)011 terms numerically, as it was done, for example, in refs. [30,55]. Nevertheless, it is usually beneficial to have analytic results available. The results presented in this article provide all integrated double-soft subtraction terms required for a description of colour-singlet decays into massive fermions. For the case of heavy-quark pair production, it is also necessary to consider integrated subtraction terms with one massless and one massive parton which are not necessarily in a back-to-back kinematics. We leave this problem for future investigations.