One-loop corrections to the double-real emission contribution to the zero-jettiness soft function at N3LO in QCD

We present an analytic calculation of the one-loop correction to the double-real emission contribution to the zero-jettiness soft function at N3LO in QCD, accounting for both gluon-gluon and quark-antiquark soft final-state partons. We explain all the relevant steps of the computation including the reduction of phase-space integrals to master integrals in the presence of Heaviside functions, and the methods we employed to compute them.


Introduction
Description of hard scattering processes at the LHC relies on QCD perturbation theory.To arrive at finite predictions for infra-red safe observables, one needs to combine all contributions to a particular cross section, that are proportional to the strong coupling constant raised to a certain power, but may differ by the number of final-state partons.This step is nontrivial as contributions with different partonic multiplicities suffer from infrared and collinear divergencies that only cancel in the sum [1][2][3].One of the several ways to simplify this step is known as the slicing method where one chooses a kinematic parameter h that distinguishes contributions with Born kinematics from those with additional final-state partons.Assuming that h = 0 corresponds to Born kinematics, one then computes the cross section for small values of h, where the kinematics is very much Born-like, and only soft or collinear final-state partons can be present in addition to the hard ones that already appear in the Born process.
Since description of soft and collinear emissions is independent of the hard process, it is not surprising that cross sections at small values of the slicing parameter are computed by integrating universal functions that describe soft-and collinear limits of squared matrix elements over phase spaces available to unresolved partons.For certain slicing variables, one can prove that these integrals factorize into so-called soft, beam and jet functions and, once these functions are known, singular dependencies of cross sections on h can be reconstructed.Furthermore, since soft, beam and jet functions separately satisfy renormalization group equations, they can be used to systematically resum logarithms of certain kinematic variables to all orders in the strong coupling constant; for the case of N -jettiness, see e.g.Refs.[4][5][6][7][8].
In this paper, we analytically compute the one-loop correction to the double-real emission contribution to the soft function for the zero-jettiness variable, accounting for both gg and q q soft final-state partons.The result for gg contribution was reported earlier in Ref. [51].In general, however, given the complexity of any three-loop computation, an independent calculation of the one-loop correction to double-real contribution to zero-jettiness soft function is clearly warranted.
The rest of the paper is organized as follows.In Section 2 we define the soft function and introduce the notations that we use throughout the paper.In Section 3 we explain how the expression for soft currents is constructed and the integration-by-parts reduction to master integrals is set up.In Sections 4 and 5 we discuss the various methods we employed to compute the master integrals.In particular, as we explain in Section 5, we found it useful to set up and solve systems of differential equations to calculate the most challenging master integrals.We discuss the computation of the boundary conditions for integrals that appear in the differential equations in Section 6.We present the result for the one-loop corrections to double-real emission contributions to the zero-jettiness soft function in Section 7. We conclude in Section 8. One-loop integrals that we used for the computations reported in this paper are collected in Appendix A. The list of all master integrals is given in Appendix B. Results for all master integrals as well as the expression for the one-loop correction to the double-real emission contribution to the zero-jettiness soft function in terms of master integrals can be found in an ancillary file provided with this paper.

