Next-to-Next-to-Leading Order $N$-Jettiness Soft Function for One Massive Colored Particle Production at Hadron Colliders

The $N$-jettiness subtraction has proven to be an efficient method to perform differential QCD next-to-next-to-leading order (NNLO) calculations in the last few years. One important ingredient of this method is the NNLO soft function. We calculate this soft function for one massive colored particle production at hadron colliders. We select the color octet and color triplet cases to present the final results. We also discuss its application in NLO and NNLO differential calculations.


Introduction
After the discovery of the Higgs boson, one of the main tasks of the Large Hadron Collider (LHC) is to discover new physics beyond the Standard Model (SM). So far there has been no sign of new physics at the LHC. In the future, it is possible that the new physics would appear in small deviations from the SM predictions. Therefore, along with the accumulating experimental data, accurate and precise theoretical predictions for the SM processes are mandatory for discovering new physics. Including higher-order effects in quantum chromodynamics (QCD) perturbation theory can help to improve the theoretical predictions. More specifically, one needs to consider fixed-order calculations for general observables, and analytical resummation or parton shower simulations to resum large logarithmic terms for certain delicately designed observables.
Soft-collinear effective theory (SCET) [1][2][3][4][5] has been widely used in analytical resummation and fixed-order calculations; see Ref. [6] for a review. One of recent progresses is the application of SCET to studying the N -jettiness event shape variable T N [7]. When T N → 0, the cross section is factorized as the convolution of a hard function, two beam functions for initial states, a soft function and jet functions [7,8], (1.1) where the N -jettiness event shape variable is defined as [7] T N = Here n i (i = a, b, 1, ..., N ) are light-like reference vectors representing the moving directions of massless external particles, and q k denotes the momentum of soft or collinear radiations. The hard function H encodes all the short distance information and can be obtained from the Wilson coefficient when matching full QCD onto SCET. The beam functions B i , (i = 1, 2), incorporate the effects of parton distribution functions as well as initialstate radiations, which are known up to next-to-next-to-leading order (NNLO) [9][10][11][12].
The jet function J n describes the final-state jet with a fixed invariant mass and has been calculated at NNLO [13,14]. The soft function S describes soft interactions between all colored particles. It has been studied up to NNLO for massless processes recently [15][16][17][18][19]. Based on these ingredients, the resummed N -jettiness distribution can be obtained [8,10,[20][21][22][23][24][25][26]. Besides the analytical resummation, the N -jettiness subtraction method, based on the expansion of Eq.(1.1) in series of the strong coupling α s , has been proposed to perform NNLO calculations [27][28][29], and been utilized to provide differential predictions for processes with final-state jets at a hadron collider [30] and processes at an electron-hadron collider [31,32]. In principle, the extension to massive external particles is possible and has been mentioned in Ref. [29]. Because there are no collinear singularities along the moving direction of the heavy particle, the definition of T N in Eq. (1.2) can control all the infrared (both soft and collinear) singularities. As a result, there is no need to introduce extra N -jettiness axes for heavy particles, such as in the case of top quark pair production [33][34][35] and single top production and decay [36,37]. The effect of including heavy external particle amounts to modifying only the soft function in Eq. (1.1). We present in this work a calculation of the soft function for one massive colored particle production up to NNLO.
We notice that the soft function for such a process in the energy threshold has already been obtained in Refs. [38][39][40], which would be useful for threshold resummation. However, the measurement function for soft function in the N -jettiness factorization is different, as indicated by Eq. (1.2). Moreover, in threshold resummation, there are no beam functions, which means the collinear-radiation effects from the initial states are not taken into account in the factorized cross section. So it can not be directly used in fixed-order calculations.
This paper is organized as follows. In the next section we briefly introduce the definition of the soft function for one massive colored particle production at hadron colliders. In section 3, we derive the structure of the soft function from its properties under renormalization group evolution. Then, in sections 4 and 5, we present the results of next-to-leading order (NLO) and NNLO soft functions, respectively. With these results at hand, we discuss its application to NLO and NNLO calculations in section 6. We conclude in section 7. In the appendixes, we show the soft function for a color-triplet fermion production and results of master integrals in virtual-real corrections.

Definition of the soft function
In this section we discuss the factorization of the cross section for one massive colored particle production and present the definition of soft function.
We consider the process where P 1 and P 2 denote incoming hadrons, Q represents the massive colored particle, and X includes any inclusive hadronic final state. The massive colored particle can be massive quarks in the SM or some new particles in extensions of the SM. For quark-antiquark initial state, the partonic process at leading order (LO) is For later convenience we introduce two light-like vectors 3) The momenta given in Eq. (2.2) can be written as where m is the mass of the particle Q. The 0-jettiness event shape variable in this process is defined as In the limit τ ≪ m, the final state contains only one massive particle along with possible soft and collinear radiations. The cross section in this limit admits a factorized form which can be derived in the frame of SCET. According to the analysis in Refs. [7][8][9], the cross section can be written as where σ 0 is the LO partonic cross section, Y is the rapidity of the partonic colliding system in the laboratory frame, the momentum fractions x a = m/ √ se Y and x b = m/ √ se −Y with √ s the colliding energy, and µ is the renormalization scale. The soft function is defined by the vacuum matrix element where T(T) is the (anti-)time-ordering operator, and Y n , Yn and Y v are the soft Wilson lines. Explicitly, they are defined as [4,41,42] with P(P) the (anti-)path-ordering operator acting on the color operator T a . The operator P k in Eq.(2.7) is defined to extract each soft gluon's momentum. The purpose of this paper is to calculate the soft function defined above up to NNLO. For comparison, the threshold soft function can be obtained by changing the measurement function in Eq.(2.7) to which has been studied before [38][39][40].

