The NNLO QCD soft function for 1-jettiness

We calculate the soft function for the global event variable 1-jettiness at next-to-next-to-leading order (NNLO) in QCD. We focus specifically on the non-Abelian contribution, which, unlike the Abelian part, is not determined by the next-to-leading order result. The calculation uses the known general forms for the emission of one and two soft partons and is performed using a sector-decomposition method that is spelled out in detail. Results are presented in the form of numerical fits to the 1-jettiness soft function for LHC kinematics (as a function of the angle between the incoming beams and the final-state jet) and for generic kinematics (as a function of three independent angles). These fits represent one of the needed ingredients for NNLO calculations that use the N-jettiness event variable to handle infrared singularities.


Introduction
The continued successful operation of the Large Hadron Collider (LHC) has led to the accumulation of a very large data set with which to study the Standard Model (SM) in unprecedented detail. The ever-increasing precision of the experimental analyses has mandated a similar increase in the precision of the corresponding theoretical predictions. Over the last few years a concerted effort has been made in the theoretical community to provide predictions accurate to next-to-next-to-leading order (NNLO) in QCD. Calculations that include colored final-state radiation are particularly challenging. Significant progress in this direction has been made recently, including the NNLO calculation of the production of dijets [1,2], V + j [3][4][5][6][7][8][9], H + j [10][11][12], H + 2 j in a e-mail: johnmc@fnal.gov b e-mail: keith.ellis@durham.ac.uk c e-mail: rmondini@buffalo.edu d e-mail: ciaranwi@buffalo.edu the vector boson fusion process [13], single top [14][15][16] and tt [17].
A critical element in the completion of a NNLO calculation is a manageable way to handle the copious infrared (IR) singularities present in component pieces of the calculation. These singularities occur in phase spaces of differing dimensionality, and cancel only when combined in a suitably-inclusive, IR-safe observable. At NLO the most widespread solutions use local subtraction terms [18][19][20]. In these approaches one subtracts a user-defined set of counterterms from a given real-emission matrix element such that the corresponding combination of matrix element plus counterterms is finite in all singly-unresolved IR limits. The counterterms are constructed in such a way as to be integrable analytically over a single unresolved parton; the integrated counterterms can then be combined with the virtual one-loop matrix elements, resulting in an analytic cancellation of IR poles. The two phase spaces are then both manifestly finite and can be integrated separately using Monte Carlo integration techniques.
The construction of a similar subtraction scheme at NNLO accuracy is a considerably more daunting task. This is primarily due to the presence of one extra unresolved parton with respect to NLO, yielding multiple overlapping singularities. Despite its difficulty, a number of subtraction schemes have been developed and successfully applied to several LHC processes [21][22][23][24][25].
Alternatives to local subtraction schemes are possible. One such method, based on a more global approach, is phasespace slicing [26]. These methods are simple to implement at NNLO, particularly if the corresponding NLO process with one extra parton in the final state is already known. In slicing methods a global parameter is used to divide the phase space into (at least) two regions. At NNLO the two regions correspond to the region which includes all of the doublyunresolved emissions, and a region which has at most one singly-unresolved parton. The IR structure in the latter region is clearly akin to that obtained in a standard NLO calculation, and hence this region is amenable to calculation using existing NLO technology. The success of the method therefore relies on the ability to calculate the region which contains doubly-unresolved partons. To this end, factorization theorems are used to calculate the cross section systematically in this region. The first slicing method applied at NNLO [27] used the transverse momentum of the final state, q T , as a slicing parameter, and the Collins-Soper-Sterman factorization theorem [28] to compute the cross section in the region of small q T . As such, this method is applicable to final states in which there is no q T associated with colored radiation, i.e. the production of color-singlet final states. A recentlydeveloped method [29,30] uses the N -jettiness event-shape variable T N [31] rather than q T , and a factorization theorem from Soft-Collinear Effective Field theory (SCET) [32][33][34][35][36]. The theorem states that the cross section in the region of small T N can be obtained from the following convolution [31,37] Here B represents the beam function, which describes initialstate collinear radiation, and J the jet function, which describes final-state collinear radiation. For these two functions, expansions accurate to O(α 2 s ) can be found in Refs. [38,39] and Refs. [40,41] respectively. The term H denotes the hard function, which is process-specific and finite. Finally, S represents the soft function, which is the main focus of this paper. The soft function is defined as the process-independent soft limit of QCD amplitudes. Of the process-independent pieces of the factorization theorem, it is by far the most complicated. For color-singlet production processes (which for the LHC correspond to zero jets in the final state), the soft function is reasonably simple and analytic expressions are known [42][43][44]. When three colored partons are present (1-jet final states for LHC kinematics), the soft function is considerably more intricate. The calculation of the 1-jettiness soft function at NLO was presented in Ref. [45]. In Ref. [46], the NNLO 0-jettiness soft function for one massive colored particle production was computed.
A method to compute the 1-jettiness soft function numerically to NNLO accuracy was presented in Ref.
[42] and the results of this calculation have since been used to compute several V + j processes at NNLO [6,7]. However, at present, there is no publicly-available computation of the 1-jettiness soft function presented in a form that can be implemented in an independent Monte Carlo code. The main focus of Ref.
[42] is to provide a methodology of computing the soft function numerically. Specific results are only presented for the qg → q configuration, and only in graphical form. Since the gg → g and qq → g configurations are absent, and due to the nature of the result presented, it is currently impossible to implement the 1-jettiness soft function at NNLO directly from the literature. The primary aim of our paper is to provide this information via an independent calculation of the 1-jettiness soft function.
Our paper proceeds as follows. In Sect. 2 we provide a general overview of the calculation. We introduce the Njettiness variable T N and present our parametrization of the phase space. In Sect. 3 we present the formulae for soft parton emission at NLO and NNLO taken from the literature, paying particular attention to their color structures. In Sect. 4 we validate the method by computing the NLO soft function for the 1-jettiness case and compare the results against known analytic formulae. In Sect. 5 we discuss our calculation for the NNLO 0-and 1-jettiness soft functions in detail. In Sect. 6 we present the obtained results and compare them against the known results in the literature. We draw our conclusions in Sect. 7.