The zero-jettiness soft function
We begin by defining the zero-jettiness soft function.The zero-jettiness variable in QCD is a limiting case of the so-called N -jettiness variable [23,24] which describes hadron collider processes without hard QCD partons (jets) in the final state, and lepton collider processes in which the number of hard final-state jets is exactly two.It is defined as follows where the list ψ i reads and P is an arbitrary normalization factor.In Eq. (2.1) n is the total number of unresolved final-state partons that need to be considered in a particular order of QCD perturbation theory, and p 1,2 are the momenta of two incoming partons that annihilate and produce a color-less final state X, or the momenta of two hard partons that are produced in a collision of an electron and a positron.
The zero-jettiness soft function is constructed from integrals of soft limits of the (various) scattering amplitudes squared, at different orders in QCD perturbation theory over phase spaces of soft partons subject to a jettiness constraint T 0 = τ.To be specific, we write the perturbative expansion of the bare soft function as where S (i) is the contribution to the soft function in i-th order in the perturbative expansion in QCD.Computation of S (i) , i = 1, 2, 3, follows same rules as the calculation of the perturbative expansion of any cross section except that i) all the matrix elements, both real and virtual, are computed using eikonal Feynman rules for fixed number of hard partons and ii) phase-space integration measure contains the zero-jettiness constraint δ(τ − T 0 ) but no energy-momentum conserving δ-function.
We also note that it is convenient to deal with the case of two hard partons in the final state, which covers the two-jet production in electron-positron collisions.In principle, computing the soft function for hadron collisions from the result for e + e − could have required a non-trivial analytic continuation.Fortunately, the zero-jettiness case is simple enough so that this does not happen [51,55,56] and the two results are actually the same.
We will now discuss the phase-space integration and fully define the contribution to the zero-jettiness soft function that we calculate in this paper.We are interested in computing the one-loop correction to the double-real contribution to the soft function.The double-real contribution appears for the first time in the calculation of S 2 ; hence, the one-loop correction to it affects S 3 .Thus, we consider the phase space for two soft particles and use n = 2 in Eq. (2.1).Since the zero-jettiness variable requires us to determine minima of various scalar products, it is convenient to split the phase space into sectors where these minima are uniquely determined, and compute contributions of these sectors separately.To this end, we introduce the Sudakov decomposition of the momenta of two final-state soft partons.We work in the center of mass frame of two hard partons p 1 and p 2 and write We choose n • n = 2, so that s = 2p 1 • p 2 = 4E 2 .We then write and find If we choose P = 2E in the definition of the zero-jettiness variable, we obtain There are two minimum conditions we need to resolve; hence, to identify the zero-jettiness uniquely, we use partition of unity Combining this with the delta function that ensures a particular value of zero-jettiness, we write For a particular final state f , the one-loop correction to the double-real contribution to the soft function is obtained by integrating the product of the tree-level (J (f ) 0 ) and the one-loop (J (f ) 1 ) eikonal currents over the corresponding phase space.Hence, we define where the sum over polarizations of the final-state partons is assumed.
The zero-jettiness constraint can be split into the same-hemisphere and the differenthemispheres contributions, given in Eqs.(2.9) and (2.10), respectively.Although it appears that four different terms need to be integrated for each of the final states f , we can simplify the calculation by making use of the symmetries between final-state partons, as well as the symmetries between two hemispheres, where appropriate.We will also use the fact that soft eikonal currents are uniform in the soft momenta; this allows us to factor out the dependence on the jettiness parameter τ .Hence, we write the one-loop correction to the double realemission contribution to the zero-jettiness soft function as where and ν f are factors that are particular to the final state f as we describe below.The two phase spaces that appear in the above equation read ) where i (2π) d−1 ).Note that, as the consequence of using just the two phase spaces, we have to compute the interferences of the one-loop and tree eikonal currents in Eq. (2.11) for different assignments of momenta k 1,2 , n and n.
We include gluons and n f quarks in our calculation.However, since we use the Feynman gauge for both virtual and real gluons, we must account for final-state ghosts to remove contributions of unphysical gluon polarizations from the final result.Hence, we write where c is the ghost field.To construct contributions of individual final states to S RRV , we use ν g = 1/2, ν c = −1 and ν q = 1.We also note that further simplifications in Eqs.(2.13) are possible if the matrix elements squared M (f ) are symmetric under the permutations of its arguments.We find that this is indeed the case for the gluon and ghost final states.However, it is not the case for the q q final state, because diagrams of the type shown in Fig. 1 are antisymmetric under the interchange of the momenta of the soft quark and the soft anti-quark.These diagrams lead to the QCD analogy of the charge asymmetry in QED, as was recently pointed out in Ref. [56].
We now summarize the steps required to compute the one-loop QCD correction to the double-real emission contribution to the zero-jettiness soft function.We elaborate on each of these steps in the remaining sections of the paper.We begin by constructing Born and one-loop two-parton currents using eikonal Feynman rules,2 compute their interferences, sum over polarizations of all final-state partons, map them onto templates in Fig. 2 and write the resulting contributions as integrals over phase spaces defined in Eq. (2.14).To compute these integrals we rewrite them as loop-like integrals using reverse unitarity [58] and use integrationby-parts technology [59,60] to express them through a small set of master integrals.However, the direct application of integration-by-parts methods to integrals with Heaviside functions is subtle [52].As follows from Ref. [52], integration-by-parts identities for integrals with θfunctions contain integrals with some or all θ-functions replaced by δ-functions.Hence, in addition to phase spaces shown in Eq. (2.14), we will also require integrals over dΦ nn δθ , dΦ nn δθ and dΦ nn δδ where the notation should be self-explanatory.We express all required integrals through master integrals, and then compute them either by performing loop-and phase-space integrations directly or by setting up and solving differential equations that these integrals satisfy.In the following sections, we elaborate on each of these steps.

