On phase-space integrals with Heaviside functions

We discuss peculiarities that arise in the computation of real-emission contributions to observables that contain Heaviside functions. A prominent example of such a case is the zero-jettiness soft function in SCET, whose calculation at next-to-next-to-next-to-leading order in perturbative QCD is an interesting problem. Since the zero-jettiness soft function distinguishes between emissions into different hemispheres, its definition involves θ-functions of light-cone components of emitted soft partons. This prevents a direct use of multi-loop methods, based on reverse unitarity, for computing the zero-jettiness soft function in high orders of perturbation theory. We propose a way to bypass this problem and illustrate its effectiveness by computing various non-trivial contributions to the zero-jettiness soft function at NNLO and N3LO in perturbative QCD.

1 Introduction Precision studies of hadron collisions moved into the focus of the particle physics community after no clear evidence for physics beyond the Standard Model has been found during the first two LHC runs.A pre-requisite for such studies is a solid theoretical framework that allows one to describe hadron collisions using quark and gluon degrees of freedom and, in this way, connect experimental data to the Standard Model Lagrangian without the need for additional modeling.The current theoretical framework is based on the concept of collinear factorization [1] that, for processes with large momentum transfer, relates hadronic cross sections to convolutions of partonic cross sections, computable in perturbation theory, with universal non-perturbative parton distribution functions.Further refinements and practical advancements of such a framework are currently among the central topics in theoretical collider physics.
Improvements in an existent framework may be provided by the soft-collinear effective theory (SCET) [2][3][4][5][6] that seeks to establish a general pattern of factorization in collider processes including both perturbative and non-perturbative physics.This effective theory defines objects that are sensitive to particular momenta "modes" such as e.g.soft, collinear etc.These objects can be calculated independently of each other and then combined to provide predictions for physical quantities such as cross sections and kinematic distributions.Key to this effort are factorization theorems that both define these objects precisely and also provide information on how physical predictions should be assembled once these objects have been computed.
There are four types of objects that appear in SCET; they are known as hard, beam, jet and soft functions.Theoretical predictions for these functions are of interest since they can be used to re-sum logarithmically-enhanced terms that appear in perturbative expansion in QCD.In addition, they can also be used to set up a slicing method for deriving fixed-order predictions for fully-differential calculations in QCD [7][8][9][10][11][12][13][14][15][16][17][18].
The slicing method requires an observable which can be used to separate the realemission phase spaces into "singular" and "regular" parts.Although, in principle, any observable can be used to do that, if an observable is chosen in a way that does not violate factorization into collinear and soft modes, a cross section, differential with respect to such an observable, should satisfy a particular factorization theorem.Such a factorization theorem would then contain soft and beam functions whose computation enables both, the construction of subtraction terms for perturbative computations and the re-summation of large logarithms that arise in theoretical predictions once the slicing variable becomes small.
In this paper we will deal with the so called zero-jettiness variable [19,20] that can be used as a slicing variable for processes where a color-less final state is produced in hadron collisions.Calculations of the corresponding soft and beam functions have been performed during the past decade.The zero-jettiness soft function has been computed through NNLO QCD in Refs.[21,22] (see also Ref. [23,24]).The zero-jettiness beam function has been calculated through NNLO QCD in Refs.[25][26][27].Studies of zero-jettiness beam functions at N3LO QCD were initiated in Refs.[24,[28][29][30] and were recently completed in Refs.[31]. 1  In this paper, we focus on certain technical aspects that arise in the computation of the zero-jettiness soft function in higher orders of perturbative QCD.Due to algebraic complexity of computations in high orders of perturbation theory, specialized tools and methods are usually employed.Chief among them is the integration-by-parts (IBP) method for loop integrals introduced in Ref. [41] and adapted to real-emission integrals in Ref. [42].However, application of the IBP method to the computation of the zero-jettiness soft function is not straightforward, as the zero-jettiness observable contains Heaviside functions 1 Jettiness soft functions for more complicated final states were studied in Refs.[32][33][34][35].Beam functions for qT factorization were calculated in Refs.[36][37][38].Automated approaches to computation of soft, beam and jet functions are discussed in Refs.[39,40].
that depend on the light-cone components of four-momenta of the emitted partons. 2 This fact, as well as the need to compute large number of complicated integrals, makes the calculation of N3LO soft function non-trivial.Our goal in this paper is to discuss possible ways to overcome these technical difficulties paving the way for the computation of the triple-real and real-virtual contributions to the N3LO soft function.As a proof of concept, we compute a number of non-trivial contributions to the N3LO soft function, that describe emission of three gluons into the same hemisphere.
The remainder of this paper is organized as follows.In Section 2 we introduce the zero-jettiness soft function.In Section 3 we explain how integration by parts can be used to simplify computation of integrals with Heaviside functions.In Section 4, we apply this method to compute the soft function at next-to-next-to-leading order in QCD.We then discuss various aspects of the N3LO calculation in Section 5. We conclude in Section 6.Some useful formulas are collected in Appendices A and B.

Heaviside functions in the zero-jettiness soft function
In this section we will briefly discuss the zero-jettiness soft function to explain how Heaviside functions appear in phase-space integrals.Consider a process where two partons with (normalized) momenta n and n collide and produce a color-neutral final state together with m QCD partons with momenta k 1,2,.,m .The jettiness variable T is defined as follows To compute the minima in Eq. (2.1), we use the Sudakov parameterization of the parton's momenta and write where k ⊥,i n = k ⊥,i n = 0, n 2 = n2 = 0 and nn = 2.It follows that To enable the choice between α i and β i in the above equation, we can partition the phase space by writing for each of the m partons. 3It is clear that the first term will contribute β i and the second one α i to the jettiness variable T .Inserting the partition Eq. (2.4) into integrals over momenta of final-state gluons we find phase-space integrals with Heaviside functions.In general, for m-emitted gluons we can obtain 2 m different terms however, using the symmetry between m gluons and the fact that the result is invariant under the simultaneous replacement of all α i 's with β i 's and vice versa, the number of different terms is dramatically reduced.For example, for m ≤ 3 which covers NLO, NNLO and N3LO cases, only two independent contributions to the soft-function need to be considered.They are 1) all m gluons are emitted into the same hemisphere and 2) all but one gluons are emitted into the same hemisphere.
These different cases can be described as integrals over the following phase spaces where and Thanks to the definition of the jettiness variable, we associate functions f j in Eq. (2.5) with Heaviside functions but, as we will see later, we will also need to consider cases where one or several of these functions f are δ-functions.
To explicitly see how these phase spaces are used, we consider an example of the realemission contribution to the soft function.Then, for NNLO and N3LO computations we require the following integrals and where Eik 3g ({k i } and Eik 2g ({k i } are properly rescaled three-gluon and two-gluon eikonal functions.We will specify these functions later.For now, suffice it to say that they depend on the scalar products of gluon four-momenta k i k j and on scalar products of gluon fourmomenta with external vectors n and n.
A standard way to simplify computation of complicated phase-space integrals is to use reverse unitarity [42] to map such integrals onto cut loop integrals for which IBP identities can be derived in a straightforward manner.This is achieved by using the formula where P ( x) is a polynomial in variables x that can be e.g certain components of gluons' momenta.Once all δ-function constraints in phase-space integrals are removed using Eq.(2.10) and, provided, that there are no other non-polynomial constraints in the integrands, one can make use of the powerful integration-by-parts technology [41] to reduce computation of a large number of phase-space integrals to a few master integrals.Unfortunately, if integrands contain Heaviside functions, this approach fails since the last condition mentioned in the previous paragraph is not fulfilled.A possible solution to this problem was pointed out by one of us in Ref. [24], where it was suggested to rewrite all θ-functions that appear in relevant integrals as follows4 (2.11) While this representation yields an integrand whose dependence on auxiliary variables z i can be computed using reverse unitarity, it also introduces one additional parametric integral per θ-function, which can become quite cumbersome.For this reason, in this paper we would like to investigate how to derive and use IBP relations for phase-space integrals with θ-functions directly, i.e. without the need to introduce additional variables.
3 Applying integration-by-parts technology to integrals with θ-functions The goal of reverse unitarity [42] is to turn phase-space integrals into loop integrals.We explained the main idea of the method in the previous section; we will now make this discussion more specific considering the zero-jettiness soft function.
To remove all δ-function constraints from the integration measure, we start with a phase-space element of a gluon i with momentum k µ i = (E i , k i ) and write it as We then re-write the δ-function as in Eq. (2.10) In addition, to deal with soft functions shown in Eqs.(2.8,2.9)we need to re-write all δ-functions that define the re-scaled jettiness as cut propagators.For example, in case of the nnn kinematic configuration, we write where Hence, if we ignore θ-functions in the integrands of S nn , S nn , S nnn , S nnn , we immediately recognize that we need to compute a collection of "loop" integrals with conventional and unconventional cut "propagators".To do that, we can apply integration-by-parts identities [41] to reduce the number of independent integrals that need to be calculated.Furthermore, there are powerful public programs such as Fire [46,47], Kira [48,49], LiteRed [50,51], and Reduze [52,53] that can perform reductions to master integrals in a highly automated and efficient fashion.
However, if relevant integrands contain a collection of θ-functions, as is indeed the case for the soft function, this procedure can not be applied.An obvious problem is that θfunctions cannot be turned into "propagators" since the mapping similar to the one shown in Eq. (3.3) does not exist.However, we would like to understand what happens if we ignore this problem and attempt to derive IBP identities for integrands with θ-functions.
To study this question, we consider the following integral where f (k) is a polynomial in momentum k and g(k) is a function that allows a standard derivation of integration-by-parts identities. 5To derive integration-by-parts identities for the integral I[θ(f ), g], we write the standard equation that is valid for dimensionally-regularized integrals.Vector v µ in Eq. (3.5) is an arbitrary vector that we do not need to specify further.Calculating the derivative, we obtain It follows that 0 The first term on the right-hand side belongs to the same class of integrals as the original one I[θ(f ), g] because it involves the same θ-function; we call this term the homogeneous part of the IBP relation.Since, by assumption, integrals of g(k) can be studied using standard integration-by-parts technology, it follows that the homogeneous terms produce a closed set of linear equations when studied on their own.The second term in Eq. (3.7) involves δ(f ); we call this term the inhomogeneous part of the IBP relation.Since we can use the generalized unitarity trick to write δ(f ) → 1/f and since f is a polynomial in k, I[δ(f ), gv µ ∂ µ f ] defines a class of integrals that can be studied on their own independent of integrals with θ-functions.In fact, obtaining integration-by-parts identities for this class of integrals can be done with standard methods.The only subtlety that we have to deal with when working with inhomogeneous terms is that the function g(k)/f (k) may contain linearly-dependent "propagators" that will have to be re-mapped onto properly-defined integral families.Although this, by itself, is not a crucial issue, it does not allow us to derive IBP relations for integrals with arbitrary powers of propagators and forces us to produce IBP relations for each of the seed integrals individually.
We thus conclude that it is possible to establish useful integration-by-parts identities for integrals with multiple θ-functions by iteratively using Eq.(3.7).It follows from that equation that the derivative of an integrand produces inhomogeneous terms, where a θfunction is replaced by a δ-function.
Therefore, by using Eq.(3.7) repeatedly, we obtain a hierarchical sequence of IBP relations containing integrals with a decreasing number of θ-functions and an increasing number of δ-functions.The IBP relations can be used to express all relevant integrals through a set of master integrals.When choosing master integrals, we try to select those that contain fewer θ-functions since they are easier to compute.We will illustrate the construction of IBP relations and their usage in the next section where we will calculate the real-emission contribution to the zero-jettiness soft function at NNLO.

IBP identities and the NNLO QCD contribution to the zero-jettiness soft function
In this Section, we show how to use reverse unitarity and modified IBP relations to compute the maximally non-abelian contribution to zero-jettiness soft function at NNLO.We define it as where [54] ω with In the next section we explicitly construct a few examples of IBP equations with θfunctions that are relevant for this case.We discuss the computation of the two terms in Eq. (4.1) after that.

An example of an IBP relation
In this section, we explain how to employ modified integration-by-parts identities discussed in Section 3. As the first step, we map all integrals that appear in Eq. (4.1) onto integral families.These families are defined by sets of linearly-independent propagators6 and additionally contain two θ-functions from phase-space measures in Eq. (4.1).Cut propagators are constructed from δ-functions δ(k 2 i ) that enforce the on-shell conditions for the emitted gluons, and also from jettiness-dependent δ-functions that appear in the corresponding phase spaces After partial fractioning, we find that we need several independent integral families in this case.For example, one integral family that is required to describe the nn-configuration reads where the subscript c denotes cut propagators.
In order to construct an explicit example of the modified IBP relations discussed in Section 3, we consider the integral Starting with this "seed integral" and following the discussion around Eq. (3.5), we derive eight different equations by computing derivatives w.r.t.k 1 and k 2 and by using vectors v ∈ {k 1 , k 2 , n, n}.For example, differentiating w.r.t.k 1 and choosing v = k 1 , we find In writing this equation, we have used the fact that homogeneous terms, which arise when the derivative does not act on the θ-function, belong to the same family of integrals as the seed integral.
As we already mentioned, inhomogeneous terms arise when derivatives act on the θfunctions.In our example the last term in Eq. (4.8) is inhomogeneous.This term vanishes because k . This is a general feature; indeed, we find that integration-by-parts identities that involve a differential operator k µ i ∂/∂k i,µ do not produce inhomogeneous terms.To prove that assertion, consider The second term on the right hand side of Eq. (4.9) vanishes which implies that differential operators of the form k µ i ∂/∂k µ i do not produce inhomogeneous contributions.As a second example, we consider the derivative w.r.t.k 1 and use v = n.We find Similar to the previous case, we have expressed homogeneous terms through integrals of the family T ex .However, in this second example, the inhomogeneous term does not vanish and requires further treatment in the IBP setup.
The reason for that is the fact that the propagators in that term are linearly dependent due to the new δ-function in the integrand δ(k 1 n − k 1 n) that becomes a rational function if reverse unitarity is used.To see this explicitly, we make the replacement δ(k , multiply the result with the partial fractioning identity and obtain We note that terms that do not contain the complete set of cut propagators were set to zero in the above equation.Furthermore, we note that it is the partial fractioning step that prevents us from writing IBP relations for arbitrary powers of propagators, as it is usually done when traditional IBP relations are derived for integrals without θ-functions.Because of this, we need to generate equations by selecting seed integrals, deriving IBP relations for them and explicitly mapping all inhomogeneous terms to simpler integral families.We further elaborate on this point in the next section.In summary, we found that modified IBP relations in the considered case relate integrals with two θ-functions to integrals with one θ-function.Written as phase-space integrals, the modified IBP relations read where the two equations arise from derivatives w.r.t.k 1 and k 2 , respectively.We can apply the same logic to integrals with a θ-function and a δ-function.Inhomogeneous terms in this case will contain two δ-functions and no θ-function.IBP relations for such integrals can be derived using conventional methods since δ-function constraints can be immediately mapped onto rational functions of parton momenta using reverse unitarity.

Reduction to master integrals
In the previous section, we presented an explicit example of a modified IBP relation that we derived starting with a seed integral with two θ-functions.We note that modified IBP relations naturally form a hierarchical structure since the smaller the number of θ-functions that a particular integral contains, the easier it is to compute it.Hence, in the course of the reduction, we try to express integrals with larger number of θ-functions through integrals with some of the θ-functions replaced by the δ-functions.Integrals that belong to the same "hierarchical level" are organized into distinct integral families.
In practice, we implement derivation of modified IBP relations in Mathematica and solve them using the "user-defined system" functionality of Kira.This requires us to define integral families with different numbers of θ-and δ-functions and to derive relations among integrals that belong to different families.We note that for conventional integrals all these steps are done automatically by publicly-available reduction programs, but we have to take care of them ourselves in the present case.
As we already mentioned, it is straightforward to derive homogeneous relations for integrals in each family for arbitrary powers of propagators.For example, for the integral family defined in Eq. (4.6), the derivative with respect to k 2 , contracted with v = n yields In writing this equation, we have defined operators î+(−) that raise (lower) the index a i of the integral T ex a 1 ,a 2 ,a 3 ,a 4 ,a 5 ,a 6 ,a 7 by one.Inhomogeneous terms, such as the one that appeared in Eq. (4.10) and discussed after that equation, are generated on an integral-by-integral basis; they require partial fraction decomposition and "on-the-fly" matching to topologies with fewer θ-functions.Although this step is straightforward, it needs to be implemented in a separate Mathematica code whose output is then fed to Kira.
Using this setup, we express the maximally non-abelian contribution to the zerojettiness soft function defined in Eq. (4.1) in terms of eleven master integrals The master integrals that appear in Eq. (4.16) read and It is interesting that the above set of master integrals is actually redundant since for i = 1, 2, 3, 4, I nn i = I nn i .This happens because in certain cases integrals depend on α i and β i (or on n and n) in a symmetric way.As an example, consider I nn 2 .We find Another interesting feature of the above set of master integrals is that only three master integrals I 5...7 contain two θ-functions; all these integrals correspond to the nn configuration.For all other master integrals, either one or both θ-functions are replaced by δ-functions; these integrals are simpler to compute than the original ones.
To understand why there are no nn master integrals with two θ-functions, we note that homogeneous parts of IBP relations are unaffected by θ-functions; hence, by solving homogeneous parts of the IBP relations we should find master integrals that would be present if all θ-functions are removed from the integrand.It is then easy to see that, in the case of nn integrals, removal of θ-functions from the integrand leads to scaleless integrals since the jettiness constraint only depends on β 1,2 in this case.

Computation of master integrals
Having discussed the reduction to master integrals, in this Section we explain how to compute the master integrals that appear in Eqs.(4.17,4.18).Since the NNLO soft function is required for the computation of N3LO soft function, we will compute S 2g N A through weight six.
We begin by calculating the four integrals that are needed to describe the nn configuration.To this end, we combine the phase-space parameterization in Eq. (2.5) and the Sudakov parametrization of the phase-space element of gluon i to write where we have used The computation of I nn 2 proceeds in a similar way.We find The integral I nn 3 involves the scalar product of the two gluon momenta.To facilitate its computation, we write 12 where dΩ and ϕ 12 is the relative angle between transverse components of k 1 and k 2 .To integrate over this angle, we introduce a new variable η defined as and write dΩ We integrate over η using the formula apply the following hypergeometric identity and obtain Following similar steps, we derive where the function X 4 reads We note that the integral in Eq. (4.29) diverges logarithmically at the integration boundaries ξ 2 = 0, 1 but the function X 4 (ξ 2 ) is regular at these points Hence, we can compute the master integral I nn 4 by subtracting divergent contributions at endpoints and adding them back.Specifically, we write The first two integrals on the r.h.s of Eq. (4.32) are trivial.The last integral is regular in the integration domain ξ 2 ∈ [0, 1] and can be computed after expansion in .We construct such an expansion using the package HypExp [56] and use the program HyperInt [57] to integrate the result over ξ 2 .We arrive at where we have discarded contributions of weight seven and higher.It remains to compute three additional integrals for the nn configuration.We begin with I nn 5 and change variables We find where We subtract the (only) logarithmic singularity at ξ 1 = 0 and obtain after the integration Following the same steps as described above, we find the result for I nn 6 Finally, for the integral I nn 7 we obtain The ξ 1 -integral is finite; we expand it in powers of and integrate using HyperInt.The result reads This concludes the computation of master integrals required for the calculation of the real-emission contribution to the zero-jettiness soft function at NNLO.

Results for the real emission contribution at NNLO
We use the master integrals computed in Section 4.  The result agrees with the one derived earlier in Ref. [24]. 9o summarize, we have shown that by constructing the integration-by-parts identities for phase-space integrals with Heaviside functions it is possible to express the real emission contribution to the zero-jettiness NNLO soft function through seven master integrals.The majority of these integrals needs to be computed by integrating over a simplified phase space with all or some θ-functions replaced by the δ-functions.This simplification is very striking in case of the nn kinematic configuration where we find that no master integrals with two θ-functions need to be computed.As we discussed earlier, this interesting feature can be readily understood if IBP technology is applied to phase-space integrals with Heaviside functions.
5 Testing the method with some N3LO contributions to the zero-jettiness soft function It appears from the discussion in the previous section that it is useful to construct IBP relations for integrals with Heaviside functions.However, given an enormous growth in computational complexity with increase in the perturbative order, it is important to check this statement by considering a more complex example.Given our interest in the N3LO QCD contribution to the zero-jettiness soft function, it is natural to check if modified IBP relations can be constructed and used to compute it.
To this end, in this section we consider the maximally non-abelian part of the realemission contribution to the soft function with all gluons emitted to the same hemisphere.We define it as follows where the function ω nn (k 1 , k 2 , k 3 ) in Eq. (5.1) reads [59] ω ) In Eq. (5.3), "permutations{k 1 , k 2 , k 3 }" describes all possible permutations of the gluon momenta k i .The four terms S (t) ik , t = a, b, c, d in Eq. ( 5.3) are contributions to the soft eikonal function that are ordered according to the structure of their collinear singularities [59].
A simple generalization of the discussion in Section 3 implies that, in order to set up IBP relations for integrals that appear in S 3g nnn , we require eight distinct classes of integrals.They correspond to integrations over the following phase spaces I θθθ ∼ dΦ nnn θθθ , I θθδ ∼ dΦ nnn θθδ , I θδδ ∼ dΦ nnn θδδ , I δδδ ∼ dΦ nnn δδδ and four more cases where the ordering of θ-and δ-functions differs from the above examples.
Similar to the two-gluon case, these classes of integrals possess hierarchical structure that we rely upon when solving the integration-by-parts identities.Indeed, a closed set of IBP relations can be derived for I δδδ using reverse unitarity.On the other hand, IBP identities for I θδδ involve I δδδ , and IBP identities for I θθδ involve I θδδ and I δδδ .Finally, integration-by-parts identities for I θθθ make use of all other integrals with θ-and δ-functions.We illustrate how this approach applies to various N3LO QCD contributions to the zero-jettiness soft function in the next sections.
where p i = n, p k = n and q i = k i , i = 1, 2, 3 in our notations.As can be easily checked, S ik is not singular when any of the two gluons become collinear to each other.This feature reduces the number of scalar products that appear in the denominators of that function, making integration over three-gluon phase space simpler.
Applying integration-by-parts identities, we find that the integral of ω (3),a nn over the phase space dΦ nnn θθθ can be expressed through six master integrals.The result of the reduction reads dΦ nnn θθθ ω (3),a (5.5) Definitions of the six master integrals are given in Appendix A, Eq. (A.1).However, it is useful to emphasize at this point that none of these integrals contains three θ-functions.
We have explained why this happens when discussing the NNLO contribution to the soft function in Section 4. Calculation of master integrals that appear in Eq. (5.5) is rather straightforward.As an example, consider a phase-space integral I 1 with three δ-functions.It reads dΦ nnn δδδ . (5.6) Since the integrand does not depend on the relative orientation of the three partons in the transverse plane, we can integrate over dΩ for i = 1, 2, 3. We then remove all δ(α i − β i )-functions by integrating over α i , i = 1, 2, 3, and find (5.7)For a less trivial example, consider the integral I 5 .It reads To compute this integral, we remove two δ-functions by integrating over α 1,2 and then integrate over β 1 to remove δ(1 − β 123 ).Then, writing β 2 = (1 − β 3 )ξ, we arrive at It is straightforward to express the integrals over ξ and α 3 through hypergeometric functions.The result reads (5.10) where we introduced a new variable u = 1 − β 3 .Although both hypergeometric functions in Eq. (5.10) are singular at u = 1, this singularity is made integrable by an explicit factor (1 − u) 1−2 in the integrand.Hence, we can directly expand the integrand in Eq. (5.10) in powers of and integrate over u using HyperInt [57].We find where we retained all contributions through weight six.
The results for all other integrals which appear in Eq. (5.5) can be obtained along similar lines; they are given in an ancillary file.Using them in Eq. (5.5), we obtain a remarkably simple result for the integral of ω (3),a nn in the nnn configuration dΦ nnn θθθ ω (3),a showing the usefulness of applying integration-by-parts technology to phase-space integrals with Heaviside functions.
where p i = n, p k = n and q i = k i , i = 1, 2, 3 in our notations.At variance with S ik contains a scalar product of two gluon momenta q 1 q 2 which causes a collinear singularity in the limit q 1 q 2 and makes it more difficult to integrate the function ω (3),b nn compared to the discussion in the previous section.However, there is no reason to expect that things may work differently for this contribution.
We have, therefore, proceeded as before and performed a reduction to master integrals.Upon checking the results numerically, we have found that reduction to master integrals produces wrong results because not all integrals that appear in the integration-by-parts identities in the course of the reduction are regulated dimensionally.
It is interesting to point out that a) integrals that appear in the function ω (3),b nn (k 1 , k 2 , k 3 ) are well-defined in dimensional regularization and b) all master integrals that appear in the expression for the amplitude obtained after the reduction do not exhibit singularities that cannot be regulated dimensionally.This means that the failure of dimensional regularization is very well hidden in the internal dynamic of the reduction process and, therefore, hard to detect.
To remedy this situation, we introduced an analytic regulator in addition to the dimensional one.The analytic regulator is introduced in such a way that the phase-space element is multiplied by a factor This modification changes the integration-by-parts equations making them ν-dependent and, therefore, significantly more complicated.However, the main steps that we have described earlier in connection with deriving IBP equations and establishing a hierarchy of integrals is not affected by the analytic regulator.
Having performed the reduction for finite values of ν and , we have found that the integral of ω (3),b nn (k 1 , k 2 , k 3 ) can be written in a remarkably simple form which, however, demonstrates very clearly why the analytic regulator is needed.We obtain i ( , ν)J i ( , ν). (5.15) As indicated in Eq. (5.15), both the reduction coefficients and the integrals J 1,..,23 are functions of ν but, studying the ν → 0 limit of Eq. (5.15), we find that for i = 1, 2, .., 22 the following holds lim (5.16) The integrals I 1,..,6 have been already discussed in the previous section and I 7,..,22 are another sixteen integrals that can be evaluated at ν = 0.These integrals are given in Eq. (A.2).The last integral, J 23 ( , ν), has a 1/ν pole but C 23 ( , ν) is proportional to ν.Therefore, C 23 ( , ν)J 23 ( , ν) produces a finite result in the ν → 0 limit that, however, is completely missed if the reduction at ν = 0 is performed.
Following this discussion and taking the ν → 0 limit where appropriate, we write the result of the reduction as follows where we renamed J 23 to J ν .As we already mentioned the integrals I 1..22 that appear in the above expression can be found in Appendix A, Eqs.(A.1,A.2).The integral that is singular in the ν → 0 limit is defined as follows (5.18) We will now explain how this "pathological" integral can be computed.
To integrate over the relative azimuthal angle between k 1,⊥ and k 3,⊥ , we use the following formula where i = 1 and j = 3.To proceed further, we change variables α 1,3 → β 1,3 /ξ 1,3 and obtain the following integral representation (5.20) We make a further change of variables, ξ 3 = rf and ξ 1 = f in the first term in the square brackets and ξ 3 = f and ξ 1 = rf in the second one.We find It is convenient to write β 1 = x(1 − y), β 2 = xy and β 3 = 1 − x, 0 < x, y < 1. Upon doing that, we find that integrations over x and y can be readily performed.We obtain where (5.23) The two hypergeometric functions in the square brackets in the last equation can be conveniently re-written using the following identity (5.24) Inserting Eq. (5.24) into Eq.(5.23), we obtain (5.25) From the above equation we readily see how the 1/ν-singularity appears; it is generated by the singularity at f = 0 in Eq. (5.25) which is not regulated dimensionally.
Although we can compute Jν by expanding it in Laurent series in and ν, we only require the 1/ν singularity of this integral.We therefore compute the residue at f = 0 and find We now use Eqs.(5.27) Calculation of the master integrals I 1,2,...,22 required for computing the integral of the function ω (3),b nn proceeds in full analogy to the above case and for this reason we do not discuss it here.Finally, if we use the reduction to master integrals Eq. (5.17) and explicit results for the master integrals that can be found in an ancillary file, we obtain the following result for the integral of the function ω (5.28) Again, we see that, apart from an unexpected (and interesting) complication related to the need to introduce an analytic regulator, it is beneficial to make use of the IBP identities to compute non-trivial contributions to the zero-jettiness soft function at N3LO.

The ω (3),c nn contribution
The third contribution to the zero-jettiness soft function is associated with the integral of the function ω where p i = n, p k = n and ik contains propagators 1/(q 1 q 3 ) and 1/(q 2 q 3 ) and, therefore, exhibits a more complicated singularity structure.However, no additional issues arise when modified IBP relations are constructed and used to reduce required integrals, so that working with the analytic regulator introduced in the previous section, see Eq. (5.14), suffices.
Performing the reduction to master integrals and taking the ν → 0 limit where appropriate, we obtain All the master integrals I 1,..,39 that appear in Eq. (5.30) are computed at ν = 0; their definition can be found in Appendix A, Eqs.(A.1,A.2,A.3).The only ν-dependent integral J ν present in Eq. (5.30) is the same integral which we have seen (and computed) in the previous section.
In principle, computation of the I-integrals that appear in the calculation of the case "c" does not differ from what we already discussed for cases "a" and "b".However, integrals with two separate scalar products of parton's four-momenta are significantly more complicated than what we have discussed so far and, to illustrate this point, we will discuss how one of them can be computed.
We consider integral I 30 defined as We can use Eq. ( 5.19) to perform integration over azimuthal angles of k 2 and k 3 .We note, however, that the momentum k 2 is special because the δ-function in the definition of I 30 implies that α 2 = β 2 .This condition removes one of the two hypergeometric functions in Eq. (5.19).Integrating over azimuthal angles and changing variables α 1,3 = β 1,3 /ξ 1,3 , we arrive at the following representation for I 30 ( As the next step, we remove δ(1 − β 123 ) by integrating over β 2 and change variables (β 1 , β 3 ) → (x, y) according to the following formula β 1 = x(1 − y), β 3 = xy.Integration over x leads to a hypergeometric function.We obtain The integral naturally splits into two parts.To proceed, we change variables ξ 1 = rξ, ξ 3 = ξ, and ξ 3 = rξ, ξ 1 = ξ, in the first and the second integral in Eq. (5.33), respectively.We also rewrite the hypergeometric functions that appear in square brackets in the above equation using the transformation We then write where The difficulty with computing these integrals is a power-like singularity at y = 1.To extract and isolate it, we transform the y-dependent hypergeometric function in the following way This representation is helpful because the y-dependence of the first term in Eq. (5.38) is simple so that integration over y can be immediately performed, and the second term in Eq. (5.38) leads to integrals with only a logarithmic singularity at y = 1.
To proceed further, we consider integral I where the two terms correspond to the two terms on the right hand side of Eq. (5.38).To compute I (a,1) 30 we integrate over y and obtain We then use this result in the expression for I (a,1) 30 and write This integral has a singularity at r = 0 and another (overlapping) singularity at r = 1, ξ = 1.The two singularities can be separated by multiplying the integrand with 1 = (1 − r) + r.
The first term in the sum removes the r = 1 singularity.Since the r = 0 singularity does not overlap with any other singularity, it can be easily extracted.On the contrary, the r = 1 singularity overlaps with ξ = 1 singularity and for this reason a slightly more sophisticated treatment is needed.To this end, we subtract the product of hypergeometric functions at r = 1 and add it back.When the difference is used in the integrand, it becomes non-singular at r = 1, so that we can expand it in and integrate.On the other hand, hypergeometric functions in the subtraction term do not depend on r so that integration over r becomes straightforward.The resulting one-dimensional integration over ξ contains a logarithmic singularity at ξ = 1 that can be easily isolated and extracted.Putting everything together, we find (5.43) The computation of I (a,2) 30 proceeds in the following way.First, we observe that in this case there are two singularities, y = 1 and r = 1.We note that the latter overlaps with the ξ = 1 singularity.To disentangle overlapping singularities, we replace in Eq. (5.36) and add the difference of the two terms back.In the difference, the r = 1 singularity is regulated, so that one only needs to extract a logarithmic y = 1 singularity.To compute the contribution of the subtraction term (the r.h.s. of Eq. (5.44)), we can integrate over y and over r to obtain yet another hypergeometric function of ξ.The resulting integral over ξ has a logarithmic singularity at ξ = 1 which can be easily extracted.We obtain Computation of the integral 30 proceeds in a similar way except that the integration over ξ can be performed right away.We find We then rewrite 2 F 1 [{−1 − 6 , −2 }, {−6 }; y] using Eq.(5.38) and integrate the two terms that appear in that equation separately.This integration is relatively straightforward since integration of the first term in Eq. (5.38) leads to yet another 3 F 2 -function and integration of the second term does not require resolution of any overlapping singularities.
Putting everything together, we arrive at the following result for the integral I 30 (5.47) Finally, using reduction to master integrals and explicit expressions for integrals given in the ancillary file, we arrive at the following result for the integral of ω  (5.48)

Differential equations and ω
(3),d nn contribution Application of integration-by-parts identities to phase-space integrals with Heaviside functions opens up an opportunity to compute complicated integrals using differential equations.We will illustrate this by considering the contribution of the function ω (3),d nn , constructed using the function S 123 p k q123 p k q1 piq2 piq3 (5.49) where p i = n, p k = n and q i = k i , i = 1, 2, 3 in our notations.A distinct feature of the function ik is that it contains a propagator 1/k 2 123 which depends on the relative orientation of all three gluons in the plane transverse to the collision axis.This feature makes it very difficult to integrate the function ω (3),d nn analytically over the phase space with θ-functions.On the other hand, the ability to write down the IBP equations and perform reduction to master integrals should allow us to use differential equations to compute even the most complicated master integrals.This is what we would like to discuss in this section.
We begin by expressing the integral of ω (3),d nn (3),d nn and derive differential equations for these integrals w.r.t. the mass parameter m 2 .We emphasize that the very possibility to use differential equations for phase-space integrals with Heaviside functions depends on our ability to set up integration-by-parts identities and reduce the derivatives of master integrals back to basis integrals.
It is quite clear that additional mass parameter makes integration-by-parts identities more complex and requires us to introduce more integrals to close them.However, once the differential equations are derived, it is in principle straightforward to compute the required integrals by solving them numerically.We do this by using the method described in Ref. [60].In what follows, we first discuss calculation of boundary conditions and then explain how to solve differential equations by considering a (relatively) simple example.

Calculation of the boundary conditions
We have seen in the previous sections that, although IBP reductions have to be performed for a non-vanishing value of the analytic regulator, the limit ν → 0 can be taken in an absolute majority of integrals after the reduction is performed.The same applies to integrals that contain propagators 1/k 2 123 that we modify by introducing the mass parameter m.We will therefore discuss calculation of boundary conditions for such integrals, setting the analytic regulator to zero.
The complexity of the boundary conditions computation strongly depends on the type of constraints that a particular integral is subject to.As we explain below, the more δfunctions a particular integral has, the easier it is to compute the boundary conditions.To understand this, consider integrals with three δ-functions In Eq. (5.50) the ellipses stand for mass-parameter-independent scalar products and i is an integer number.In all integrals that refer to the nnn configuration, all β-variables are restricted because β i > 0, i = 1, 2, 3 and Moreover, for integrals with three δ-functions, all α-variables are equal to β-variables and k 2 i,⊥ = α i β i = β 2 i , i = 1, 2, 3. Hence, integration in Eq. (5.50) is performed over a finite region of the three-particle phase space, so that the approximate form of the integral in the m 2 → ∞ limit is simply obtained by expanding the propagator 1/(k (5.51) It is obvious that a Taylor expansion of I δδδ -integrals in k 2 123 /m 2 produces integrals where a "massive" propagator is absent.As the result, once the expansion is performed, we can use integration by parts for k 2 123 -independent integrals to express any expanded I δδδ integral through master integrals computed in the previous sections.
Consider now an integral that contains two δ-functions and one θ-function.We choose four-momenta k 1,2,3 in such a way that the θ-function depends on k 1 (more precisely on α 1 and β 1 ) and write At variance with I δδδ integrals, the integration over α 1 in I θδδ integrals is not restricted from above.This implies that, in the limit m → ∞, there are contributions from the region α 1 ∼ m as well as α 1 ∼ 1.Since the phase-space element scales as α − 1 and since in the limit α , we find that all integrals I θδδ have the following asymptotic dependence on the mass parameter m in the m → ∞ limit (5.53) In this formula i 1 , i are two integers that are particular to the integral under consideration, and A 1,2 are independent of the mass.To compute A 1 we need to Taylor-expand the integrand in k 2 123 /m 2 and then use integration by parts to reduce the resulting integrals to massless θδδ-integrals.To compute A 2 , we need to drop the θ(α 1 − β 1 ) constraint, as it is only relevant for α 1 ∼ β 1 ∼ 1, and then simplify a particular integral under the assumption that α 1 ∼ m 2 α 2 , α 3 , β 1 , β 2 , β 3 .We note that in this case all the relevant integrals can be straightforwardly computed in a closed form in terms of hypergeometric functions.
To illustrate how contributions of different regions can be computed, we consider one of the integrals with 1/(k 2 123 + m 2 ) propagator We would like to compute the leading contribution to O(m −2 ) branch of this integrals in the limit m → ∞.We consider two distinct contributions α 2 ∼ m 2 1 and (5.57) However, in the case α 3 ∼ m 2 , there are additional dependences of the integrand on k (5.58) We conclude that the leading m → ∞ contribution to the O(m −2 ) branch of the integral B 1,δθθ in Eq. (5.56) arises from the region where α 2 ∼ m 2 and all other variables are of order one.
To compute this contribution, we notice that, upon taking the limit 123 + m 2 , the only dependence of the integrand on relative azimuthal angles of massless partons resides in the simple scalar products, e.g.k 1 k 3 in B 1 .This integration can be easily performed following the discussion in the previous sections.We find (5.59) 11 As we explained earlier, there are no master integrals with three theta-functions in case of the nnn configuration.
Therefore, to determine the O(m −4 ) branch of any I θθδ integral, we need to compute At this point, we can again use the integration-by-parts method and express all relevant integrals F (q 2 , qn, qn) as linear combinations of five master integrals.They read (5.66) Note that these integrals do not contain scalar products of the gluon four-momenta.This happens because such scalar products are either simple as e.g. in case of k 1 k 2 = q 2 /2 or can be simplified for large α 1,2 as e.g. in It is convenient to compute the integrals B 1,..,5 in the rest frame of q where k 1,2 are back-to-back.Then, using the following result for the angular integral dΩ we easily find where

.69)
A typical function F (q 2 , qn, qn) is given by a linear combination of the integrals B 1,..,5 .For example, for one of the boundary integrals that we will refer to as B 2 this function reads (5.70) Hence, to determine the m −4 branch of the corresponding integral, we need to compute qn, qn). (5.71) To calculate the integral in Eq. ( 5.71), we need to choose a convenient parameterization to integrate over q.We do this in the following way.We introduce the Sudakov decomposition for the vector q and write q = 1 2 α q n + 1 2 β q n + q ⊥ .It follows that (5.72) In addition, the four-vector q needs to be time-like, q 2 > 0. This implies q 2 = α q β q − q 2 ⊥ > 0. (5.73) We would like to simplify the mass-dependent denominator in Eq. ( 5.71).To do that, we first re-write it using Sudakov variables where we introduce a new variable t = α q − q 2 ⊥ .Positivity of q 2 requires α q > q 2 ⊥ β q . (5.75) (5.76) Hence, we can write where the integration boundaries are 0 < β q < 1, 0 < α q < ∞ and α q (1 − β q ) < t < α q .Using the fact that the integrand in Eq. (5.71) does not depend on the azimuthal angle, we integrate over dΩ (d−2) and write (5.78) To integrate further, we change variables α q → ξ with α q = t/ξ and 0 < ξ < 1.Since q 2 = t − α q (1 − β q ) = t ξ (ξ − (1 − β q )) > 0, (5.79) the integration boundary for β q becomes 1−ξ < β q < 1.To accommodate these boundaries in a natural way, we change variables β → r, with β = 1 − ξr, 0 < r < 1. Upon changing variables and using explicit expressions for integrals B 1,2,4 , we note that integration over t factorizes and can be performed easily.We obtain  (5.81) Calculation of O(m −4 ) branches for other integrals that are needed to determine boundary conditions proceeds along similar lines.Integration over t can always be performed exactly and the subsequent integration over ξ and r is then completed by constructing subtraction terms of end-point singularities to facilitate expansion of integrands in .

Numerical solution of the differential equations
Having discussed the computation of the boundary conditions, we need to explicitly write down and solve the differential equations.We note that we need to consider about two hundred mass-dependent master integrals in total to close the system of differential equations.Since it is impossible to discuss such a huge system of differential equations in any detail, we decided to choose a small eleven-by-eleven sub-system and discuss it in a comprehensive way.
We consider eleven integrals

Conclusion
In this paper we discussed computation of real-emission integrals for observables that contain Heaviside functions.This is an interesting problem because reverse unitarity [42] cannot be immediately applied to map such integrals onto multi-loop integrals, preventing straightforward use of integration-by-parts identities in such cases.We discussed a way to re-introduce integration-by-parts technology into the computation of such integrals and showed that the resulting IBP relations have a clear hierarchical structure since, in addition to original integrals, there appear integrals with Heaviside functions replaced with δ-functions.Integrals with δ-functions are, however, much simpler since they can be dealt with using reverse unitarity and, thus, IBP equations for them are self-contained.In addition, IBP relations provide a foundation for deriving differential equations for real-emission master integrals with Heaviside functions that can be solved numerically even if analytic integration becomes too difficult.

1 Introduction 1 2 3 3 4 . 1 5 38 A 38 B
Heaviside functions in the zero-jettiness soft function Applying integration-by-parts technology to integrals with θ-functions 5 4 IBP identities and the NNLO QCD contribution to the zero-jettiness soft function 7 An example of an IBP relation 7 4.2 Reduction to master integrals 10 4.3 Computation of master integrals 11 4.4 Results for the real emission contribution at NNLO 14 Testing the method with some N3LO contributions to the zero-jettiness soft Master integrals for nnn contribution to the soft function Matrices for the differential equations in Section 5.4.2 41

5. 1 10 S
The ω (3),a nn contribution The simplest contribution to consider is an integral of ω (3),a nn .According to Eq. (5.3), this function is constructed from the function S (a)ik , which reads[59]

5. 2
The ω (3),b nn contribution We consider the second contribution to the zero-jettiness soft function given by an integral of the function ω (3),b nn .This function is constructed from the function S piq1 piq3 p k q2 p k q3 (−pip k piq12)

2 , 3
in our notations.At variance with the functions S (a,b) ik , the function S (c) ik defined in Ref.[59], to the zero-jettiness soft function.The function S
2 123 + m 2 ) in k 2 123 /m 2 .It follows that all I δδδ integrals have a particularly simple asymptotic mass dependence lim m→∞