N -jettiness variable T N
For a parton scattering event the N -jettiness variable T N [31] is defined as In units whereh = c = 1, T N therefore has the units of mass.

Sudakov decomposition
A convenient way of parametrizing the momenta appearing in the phase-space integrals that enter the calculation of the soft function is to use a Sudakov decomposition of the momenta in terms of two of the momentap i that appear in Eq. (3). We first define a shorthand notation for the quantity that appears in Eq. (3), namely the projection of a vector q along the direction ofp i , The parton momentum q μ can then be expanded as where y i j = 2p i ·p j and q i j⊥ is transverse to the plane spanned byp i andp j . The Sudakov expansion forp k , which is not one of the Sudakov base vectors, iŝ We can calculate q k , the projection of q on a non-Sudakov base vector, and obtain where φ qk is the angle in the transverse plane between q and p k , and The ratio of the projection along a non-Sudakov direction k to the projection along a Sudakov direction i or j is given by, where x i j = q i /q j . For the case of hadronic collisions where two of the directions,p 1 andp 2 , are those of the beams we have, so that y 12 = 1. We note that this notation follows that of the calculation of the NLO soft function [45]. It differs from that of Ref.
[42], in which relevant quantities are expressed in terms of pure directions, n i . Equivalent expressions can be obtained by making the replacement,

Measurement function
Written in terms of the projected momenta, the definition of the N -jettiness is The minimum over i in this equation provides a natural division of the phase space for extra emission into regions where each projection is smallest. For the case of 1-jettiness we will label the three hard directions as i, j and k. The singleemission phase space is then partitioned by inserting a measurement function F where, and F i corresponds to q 1 being closest to direction i so that, for instance, We note that in the original calculation of the NLO soft function [45], a further hemisphere decomposition of the measurement function was used to separate the divergent and finite parts of the calculation. However, we will not pursue that method for our NNLO calculation. We can extend the decomposition of F to the double-emission case by writing the measurement function as, where the sum runs over the nine combinations of a, b ∈ {i, j, k}. The notation F ab implies that q 1 is closest to direction a and q 2 is closest to b with, for example, As it will be explicitly shown in Sect. 3, the integrals that we have to evaluate in order to calculate the soft function have an eikonal form in which there are two emitting directions i and j. We attach the labels of the emitters as superscripts to the measurement functions. In this notation we can write out the decomposition in Eq. (16) explicitly as Each term corresponds to one of four cases [42], as indicated: 1. Both q 1 and q 2 closest to the same emitting direction; 2. q 1 and q 2 closest to different emitting directions; 3. One of q 1 and q 2 closest to an emitter, while the other is closest to the non-emitting direction k; 4. Both q 1 and q 2 closest to the non-emitting direction k.
In general, there is one more case where q 1 and q 2 are closest to different non-emitting directions k, l, but this is absent in the 1-jettiness case.

Phase space -single emission
Let us start by considering the phase space for a single emission: Using the Sudakov variables q i and q j defined in Sect. 2.2 above, we can rewrite the integration measure as The phase space then becomes where dΩ (d−2) is the (d − 2)-dimensional angular measure. Setting d = 4 − 2 and using the standard expression for the angular measure after integrating over unconstrained angles given in Eq. (A.12), we obtain P S (1) It is convenient to normalize the remaining angular integration so that it integrates to one using, The final expression for the single-emission phase space is then,

Phase space -double emission
For the double-emission phase space we employ the same Sudakov decomposition as the single-emission case and, following the same steps as above, we have P S (2) The integral over the transverse space for q 1 can be performed just as in the single-emission case, with the result given in Eq. (A.12). The integral over the transverse space for q 2 is more complicated since one of the angles cannot be integrated out; the form of the integral is given in Eq. (A.10). Combining these expressions we arrive at the final form for the phase space, where 3 Soft function and soft radiation

