Two-loop $N$-jettiness soft function for $pp \to 2j$ production

We calculate the two-loop soft function for $pp \to jj$, $e^+e^- \to 4j$ or $e p \to 3j$ when the $N$-jettiness observable is measured. The result presented here makes up the necessary piece for realizing the full next-to-next-to-leading order predictions for these processes using the $N$-jettiness subtraction scheme.


I. INTRODUCTION
The accumulation of the data set and the continuously increasing precision of the experimental analyses bring the physics program at the Large Hadron Collider (LHC) and other experiments the ability to search for new physics from the small deviations of the data from the Standard Model predictions. In order to fulfill such searches, the precision of the theoretical predictions from Quantum Chromodynamics (QCD) has to match the small experimental errors to reliably interpret the collider data. Currently, the only theoretical tool, with no known substitutes, enables us to make predictions for collider phenomenologies out of the first QCD principle is the perturbative calculation based on the systematic expansion of the strong coupling constant α s or α s L if large logarithmic correction L is present. In recent years, tremendous efforts have been made in realizing the theoretical predictions for the physical processes at the colliders with a full control of the final state kinematics beyond the next-to-leading order (NLO) accuracy. On the fixed order side, significant progress has been made in the past few years, which includes achieving the fully differential gluon-gluon fusion Higgs (H) production [1,2] and the single jet production in deep-inelastic scattering [3] at next-to-next-to-next-to-leading order (N 3 LO) in α s and the next-to-next-to-leading order (NNLO) calculations of the single inclusive jet or dijet productions [4], V /H + j processes [5][6][7][8][9][10][11], Higgs production through the vector boson fusion [12,13], single top [14,15] and tt productions [16,17] at the LHC, the jet productions in deep-inelastic scattering [18,19], and the heavy flavor production in deep-inelastic neutrino scattering [20].
The rapid growth in the higher order predictions benefit from either the new higher loop calculations achieved or the theoretical schemes developed for dealing with the singularities in the real emissions. Recent achievements of the analytic two-loop five-point amplitude [21][22][23] highlight the current status of the loop calculations for the collider processes, and meanwhile new ideas are keeping on popping out [24][25][26][27]. To handle the infrared (IR) singularities in the real corrections, various schemes have been proposed and successfully applied to the LHC phenomenology, such as the local subtraction schemes including the antenna subtraction [33], the sector improved residue subtraction [34] and the nested soft-collinear subtraction [35], as well as the physical observable based global schemes like the q T -subtraction [36], inclusive jet mass subtraction [37] and the N -jettiness subtraction [5,38]. Although the global schemes are always criticized for the numerical instability, they have realized many interesting higher order predictions for the physical processes at the hadron colliders including the numerically most challenging Higgs production at N 3 LO [1].
Among the local schemes, N -jettiness scheme is the known scheme being able to handle the IR singularities for generic jet processes at the hadron colliders. The scheme utilizes a threshold cut-off τ cut in the N -jettiness event shape τ [39] to distinguish between the NNLO N jet (when τ < τ cut ) and the NLO N + 1 jet configurations (τ > τ cut ). Below the τ cut , a Soft-Collinear-Effective-Theory [40][41][42][43][44] approach based on the factorization theorem [39] is used to approximate the full QCD contribution. Here B stands for the beam function due to the initial state energetic radiations while J's are the jet functions for the final state collinear radiations. The hard function H encodes the hard loop information. S is the soft function for soft radiations. The trace is over the relevant color space. All the ingredients in the factorization theorem are universal except for the hard function H which depends on the specific process under consideration. The difference between the factorization theorem and the full QCD is suppressed by power corrections in τ cut , which have been studied extensively recently for the 0-jettiness case [45][46][47][48][49]. The universal beam and the jet functions are all known analytically to two-loops [50][51][52][53] and recently some of them are even calculated at three-loop order [54][55][56]. As for the quark beam function, its NNLO longitudinal polarized counterpart is also known in the literature [57]. The calculation of the soft function is complicated by its explicit dependence on the full N -jettiness measure, therefore its analytic NNLO calculation is by far out of reach and people turn to numerical solutions. A generic framework was proposed in [28] for calculating such soft functions with arbitrary numbers of jets and the 1-jettiness soft functions in pp-collision were calculated using this scheme or similar approach for both the massless jet [28,29] and the massive heavy flavor productions [30,31]. However, at NNLO, the most general calculations can only be covered when one considers the N -jettiness soft function with four external hard legs, which is missing in the original calculations. In this manuscript, we present such calculation by extending the previously developed numerical method and show sufficient calculating details. We noticed a recent conference proceeding [32] on the 2-jettiness soft function using pySecDec. Our work uses a different and independent approach which can serve as a cross check of the calculations. The rest of this work is organized as follows. We discuss the numerical setups and the phase space parameterization in detail in Section II. In Section III, we present our numerical results for the 2-jettiness observable in hardon-hadron collisions. Last we conclude in Section IV