Renormalization
Before showing the details of our calculation, it is instructive to study the structure of the soft function first, which can guide our calculation and serve as a check of the final results. The soft function in Eq.(2.7) is defined in terms of bare parameters and contains ultra-violate divergences when calculated in QCD perturbation theory. These divergences can be absorbed into counterterms, leaving the renormalized soft function finite. Following the convention, we make use of dimensional regularization, i.e., generalizing the dimension of space-time to d = 4 − 2ǫ. After renormalization, the soft function depends on the renormalization scale µ in Eq.(2.6). This intermediate scale can appear only in logarithmic forms or in α s (µ), the former completely determined by the corresponding counterterm. The similar situation happens for the other components in Eq.(2.6). Since the all-order cross section does not depend on this intermediate scale, the dependence on the scale µ will cancel against each other order by order in α s . As a consequence, we can derive the logarithmic terms of the scale µ in the soft function from the knowledge of the NNLO hard and beam functions. According to dimensional analysis, the bare soft function can be written in QCD perturbation theory as where we use renormalized strong coupling α s and its renormalization factor Z αs = 1 − β 0 α s /(4πǫ) + O(α 2 s ). The coefficients s (n) do not depend on τ . It is convenient to discuss the renormalization group equation by using Laplace transformed soft function, given bỹ Then the corresponding renormalized soft functions is defined as where the renormalization factor Z s satisfies the differential equation with γ s the anomalous dimension of the soft function.
If we know the result of the anomalous dimension γ s , integrating the above equation as in Refs. [43,44], we obtain a closed expression for Z s Here we assume the expansion series Now expanding Eq. (3.3), we get the renormalized NLO and NNLO soft functions in Laplace spaces Since the renormalized quantities are finite, we can infer the coefficients of ǫ −i , i = 1, 2, ..., in the bare soft functionS.
In the above derivation, we have assumed the knowledge of γ s . Actually, we do know it from the independence of the cross section on the renormalization scale µ. Performing Laplace transformation of Eq.(2.6), we find whereB i is the beam function after Laplace transformation. The NLO and NNLO beam functions have already been calculated in Refs. [9][10][11][12]. It is found that the beam function satisfies the same renormalization group evolution as the jet function to all orders [9]. We write the renormalization group equation for the beam function as where T i is the color generator associated with the i-th parton [45,46], γ cusp and γ i B are anomalous dimensions controlling the double and single logarithms, respectively [11,12]. The subscript or superscript i can be 1 or 2, denoting the beam functions of the two initial-state particles, respectively.
The renormalization group equation for the hard function was studied at NNLO in Refs. [47,48]. In our case, there is only one color basis. It is straightforward to organise the renormalization group equation for the hard function as The anomalous dimensions γ 1,2 and γ Q , associated with the initial-and final-state particles, can be found in Refs. [47,48] and references therein. Inserting Eqs. (3.8-3.9) to Eq.(3.10), we obtain the anomalous dimension of the soft function, of which the explicit expression is available up to NNLO.