Input generation and the reduction to master integrals
Although it is straightforward to construct the input expression for the interference of tree and one-loop soft currents, it requires significant amount of bookkeeping, especially at N3LO QCD.Because of that, we developed a tool chain based on the generic representation of Feynman diagrams that we have to compute.We start with the template shown in Fig. 2 where we always draw a one-loop amplitude to the left of the cut and the Born amplitude to the right of the cut.The amplitude to the right of the cut is supposed to be complexconjugated but it is simple to perform this operation for Born amplitudes.
To better control the input, we generate the one-loop amplitude and the tree amplitude separately.We work with Feynman rules adapted to describe the soft limit of QCD.To this end, in addition to conventional QCD Feynman rules, we include eikonal rules that parameterize the propagation of a hard parton and its interactions with soft gluons. 3o generate Feynman diagrams from QCD Feynman rules supplemented with additional eikonal lines and their interactions, we developed a custom model file for the package DIANA [61] which calls QGRAF [62] and creates a suitable input for Feynman diagrams that, together, provide soft currents, both tree and one-loop.Furthermore, as a check, we computed one-loop and tree two-parton currents starting from diagrams constructed from conventional QCD Feynman rules and then treating the emitted gluons and quarks, both real and virtual, in the soft approximation.The obtained expressions for the one-loop and tree amplitudes are simple enough to construct their interferences and map them onto "loop" diagrams shown in Fig. 2.This latter step is useful to simplify further calculations that we now describe.
The generated expression for the interference of the one-loop and tree eikonal currents is clearly the same for nn and nn kinematic configurations.However, since the expression for the jettiness variable is different in the two cases, delta-functions that ensure that the phasespace integration is performed at fixed τ are also different.In the context of reverse unitarity, these delta-functions are mapped onto propagator-like objects which then induce different relations between propagators in each of the configurations.Since these linear dependencies have to be resolved to define unique "loop-integral" families, these families are necessarily different for nn and nn configurations.
To deal with a large number of integrals, where loop and phase-space integrations are intertwined in a complex way, it is customary to apply integration-by-parts (IBP) reduction in the momentum space.In the context of the computation of the zero-jettiness soft function, the application of the IBP method requires non-standard modifications, as was first pointed out in Ref. [52].This happens because the IBP reduction of integrals with theta-functions forces us to consider inhomogeneous linear relations between integrals.The inhomogeneous nature of the system of equations that needs to be solved makes it impossible to use any of the well-tested public codes that are used in large-scale conventional multi-loop computations for the integral reduction.Because of that, we discuss several reduction-related points below that are particular to the problem at hand and which were essential for the successful application of the IBP technology to the computation of one-loop virtual corrections to the double-real emission contributions to the zero-jettines soft function.
Our starting point is the interference of the one-loop and tree eikonal currents expressed through a variety of integrals with two theta-functions.We can write integration-by-parts equations for these integrals, but these equations will contain inhomogeneous terms which are generated when a derivative with respect to the momentum of one of the soft partons acts on a θ-function.Since our goal is to use this system of equations to express a large number of the original integrals in terms of a small(er) number of simpler integrals, we need to introduce criteria to order integrals by their "simplicity".The most important criterion that we use for this purpose is the total number of theta-functions in an integral.We will refer to this parameter as level.We then notice that, since the IBP equations are obtained by computing exactly one derivative of the integrand, any IBP equation that appears in the calculation, can be cast into the following form where integrals on the left-hand side have level L and integrals on the right-hand side have level L − 1.This observation has to be contrasted with the result of the full reduction which implies that an L-level integral will be expressed through master integrals with levels from L to zero.In writing Eq. (3.3) we have assumed that there are n L master integrals M L i at level L. To order integrals further, we need to elaborate on the relative complexity of integrals with the same level value.
As we have discussed, in the course of deriving the IBP relations, an integral can move from an original level to a lower level.When this happens, an additional partial fractioning is required since delta functions, that appear in the integrands once θ-functions are differentiated, are treated as new propagators in the context of reverse unitarity.These new propagators often depend linearly on the propagators that appear in original integrals and such linear dependencies have to be resolved before one can proceed with the reduction.
The need to perform partial fractioning for each new IBP equation, except for those that are generated at the lowest level, is one of the major complications to efficiently setting up and solving the system of equations in this case, as compared to more conventional cases.Our strategy was to construct templates to generate IBP equations, perform partial fractioning, and systematically map integrals that only contain linearly-independent propagators on integral families defined at each level.In spite of the fact that the application of templates to large intermediate expressions with the goal to identify and classify integrals was crucial for the success of the reduction procedure, it still proved to be time-and resources-consuming process.
Similar to conventional cases, there are symmetry relations between integrals, and the actual reduction can be performed for smaller sets of truly independent integral families.Symmetries within a given family are based on the invariance of integrals under linear transformations of loop and/or final-state parton momenta but these transformations must be consistent with constraints imposed by an observable.In the context of loop computations, such linear transformations are applied to an integral described by a vector of powers of the propagators I[n 1 , n 2 , . . .] where negative values are used to describe scalar products in numerators.A momentum transformation then leads to a re-shuffling of positive indices and to momenta shifts in the numerators.The possible shift-identities are universal for all integrals characterized by the same set propagators in positive and negative powers.Following the standard terminology, these "sub-families" are referred to as "sectors".Mapping relations can be constructed between sectors; these relations have to be ordered to uniquely specify which sectors are kept and which ones have to mapped.
Mapping rules between integrals can be constructed by using the method first introduced in Ref. [63].This method relies on identifying relations between different loop integrals by bringing their Feynman parameter representations into a unique form so that their similarities become manifest.To apply it to integrals with delta and theta functions, we first map all the relevant integrals on to auxiliary loop integrals that contain all propagators of the original integrals and additional propagators constructed from arguments of each delta and theta function.To ensure that propagators that originate from delta and theta functions are not confused with conventional propagators, we add auxiliary masses to each type of new propagators (3.4) We then search for symmetry relations between different integrals following the procedure described in Ref. [63].
Following the above remarks, we introduce the global ordering of all sectors that appear in our problem using the L, T L , S T L ,L tuples, where L is a level, T L is a unique family label defined for the level L and S T L ,L is a sector identifier for a given family.Defining the mapping rules at each level, we construct the final list of tuples which contains unique sectors for which IBP relations must be constructed and solved.Ordering of tuples in complexity is fixed once and for all, and it tells us the "correct" order of IBP reductions in each of the unique sectors.
To generate IBP equations we use lists of seed integrals.It is difficult to accurately predict the maximal powers of the numerators and denominators of these seed integrals that are required for a successful reduction.We estimate these powers from integrals that appear in the expression for the soft function and we extend it to include all possible integrals that appear at the intermediate stages.The IBP reduction starts at the maximal level (L=2 in our case), where IBP identities are generated starting from the integrals in the seed list.All integrals that appear when the IBP identities are derived are mapped onto unique sectors.The same procedure is repeated for all simpler sectors and lower levels, until only L = 0 integrals are added when new IBPs are derived.Level-zero integrals do not contain theta functions; hence, they can be dealt with using the standard reduction methods implemented in Kira [64,65].To this end, we generate large look-up tables which contain all possible delta-only integrals that appear as inhomogeneous terms during the reduction.
To solve the constructed system of IBP relations, we start with the simplest sector with L = 1, where for the first time inhomogeneous L = 0 contributions appear.To deal with these inhomogeneous equations in a selected sector, we use reduction tables for integrals from all simpler sectors and construct a linear system readable by Kira where all integrals from simpler sectors are replaced with the corresponding master integrals.
We note that Kira's ability to deal with user-defined systems tremendously simplifies and speeds up the solution of the system of linear equations.Indeed, once the reduction tables that are needed to replace the integrals from simpler sectors through corresponding master integrals are provided, Kira can reduce the system to the minimal set of master integrals in an efficient way using modern approaches to the IBP reduction problem, such as finite field sampling [66].Since this strategy assumes that reduction results in a given sector are re-used in more complex sectors, we prepare reduction tables for all possible integrals in a particular sector which are limited by some powers of propagators in the denominators and some powers of irreducible scalar products in the numerators.It is important to note that when working with sectors at the same level, reductions of integrals in sectors with the same number of propagators can be done in parallel since they only depend on the reductions in simpler sectors which are assumed to be already available.