II. COMPUTATIONAL SETUP
In this section, we summarize the setups used in our calculation. We follow the general numerical framework proposed in [28], and show its capability in dealing with N -jettiness with four hard reference vectors in pp → jj, ep → 3j and 0 → 4j. We note that the calculations with four reference vectors go through all possible configurations in calculating any N -jettiness soft functions at NNLO. We clarify some of the computational details of this framework.
For simplicity, we stick to the hadronic N -jettiness definition which reads where n i 's are the N light-like hard reference vectors along the initial state hadronic beams or the final state jet axes and q i 's are the four momentum of the soft radiations. We note that here all the hard vectors n i 's are physical objects and are defined in 4-dimension while the soft momentum q i 's can potentially be un-resolvable and are therefore d-dimensional quantities, with d = 4 − 2 in the dimensional regularization. The soft function with n real soft emissions then can be calculated by performing the phase space integration where we have introduced the N -jettiness measurement function and S is the soft current matrix element for the soft radiations which, up to NNLO, has been known for a long time [58,59]. For N -jettiness with 3 reference vectors , to NNLO, the soft current can be expressed in terms of a sum of dipole contributions, while starting from four reference vectors, like the case we study here, new triple-pole structure involving 3 eikonal lines arises which is induced by the one-loop soft current. The triple-pole contribution, which will be calculated later in this work, is the only remaining piece for calculating the N -jettiness soft function with arbitrary external eikonal lines at O(α 2 s ). For completeness, we summarize the explicit form of the soft current S to O(α 2 s ) in the Appendix V A. In order to compute the soft function numerically, the key step is to find a suitable momentum parameterization to allow us to write the phase space integration in the form of which is ready for the Laurent expansion so as for us to identify all the -poles and evaluate the coefficients of the poles numerically. We thus now turn to the discussion of the phase space parameterization for single and double real emissions with four reference directions, separately. We note that though we are focusing on the four reference vector case, the parameterization is general enough to be applicable directly to any N -jettiness computation at the NNLO. In all cases below, we utilize the Sudakov decomposition to write an arbitrary light-like four vector k µ with respect to two of the reference vectors n i and n j as with k + = n i · k, k − = n j · k and k ⊥ · n i = k ⊥ · n j = 0 and |k ⊥ | 2 = 2 n ij k + k − . Here we have introduced the notation n ij ≡ n i · n j .
A. single real emission phase space parameterization With 4 reference vectors, n i , n j , n k and n l , at hand, for one real emission with momentum q µ 1 for computing the NLO soft function or the NNLO real-virtual contribution, we can parameterize all the involved vectors with respect to n i and n j , which are where we have used the freedom to choose the azimuthal angle of the vector n k and the additional -angles other than α 1 to be 0. Here α 1 , φ l , φ 1 ∈ [0, π] with c φ l = n ik n jl + n il n jk − n ij n kl 2 √ n ik n il n jk n jl .
Here the non-zero φ l makes n k and n l span the azimuthal plain and therefore α 1 is necessary for parameterizing the d-dimensional momentum q 1 , since q 1 in general lies outside the pain n k and n l belong to. We note that since the reference vectors n i 's are physical, for arbitrary numbers of the reference vectors, their "⊥"-components will all be lying in one azimuthal plain, and thus Eq. (7) is the most general parameterization for computing the N -jettiness function with arbitrary N . For the processes similar to what we are considering here, the transverse components n ⊥ k and n ⊥ l have to align with each other (φ l = 0 or π) due to the momentum conservation, therefore one can simplify the parameterization by setting α 1 = 0 which reduces to the parameterization used in [28,29] for computing 1-jettiness.
Although in general n i and n j can be chosen arbitrarily for performing the Sudakov decomposition, in order to correctly isolate the -poles, it is more useful to choose n i so that n i · q contributes to τ as set by Θ N,i while n j · q appears as one of the singular poles in the denominator of the soft current matrix S. With this parameterization, for non-zero α, the phase space integration for single emission can be modified to where we have made the variable changes It can also be shown that for α 1 = 0, one can further modify the phase space integration to make φ 1 ∈ [0, 2π] when taking into account the properties of the soft current and the N -jettiness observable, as has been done in [28]. Eq. (9) is our starting point for calculating the single emission contributions to the soft function numerically and we list in the Appendix V B all the final forms of the phase space integrations suitable for direct numerical evaluation of the one emission contributions.