NLO soft function
The LO soft function is trivial and has been given explicitly in Eq.(3.1). In this section, we present its NLO result. Expanding the soft Wilson lines in Eq.(2.7) in series of the strong coupling, we write the NLO soft function as where e γ E ǫ is inserted because we use MS renormalization scheme. d µν (q) is the polarization tensor, defined as with ε µ (q) the polarization vector of the soft gluon. The factor J µ(0) a (q) is the LO one-gluon soft current, or the eikonal current,  with a the color index. The color conservation implies q · J (0) a (q) = 0 so that we can choose The factor F (n,n, q) is a measurement function, embodying the constraint from the 0jettiness variable. In the center-of-mass frame, it is defined as where we use the notations, q + = q · n and q − = q ·n.
After performing the phase space integration, we obtain the NLO bare soft function where ψ 0 (x) is the digamma function which is defined as ψ 0 (x) = d ln Γ(x)/dx .

NNLO soft function
Due to the interest of practical application, we consider two kinds of color structures. One is 3 ⊗3 → 8, denoting a massive color octet production from a triplet and an anti-triplet initial states. The other is 3 ⊗ 8 → 3, representing a massive color triplet production from a triplet and an octet initial states. We will present the result of a color octet production in the main text while leave the result of a color triplet production in the appendix. The extension to other color structures can be obtained immediately.
The NNLO contribution consists of two parts, i.e., The first part is the virtual-real correction, i.e., the one-loop virtual corrections to LO soft gluon current; the second part is the double-real correction, i.e., the corrections with a double-gluon soft current or a massless quark-pair emission. Because they have different final states, they can be calculated independently.

Virtual-real corrections
The virtual-real contribution to the NNLO soft function comes from the situation when two emitted soft gluons contact together to form a loop, as shown in Fig. 1. In principle, one needs to calculate this kind of contribution by expanding the soft Wilson lines in Eq.(2.7) and selecting the diagrams with a loop. In practice, one can utilize the result of the soft limit of one-loop QCD amplitudes which has been studied in Refs. [49][50][51] and Ref. [52] for massless and massive external particles, respectively. The unrenormalized one-loop soft current can be written as [52] J µ(1) where f abc are the structure constants of the group SU c (3), and the factor g ij (ǫ, q, p i , p j ) is given by  The virtual-real contribution to the soft function is given by

