N -jettiness soft function at next-to-next-to-leading order in perturbative QCD

: We derive a compact representation of the renormalized N -jettiness soft function that is free of infrared and collinear divergences through next-to-next-to-leading order in perturbative QCD. The number of hard partons N is a parameter in the formula for the finite remainder. Cancellation of all infrared and collinear singularities between the bare soft function and its renormalization matrix in color space is demonstrated analytically.


Introduction
The need for precise description of high-multiplicity final states produced in hard processes at the Large Hadron Collider (LHC) led to many studies of the so-called subtraction and slicing schemes for perturbative QCD computations .The role of these schemes is to enable cancellation of infrared and collinear divergences between virtual and real contributions to cross sections without compromising the fully-differential nature of theoretical predictions.
In spite of the fact that slicing and subtraction schemes aim at achieving identical goals, their discussions in the literature proceed quite differently.Indeed, a typical pathway to the construction of a subtraction scheme involves studies of how factorization of matrix elements in the soft and collinear limits, and factorization of phase spaces for final-state particles can be combined to make the subtraction terms both observable-independent and analytically (or numerically) integrable in d-dimensions. 1On the contrary, advanced slicing schemes start from the choice of an observable for which physical cross sections factorize into universal quantities, in the limit where the selected observable is small.
In spite of a somewhat different philosophy, it is clear that a very tight connection between slicing and subtraction approaches must exist since variables employed for slicing calculations are infrared safe.It follows that a differential cross section for the small value of the selected slicing parameter -which can be considered as the major quantity to understand in order to enable the construction of a slicing scheme -should be computable with subtractions.
We believe that existing computations in the context of existing slicing schemes do not make this connection sufficiently clear and do not exploit it.To see this, consider available computations of the N -jettiness soft function at next-to-leading (NLO) [16] and at next-tonext-to-leading order (NNLO) [17,18,[31][32][33][34][35][36][37] in perturbative QCD.These computations are performed by mapping the available phase space of soft-gluon emissions onto hemispheres along the different hard directions, computing many of the required integrals numerically, and studying cancellation of infrared and collinear divergences against the soft-function renormalization matrix a posteriori.Furthermore, the soft-function renormalization matrix is constructed from the consistency requirement that soft and collinear divergences cancel when different elements of the factorization theorem (the hard function, the beam function, the jet function and the soft function) are combined to produce a physical cross section.While this computational method works fine in practice, it is in stark contrast with the current developments in understanding subtraction schemes where one attempts to find suitable subtraction terms for a generic N -jet problem, integrate them over unresolved phase spaces, show analytically the cancellation of all the 1/ϵ poles and derive a compact representation for the finite remainder [24,38].
The goal of this paper is to show that borrowing certain ideas from the recent developments in constructing generic NNLO QCD subtraction schemes [38] is beneficial for computing essential ingredients for modern slicing computations.Our prime example is the N -jettiness soft function [16] at NNLO QCD that we study for a generic value of N .We analytically demonstrate the cancellation of all soft and collinear 1/ϵ divergences that appear in the calculation of the N -jettiness soft function against the corresponding renormalization constant, and derive a simple formula for the finite remainder that is valid in four space-time dimensions.We note that very recently a calculation of N -jettiness soft function for arbitrary N was reported in Ref. [37].The computational method used in that reference is, however, very different from ours since it relies on the extension of numerical methods developed in Refs.[35,[39][40][41].
The remainder of the paper is organized as follows.In the next Section we define the N -jettiness variable [42] and the soft function that we study in this paper.In Section 3 we discuss the renormalization of the soft function.We continue with the calculation of the N -jettiness soft function at NLO QCD in Section 4. We then proceed to the analysis of the NNLO case.We start by discussing in Section 5 how the NNLO QCD computation of the N -jettiness soft function can be organized.We first focus on the two-gluon final state.We discuss the contribution of the uncorrelated emissions of two soft gluons in Section 6.We analyze the triple-color correlated terms in N -jettiness soft function in Section 7. Correlated emissions of two gluons are discussed in Sections 8, 9, 10 and 11.We study the contribution of the soft q q pair to N -jettiness soft function in Section 12.The final result for the soft function can be found in Section 13.Numerical implementation and comparison with results available in the literature is described in Section 14.We conclude in Section 15.Some useful formulas are collected in appendices.

Definitions of basic quantities
We consider partonic processes with an arbitrary number of hard color-charged partons that we denote by N The first process refers to e + e − annihilation, the second to deep-inelastic scattering and the third to hadron collisions.Each parton in Eq. (2.1) carries a four-momentum p i and a color charge T i .All partons are assumed to be massless.
The cross sections of processes shown in Eq. (2.1) are affected by the radiation of soft quarks and gluons, both real and virtual.To quantify their impact, an observable that controls the softness of the radiation is needed.In this paper, we will use the N -jettiness T [42] for this purpose.
To define T , we proceed as follows.For a generic process, we split all final-state partons that we have to consider in a particular order of perturbation theory into resolved R and unresolved U ones.We define a list where x ∈ U and h i ∈ R. Quantities P h i are arbitrary normalization constants.We then define the N -jettiness variable as [42] T (R, U) To compute radiative corrections to processes in Eq. (2.1), we need to consider virtual and real-emission corrections.Since we work through next-to-next-to-leading order, we need to include virtual corrections up to two loops, one-loop corrections with and without one soft parton and the real emission processes with up to two soft partons.For virtual corrections to processes without additional soft radiation, U is an empty set and so T is zero.The same is true for leading order processes shown in Eq. (2.1).
We note that virtual corrections do not need to be computed explicitly.Indeed, since infrared and collinear singularities of a generic two-loop amplitude are known [43][44][45], it is customary to separate them from the short-distance contribution to the two-loop amplitude known as the "hard function" and include them into the renormalization constant of the soft function, that we discuss in the next section.Hence, to compute the soft function at NNLO we need to study contributions with one or two soft partons to cross sections of processes in Eq. (2.1), subject to N -jettiness constraint.
The perturbative expansion of the bare soft function can be written as where τ is the value of N -jettiness and a s = α s (µ)/(2π) is the renormalized QCD coupling constant.The quantities S τ,i (τ ), i = 1, 2, can be computed starting from universal functions that describe the behaviour of QCD matrix elements when some of the radiated partons are soft.Before proceeding, it is useful to point out that it is very convenient to work with the Laplace transform of the soft function [37].The Laplace transform is defined as follows (2.5) As we will see from the explicit computations described below, the dependence of functions S τ,i , i = 1, 2, on τ is very particular and is given by simple powers τ −1−nϵ , where n = 2, 4 etc.For these functions, the Laplace transform evaluates to Note that in the last step we introduced a variable ū = ue γ E to simplify the expansion of the corresponding Γ-functions in ϵ.In what follows we will use Eq.(2.6) to compute the Laplace transform of the perturbative expansion of the bare soft function.

Renormalization of the soft function
Perturbative computations of the soft function lead to expressions that contain soft and collinear divergences that manifest themselves through 1/ϵ poles.These divergences are removed by a dedicated renormalization [16].
The renormalization of the soft function is a multiplicative (matrix) renormalization in the Laplace space [37].We denote the Laplace transforms of the bare and renormalized soft functions by S and S, respectively.The two functions are related by the following equation where Z is a matrix in color space, see Appendix A. While it is entirely possible to compute the unrenormalized soft function S and then remove the divergences by using Eq.(3.1), it is beneficial to study suitable combinations of Z and S that are actually needed for computing S.
To find these combinations, we write expansions of Z, S and S in powers of α s substitute them into Eq.(3.1) and solve for S1,2 .We find It was recently pointed out in Ref. [38], in the context of the application of the nested soft-collinear subtraction scheme to arbitrary multi-jet processes, that it is beneficial to separate iterations of O(α s ) soft, soft-collinear and virtual contributions to cross sections from the rest of NNLO contributions.Following the same logic, we write It is useful to do the same for the NNLO contribution to the renormalization matrix It is then easy to see that S2 can be written as We will show below that this representation is beneficial for computing the renormalized N -jettiness soft function.Indeed, when representation in Eq. (3.6) is used, the required effort is minimized because cancellations of 1/ϵ poles between the different quantities are identified relatively early in the course of the computation.However, before discussing how the NNLO computation is performed, we will illustrate our approach by calculating the NLO contribution to the N -jettiness soft function.

N -jettiness soft function at NLO
We will explain main ideas of our approach by considering N -jettiness soft function at NLO.Since one-loop virtual corrections to N -jettiness soft function do not need to be considered, we focus on the single real-emission contribution below.We use label m to denote a soft gluon and write its momentum as p m .The gluon m is the only unresolved parton, therefore U = {m}.We choose , where E i is the energy of the parton h i in the chosen reference frame, and obtain (c.f.Eqs.(2.2,2.3)) In Eq. (4.1) E m is the energy of the unresolved gluon m and ψ m is a function defined as where ρ ij = 1 − ⃗ n i • ⃗ n j and ⃗ n i is a unit vector which points in the direction of the threemomentum of parton i.
The real-emission single-gluon soft contribution at fixed N -jettiness τ reads S,ij (τ ), ( where we use the notation (ij) to indicate that the sum runs over all i, j ∈ R with the constraint i ̸ = j.Furthermore, T i is the color charge operator of parton i, I g s is the bare QCD coupling constant and is the soft eikonal function.It is convenient to write the phase space element of the gluon m as follows We then combine the coupling constant squared with the remnant of phase space normalization factor and note that where a s (µ) = α s (µ)/(2π) and α s (µ) is the QCD coupling constant renormalized in the MS scheme.We work to first order in a s and use Eq.(4.6) to integrate over energy of the gluon E m in Eq. (4.4).We find [dΩ where ⟨..⟩ m indicates integration over directions of the vector ⃗ n m .
To proceed further, we compute the Laplace transform of the soft function and obtain Integration over directions of the vector ⃗ n m induces collinear singularities when ⃗ n m ||⃗ n i and ⃗ n m ||⃗ n j .To isolate them, we make use of the fact that collinear singularities of the integrand are fully controlled by the remnants of the eikonal function in Eq. (4.10) and that the collinear limits of the function ψ m can be predicted for an arbitrary hard process.Indeed, consider the collinear limit m||i as an example.In that limit, ρ im is the smallest ρ-value so that lim It turns out useful to rewrite the integrand in Eq. (4.9) as follows We then define where the function g ij,m contains higher powers of ϵ in addition to the first power shown explicitly.Hence, we write Since the function ϵg ij,m vanishes in all collinear limits and is O(ϵ), it only provides a finite contribution to the N -jettiness soft function that will have to be computed numerically for a given configuration of hard partons.On the contrary, the first integral in the right-hand side of Eq. (4.14) can be computed explicitly, see Eq. (B.5).We write where Putting everything together, we find It follows from the above equation that all 1/ϵ poles have been separated from the finite remainder and that the remnant of the complicated N -jettiness function appears in the finite remainder only.
It remains to combine S 1 with the renormalization matrices Z 1 and Z + 1 , c.f. Eq. (3.3).The expression for Z 1 can be found in Appendix A; similar to the NLO soft function shown in Eq. (4.17), it is given by the sum of products of the color-charge operators with some coefficients.Combining the different contributions and discarding all terms beyond O(ϵ 0 ), we find where This result can be easily evaluated numerically for an arbitrary number of hard partons N .The logarithm ln(ψ m ρ ij /(ρ im ρ jm )) provides an infrared regulator for an integral that would have exhibited collinear divergencies without it.

N -jettiness soft function at NNLO
We will now explain how to extend the above approach to compute the NNLO contribution to the N -jettiness soft function.It is clear that the NNLO computation is significantly more complex than the NLO one.However, we will show that it is possible to split the calculation into several independent parts making the entire problem simpler.
The NNLO contribution to the bare soft function reads where S 2,RR is the double real-emission contribution, S 2,RV is the real-virtual contribution and the last term arises because of the renormalization of the strong coupling constant in the NLO soft function.The leading order QCD β-function is defined as where C A = N c is the Casimir of the group SU (N c ), T R = 1/2 and n f is the number of massless quarks.The double-real emission contribution is obtained by integrating the double-soft eikonal current [46] over the phase space of two soft gluons subject to the N -jettiness constraint.Following Ref. [46], we further split the double-real contribution into correlated and uncorrelated pieces S 2,RR,τ = S 2,RR,τ,T 4 + S 2,RR,τ,T 2 . (5.3) We note that the subscript τ indicates that the Laplace transform of the corresponding quantities has not yet been computed.The uncorrelated contribution reads where and the eikonal functions S ij (m) and S kl (n) are defined in Eq. (4.5).The correlated contribution S 2,RR,T 2 reads where The eikonal function Sgg ij (m, n) can be found in Appendix C.
The real-virtual contribution reads [46] where and λ ij = 1 if both i and j refer to incoming or outgoing partons and zero otherwise.Furthermore, the normalization factors in Eq. (5.8) read and we have defined F kij = f abc T a k T b i T c j to describe anti-symmetrized product of three color charge operators.Finally, the functions I RV,ij and I kij read and (5.12) In the next section we will show that (5.13) Therefore, the quantity S 2,r introduced in Eq. (3.4) reads (5.14)

Uncorrelated emissions of two soft gluons
We have argued that S 2 contains an iterated contribution of the next-to-leading order soft function S 1 and that it is beneficial to separate it from the rest of S 2 , c.f. Eq. (3.4).We will now show that this iteration can be identified with the Laplace transform of the term S 2,RR,τ,T 4 in Eq. (5.4).To this end, we write It is straightforward to integrate over energies of two gluons E m,n .Indeed, the first integration (say over E n ) is elementary since it just removes the δ-function.The second integration over E m extends from 0 to E max m = τ /ψ m , i.e. the energy of the gluon m which corresponds to E n = 0.The integral turns out to be of a Beta-function type, and we obtain Using this result, and performing the Laplace transform, we derive (ij),(kl) where the last step follows from the symmetry between (ij) and (kl) summation indices, and from the comparison with the results for the NLO N -jettiness soft function derived in Section 4. As we already explained in Section 3, it is straightforward to combine the above result with the iterated contribution of the renormalization constants to arrive at the relevant contribution to the renormalized soft function.

Terms with three color charges
We will next discuss the computation of the contribution that depends on triple products of color charges.It originates from the following terms in the expression for S2 where the real-virtual triple-color correlated contribution is given by the Laplace transform of the last term in Eq. (5.8).
The commutators are easy to compute using a general formula in Ref. [38]. 2 We find where the last step follows from the definition of L ij in Eq. (4.19) and the (anti)-symmetry of F kij with respect to permutations of its indices.
The second commutator reads where in the last step we replaced k ↔ j and used the symmetry of λ kj and the antisymmetry of F kij .Furthermore, it is convenient to replace λ kj in the above expression with κ kj defined in Eq. (5.9).This is possible because The first of these equation holds because the summation over i and k involves a symmetric and an anti-symmetric tensor.The second is fulfilled because one can sum over j using conservation of color charges and use the antisymmetry of F kij to set the resulting terms to zero.Hence, we can write The real-virtual triple-color correlated contribution evaluates to We need to combine the three contributions to Stc 2 shown above, demonstrate that all 1/ϵ poles cancel and derive a representation for the finite remainder.We will do this in steps starting from contributions that exhibit dependence on the N -jettiness function.
There are two such contributions -the commutator in Eq. (7.5) and S RV,tc .We define Both terms in Eq. (7.7) involve integrals of various combinations of ρ ij 's and ψ m over the directions of the momentum of the gluon m.Following the discussion of the NLO contribution to the soft function, we write ki,m ) Using Eq. (7.8), and focusing on the divergent contributions that contain g (2) and g (4) , we find ki,m − g where the last step follows from the equality g (2) = g (4) through order O(ϵ 0 ).Extending the above computation in a straightforward way, we derive the finite Njettiness-dependent triple-color correlated remainder.We obtain Next, we discuss triple-color correlated contributions that do not involve jettiness functions.The first one is given in Eq. (7.2); it contains explicit 1/ϵ 2 divergence and does not require further discussion.The second contribution is extracted from Eq. (7.5) using Eq.(7.8).The required angular integral is straightforward to compute and we do not discuss it further.
It is slightly more challenging to calculate a similar contribution to S RV,tc .However, we can cast it into a form that allows us to expand it in powers of ϵ.The corresponding expression reads The first two integrals on the right-hand side are straightforward to compute but the last one is challenging.Since this integral is O(ϵ), it contributes to the divergence of the triple-color correlated term so that we need to analytically compute it.
The simplest way to do this is to use the results of Ref. [38] where the analytic expression for the following integral was derived.The idea is to rewrite the integral in Eq. (7.12) in the same way as the integral in Eq. (7.11) and then take the difference of the two results.Upon doing that, one finds that the complicated finite integral appears at order O(ϵ 2 ) only; this integral can be computed numerically alongside with N -jettiness dependent contributions.
Putting everything together, we find that all divergences in triple-color correlated contributions cancel, and a finite result is obtained.It reads where We note that the function L ki is given in Eq. (4.19), the function W kij reads and the function Ḡikj r,fin can be found in Eq. (H.16) of Ref. [38].

The remaining contributions
We are left with the discussion of the contributions to the soft function that are not iterations of next-to-leading order terms, and do not contain triple-color correlators.We write them as where S 2,RR,T 2 is defined in Eq. (5.6) and is the part of the real-virtual contribution that depends on the products of two color-charge operators.
The calculation of I RV,ij is straightforward.We obtain Next, we perform the Laplace transform and find This contribution contains both explicit 1/ϵ terms and collinear divergences that arise after integration over directions of the gluon m.The latter can be isolated by following what has been done for computing the soft function at next-to-leading order.We also note that terms that involve the jettiness function ψ m start contributing at order 1/ϵ 2 .We will discuss how the jettiness-dependent divergent contributions cancel out once we isolate similar terms in S 2,RR,T 2 in the next Section.9 Correlated emissions of two soft gluons: the strongly ordered limit We consider the remaining contribution to S 2,RR,T 2 .It reads where Since the function Sgg ij is symmetric under the permutation of m and n, we can introduce the energy-ordering E n < E m (E n = ωE m ) and use the fact that the scaling of Sgg ij with the common gluon energy "scale" E m is uniform.We then write where ψ mn = ψ m + ωψ n .We note that from now on, for simplicity of notation, we will always assume that in the evaluation of Sgg (m, n) the energy of the gluon m is set to 1 and the energy of the gluon n is set to ω.Finally, we perform the Laplace transform and find where The representation of function I ij in Eq. (9.4) is the starting point of our analysis of the correlated emissions of two gluons.To compute I ij we identify and subtract various singular limits of the integrand in Eq. (9.4).We begin with the ω → 0 limit which we will refer to as the "strongly-ordered" one, in the sense of the energy ordering.The stronglyordered limit of Sgg ij evaluates to We then need to compute To integrate over ω, we change the integration variable ω → 1/ω and find We then rename m ↔ n and obtain we finally find The last term in Eq. (9.10) is fully factorized and can be computed in the same way as the uncorrelated contribution.We only need to understand how to deal with the first term in Eq. (9.10).
To this end, we note that the angular integral is singular in the three collinear limits m||i, n||j and m||n.To regulate the first two, we define the following functions Then, we split the integral into contributions with and without functions f im and f jn .This splitting introduces unregulated singularities in some integrals; to address this issue, we modify the integral in Eq. (9.11) by incorporating an additional (analytic) regulator, and write it as Using the functions f im , f jn defined in Eq. ( 9.12), we rewrite it as follows We need to compute the above integral and take the limit ν → 0 before the expansion around ϵ = 0 is performed.It is clear that in the second, the third and the forth terms in the right-hand side in Eq. (9.14), we can integrate over directions of n or directions of m or directions of both.We begin with the computation of the last integral in the right-hand side of Eq. (9.14) where we first integrate over n.The integration is straightforward [47] and we obtain , this integral appears to be proportional to the analytic regulator.Because we are supposed to take the ν → 0 limit at fixed ϵ, it is tempting to conclude that the whole integral vanishes.However, this conclusion is wrong because the remaining integral over directions of ⃗ n m is ill-defined for ν = 0.
To extract the 1/ν singularity from the remaining integration over ⃗ n m in Eq. ( 9.15) we note that it originates from a kinematic configuration where m||j.Hence, to extract it, we take the collinear limit in all quantities of the integrand except 1/ρ 1−ϵ+ν jm , whose integration produces the 1/ν singularity.Taking the ν → 0 limit, we find It is easy to convince oneself that the second and the third integrals in Eq. (9.14) cannot be rescued using this mechanism and just vanish, and the first integral does not need the analytic regulator since all divergences are dimensionally regularized.
Hence, we conclude that the angular integral needed to describe the strongly-ordered soft limit is given by the following equation The usefulness of this representation stems from the fact that in the remaining integral on the right-hand side m||i and n||j collinear singularities are regulated since f im and f jn vanish in these limits.The singular limit which remains is the m||n one.We therefore subtract the m||n singularity and write3 where Cmn = I − C mn .We then simplify the last term.We write The final result for the integral I so ij defined in Eq. (9.10) reads We use this result to write the strongly-ordered contribution to the uncorrelated realemission part of the N -jettiness soft function as follows We continue with the discussion of how to extract divergences from S so 2,RR,T 2 .We note that the first term in the sum on the r.h.s. in Eq. (9.21) is already finite because functions f mi ∼ f nj ∼ O(1) vanish in the collinear limits m||i, n||j, respectively and the singularity at m||n is regulated by Cmn operator.Since the ϵ-dependence of the second term on the right-hand side in Eq. (9.21) is explicit, we need to focus on the remaining two terms.We rewrite them to make their dependencies on the N -jettiness constraint explicit.We find and When these equations are used in Eq. (9.21), higher order terms in the expansion of functions g (4) and g (2) in ϵ become necessary.Fortunately, it is easy to compute them since Furthermore, since these functions regulate all divergences in the angular integrals where they appear, Eq. (9.24) can be expanded in powers of ϵ before integration over directions of the momentum of gluon m is performed.Comparing the divergent contributions of all ψ-dependent terms in Eq. (9.21) with similar contributions to the real-virtual soft limit in Eq. (8.4), we find that they cancel.The remaining strongly-ordered contributions are accounted for when the final result for the soft function is assembled.

Collinear subtractions
We continue with the discussion of the uncorrelated contribution to N -jettiness soft function.Our starting point is Eq.(9.4).We have computed the strongly-ordered contribution to that equation and we will consider the remaining contribution in this section.Hence, we write where N u can be found in Eq. (9.5), Sω = I − S ω and S ω is the operator that enforces the strongly-ordered limit.Furthermore, the eikonal function in Eq. ( 10.1) should be evaluated assuming that the four-momentum of the gluon m is (1, ⃗ n m ) and the four-momentum of the gluon n is ω(1, ⃗ n n ).
The integrand in Eq. ( 10.1) possesses a complicated singularity structure.We separate double-collinear and triple-collinear singularities introducing partition functions (see e.g.Ref. [38]) and write Sω [ where the first term describes the contribution of the double-collinear and the second one -of the triple-collinear partition.We will discuss them in turn.
We begin with the double-collinear contribution.It reads By construction, the partitions ensure that the only allowed singularities are m||i, n||j, and m||j and n||i, respectively.However, it is easy to see that the integrand is not singular in these limits.Hence, we can expand ψ 4ϵ mn in powers of ϵ in Eq. ( 10.3).The first term of this expansion corresponds to the contribution of double-collinear partitions to the soft-subtracted integral of ω 2 Sgg ij without the N -jettiness constraint.Since this term is multiplied by 1/ϵ, we need to compute it analytically if we want to show the cancellation of all ϵ-singular terms without relying on numerical calculations.We will show that it is possible to extract this contribution from the result of Ref. [48], so that its computation is not required.The jettiness-dependent term that arises in the expansion of Eq. ( 10.3) in powers of ϵ, contains a factor 4ϵ ln ψ mn and, thus, only provides an ϵ-finite contribution that we evaluate numerically.
The term I tc ij is a triple-collinear contribution.In this case we write where w tc = w mi,ni + w mj,nj .( To facilitate the computation of Sω [I tc ij ], we introduce the standard sectors to further partition angular integrations, see e.g.Refs.[1,9,10,38].We also assume that when we subtract m||n singularity, we need to simplify the phase space but a similar simplification is not needed for the triple-collinear limit.Finally, we note that there are no double-collinear m||i, n||i and similar singularities.Hence, we only need to subtract triple-collinear and m||n double-collinear singularities in the integrand in Eq. (10.4).With this in mind, we write4  Since the 1/ϵ divergence in the last term in Eq. (10.6) arises from a prefactor, we can expand the integrand in powers of ϵ.In particular, we apply the expansion to the function ψ 4ϵ

Sω [I
mn .The term with ln ψ mn is O(ϵ 0 ) and the term without it provides an O(1/ϵ) divergent contribution.However, this contribution is the same as in the case when no N -jettiness constraint is imposed.To make this more explicit, we combine double-and triple-collinear partitions, and write In the above formula, the first, the second, the third and the fourth terms are divergent.Our strategy will be to compute the first and the second terms explicitly, and extract the third and the fourth terms from the result in Ref. [48] making use of the fact that they do not depend on the N -jettiness function ψ mn .
To make this connection explicit, we write a formula that describes the extraction of singularities of the integral of the double-eikonal function in a situation when no Njettiness constraint is imposed.This quantity, computed in Ref. [48], can be cast into the following form where We can now write an expression for J ij , where the various singularities are extracted.We obtain (10.10) Since the left-hand side of the above equation was computed through O(ϵ 0 ) in Ref. [48], we can use this result to determine the last two terms that appear on the right-hand side in Eq. (10.10) and also enter Eq. (10.7).This is a sensible thing to do because their computation appears to be particularly challenging.Clearly, to make this happen, we have to compute all the other terms that appear in Eqs.(10.7,10.10)through finite terms in the expansion in ϵ.We describe details of the required computations in what follows.

The double-collinear m||n contribution with N -jettiness constraint
To proceed, we need to compute the m||n limit of the soft-subtracted eikonal function Sij , given by the first term in Eq. (10.7).In principle, this is straightforward, except for the fact that the m||n singularity appears to be stronger than logarithmic.Although this is not the case, it implies that the corresponding contribution needs to be carefully extracted.We find where the three-vector κ n is defined through the following equation with the additional constraints We also need the collinear limit of the phase space in b-and d-sectors.The corresponding parametrization can be found in Refs.[1,38].We obtain [dΩ where

.15)
As the next step, we need to integrate over azimuthal angles parameterized by dΩ and dΛ in Eq. (10.14).To facilitate this, we write where v µ = (n µ i /ρ im − n µ j /ρ jm ), and the four-vector κ satisfies where t µ = (1, ⃗ 0) and e µ m = (0, ⃗ n m ).The averaging of κ µ κ ν over [dΩ The derivation of this formula as well as the definition of vector r µ i are given in Appendix F in Ref. [38].
To further simplify this expression, we write the metric tensor of the transverse space as follows Since we find t µ v µ = −e m,µ v µ .(10.21)This allows us to write A straightforward computation gives and We note that the function W m ij develops a singularity in the limit m||j but this singularity is regulated by the partition function.
Finally, we note that the contributions of sectors b and d are equal.Integrating over x 4 and adding the two sectors, we obtain In Eq. (10.25), we introduced the following quantities and the analogs of the cusp anomalous dimension It is convenient to write the result for the sum of i-and j-partitions.To do this in an optimal way, we need to rearrange contributions that appear in Eq. (10.25).The key point in this rearrangement is to make the divergent contribution symmetric with respect to the replacement i ↔ j.We illustrate how this is done by considering the term proportional to γ g in Eq. (10.25).We write and observe that the first term in the right-hand side of Eq. (10.28) can be easily combined with the m||j, n||j partition and the last term is finite.Therefore, it can be computed by expanding in ϵ before integration.Note that in this term the ψ m -dependent contribution appears at order O(ϵ 2 ) which, given the 1/ϵ 2 prefactor in Eq. (10.25) is the last (finite) order in the ϵ-expansion that we care about.Using these considerations, we combine contributions of the m||i, n||i and m||j, n||j partitions and obtain We note that we have used w mi,ni m||n + w mj,nj m||n = 1. (10.30) The term A ij reads As we already mentioned, the second and third terms on the right-hand side of Eq. (10.31) can be expanded in ϵ because expressions in square brackets protect the integrands from developing collinear singularities.Furthermore, since the expressions in square brackets start contributing at O(ϵ), the N -jettiness dependence that arises from factors ψ 4ϵ m , only appears in O(ϵ 0 ) contributions to the soft function.
Next, we discuss the B (ij) term in Eq. (10.29).In this case terms with W m ij provide integrable contributions.Taking this into account and repeating the steps we employed to rewrite A ij , we find (10.32) It follows from Eq. (10.29) that B ij is needed through O(ϵ) to compute the relevant contributions to the soft function.Therefore, for the analysis of the divergent contributions to the soft function, in Eq. ( 10.32) we require the first term without g (4) and the last term at ϵ = 0.

The triple-collinear subtraction term with N -jettiness constraint
The second term we require is the triple-collinear limit of the eikonal function subject to the N -jettiness constraint.The corresponding equation reads There are four contributions to the above equation -the i-and j-triple-collinear partitions and the m||n subtractions terms applied to the triple-collinear limit to remove subdivergences.The calculation of all these terms is similar, so that we take the i-triplecollinear partition and discuss first the triple-collinear limit.
For further reference, we define where we have used It follows from this equation that in the triple-collinear limits the dependence on the Njettiness function disappears since we only need its limiting expression.
To proceed, we require the triple-collinear limit of the soft-subtracted eikonal function.It is straightforward to obtain it; we find . (10.36) We then write where To simplify the calculation of this quantity we use the (partial) symmetry of the integrand under ω → 1/ω transformation as well as the invariance of the integration measure under m ↔ n, to extend the integration over ω in Eq. (10.38) to ω = ∞ and to simplify the integrand.We find It turns out that it is relatively straightforward to compute B (i,tc) gg .To do that, it is convenient to employ a Mellin-Barnes representation [49] for 1/(ρ im + ωρ in ) 1−4ϵ and then integrate over ω and the directions of m and n.We note that this requires the introduction of a second Mellin-Barnes representation for integrals that involve ρ mn .This is typically done using a particular Mellin-Barnes representation for the hypergeometric function (10.40) The two-dimensional Mellin-Barnes integrals that one obtains are relatively simple as one integration can always be performed using the second Barnes lemma.The remaining integration is either performed using the first Barnes lemma or by mapping appropriate products of Gamma-functions onto parametric integrals using Euler's representation of a Beta function [49].To perform the parametric integrals we use Hyperint [50].We find We continue with the discussion of the double-collinear subtraction of the triplecollinear contribution.Focusing on the i-collinear partition, the corresponding expression reads .42)We now discuss how to compute the required limits.Our starting point is Eq. ( 10.36) where we are supposed to take m||n limit; this limit is made specific through the parametrization of scalar products in b and d sectors.The challenge in taking the limit is related to the fact that 1/ρ 2 mn appears in the first term in Eq. ( 10.36) which is too strong a singularity and so an expansion around the limit is needed.Fortunately, this can be avoided if we recall that in b and d sectors the following equation holds (c.f.Refs.[1,38]) Here N is the function that appears in the phase-space parametrization of b and d sectors [1,38].For the purposes of computing the double-collinear limit of a triple-collinear limit, we need to take x 4 → 0 in Eq. (10.43).Upon doing that, we find [1,38] lim If we combined the above expression together with Eq. ( 10.14), we easily obtain (10.45) Note that we have used dΛ λ = (1 + 2ϵ)/2, see Refs.[1,38] for details.We also note that the contribution of the sector d is identical.The rest of the calculation is straightforward.Integrating over directions of m and over the energy parameter ω, and combining contributions of sectors b and d we find is defined in Eq. (10.15).
Combining the different results, we obtain As expected, we observe the cancellation of the leading O(1/ϵ 3 ) singularities between the triple-collinear contribution and its double-collinear limit, indicating that the subtraction of m||n subdivergences was properly performed.

Correlated emissions without N -jettiness constraint
As we mentioned earlier, we would like to use the result of Ref. [48] to deduce certain contributions to the N -jettiness soft function.The corresponding expression is given in Eq. (10.10); our goal is to determine the sum of the fourth and the fifth terms in the righthand side of that equation from the known result for J ij .It is clear that this is only possible if all other terms on the right-hand side of Eq. ( 10.10) are known.Below we describe the calculation of the strongly-ordered, double-collinear m||n and the triple-collinear limits of the correlated emissions without the N -jettiness constraint.

Strongly-ordered limit without N -jettiness constraint
We begin by considering the following integral Using the result for the soft limit of the eikonal function Sgg ij , we find The two terms that need to be integrated in Eq. ( 11.2) are quite different.The second one is easy to compute because integrations over directions of m and n factorize.For each of the integrals we use where Integration of the first term on the right-hand side in Eq. ( 11.2) requires more effort.We first use the result in Eq. (B.5) to integrate over n, then write the obtained hypergeometric function using the representation in Eq. (10.40) and integrate over directions of m.This gives another hypergeometric function for which we introduce Mellin-Barnes representation again.Finally, we find (11.5) The integration over z 1 can be performed with the help of the Barnes' second lemma.Once this is done, the resulting integral can be transformed to a parametric integral using Euler's integral representation of Beta-functions.This gives where η ij = ρ ij /2.The resulting integral over x is straightforward to compute using the program Hyperint [50].
Performing the integration, combining the results for the two integrals that appear in Eq. ( 11.2) and expanding in powers of ϵ, we find In addition to standard polylogarithms, also the Nielsen polylogarithm S 2,2 (x) appears in the above formula.

Double-collinear limit without N -jettiness constraint
The easiest contribution to compute is the double-collinear m||n subtraction term since we can borrow the calculation from Section 10.1 nearly verbatim.We then find The various quantities in this expression differ from similar quantities when the N -jettiness constraint is imposed.In particular, and Furthermore, the anomalous dimensions read in this case (11.11)

Triple-collinear limit without N -jettiness constraint
The next contribution that we require is the triple-collinear one, from which the doublecollinear limit is to be subtracted.Similar to the calculation with the N -jettiness constraint, we find it to be convenient to calculate the two terms separately and then subtract them from each other.Furthermore, the two triple collinear configurations, m||n||i and m||n||j, are the same.Hence, we consider one of them and multiply the result by a factor two to account for the other one.
An important difference with respect to a similar computation with the N -jettiness constraint is the fact that the ω → 1/ω transformation is not a symmetry anymore.For this reason the closed-form integration over ω becomes impossible.Hence, we proceed as follows.We first employ the Mellin-Barnes representations to integrate over angles of m and n.Similar to the N -jettiness case, we find that one of the Mellin-Barnes integrations can be performed using the second Barnes lemma.We then map the remaining Mellin-Barnes integrals onto multi-variable integrals using Euler's integral representation of Betafunctions and then integrate over the new (Euler) variables as well as over ω.We find (11.12) The next step is to compute the double-collinear limit of the triple-collinear subtraction term without the N -jettiness constraint.The calculation follows the same steps as what is discussed in the previous Section where the N -jettiness constraint was taken into account.We find the following representation of the corresponding quantity The remaining integration over the directions of parton m is straightforward.We find (11.14)

Extracting the required combinations of the various contributions
As we already mentioned, J ij defined in Eq. (10.10) was computed in Ref. [48].We can use this result, supplemented with our computation of the various limits without the N -jettiness constraint, to derive analytic results for the quantity that is needed for computations with the N -jettiness constraint.
Specifically, we use Eq.(10.10) to write Using the results of the calculation discussed earlier, we obtain where To save space, we do not show O(ϵ 0 ) terms in Eq. (11.16) although they are needed for the final result.It is, however, straightforward to compute them following the preceding discussion.
We can use the result shown in Eq. (11.16) to compute the remaining contribution to the N -jettiness soft functions that originates from correlated emissions.Upon doing that, we find that the divergent contributions cancel and a finite remainder is obtained.We present the corresponding results in Section 13.

Final-state fermions
We continue with the discussion of the contribution of the soft q q pair to the N -jettiness soft function.In this case, the entire contribution comes from correlated emissions.It reads where The eikonal function Sqq ij is symmetric under the interchange of q and q momenta; it can be found in Appendix C. Hence, we proceed in exactly the same way as with the gluons in that we introduce the energy ordering E n < E m , write E n = ωE m , with 0 < ω < 1, and integrate over E m .After the Laplace transform, we obtain where we have to set E m to one and E n to ω when computing Sqq ij (m, n) in the last equation.We continue with the investigation of the various limits of the above equation.In principle, these limits are identical to the ones in the gluon case except that the q q eikonal function does not possess a strongly-ordered singular limit.Accounting for this and following what has been done for the two-gluon final state, we write We continue with the discussion of the individual terms.The first one is the double-collinear m||n limit.It reads The vectors v µ and κ ν n are defined in Eq. (10.16).We compute this contribution following the discussion of the two-gluon case.Borrowing most of the calculation from there, we write the result for the sum of b and d sectors with and Structurally, this result is identical to the m||n contribution in the two-gluon case except that different collinear anomalous dimensions appear.The analysis, therefore, proceeds along identical lines and the result can be written in the same way as in the two-gluon case, c.f. Eqs.(10.29, 10.31, 10.32) up to the change in the anomalous dimensions and the overall normalization factor.
As the next step, we need to compute the triple-collinear limit of the q q eikonal function.The calculation follows what has already been done for the two-gluon case.We obtain Since this expression is also very similar to the two-gluon case, we proceed accordingly.As the first step, we change the integration variable ω → 1/ω and extend the integration over ω to ω = ∞ using the symmetry of the integrand.We then write Σ i,tc q q = N u ϵ B i q q, (12.10) where Using the m ↔ n and ω → 1/ω transformations, we rewrite the above expression as 12) The calculation proceeds along the same lines as in the two-gluon case.We obtain The double-collinear subtraction term of the triple-collinear contribution reads Following the very same procedure as in the two-gluon case, we arrive at Combining the above results and accounting for the contribution of the second partition function, we obtain the result for the triple-collinear limit of the q q contribution to the integrated eikonal function subject to N -jettiness constraint To determine the collinear-regulated contribution in Eq. (12.4), we proceed in the same way as in the two-gluon case.Namely, we write a similar representation for the "energyregulated" integral of the q q contribution to the soft function computed in Ref. [48].It reads (12.17) We observe that, up to a normalization factor, the last two contributions are the same as the missing divergent ones in Eq. (12.4).Since the result for J q q ij is known [48], we can extract the contribution of the last two terms in Eq. (12.17) if the first two terms are known.We note that these terms are similar but not identical to the ones that we already computed because the N -jettiness function ψ mn does not appear in Eq. (12.17).However, their computation is nearly identical to what has been discussed in the two-gluon case and for this reason we do not repeat it.We find To obtain the above formula, we used the following result for the triple-collinear contribution in the case without the N -jettiness constraint.13 The renormalized N -jettiness soft function We are now in a position to present the renormalized N -jettiness soft function through NNLO in the perturbative expansion in QCD.We write It is convenient to introduce the following short-hand notations The NLO contribution reads The NNLO contribution is constructed from different pieces.We write The function G ij reads In the above formula δ ij = θ ij /2 and θ ij is the relative angle between partons i and j.The Clausen functions are defined through their relation to complex-valued polylogarithms Li n (e iz ) − Li n (e −iz ) . (13.6) The functions S a 1 ,a 2 and G a 1 ,a 2 ,...,am (x) are the Nielsen and Goncharov [51] polylogarithms, respectively.The auxiliary functions that we employed while writing the above formula read (13.7) The function Q ij evaluates to The function G triple kij , defined in Eq. (7.14), corresponds to the finite remainder of the triplecolor correlation terms.

Numerical implementations and checks
We have implemented the above formulas into a Fortran code.The implementation is, in principle, straightforward.To compute terms that require integration over directions of gluons' momenta, we use the phase space parametrization described in Ref. [1,9,10].A possible choice of partition functions can be found in Ref. [52].
As we already mentioned, the N -jettiness soft function was recently discussed in Ref. [37] where a plethora of numerical results for various kinematic configurations have been presented.We have checked many of the presented results for various values of the parameter N and found excellent agreement.We present some examples of such a comparison in Appendix D.
We note that to obtain the high-precision numbers shown in Appendix D, we have used a very large number of sample points which results in runtimes that can be as large as a few minutes per dipole per kinematic point.However, to obtain the same numbers with a better-than-percent precision, we need a relatively small number of sampling points which results in runtimes of an order of a few seconds per dipole per phase space point.We note that runtimes for individual dipoles are largely independent of N .Since (N + 2)(N + 1)/2 dipoles need to be calculated to determine the N -jettiness soft function for a given phase space point, O(1 − 2) minutes will be required to do that in case of three-or four-jet production at a hadron collider, or six-jet production at an e + e − collider.
The situation with triple-color correlated terms is quite similar.Since in this case integration over the direction of a single gluon is required, the integration converges rather fast.In general, for an N -jettiness case, the evaluation of (N + 2)(N + 1)N independent functions G triple kij is required from which a smaller number of independent color-correlated contributions can be constructed.We need about twenty seconds to compute all the required G triple kij functions for N = 2 and, therefore, we will need about a minute to compute them for N = 3 and about two minutes for N = 4.

Conclusions
We have described the computation of the N -jettiness soft function through NNLO in QCD.Keeping N as the parameter, we demonstrated the cancellation of all 1/ϵ poles analytically against the soft-function renormalization matrix.Furthermore, we derived a simple representation for the finite, jettiness-dependent remainder valid for an arbitrary number of hard partons N .We compared our numerical results for N = 1, 2 and N = 3 with the ones recently presented in Ref. [37] and found excellent agreement.We have also found that the representation of the finite remainder that is derived in this paper leads to fast and rapidly convergent integration which is important if the results for N -jettiness soft function are to be used for the computation of higher-order corrections for multi-jet production at colliders.
From a more general perspective, this study was motivated by a question of how advances in developing NNLO subtraction schemes for generic processes at the LHC can be used to derive representations for the building blocks of modern slicing methods that rely on choosing observables to define suitable slicing variables.We have shown that, at least for N -jettiness, the benefits of applying subtraction-inspired methods are significant.We note that the application of the subtraction-inspired methods for computing the N -jettiness soft function becomes possible if one departs from the (by now) standard approach [16,18,34] to such computations where hemisphere soft functions are used as elementary building blocks and, instead, interprets N -jettiness as one of the many infrared safe observables that can be studied using available subtraction schemes.

C Eikonal functions
For completeness, we provide the double real-emission eikonal functions that we use to compute the relevant contributions to the N -jettiness soft function.Following Ref. [46], for two soft gluons, we define The function S gg ij (m, n) reads (C.3) In case the soft partons are a quark and an anti-quark, we write Sqq ij (m, n) = 2S q q ij (m, n) − S q q ii (m, n) − S q q jj (m, n), (C.4) where [46] S q

D Numerical Checks of N-Jettiness Soft Functions
In this section, we compare results for the N -Jettiness soft function obtained in this paper with those in Ref. [37].For this comparison, we define two functions G nl ij and Q nl ij that are obtained from G ij and Q ij by setting all terms with L ij to zero We also note that for the comparison of the triple-color correlated contributions, we do not set these logarithms to zero but take ū = µ = 1 instead.As the result, we should set L ki → 1/2 ln η ki in Eq. (7.14).

D.1 1-Jettiness
In order to compare our numerical results, we parameterize our phase space in a similar way as done in Ref. [37].In the case where there are two back-to-back beams and one jet confined in a plane, we parameterize the scattering by the angle θ 13 using the following momenta, n 1 = (0, 0, 1), n 2 = (0, 0, −1), n 3 = (sin θ 13 , 0, cos θ 13 ).(D.3) In order to compare a phase space point where the jet is separated from the beams, we take θ 13 = 12π 25 . (D.4) In Table 1 we present the results for the non-logarithmic coefficient of the renormalized 1-jettiness soft function, which are in agreement with the ones in Ref. [37].

D.2 2-jettiness
We again consider the configuration where there are two back-to-back beams.We parameterize the scattering of the two jets by θ 13 , θ 14 and ϕ 4 as n 1 = (0, 0, 1), n 2 = (0, 0, −1), n 3 = (sin θ 13 , 0, cos θ 13 ), n 4 = (sin θ 14 cos ϕ 4 , sin θ 14 sin ϕ 4 , cos θ 14 ). (D.5) Considering a generic phase space point where the beams and jets are separated from each other, we take Results for all possible dipole coefficients are presented on Table 2.In the case of 2-jettiness, the tripole contribution can be reduced to a single color structure, f ABC T A 1 T B 2 T C 3 .The renormalized sum of all tripole contributions in this configuration is 1064.77± 0.08, which is in agreement with the value (1064.6 ± 0.1) reported in [37].

D.3 3-jettiness
We parameterize the scattering of the additional jet by θ 15 and ϕ 5 n 5 = (sin θ 15 cos ϕ 5 , sin θ 15 sin ϕ 5 , cos θ 15 ). (D.7) We compare in Table 3 our 3-jettiness dipole contributions to the benchmark result in [37] by taking the following phase space point For the tripole contribution in this configuration, there are four independent color structures.In Table 4 we present our results for each of these configurations, as defined in [37].