The master integrals
Once the reduction problem is solved, master integrals have to be computed.It follows from the above discussion that their calculation requires integration of one-loop integrals and various factors that depend on momenta of final state partons, over different phase spaces.We find that, in total, 61 master integrals need to be computed; they are listed in Appendix B. All one-loop integrals, needed to compute the master integrals for the soft function, contain two-, one-or no eikonal lines.All of them, except the five-point integrals, can be easily expressed in terms of hypergeometric functions; the results can be found in Appendix A. Furthermore, there is a large number of master integrals for which integration over phase spaces of two soft partons can also be performed with a relative ease.To show an example of how such computations are performed, consider master integral defined in Eq. (B.56).It reads where k 12 = k 1 + k 2 and 2,1 10111 can be found in Appendix A. We write it here one more time for the ease of reference where Ω (n) = 2π n/2 /Γ(n/2) denotes the surface of a unit sphere embedded in n dimensions.
To proceed with the calculation of the integral I 56 , we begin by integrating over the relative azimuthal angle between transverse momenta of the two gluons.Since α 1 β 2 > α 2 β 1 for nn integrals, we find 3) The required integral becomes We change variables where and use the δ-function to eliminate β 1 .We find To isolate singularities in this integral, we note that hypergeometric functions in Eq. (4.6) do not contain isolated singularities.In fact, ) has a logarithmic singularity at α 2 = 0. To isolate non-integrable singularities, we perform partial fractioning of the prefactor in Eq. (4.6) and find . (4.7) If the last two terms in the right hand side of Eq. (4.7) are used instead of the expression in the left hand side of this equation in Eq. (4.6), integrations over α 2 and r 1,2 can be performed upon expanding the integrand in ε since it does not contain non-integrable singularities.We use HyperInt [67] to perform the required integrations.
The second term in the right hand side of Eq. (4.7) has a singularity at α 2 = 1.However, α 2 = 1 is a regular point of the hypergeometric function and of an ε-dependent term (α 2 + (1 − α 2 )r 2 ) −ε in Eq. (4.6).Hence, focusing on the second term in Eq. (4.7) and the α 2 -integration we write The first term in the curly brackets on the right hand side of the above equation can be immediately integrated, and the term in the square brackets expanded in ε as it is not singular at α 2 = 1.The following integrations over r 1,2 are also not singular and can be easily performed using HyperInt.
We continue with the discussion of the first term in (4.7).Substituting it into the rest of Eq. (4.6), changing the integration variable α 2 → u, where and integrating over r 1 , we obtain the following integral to compute This integral is divergent at u → 1.To extract this divergence, we note that the following equation holds (4.12) To prove it, we first use the standard identity that relates hypergeometric functions with arguments u and 1 − u and, then employ one of the Gauss' relations to shift parameters of one of the hypergeometric functions.When we use Eq.(4.12) in Eq. (4.11) the first two terms on the r.h.s. of Eq. (4.12) can be integrated over u in terms of hypergeometric functions and, since the remaining integration over r 2 is non-singular, it can be performed expanding the integrand in powers of ε and integrating using HyperInt.The third term on the r.h.s. of Eq. (4.12) does not generate singularities when integrated over u and r 2 in Eq. (4.11).Thus, we integrate it with the help of HyperInt after expanding in ε.
Finally, putting everything together, we find the following result for the integral defined in Eq. (4.1) where we performed an expansion in ε up to O(ε) when weight-six contributions appeared for the first time.Although the absolute majority of the required integrals can be computed following similar steps, there are several integrals for which doing this appears to be difficult.As an example, we consider the integral in Eq. (B.49) which we also display here for convenience In Eq. (4.15) 1,2 11011 is a one-loop four-point function that can be found in Appendix A. It reads The difficulty with evaluating the integral in Eq. (4.15) is related to the appearance of the invariant mass of the two gluons k 2 12 in the argument of the hypergeometric function in Eq. (4.16).Since k 2  12 depends on the relative angle between transverse momenta of the two gluons, the integration over this azimuthal angle becomes very complicated.It is quite possible that, with some effort, one can find a way to expand this integral in powers of ε or derive a suitable Mellin-Barns representation for this integral, but we decided not to pursue this effort.Instead, we have opted for constructing a system of differential equations for the remaining integrals by adding auxiliary parameters to the problem.Below we discuss how this can be done.