B. double real emission phase space parameterization
Now we sketch the parameterization for the double real emission with momenta q 1 and q 2 at NNLO. We want to add on top of Eq. (7) the momentum q 2 parameterization, meanwhile require that the dot product q ⊥ 1 · q ⊥ 2 only relies on the azimuthal angle between these two vectors, so as to use the non-linear transformation [28] in a later stage for the purpose of pole isolation. To find the particular parameterization of q 2 , we first rotate Eq. (7) around the n k ⊥ axis by α 1 to eliminate the α 1 dependence in q µ 1 , then we rotate in the azimuthal plain by φ 1 to align the x-axis with q ⊥ 1 to further remove the dependence of φ 1 in the parameterization of q µ 1 . These rotational operations lead to the momentum representation where we introduced the notation c α 1 It is now straightforward to see that to satisfy our requirements, q µ 2 can be parameterized (in this frame after rotation) as The β angle here is necessary to guarantee that thecomponents of q 1 and q 2 are not always aligned with each other. Due to the Lorenz invariance of the individual emission phase space measure, the rotational operations used to derive the parameterization of q 2 do not introduce any complexity to the Jacobian of the phase space integration. We also note that the parameterization derived here is the most general for computing the double real contributions to any N -jettiness soft functions. As for the processes we are interested in in this work, using the same momentum conservation for the single real emission case, one can set α 1 and β to 0 to simplify the phase space integration.
Within this parameterization, the dot product q 1 · q 2 takes the form which generates a linear singularity in the double real soft current if q 1 be collinear to q 2 is possible, which prevents us from extracting the -poles using the Laurent expansion in Eq. (5). To extract the poles, the angle c φ 2 is further parameterized as c φ 2 = 1 − 2η, with the non-linear parameterization if the soft momentum q 1 and q 2 can be collinear to each other, and if q 1 and q 2 can no way be collinear as constrained by the jettiness measurement, for instance when τ = n i · q 1 + n j · q 2 .
Last we comment on the choice of n i and n j in the Sudakov decomposition. The same logic follows the single emission case. The reference vectors n i and n j are chosen if n i · q 1 and n j · q 2 contribute to τ , i.e. τ = n i · q 1 + n j · q 2 . Otherwise, if both n i · q 1 and n i · q 2 contribute to the jettiness observable, then n i will be naturally picked as one reference axis for the decomposition, and for convenience, n j is chosen in the way that the double real soft current is singular as n j · q 1 → 0 and/or n j · q 2 → 0. Additional variable changes are needed to map all variables onto the regime [0, 1] suitable for doing the Laurent expansion. The mapping is standard in the sector decomposition calculations and has been detailed in [29,31]. After the mapping, in the cases where one encounters the fractional power of the form s − 1 2 + with s ∈ [0, 1], a further variable transform s = sin 2 π 2 x , with x ∈ [0, 1], is then used.