Soft function
We define the unrenormalized N -jettiness soft functioñ S(T N ) as a perturbative series in powers of the bare strong coupling α s : We renormalize the coupling constant by performing the replacement The renormalization factor is given by, (30) with C A = 3, T R = 1 2 , N F = 5, and α s ≡ α s (μ) at the renormalization scale μ. The coefficients of the perturbation series of the renormalized soft function S(T N ) in terms of the unrenormalized ones then read: The leading-order contribution is simply S (0) (T N ) = δ(T N ), since at leading order there is no emitted radiation.

Soft radiation at NLO
We start by computing the soft function at NLO, which allows us to illustrate the main features of the method as well. The result for the NLO soft function is known analytically [45] and can be used to validate our numerical evaluation. The NLO corrections to the leading-order soft function are made up of two different contributions: Born-type processes with one-loop corrections ("virtual" corrections) and tree-level processes with the emission of one additional parton ("real" corrections). The former only contribute at T N = 0. Since we are considering corrections on massless eikonal lines, there is no dimensionful quantity to carry the dimension of the one-loop integrals; their contribution therefore vanishes in dimensional regularization. The only contribution at NLO is therefore the one from real radiation. The form of the squared amplitude representing the emission of a single soft gluon is well known. Using the same notation that will be employed at NNLO we can write the factorization at O(g 2 ) as, c.f. Eq. (12) of Ref. [47]. The eikonal function S i j (q) is given by Note that here we have introduced the normal M S factor S , where γ E is the Euler-Mascheroni constant. The emission of a soft gluon produces color correlations that are indicated in Eq. (32) by the subscripts i and j in M Factoring out the leading-order amplitude squared we thus have a simple expression for the soft-gluon approximation, where we have limited the sum that appears in Eq. (32) to i < j and added the consequent factor of two. Introducing the measurement function of Eq. (14), the total result for the soft function at NLO is then, 123 where, in the last line, we have explicitly indicated the choice of momenta in the Sudakov decomposition in the phasespace. The explicit expressions for the quantities in curly braces in Eq. (37) are given in Eqs. (15) and (24)

Color conservation
The eikonal expressions given above are written using colorspace notation [20]. Using color conservation we have, Thus for the case of 0-jettiness ( j ∈ {1, 2}) we find that 3 and cyclic permutations. For the 1-jettiness case we can write, All products T i · T j can therefore be expressed in terms of sums of Casimirs T 2 i (and vice versa) with

Soft radiation at NNLO: real-virtual
As earlier discussed, the virtual diagrams do not contribute because of scaling arguments, so that at NNLO we are left with the real-virtual contribution and the double-real contribution (to be considered in the following two sections). Since the real-virtual corrections involve only one real emission, the calculation follows closely the NLO case. The one-loop contribution to the soft-gluon current has been given in Ref. [47]. Combining Eqs. (23) and (26) of Ref. [47], the O(g 4 ) real-virtual contribution is, where S i j (q) is given in Eq. (33). The notation represents a sum over values of the indices that are distinct (for instance for the second term this explicitly means i = j, j = k, k = i). The calculation of 1-jettiness does not permit such a contribution from the second term and therefore we do not consider it further.
The real-virtual contribution is thus, where we have used and extracted an overall factor We have also identified a factor (shown in square brackets in the third line of Eq. (42)) that will be naturally cancelled by a corresponding one in the phase space, c.f. Eq. (24). Finally, we note that the method for calculating this contribution will be very similar to the one used for the NLO soft function, after the replacement of the integrand factor y i j /q i q j by [y i j /q i q j ] 1+ . For the real-virtual contribution to the NNLO soft function we therefore have 3.5 Soft radiation at NNLO: qq emission Turning now to the double-real emission, we first consider the radiation of a soft qq pair. The factorization of QCD amplitudes in the double-soft limit has been given in Refs. [48,49]. From Eq. (95) of Ref. [49] we have, Using the color identities of Eq. (39), we can rewrite the result for the matrix element as, where the quark eikonal result always appears in the special combination [42] From Eq. (96) of Ref. [49] the soft quark pair production result for I i j is Inserting this result into Eq. (48), the result for quark pair emission can be written in terms of two functions where Using Eq. (47), we see that the contribution to the soft function due to the emission of N F flavors of light quarks is given by, where F i j is the measurement function, Eq. (18).

Soft radiation at NNLO: gg emission
The case of two-gluon emission gives rise to both Abelian and non-Abelian contributions. The Abelian two-gluon matrix element squared is given by the product of two single-gluon currents weighted by a factor of 1 2 . The integrations over the two emitted momenta factorize, so that the Abelian twogluon emission result is determined by the NLO result, and we will not consider it further.
The result for non-Abelian soft radiation has been given in Eq. (108) of Ref. [49] and is proportional to The two-gluon soft function is given in Eq. (109) of Ref. [49], Using color conservation, Eq. (39), it is clear that we only require the combination, This result can further be decomposed as [42], where J I i j , J I I i j are given in Eqs. (51) and (52) and 123 Thus the final result for this contribution to the soft function is where S f = 1/2 is the statistical factor for the two final-state gluons, and F i j is the measurement function, Eq. (18).