Differential equations for master integrals
In this section we explain how to set up differential equations for master integrals that appear in the calculation of the one-loop correction to the double-real emission contribution to the zero-jettiness soft function.The key observation is that we can replace θ-functions in all master integrals that we need to compute with integrals of δ-functions that depend on the argument of the original theta function and an auxiliary parameter.We then use these parameters to derive and solve the differential equations for descendants of the original master integrals that contain these modified delta functions.
To illustrate this point further, we consider a θδ integral Since this is an nn-integral, it contains a theta function θ(α 1 − β 1 ).We then use the integral representation [49] 4 and write where J θδ (z) = dΦ nn δzδ α 1 f (k 1 , k 2 ). (5.4) The modified phase space in the above equation is defined as Since there is a parameter z in an integral with the δ-functions, we can use reverse unitarity to derive differential equations for J θδ (z) integrals.Exactly the same can be done for θθ-integrals where the only difference is that we will have to rewrite each θ-function as an integral of a δ-function so that resulting auxiliary integrals depend on two parameters z 1,2 .
Although δδ integrals are, in principle, simpler, not all of them are simple enough to allow seamless direct computation.For this reason, it is useful to set up a system of differential equations also for δδ integrals.Since the trick that we used for θδ and θθ integrals cannot be used for δδ integrals, we need to introduce an auxiliary parameter in a different way.Specifically, we consider a δδ-integral and write it as where (5.8) The parameter x fixes the invariant mass of the two soft gluons [51].We note that the integration boundaries in Eq. (5.7) are easily established using the Sudakov parametrization for the final-state gluons' momenta.Similarly to the case of z-dependent auxiliary integrals discussed earlier, it is straightforward to derive differential equations for J δδ (x) using reverse unitarity.
To summarize, following the above discussion we conclude that all master integrals that we require to calculate the one-loop corrections to the double-real emission contribution to the zero-jettiness soft function can be written in the following way where ⃗ J δδ,θδ,θθ are the parameter-dependent descendants of the original master integrals whereas ⃗ f δδ,θδ,θθ are x-, z-and z 1,2 -dependent master integrals.The functions R δδ,θδ,θθ are the reduction coefficients which arise when ⃗ J integrals are expressed in terms of ⃗ f integrals; they are rational functions of ε and x, z or z 1,2 .Since ⃗ J integrals do not contain θ-functions, the reduction of ⃗ J to ⃗ f , is performed using the standard Kira interface.As we already mentioned, ⃗ f integrals satisfy differential equations which can be constructed by differentiating w.r.t x, or z, or z 1,2 , and then reducing new ⃗ f integrals back to the original set.We note that the set of ⃗ f integrals should be large enough so that the system of differential equations closes.
Once the integral representations of ⃗ I integrals shown in Eq. (5.9) and the differential equations for ⃗ f integrals are derived, we solve these differential equations and determine ⃗ f integrals provided that an appropriate set of boundary conditions is known.Since these integrals are sufficiently complex, we can only compute them expanding in ε.We emphasize that such results for master integrals are valid in the bulk of the corresponding integration intervals but not at the boundaries where the integrals are singular.For this reason, it is in general not possible to use the ε-expanded results for the master integrals to integrate over x, z and z 1,2 because the corresponding integrals diverge at the integration boundaries.
To explain how, in spite of this problem, the ⃗ f integrals are computed and then integrated over auxiliary parameters, we note that a basis for these integrals can be chosen that brings the system of differential equations to the so-called ε-form [69].We will refer to this new basis of ⃗ f integrals as ⃗ g integrals.For integrals that depend on a single variables, Eq. (5.9a) and Eq.(5.9b), the differential equations in the ε-form read ∂⃗ g ∂y = εA(y)⃗ g, ⃗ f = T⃗ g. (5.10) For the case of two variables Eq. ( 5.9c) we obtain a system of coupled partial differential equations ∂⃗ g(y 1 , y 2 ) We note that we changed variables x = y 2 and z = y 2 for single-variable integrals and z i = y 2 i for integrals in (5.9c) to write the differential equations Eqs.(5.10,5.11).For all the required cases the transformation matrices T were found using CANONICA [70].
It turns out that the single-variable differential equations Eq. (5.10) in our case are such that the general solution expanded in ε can be immediately written in terms of harmonic polylogarithms [71].For the two-variable case Eq. (5.11), the general solution of the system is expressed through the generalized polylogarithms [72]. 5To construct particular solutions we need to determine boundary conditions.We do this by calculating selected integrals around chosen singular points of the differential equations.To this end, consider a singular point y 0 .Since this is a regular-singular point of the relevant systems of differential equations, the expansion of the integrals around it can be written as follows where ⃗ C y 0 is the vector of constants and the matrix U y→y 0 reads6 All possible values of parameters λ as well as recursive relations between matrix coefficients c λ,n can be obtained from the differential equations Eq. (5.10) which the matrix U y→y 0 satisfies.We note that the rational matrix A in that equation should also be expanded in series around y = y 0 .The resulting solutions contain a few constants that need to be determined separately.To this end, one computes the expansion of the original ⃗ f integrals at y = y 0 by calculating ⃗ f = T⃗ g.The constants ⃗ C y 0 can then be determined from the explicit computation of a few "branches" of ⃗ f integrals in the vicinity of y = y 0 .This step completes the determination of particular solutions for single-variable integrals ⃗ f .
Once integrals ⃗ f are known as an expansion in ε, their linear combinations with rational coefficients need to be integrated over x or z from zero to one, to recover the original master integrals, c.f. Eq. (5.9).As it turns out, these integrations, are singular at the boundaries.The singularities can appear both in matrices R and in the integrals ⃗ f themselves.To deal with these singularities, we need to construct ε-exact approximations to integrands in Eqs.(5.9a,5.9b)around singular points such that differences between actual and approximate integrands are integrable in the whole interval [0, 1] for ε = 0. Below we explain how this can be done using the ε-expanded solutions in the bulk of the interval provided that they are known to a sufficiently high ε-order.
The challenge is to construct an ε-exact solution of the differential equation in the vicinity of the singular point y = y 1 from the available ε-expanded solution valid at y ̸ = y 1 .This is an "inverse" problem to what we did when we constructed particular solutions from the approximate computation of integrals at y = y 0 .Hence, if we need to calculate the εresummed approximation to the solution of the differential equations at a point y = y 1 , we write an Ansatz similar to Eqs. (5.12, 5.13) with y 0 replaced by y 1 , expand the approximate result in ε and the exact one in y − y 1 , and determine the vector of constants ⃗ C y 1 .Once this is accomplished, we can compute the expansion of ε-exact integrals in y − y 1 to any order using recurrence relations derived from the differential equations, thereby determining the required subtraction terms.
We now explain how to do this in an efficient way.We start by writing the vector of unknown constants ⃗ C y 1 as where ⃗ g is known as an expansion in ε.In principle, we can compute the matrix U y→y 1 following the steps discussed in connection with the evaluation of the matrix U y→y 0 .However, it turns out to be more efficient to derive and solve the differential equation satisfied by U −1 y→y 1 instead of solving an equation for U y→y 1 and then inverting the matrix.We note that the equation for the inverse matrix is easily derived from an equation that the matrix U satisfies.We use these results to enable the calculation of the original master integrals ⃗ I. Accounting for the change of variables x, z → y 2 , we write (5.15) The subtraction term ⃗ S removes all singular terms of the expression 2yRT⃗ g in both y → 0 and y → 1 limits; it is obtained from expansions of the corresponding U -matrices around singular points using representation shown in Eq. (5.13).It reads By definition the integration of the first term of Eq. (5.15) converges on the interval y ∈ [0, 1] and we perform it using the package HyperInt [67].
The integration of the subtraction term is elementary since it only contain powers of y and (1 − y).The only difficulty here is a need for a deep expansion of vectors ⃗ C y→0 and ⃗ C y→1 in ε.However, using procedures explained above it is possible to construct such expansions provided that solutions to differential equations in an ε-expanded form are known to sufficiently high orders in ε.
We treat integrals that depend on two integration variables in a similar way, except that the boundary conditions now are series in two variables.This leads to certain changes in the described algorithm.We start from the system of differential equations in a canonical form, where and solve the first differential equation.We then find where ⃗ b is now a function of the second variable.Substituting Eq. ( 5.19) into the differential equation for y 2 , c.f. Eq. (5.17), we obtain By construction, this equation does not depend on y 1 and, therefore, can be solved following what has been discussed earlier.Hence, we bring this equation to a canonical form with the help of an additional transformation ⃗ b = T ′ ⃗ b ′ and write its solution as where ⃗ d is a vector of constants.Finally, the solution of the original two-variable system is written as follows The constants ⃗ d can then be determined by comparing the expansion of the integrals ⃗ f around chosen singular points with the expansion of matrices on the right hand side of the above equation.We note that, in variance with two single-variable cases, matrix U (y 1 , y 2 ) contains generalized Goncharov polylogarithms, whereas matrix V depends on harmonic polylogarithms only.
Finally, to construct the subtraction terms which are required to enable integration over y 1,2 , we proceed in the same way as in the single-variable case.The difference, however, is that we need to construct generalized expansions in two variables.Since we need ⃗ f (y 1 , y 2 ) integrals to compute θθ integrals and since these integrals only appear in an nn configuration which is less singular in comparison to nn, it appears that a fairly small number of such subtraction terms needs to be constructed.To manipulate and integrate generalized polylogarithms that depend on two variables, we use their implementation in HyperInt.