III. NUMERICAL RESULTS
In this section, we present our main results. We highlight our results using the 2-jettiness soft function to O(α 2 s ) in pp → jj. Other processes such as e + e − → 4j can be easily obtained by changing properly the values of λ ij in Eq. (21), without doing any additional integrations. The 2-jettiness soft function depends on four hard directions, which we denote as n 1 , n 2 , n 3 and n 4 . We align n 1 and n 2 with the incoming beam axes in the ±z-directions, and let n 3 and n 4 lie along the outgoing jet directions. For the purpose of presenting, though not at all required in our calculation, we let the incoming partons carry exactly the same energy in the collision and show the 2-jettiness soft function as a function of n 13 and we present the results in the Laplacian space ∞ 0 dτ e −τ z σ(τ ). We start from the NLO soft function. At NLO, the renormalized 2-jettiness soft function can be written as a sum of dipole terms where C (1,i) ij are the coefficients of the i -poles returned by our calculation using the method for single emission as sketched before. The coefficients of the logarithmic termsL 2 andL can be predicted by solving the renormalization group equation (RGE) of the soft function and therefore serve as checks of our numerical approach.
In fig. 1, we show the NLO coefficients of theL n (n = 2, 1, 0) terms for dipole S dots represent the coefficients forL 2 ,L 1 andL 0 from our numerical set-ups, respectively. The solid lines are obtained by solving the RGE analytically. We see perfect agreements between the analytic results and our numerical predictions, which implies the validity of the computational framework at this order. At NNLO, other than the abelian terms which can be obtained by exponentiate the NLO results, the soft function receives non-abelian contributions from the qq-emission, the ggemission, the real-virtual corrections and the renormalization. The real-virtual corrections can be further broke down to the dipole and the triple-pole contributions, and the rest are of the dipole forms. The summation of all the dipole contributions can be organized by the color factors C A and N F T R , while the triple-pole real-virtual corrections can be massaged into one single color term f abc T a 1 T b 2 T c 3 with the help of the color charge conservation T i = 0. The final results of the NNLO soft function can therefore be written as where all the NNLO dipoles S (2) ij and triple-pole S (2) 123 in the non-abelian terms can be written as a series in terms of the logarithmsL. And once again, all the coefficients of those logarithmic termsL n with n = 3, 2, 1 can be predicted by RGE with the knowledge of the full NLO results obtained before.
The coefficients of the N F T R terms for the dipoles S 12 and S (2) 13 are shown in fig. 2, which include the contributions from the qq emission and the N F term in the renormalization. Perfect agreements are observed between the numerical calculations of theL n coefficients for n = 3, 2, 1 as shown in colored dots and the RGE predictions in solid lines, which validates not only the NNLO computations of the N F contributions but also the correctness of the non-logarithmic term in the NLO calculations. The NNLO non-logarithmic prediction (in black dots) displayed here is not presented in the original NNLO N -jettiness soft paper [28] and is one of our main results of this work. We note that the NNLO results here and below are also presented in a conference proceeding [32] in a different way with a different approach.
Same situation happens to the C A term, as displayed in fig. 3, in which we also compared L n coefficients from the numerical calculation against the results from the RGE to find no difference. The non-logarithmic coefficients are represented by the black dots.
Last we give the prediction for the triple-pole term S (2) 123 , which arises firstly for four hard reference directions. Just like the dipole contributions, the logarithms in the triplepole are still predictable from the RGE. The analytic predictions are once again found to . agree with our numerical results, as clearly seen in the left panel of fig. 4. Other than that, we also compare the -pole terms without theL's from our numerical calculations with the ones predicted by the RGEs and find satisfactory agreements. The result for the nonlogarithmic term (left panel in fig. 4) as well as the term entirely unpredictable from RGE (right panel, normalized further to −2π) are indicated in black dots, which finalizes all the NNLO contributions to the 2-jettiness for the pp → jj process. .