Soft function at NLO
In this section we will make the contributions shown in Eq. (37) explicit and then assemble the complete NLO soft function.
To evaluate this contribution we first perform the change of variables, The measurement function defined in Eq. (15) then becomes, Parameterizing the phase space given in Eq. (24) using the same set of variables, Eq. (61), and combining the two gives, The matrix element from Eq. (36) can be written as, This contribution to the NLO soft function is then, For this contribution we use the change of variables, and note that the projection that enters the matrix element (q j ) can be related to these through The combination of measurement function and phase space is trivially related to the expression in Eq. (63) by cyclic permutation of the labels i, j and k. The matrix element is, which yields the expression for this contribution,

0-jettiness
The NLO soft function for 0-jettiness is straightforward to compute analytically and we need not resort to the numerical methods that we will employ throughout the rest of this paper. A short description of this calculation is given in Appendix E for completeness.

1-jettiness
In the 1-jettiness case, the calculation of the soft function at NLO is more involved since one of the θ -functions depends on an angle and therefore we have to resort to numerical integrations. In the 1-jettiness case we can have three different leading-order configurations, gg → g, qq → g, and qg → q, where the configuration determines the color factors that appear in Eq. (36) through the relations presented in Sect. 3.3. For the LHC kinematics we definep 1 andp 2 as the directions of the initial-state partons, as in Eq. (11), and p 3 as the direction of the final-state parton. The soft function is given by Eq. (37), where each of the contributions is evaluated using either Eqs. (65) or (68). The integrals are evaluated using the sector decomposition approach [23,50,51]. The term x −1+ i j in Eq. (65) is expanded in terms of delta and plus distributions by means of the relation where, given a sufficiently-smooth function f (x), the plus distributions are defined as After using the expansion and expanding the integrals as Laurent series in , the coefficients of the series are obtained by numerical integration. This integration is straightforward and has been performed with a Fortran code using double precision accuracy without, in this case, any additional cuts for numerical safety. 1 The final expression for the soft function is obtained by expanding the overall factor T −1−2 1 in terms of delta and plus distributions, using Eq. (71) starts at order −1 , we note that it is necessary to compute all integrals up to O( ). For our purposes, namely the factorized cross section defined in Eq. (1), we require only the O( 0 ) part of the renormalized soft function. The result is then given as the coefficients of δ(T 1 ) and L n (T 1 ) with n = 0, 1: We now present the numerical results for the 1-jettiness NLO soft function. We specialize our calculation for LHC kinematics, where we have y 12 = 1 and momentum conservation gives y 23 = 1 − y 13 so that the result is a function of y 13 alone. We compute the soft function for the three channels (gg → g, qq → g, and qg → q) with 15 different values of For each channel we compare the coefficients of δ(T 1 ), L 0 (T 1 ), and L 1 (T 1 ) (denoted by C i with i = −1, 0, 1 in our notation) with the known analytic results in the literature [45]. We find excellent agreement with the known results for all channels and all coefficients, as demonstrated in Fig. 1. The agreement is at the level of 0.2% or better. In Fig. 2 we plot the O( 2 ) coefficient of the series expansion in of the soft function before the term T −1−2 1 has been written out using Eq. (71). The O( 2 ) term does not contribute at NLO once T −1−2 1 is expanded, but instead enters the coefficient of δ(T 1 ) in the renormalization contribution to the NNLO soft function, as shown in Eq. (31) and explained in detail in the next section. Since this contribution is not available in the literature, we perform a fit of our numerical results and present the fit here for completeness: The coefficients k (m,n) are collected in Table 1. For the channels gg → g and qq → g we set k (m,n) = k (n,m) (with m = n) since these channels are symmetrical in the initial state.

Soft function at NNLO
From Eq. (31) the renormalized N -jettiness soft function at NNLO, S (2) We have already discussed the calculation of the first term in Eq. (75) in the previous section, so we now study the second term in more detail. At NNLO three types of corrections arise: two-loop virtual corrections, one-loop corrections with one real gluon emission ("real-virtual"), and double-real corrections, consisting of the emission of qq and gg pairs. In the latter case, the emission of two gluons is made up of two diagrammatic contributions: an Abelian contribution, i.e. two single-gluon currents, and a non-Abelian contribution, i.e. emission of a gg pair through a three-gluon vertex. In light of this, and since the virtual two-loop corrections vanish in dimensional regularization, we can rewrite the term S (2) (T N ) as the sum of four non-vanishing contributions: withS (2) gg (T N ) indicating the non-Abelian part of the doublereal gg corrections. The Abelian termS (2) ab (T N ) can be obtained directly from the NLO soft function thanks to wellknown exponentiation theorems [52,53]. The calculation of the real-virtual contributionS (2) RV (T N ) is conceptually identical to the calculation of the NLO soft function and has been 123 Fig. 1 Numerical results for the 1-jettiness soft function at NLO for the partonic channels gg → g, qq → g, and qg → q. We plot the coefficients of L 1 (T 1 ), L 0 (T 1 ), and δ(T 1 ) as functions of y 13 ∈ [0, 1]. The known analytic results are taken from Ref. [45] described in Sect. 3.4. We therefore spend the remainder of this section discussing the calculation of the double-real corrections represented byS (2) gg (T N ) andS (2) qq (T N ). The decomposition of the relevant eikonal matrix elements into a set of basis integrals is given in Sects. 3.5 and 3.6 for the qq and gg cases respectively. We will detail the transformations necessary to render the basis integrals, obtained by combining the double-emission phase space with the eikonal integrands, amenable to numerical evaluation. We will provide a number of illustrative examples of integrals, leaving a complete enumeration of all contributions to Appendix C.