Boundary conditions
To compute integrals by solving the differential equations, we require boundary conditions.The point for which the boundary conditions are calculated is chosen in such a way that computation of integrals simplifies.To understand for which values of x, z and z 1,2 such simplifications happen, we need to consider the dependence of phase spaces and various scalar products on these parameters.The most complicated scalar product is 2k 1 k 2 because of its dependence on the azimuthal angle between momenta of the two gluons k 1,2 .Hence, an appropriate choice of the kinematic point should simplify this dependence.
We begin with the discussion of δδ integrals.In this case, the phase-space reads Since for δδ integrals we can trade the integration over cos φ 12 for the integration over x.We find where is the normalization factor.A convenient point to compute the boundary conditions is x → 1.Indeed, since the argument of the square root in the above equation needs to be positive, the integration region over β 1 is given by the following equation It follows, that if x is close to one, β 1 is approximately 1/2.Since β 1 = 1/2 is a non-singular point, computation of master integrals simplifies.The θδ and θθ integrals are similar.We consider the latter case, as it is more general.We assume that α 1 = β 1 /z 1 and β 2 = α 2 /z 2 .Using this parameterization, we find the following expression for 2k This expression simplifies if we consider the limit where both z 1 and z 2 become vanishingly small since the dependence on the azimuthal angle disappears.This simplification allows us to compute the boundary conditions for master integrals with a relative ease.We also note that a similar argument applies to θδ-integrals since in that case we obtain 2k 1 k 2 from the above equation by setting either z 1 or z 2 to 1.
We will now discuss two examples of the boundary integrals.One of the δδ-integrals that we need reads where the one-loop integral can be found in Appendix A. Using the expression for 2,1 11101 , and the fact that dΦ nn δδ δ(2k 1 k 2 − x) equals to dΦ in Eq. ( 6.1), we find where the integration over β is performed on the interval shown in Eq. (6.5).
We then consider the limit x → 1 in which case, as we explained earlier, the integration over β is restricted to β ≈ 1/2.It follows that where we used lim Finally, we note that all other x-dependent integrals can be studied in a similar fashion.
We continue with the discussion of the boundary conditions for θδ-and θθ-integrals.We will only consider one of the θθ integrals since the analysis of the θδ ones is rather similar.To this end, consider the following integral Using the expression for 2,1 11101 from Appendix A, we write (6.12) We observe that at small values of z 1,2 the dependence on cos φ 12 disappears from the integrand.Hence, the leading asymptotic behavior of J θθ (z 1 , z 2 ) in the limit z 1,2 → 0, follows from the simplified integral lim To integrate over β 1 in the limit z 1,2 → 0, we note that three integration region lead to particular dependencies on z 1,2 ; they are Contributions of each of these regions can be calculated independently of each other by expanding the integral in variables that are small in a particular region and extending the integration region to ensure that homogeneous scaling of the remaining integral in the small variable is achieved [73].Since we are interested in the leading asymptotic in each region, this is a relatively simple thing to do.Hence, we write lim The contribution of the "Taylor-expansion" region F T is computed by setting z 1 = z 2 = 0 everywhere in the integrand of the function F θθ .We then find In the region where β 1 ∼ z 2 the argument of the hypergeometric function is small, O(z 2 z 1 ), so that the hypergeometric function can be Taylor-expanded.Furthermore, to leading order in β 1 ∼ z 2 ≪ 1, we can replace To ensure that the resulting integral scales uniformly with z 2 , we extend the integration over β 1 to infinity and find Finally, we perform a similar analysis for β 1 ∼ 1 − z 1 .The only difference with respect to the previous case is that the hypergeometric function cannot be simplified since 1) in this region.Hence, changing variables β 1 = 1 − ξ, we write we present the result of the calculation of the one-loop correction to the emission of two partons where color charges of hard emitters are not specified.We find RRV,q q = C R C A (n f T F ) The label R = F, A in the above formula refers to the SU (3) representation of the hard emitters; furthermore, C F = (N 2 c − 1)/2N c , C A = N c and T F = 1/2 are the usual Casimir invariants of the gauge group SU (N c ).
We have checked the above result in several ways.First, we observe that there is no RRV,gg and RRV,q q.This is the consequence of the exponentiation of soft emissions in Abelian gauge theories as follows from considering R = F case.Next, the complete input for computing tree-and one-loop eikonal currents was assembled using two different methods.All master integrals as well as some integrals that appear in the soft function before the reduction to master integrals were computed numerically and compared with the analytic results.For the numerical evaluation, we used pySecDec [74] and MB [75] packages.Finally, our results for gg and cc final states agree with their recent evaluation reported in Ref. [51].We note that the result for S gg RRV reported in the ancillary file of Ref. [51] appears to contain irrational numbers of weight seven which are absent in our result.However, all these weight seven quantities disappear once the result reported in Ref. [51] is carefully simplified.