Double-real corrections
We now consider the double-real contribution to the NNLO soft function. The selected diagrams are shown in Fig. 2. This part can be consider as the real correction to the LO soft current. Instead of using SCET to write down the amplitudes for each diagram, we make use of the results in Refs. [55,56] where the infrared behaviour of tree-level QCD amplitudes at NNLO has been analyzed. There are two types of contributions in this part, i.e., the emission of two soft gluons and the emission of a soft quark-antiquark pair. Schematically, S We will discuss them separately. The two-gluon soft current is given by [55,56] J µν(0) a 1 a 2 (q 1 , q 2 ) = where we have denoted the momenta of the two gluons as q 1 and q 2 , respectively. Then we can get its contribution to the NNLO soft function, where F (n,n, q 1 , q 2 ) is the measurement function for double-real emission, One can see that the whole phase space is partitioned to four pieces. Actually, since the two-gluon soft current is symmetric under the exchange of the two gluons, one needs to calculate only two of them.
In the case of the emission of a soft quark-antiquark pair, the tree-level amplitude square is factorized as [55] |M(q 1 , q 2 , p 1 , with the color factor T R = 1/2 in the SU c (3) group. Therefore we can obtain the corresponding contribution to the NNLO soft function as T ij (q 1 , q 2 )F (n,n, q 1 , q 2 ) .

(5.12)
It is convenient to perform the phase space integration in the light-cone coordinates due to the constraint from the measurement function. In particular, we take the parametrization of the phase space as 13) and express any scalar production of two momenta in terms of their components in the light-cone coordinates. Then we insert two identities to extract the contributions from the two hemispheres. The choice of q ± i , i.e., q + i or q − i , depends on which piece in Eq.(5.9) is selected. It is easy to calculate the integration involving the δ function. After that, we parametrize the variables τ 1 and τ 2 as Finally, the integrals we need to calculate boil down to four-fold integrals over a unit hypercube. In fact, there are in total over two hundreds of integrals in our calculation because of the different kinematic structure; e.g., the momentum p i in Eqs. (5.7,5.11) can be massless or massive. Another reason is that the constraints imposed by the measurement function in Eq.(5.9) hinder the application of the normal reduction method to present a large quantity of scalar two-loop integrals in terms of a few master integrals. For illustration, we show explicitly one example of such integrals, For our purpose to obtain an NNLO cross section, we only need the exact results up to O(ǫ) (notice the known factor τ −1−2nǫ in the final result). As a result, we decide to calculate these integrals numerically up to the desired precision, as in the calculations of the dijet soft functions for e + e − collisions [57] and the one-jettiness soft function for proton-proton collisions [18]. We have adopted two different methods to deal with these four-fold integrals so that they can provide a cross-check. In the first method, we apply the Mellin-Barnes representation and then evaluate the Mellin-Barnes integrals numerically using MB packages [58,59]. In the other method, we make use of the SecDec [60,61] program, which is based on the sector decomposition method [62], to calculate these integrals numerically. For example, the integral in Eq. ( We find that the results from both methods agree well with each other. The final result of the double-real contribution is As mentioned below Eq. (3.7), the coefficients of the poles in the NNLO bare soft function can be determined from the anomalous dimensions. Meanwhile, we have obtained the analytical result of s (2) VR up to O(ǫ) in the last subsection. As a result, we can derive the result of s (2) DR up to O(ǫ 0 ) without any hard calculation. The comparison between this result and Eq.(5.18) is presented in Tab. 1. We can see that the two results agree very well, which can be considered as a nontrivial check of our calculation.

Renormalized soft function
In the last two subsections, we have already obtained the bare soft function at NNLO. It is straightforward to perform the renormalization according to Eq. (3.3). In the end, we obtain the following renormalized soft function in Laplace space for a color octet productioñ DR in two different methods. In the first method, denoted by SCET prediction, the results are obtained from the anomalous dimensions and the analytical expression of s (2) VR . In the second method, denoted by real calculation, the results are what we calculate by Eq. (5.18). Their difference is defined as the real calculation minus the SCET prediction.
with 20) and A 1 is the coefficient of ǫ in s DR + s (2) VR . Notice that the scale dependence is encoded in the variable L. The above equation is the main result in this work. It is ready to extend to other color structures.