Double-real corrections
We first identify a common overall factor that is associated with the total angular volume, in the phase-space parametrization of Eq. (26), where and S is given in Eq. (34). We will always be able to trivially rescale one of the angular integrals so that the range of integration is between zero and one, with φ = π x φ , and the integral over β can be recast in similar fashion, by performing the change of variable cos β = 1 − 2x β . We thus rewrite Eq. (26) as, where we have written the remaining φ 2 integral in a form that anticipates the Jacobian factor that will eventually emerge as in Eq. (79). Finally, in our eventual evaluation of the x β integral we will replace this form with a simpler one in which singularities can only appear at x β = 0, 2 As before, we compute the integrals numerically using sector decomposition and the plus-distribution expansion given in Eq. (69). As in the NLO case, the choice of the phase-space parametrization will depend upon which case is picked out by the measurement function. We will discuss each of the cases described in Sect. 2.3 in turn. The combination of the matrix elements, phase space, and measurement function for each case will be expressed as, where the coupling-associated factor F is defined in Eq. (C.25). In this equation X labels the division of the eikonal approximations, Eqs. (51), (52), and (58), and a, b the different cases. We have extracted an overall factor that will be universal across all contributions, (84)

Case 1: F i j ii
The appropriate change of variables for this case is, We then have, An alternative procedure [54] is to split the integration range in half, into [0, 1 2 ) and [ 1 so that the combination of the phase space in Eq. (81) with the measurement function F i j ii yields, where φ 1 = π x φ 1 . In the simplest case we will be able to perform a rescaling of φ 2 as in Eq. (79). However, in general this is not true. The reason for the additional complication is the presence of the factor q 1 · q 2 in the eikonal factors T i j and U i j . Although the Sudakov decomposition is very convenient for describing the other dot products, the expression for this one is considerably more complicated, where z 12 = 1 2 (1 − cos φ 12 ) and φ 12 is the angle between q 1 and q 2 . A method for handling this denominator has been outlined in Refs. [42,55] and we follow this strategy here. We map z 12 to a new variable λ through the relation, or the inverse, The Jacobian associated with the transformation from z 12 to λ is, In terms of the new variable λ we have, so that Further, from Eq. (89) we have that, We observe that by means of this change of variable the collinear singularity q 1 · q 2 → 0 has been mapped to s → t.
In the presence of this denominator (and associated factor of z 12 ) we must be careful to ensure that the relevant angles are handled appropriately. It is convenient to perform the integration in a rotated frame, c.f. Appendix B. Following that logic, we replace the measure dβ dφ 2 with dβ 12 dφ 12 and, from Eq. (B.23), relate the angle φ 2 to φ 1 , β 12 , and φ 12 through, cos φ 2 = cos φ 1 cos φ 12 − sin φ 1 sin φ 12 cos β 12 .
The final integral in Eq. (87) then becomes, where λ = sin 2 (π x φ 2 /2). Thus the final form for the phase space is, where φ 1 = π x φ 1 , cos β 12 = 1 − 2x β 12 , φ 12 is defined through Eq. (94), and φ 2 can be obtained from Eq. (95). The matrix elements that must be evaluated to compute the double-real contributions for this case are given in Appendix C.1. They are expressed in terms of the new variables ξ , s, and t and to aid in the evaluation of contribution I I I they have been further subdivided into integrals that are simpler to compute.
To illustrate a further complication that arises for these integrals, consider the evaluation of contribution I . Combining the phase space from Eq. (97) with the matrix element in Eq. (C.26) we obtain, where B R R is given in Eq. (78) and we have removed an overall factor of N i j in the definition of I (2),i j,I ii , c.f. Eq. (83). In order to be able to perform the integration, we have to handle the line singularity associated with the |s − t| −1−2 term. Following Ref.
[42], we do so by partitioning the integral into two contributions by means of the identity, In the s > t sector we then perform the change of variables while in the t > s sector we have such that, in general, all singularities are located at x 2 = 0 and x 3 = 0. The only complication is that, for the integrand J I I I b i j , this procedure also yields singularities at x 3 = 1. However, these are simple to handle by using a simplification similar to the one in Eq. (82), but carried out for x 3 . These transformations are sufficient to treat all of the singularities in this case.

Case 2: F i j i j
The change of variables for this case is, In combination with the measurement function, this results in a phase-space parametrization that is very similar to the simplest one for case 1, where φ 1 = π x φ 1 and φ 2 = π x φ 2 . The remaining invariant that enters the matrix elements (given in Appendix C.2) becomes, where z 12 = 1 2 (1 − cos φ 12 ) and φ 12 is obtained from, c.f. Eq. (B.19), This does not require any further reparametrization of the phase space because factors of 1/q 1 · q 2 do not lead to singularities as s, t → 0.

Case 3: F i j ik
This case uses the following change of variables, based on a Sudakov expansion with respect to i and k, The phase-space parametrization becomes, where φ 1 = π x φ 1 and φ 2 = π x φ 2 . Since one of the emitting lines is no longer one of the vectors in the Sudakov expansion, more invariants must be defined. The remaining quantities are, where z 12 = 1 2 (1 − cos φ 12 ) and φ 12 is again defined through Eq. (B.19). There are no further singularities to disentangle in this case, for which all matrix elements are specified in Appendix C.3.

Case 4: F i j kk
We perform the Sudakov expansion with respect to base vectors i and k and then the following change of variables, After this transformation the phase space parametrization is, The other invariants that enter the matrix elements are, Since this denominator takes the same form as the one considered in case 1, it must be handled in a similar manner. Explicitly we have, where the definitions of all the angles are taken over from case 1. The matrix elements for this case are given in Appendix C.3. In the calculation of the contribution J I i j an additional subtlety arises. Consider the evaluation of the s > t sector that appears after the partitioning of Eq. (99). Following the change of variables in Eq. (100) we have, The penultimate line of this equation appears to indicate that a plus-distribution expansion is required for x 2 and that the x 3 integration is too singular as x 3 → 0 (i.e. s → t in the original variables). However, a careful analysis reveals that neither of these is true. Instead, these additional denominator factors are actually regulated by the numerator, To see that this is the case we write out the expressions for A ki, j explicitly and use the relation between the angles in Eq. (95) to find, This expression makes it clear that the x 2 denominator factor is harmless. Moreover, we observe that in the limit x 3 → 0, so that, The use of this limit is essential in order to obtain the correct subtraction of the singularity at x 3 = 0.

Results
We first make a few observations regarding our numerical integration procedure. Each contributing integral is computed using VEGAS in double precision. Since many of the integrals contain square-root singularities, we routinely perform remappings to remove such factors and improve the convergence of the numerical integration, Moreover, in order to avoid numerical instability when evaluating the double-real integrands extremely close to singularities, we impose a tiny cut on every integration range: We do not observe any sensitivity of our results to reasonable variations of δ, within Monte Carlo uncertainties, and choose δ = 10 −12 for the final results presented below. We have additionally checked that running the code in quadruple precision with a cutoff reduced to δ = 10 −22 does not alter the results. In order to provide the reader with an example of our raw results, and a point of comparison for an independent implementation, we present the numerical value of all double-real integrals at a single phase-space point in Appendix D. Each integral is typically evaluated with an uncertainty that is far smaller than one percent, but which can be at the percent level for a few contributions where the absolute value is very small. This accuracy is more than sufficient to obtain the NNLO soft function at the level necessary for phenomenology.
Having described all of the necessary ingredients to perform the calculations, we can now present our numerical results for the 0-and 1-jettiness soft functions at NNLO. In particular, we focus on the non-Abelian contribution to the soft functions. Following the notation of Sect. 5, we define the non-Abelian part of the NNLO soft function as the sum of four different contributions: Each individual contribution is computed as explained in the previous sections. After performing the sum, the non-Abelian soft function corresponds to the O( 0 ) contribution to the total. This is written as where C n with n = −1, . . . , 3 are numerical coefficients, functions of the invariants y i j .
The gg channel is obtained by simply rescaling by a factor C A /C F . A comparison of our numerical evaluation of the NNLO soft function with the analytic results, for the values of y 12 above, is shown in Fig. 3. Apart from regions where the soft functions are very close to zero, the numerical and analytic results agree perfectly, at the level of a few per-mille or better.

1-jettiness
We present the non-Abelian contribution to the 1-jettiness soft function. As in the NLO case, we specialize our calculation for LHC kinematics. The coefficients C n are functions of y 13 only since y 12 = 1 and y 23 = 1− y 13 . We compute the coefficients C n for the three different partonic configurations (gg → g, qq → g, and qg → q) and for 21 different values of y 13 : The coefficients C 0 -C 3 are known analytically and can be derived, for instance, from renormalization group constraints [30,42,45]. They therefore provide a useful check of our calculation, particularly in the case of C 0 , which receives contributions from all the basis integrals and expansions that must be performed for C −1 . Our results are shown in Fig. 4 and indicate that, as at NLO, the numerical integration is accurate to the per-mille level when compared with the analytic results for C 0 -C 3 . The calculation of the endpoint contribution C −1 is the central result of this paper. Our results for this contribution are shown separately in Fig. 5, for each of the three configurations. We have performed fits to our results using the functional form, and these fits are also shown in the figure. Uncertainties on the fits are estimated at the 95% confidence level, and are plotted in the figure (they essentially correspond to the line thickness). The fits provide a very good description of the endpoint coefficient and are far more efficient for subsequent evaluation of the soft function. The values of the fit coefficients of Eq. (123) for the three cases are collected in Table 2. For the channels gg → g and qq → g we set c (m,n) = c (n,m) (with m = n). We have tested the fits using randomly-generated phase-space points for comparison, i.e. with randomly-generated values for y 13 ∈ [0, 1] (obtained in our case using Mathematica [56]). Within their respective MC uncertainties, every event agrees with the prediction from the fits. When the MC uncertainties are small the agreement is within 0.5%, while for points for which the coefficients are close to zero the MC uncertainties are significantly larger and the agreement is at the few percent level.
To assess the overall impact of the corrections to the soft function, in Fig. 6 we show the NLO contribution along with the Abelian and non-Abelian pieces of the NNLO contribution. It is clear that the Abelian pieces are significantly larger than the non-Abelian pieces over the entire phase space for each partonic configuration. This can be readily understood from the analytic structures of the two pieces. The δ(T 1 ) coefficient of the Abelian part arises from the convolution of the NLO result with itself, and therefore it has two contributing parts. Firstly, there is the "pure" δ(T 1 ) piece, which arises from the δ(T 1 ) coefficient of the NLO soft function squared. This term is comparable to the non-Abelian part Fig. 4 Numerical results for the non-Abelian part of the 1-jettiness soft function at NNLO for the partonic channels gg → g, qq → g, and qg → q. We plot the coefficients of δ(T 1 ) and L n (T 1 ) with n = 0, 3 as functions of y 13 ∈ [0, 1]. The analytic results are taken from Refs. [30,42,45] in size. Secondly, there are mixed contributions which arise from convolutions of the form where the coefficients V mn k can be found, for instance, in Table 1 of Ref. [30], and are roughly O (1). After both contributions are combined together the dominant part of the total coefficient is determined by the V 01 −1 C 0 C 1 term. Such a term is not present in the non-Abelian calculation.

Generic kinematics
Finally, we compute the non-Abelian part of the 1-jettiness soft function at NNLO with generic kinematics, i.e. for the case where the scattering process has three colored partons at Born level, but where two of them do not represent the directions of incoming beams (as is the case for LHC kinematics). We parametrize the three directionsp 1 ,p 2 , andp 3 as 123 (1, 0, 0, 1) with θ 2 , θ 3 ∈ [0, π] and φ ∈ [0, 2π ]. For the invariants y i j (with i, j ∈ {1, 2, 3}) we therefore explicitly have We observe that with this parametrization the LHC limit is recovered by setting φ = 0 and θ 2 = π . We compute the coefficients C n of the non-Abelian part of the soft function for the three different partonic configurations (ggg, qqg, and qgq) by choosing 200 random values for θ 2 , θ 3 , and φ. We then perform numerical fits to our results for the coefficient of δ(T 1 ), C −1 . The functional form of the fits is taken to be In order to obtain accurate fits, we retain all 64 coefficients in Eq. (127) for each partonic channel. The values of the coefficients c (k,m,n) for the three cases are collected in ancillary files that we include in the arXiv submission.
We can test the validity of the generic fits by ensuring that they correctly reproduce the dedicated LHC fits obtained in Sect. 6.2 when choosing y 12 = 1 and y 23 = 1 − y 13 in Eq. (127). We therefore define the following ratio In the above equation C ab −1 (y 13 ) represents the δ(T 1 ) coefficient of the Abelian part of the soft function (evaluated for LHC kinematics). We have added it to the numerator and denominator such that R fit compares the total NNLO δ(T 1 ) coefficients using the two different fits. This combination is the only one that is relevant for phenomenological applications. In addition, since C −1 vanishes for certain values of y 13 , the sum of Abelian and non-Abelian pieces helps to ensure that the ratio is not dominated by the regions in which one of the fits is close to zero. The ratios of the two fits are shown in the left panel of Fig. 7. By inspecting Fig. 7 we see that for the ggg and qqg channels the generic fit reproduces the dedicated LHC fit to better than 1%. On the other hand, the qgq channel is poorly behaved in the region y 13 ∼ 0.15. This is due to the vanishing of both the Abelian and the non-Abelian pieces in this region, which causes large sensitivity to fitting uncertainties. Away from this region the general fit does a good job (within 1%) at reproducing the dedicated LHC fit.
Additionally, as in the LHC case, we can test the fits by generating random phase-space points and comparing the numerical results of the non-Abelian piece with the predicted values from the corresponding fits. In this case we perform the comparison by generating 450 phase-space points (i.e. 450 random values for θ 2 , θ 3 , and φ) and plotting the ratios between fit and numerical values as histograms. The results are shown in the right panel of Fig. 7. We observe that for all three channels the histograms are well-centered around the value of 1. In particular, the values obtained with the fit lie within 2% of a dedicated calculation of a configuration for 74, 75, and 74% of the total number of random phasespace points respectively (ggg, qqg, and qgq), and within 5% for 87, 88, and 88% of the phase-space points. Similarly to the LHC case, the agreement between fitted and numerical values is better for points with smaller MC uncertainties and coefficients that are not close to zero. We therefore believe that the fits of Eq. (127), which are valid for non-LHC kinematics, could be successfully used, for instance, in ep and e + e − applications.

Conclusions
In this paper we have presented a calculation of the nextto-next-to leading order (NNLO) 1-jettiness soft function. The soft function is a component part of NNLO calculations which employ slicing methods based on the N -jettiness global event shape variable. In particular, this function is a required piece of the calculation of differential pp → X + j type processes at NNLO, in which X represents a color singlet, for instance a single vector boson.
Our calculation bears the traditional hallmarks of NNLO calculations in regards to its complexity in unresolved limits. In order to deal with these issues we have employed a numerical approach which uses sector decomposition of the relevant phase-space integrals to disentangle overlapping singularities present at this order. In order to ensure the correctness of our results we have implemented two completely independent computational codes, and validated them against one another. We have further validated our results by recomputing the known 0-jettiness and 1-jettiness results at NNLO and NLO accuracy respectively. As a final validation of our results we have checked the known analytic pieces of the 1-jettiness soft function at NNLO that can be derived from renormalization group arguments. The primary result of our calculation is a numerical determination of the δ(T 1 ) endpoint contribution, which cannot be deduced from the RGE's alone. We have computed this contribution for LHC kinematics (scattering processes with two back-to-back partons and one final-state jet at Born level) and for generic kinematics (three colored partons at leading order). In order to disseminate our results we have produced polynomial fits to the results of our numerical integration for both configurations. We have additionally checked both fits using randomlygenerated phase-space points, and find excellent agreement with our Monte Carlo output.
Our fits represent the first such results presented in the literature. A previous calculation of the 1-jettiness soft function at NNLO for LHC kinematics has been published [42]. It does not contain the information required to implement the results in a standalone Monte Carlo program. We therefore believe our calculation, and the corresponding numerical fits, will be useful for those interested in applying jettiness-slicing and subtraction methods to NNLO calculations where three colored particles are present. Further applications of this method are certainly possible in future, for instance to analyze a vari-ety of global event shape definitions or to increase the number of partonic scatters under consideration, such as in the calculation of the 2-jettiness soft function at NNLO. We leave such applications to a future study.
A (d − 2)-dimensional vector can be written as, where the components before the semi-colon are the extradimensional pieces and the vector n is a unit vector in the first extra dimension. We arrive at this frame by setting θ i = π/2 for i = 1, d −5. Let us choose vectors in the transverse plane given by so that for q 1 we have set θ d−3 = π/2, θ d−4 = φ 1 and for q 2 we set θ d−3 = φ 2 , θ d−4 = β. The integral over unconstrained angles is, so that, after integrating over them, we get the following expression By performing the integration over β we get Using the standard result, , (A. 13) the total volume of angular integration is (A.14)

Appendix B: Rotational invariance of the solid angle integral measure
Perform an angle-radius decomposition of the coordinates where x is a coordinate in the transverse plane beyond the usual two. Hence Now perform a rotation about the -axis by an angle φ 1 : These coordinates can also be parametrized by introducing new angles, . (C. 31) We note that since J I I I a i j does not depend on λ it may be treated using the phase-space measure in Eq. (87). All other contributions require the use of Eq. (97) and the further partitioning indicated in Eq. (99).

Appendix C.2: Case 2
Using the relations in Eq. (102), the matrix elements are given by, (C.35) All contributions may be evaluated using the phase-space parametrization given in Eq. (103).

Appendix C.3: Case 3
Using the relations in Eq. (106), the matrix elements are given by, (C.37) (C.39) All contributions may be evaluated using the phase-space parametrization given in Eq. (107).

Appendix D: Results for double-real integrals
Results for the O( ) coefficient of each of the basic integrals entering the double-real emission calculation, at a sample phase-space point, are given in Table 3. These coefficients are the ones that enter the calculation of the endpoint contribution to the NNLO soft function and may be useful to the reader interested in reproducing the results of our calculation.

Appendix E: Analytic result for 0-jettiness at NLO
For LHC kinematics, the 0-jettiness case refers to processes with two initial-state colored partons and no final-state jets (plus any non-colored final-state particle). The only allowed leading-order configurations are therefore either a pair of gluons or a quark-antiquark pair. For the purposes of our calculation, this means that i, j ∈ {1, 2} and that we need not consider the phase-space sector F i j k . Moreover, there is no θfunction involving an angle left in the problem. As discussed in Sect. 3.3, color conservation means that the overall color factor, C = −T 1 · T 2 , is given by either C = C F for a quarkantiquark pair or C = C A for a pair of gluons. The 0-jettiness soft function at NLO is thus given by Eq. (37) after setting For color-singlet production at the LHC the two initial-state partons are back-to-back so that y 12 = 1 (L 12 = 0) and the soft function simplifies to, (E. 48) 123 Table 3 Results for all the double-real integrals at the point y 13 = 0.9 with i = 1, j = 2, and k = 3. The