Conclusions
In this article, we described an analytic computation of the one-loop correction to the doublereal emission contribution to the zero-jettiness soft function.Our calculation is based on reverse unitarity and its modification [52] required to deal with phase-space integrals that contain Heaviside functions.Our result for gg final state agrees with the recent computation of the same quantity reported in Ref. [51].In addition, we also presented the one-loop correction to the contribution of the q q final state to the zero-jettiness soft function.
We explained in detail the way to construct the input for tree-and one-loop two-parton currents and the organization of the integration-by-parts reduction since publicly available tools cannot be used for this purpose when Heaviside functions are present.We discussed how master integrals are calculated; we present results for all of them in an ancillary file provided with this submission.
Finally, we briefly comment on what remains to be done to complete the analytic calculation of the zero-jettiness soft function.The most non-trivial contribution that still needs to be calculated is the part of the triple-real emission contribution that describes the kinematic configuration where two gluons are radiated into one hemisphere and the third gluoninto the opposite hemisphere.We are quite confident that methods employed in Ref. [53] to compute the same-hemisphere triple-real contribution to the soft function can be successfully adapted to this case as well, and we look forward to addressing this challenge.work of DB is supported in part by the Swiss National Science Foundation (SNSF) under contracts 200020 188464 and 200020 219367.

A One-loop integrals
where For the calculations described in this paper, we need a number of loop integrals.The results that we have relied upon are given below.

B List of master integrals
To present a list of master integrals, we define

Figure 1 .
Figure 1.Example diagram which are anti-symmetric under the interchange of quark and anti-quark momenta.

Figure 2 .
Figure 2. Different contributions to soft amplitude squared