Application to fixed-order calculations
Given that the N -jettiness subtraction has proven to be successful in performing NNLO differential calculations, we anticipate that the soft function we just obtained can be used in fixed-order calculations. In this section, we take the massive color-octet vector boson production at the LHC as an example, but the method is general.
In fixed-order calculations, the cross section is divided to two parts where τ cut is an intermediate cutoff parameter. When τ cut is small compared to the typical scale of the process, such as the mass of the massive color octet, the first part can be approximated very well by Eq. (1.1). The beam functions have been already known up to NNLO [9,11,12] and the soft function is obtained up to NNLO in this work. The hard function is available only up to NLO [63]. Therefore, we can calculate this part using these functions in effective field theory. The second part, due to the constraint τ > τ cut , corresponds to the process of a color octet associated with a jet production. We only need its cross section at (N)LO for the calculation of dσ in Eq.(6.1) at (N)NLO. As a result, it can be computed with various existing packages, e.g., MadGraph5 aMC@NLO [64]. The detailed discussion of the higher-order effects in a colored particle production using the N -jettiness subtraction deserves another paper, which we leave to future work. Instead, we now show a spectrum of the τ distribution when τ is small in Fig. 3. We consider a color octet vector boson production at the 13 TeV LHC. It has a mass of 1 TeV and couplings r L = r R = 1, which are defined in Ref. [63]. The renormalization and factorization scales are both chosen at 1 TeV, and the CT14NLO [65] parton distribution function is used. We provide the results calculated by Eq. (2.6) in SCET and the automatic package MadGraph5 aMC@NLO with the model generated by FeynRules [66]. From Fig. 3, we observe that they are coincident with each other if τ < 1 GeV, which means that the prediction of the effective field theory reproduces the full QCD fixed-order calculation in the limit τ → 0, as we just discussed.
Notice that though the spectrum shows that the differential cross section becomes very large when τ decreases, it does not mean the prediction for an observable is infinite. Actually, there is another large contribution, which has an opposite sign, from the term with δ(τ ) in each order of α s , as shown in Eq.(3.1) after the expansion The plus distribution is defined as If an observable is infrared safe, i.e., not sensitive to the soft and collinear radiations, then one needs to integrate over the region near τ = 0, e.g., from 0 to ∆(> 0), and thus the final result is finite.

Conclusions
Precise theoretical predictions are important for both exact estimate of the backgrounds and the search for new physics. To obtain reliable predictions, QCD higher-order effects need to be taken into account. N -jettiness subtraction is one of the efficient methods to compute the fully differential cross section at NNLO. In this work, we calculate one of the indispensable ingredients in N -jettiness subtraction method, i.e., the soft function, for one massive colored particle production at hadron colliders. The NLO soft function is obtained to all orders in ǫ while the NNLO soft function is only calculated up to O(ǫ 0 ) for our purpose. We have used Mellin-Barnes integral representation and sector decomposition method to derive the above results. The structure of the soft functions is in coincidence with the expectation from analysing the scale independence of the cross section. We present the final results for the color octet in the main text and color triplet in the appendix, respectively. After that, we briefly discuss its application to the fixed-order calculations. For illustration, using this soft function along with other functions already known before, we reproduce the behavior of the NLO cross section for a massive color-octet vector boson production when the 0-jettiness variable τ is small. If the proper hard function is available, we can calculate the differential cross sections at higher orders. Our method can also be used to calculate the soft function for single top or top quark pair productions.
A Soft functions for a color triplet particle production In this appendix, we collect the renormalized soft functions for a direct top production from quark-gluon initial states via a flavor-changing neutral current. Since all the master integrals are the same as in the case of a color octet production, we omit the intermediate details and show the final results directly. The NLO soft function is s (1) = −(C A + C F )L 2 + 2C F L − 1 2 (C A + C F )π 2 + 8 ln 2 C F . (A.1) The NNLO soft function is where δ → 0 + is introduced to indicate the way of continuation. For p 2 i = m 2 i and p 2 j = 0, they are given by For p 2 j = m 2 j and p 2 i = 0, they are given by