IV. CONCLUSIONS
In this work, we present the calculation of the N -jettiness soft function with four hard reference directions using the computational framework proposed in [28] using the approach of sector decomposition. We generalize the parameterizations used for the 1-jettiness calculation to the arbitrary N -jettiness case. We specifically calculated the NNLO N -jettiness soft function with four external hard legs, where the last new piece at this order, the triple-pole configuration, first arises. We managed to isolate all the singularities from the soft function integrands and reduce the computation to a set of numerical integrals which are ready to evaluate. We check the numerically computed logarithmic terms in the soft function against the predictions from the RGEs and found perfect agreements. The non-logarithmic terms which can not be obtained through the known RGEs are calculated directly in this work. We expect the result obtained in this manuscript to find its quick applications in both the fixed NNLO calculations and achieving the parton shower matching [60] for the relevant processes. For completeness, we list all the soft current matrix elements here, which were first derived in the seminal works [58] and [59] by Catani and Grazzini. The NLO soft current with one soft emission q is given by with the NLO dipole given by The NNLO real-virtual correction can be written as the sum of the dipole S RV,dipole and triple-pole S (2) RV,tri. contributions, with the dipole contribution be , (20) and the triple-pole term where λ ij = 1 if i and j are both incoming or outgoing, otherwise λ ij = 0. The NNLO real-real emission is given by from the qq emission. To get this form we have used the color charge conservation i T i = 0. Here and where q 1 and q 2 are the momentum carried by the real emissions. We denote the term proportional to 2 n ij q 1 · q 2 as I II ij and the rest I I ij . And the gg emission gives the non-abelian contribution with where and the current in the strongly ordering limit The abelian piece for the gg double real emission is

B. numerical integrations for single emission
Here we list all final numerical integrations for evaluate single emission contributions. We note again that though we only show the N -jettiness with four hard directions n i , n j , n k and n l , the parameterization here presented here is general enough for any N -jettiness soft functions. To extend to the case with more reference vectors, one simply insert more θ functions originated from the N -jettiness measurements. We use the notation q||n i to denote the separation of q and n i is the smallest and thus τ = n i · q.

final forms for numerical integration at NLO
At NLO, the contribution can be written as a sum of dipoles. We assume that the soft momentum q is emitted from dipole ij. For the integrations shown below, we have extracted out an overall factor −16π 2 αs 2π e γ E 4π T i · T j .
For q 1 n i , the final form suitable for numerical integration is where we have defined and As for the case in which q 1 n j , the result can be obtained by switching i and j in Eq. (30).
While for q 1 n k , we choose n k and n i to be the reference vectors, to have q + 1 = n k · q 1 , q − 1 = n i · q 1 . Following the same variable changes shown in Eq. (10), we get the final form for the numerical evaluation where Again, we can switch k and l to obtain the q 1 n l contribution.

final forms for 2-parton correlated real-virtual contribution
For the 2-parton correlated real-virtual contribution, we have for our numerical evaluation • case 1: q 1 n i .
• case 2: q 1 n k . I k, (2) ij,RV = π 4 n ij n ik n ij n ik 4 1 2π Here we have normalized to The q||n j and q||n l cases can be easily obtained by switching properly the indices.

final forms for 3-parton correlated real-virtual contribution
Now we turn to the triple-pole case in the real-virtual correction, in which 3-parton correlated emission contributes. This configuration first arises in this case with four external legs. We assume the soft emission is from triple-pole ijk. And we normalize our results to where F ij i has been defined before.
• case 2: q 1 n j We let q + 1 = n j · q 1 , q − 1 = n i · q 1 to find • case 3: q 1 n k We let q + 1 = n k · q 1 and q − 1 = n l · q 1 to get I k,(2) ijk = π 4 n ik n kl n ij n kl 4 1 2π Here we redefine F ij k (x 3 ) as • case 4: q 1 n l We let q + 1 = n l · q 1 and q − 1 = n k · q 1 to have I l,(2) ijk = π 4 n ik n kl n ij n kl 4 1 2π where