The analytic two-loop soft function for leading-jet $p_T$

We present the calculation of the two-loop soft function for the transverse momentum distribution of the leading jet produced in association with any colour-singlet system (e.g.~a Higgs or a $Z$ boson). This constitutes a central ingredient for the resummation of the above distribution as well as the jet-vetoed cross section at the next-to-next-to-next-to-leading logarithmic order, both of which play an important role in the precision physics programme at the Large Hadron Collider. The calculation is performed in soft-collinear effective theory with an appropriate regularisation of the rapidity divergences that occur in the phase-space integrals. We obtain analytic results by employing an exponential regulator and by taking a Laurent expansion in the jet radius $R$. All expressions are presented as ancillary files with this article.


Introduction
Precise tests of the Standard Model at the LHC involving the production of colourless final states require accurate control over QCD radiation. Jet vetoes are commonly used to reduce the unwanted QCD background by rejecting, or vetoing, events containing an energetic jet, i.e., with transverse momentum exceeding some cut-off value p veto T . The presence of the cut-off, typically much smaller than the hard scale for the process Q, leads to large logarithmic terms of the form α n s ln k (p veto T /Q), with k ≤ 2n, in the perturbative series of the jet-vetoed cross section. Such corrections can spoil the convergence of the perturbative expansion and thus require all-order resummation. The first resummation of such logarithms was performed in Ref. [1] at the next-to-leading logarithmic (NLL) accuracy, and subsequently at next-to-next-to-leading logarithmic (NNLL) accuracy in Ref. [2] using a direct calculation in QCD, as well as in Refs. [3,4] in the framework of softcollinear effective field theory (SCET) [5][6][7][8][9] by means of renormalisation group methods. The above NNLL calculations also include a numerical extraction of the O(α 2 s ) constant terms relative to the Born (often referred to as NNLL ), whose analytic calculation will be the focus of the study initiated by the present paper.
A great variety of phenomenological applications for various colour-singlet processes at the LHC has followed these results. In the important area of Higgs physics, state of the art predictions for the jet-vetoed (zero-jet) cross section based on the combination of NNLL resummation with a next-to-next-to-next-to-leading-order (N 3 LO) fixed-order calculation have been presented in Ref. [10], including an account of heavy-quark-mass effects [11] and the resummation of logarithms of the jet radius of the type α n s ln n R 2 following Ref. [12].
Similarly, predictions for the exclusive one-jet cross section have been presented in Ref. [13], and combined with the zero-jet bin in Ref. [14]. More recently, further technology has been developed for the calculation of differential distributions in the zero-jet bin [15], as well as to incorporate the effect of rapidity cuts on jets [16]. An alternative definition of jet vetoes encoding a dependence on the rapidity of the vetoed jets has been presented and studied in Refs. [17][18][19]. Finally, applications to other colour-singlet production processes, such as the production of multi-boson electroweak final states, have been carried out in Refs. [20][21][22][23][24][25].
Going beyond NNLL accuracy is crucial for keeping up with the foreseen performance of the LHC experiments. Previous considerations beyond NNLL were made in Ref. [3]. A numerical calculation of the two-loop soft function was presented in Ref. [26], albeit using a different regularisation scheme to what we consider here. A full N 3 LL calculation remains, however, unavailable. This paper initiates a systematic study of the higher-order terms for arbitrary colour-singlet processes. The separation of scales achieved within SCET allows the formulation of a factorisation theorem that we take as a starting point. The jet-veto cross section is obtained as the integral of the leading-jet transverse momentum distribution, which belongs to the class of SCET II problems. These are affected by the so-called factorisation (or collinear) anomaly [27,28], which is connected to the presence of rapidity divergences in the ingredients of the factorisation theorem. Such divergences are not regulated by the standard dimensional regularisation scheme and therefore an additional (rapidity) regulator must be introduced. The introduction of this regulator leads to a more involved structure of the renormalisation group equations that govern the evolution of the factorisation's building blocks. Various approaches have been proposed in the literature to restore factorisation and achieve a consistent resummation. In this article, we adopt the rapidity renormalisation group method [29,30], which relies on the introduction of a new scale ν, whose renormalisation flow allows one to consistently resum the logarithms associated with rapidity divergences.
This article presents the first analytic two-loop computation of one of the critical ingredients for N 3 LL resummation of the logarithms ln(p veto T /Q), namely the soft function. The article is organised as follows: In section 2, we review the factorisation theorem for the production of a colour-singlet with a veto on the transverse momentum of the leading jet and discuss the definition of the soft function in the presence of a rapidity regulator. Section 3 discusses our computation of the two-loop soft function. The results and various checks are presented in section 4, and the conclusions are given in section 5.

Factorisation of leading-jet transverse momentum in SCET
In this section we describe the formalism for the jet-vetoed cross section in colour-singlet production, which is defined as the integral of the leading-jet transverse momentum distribution up to a cut p veto T . As mentioned in the introduction, we adopt the formalism of the rapidity renormalisation group [29,30] to deal with the rapidity divergences, as opposed to the original treatment of Ref. [28] in which the exponentiation of the collinear anomaly follows from the consistency of the EFT formulation.
Our starting point is the factorisation theorem for this observable in SCET [3,4,31]. We consider the production of a generic colour-singlet system of squared invariant mass Q 2 , with a veto on the transverse momentum of the leading jet, p jet T < p veto T . In the limit p veto T Q, the cross section differential in the kinematics of the colour-singlet system is given by where dΦ Born denotes the full Born phase-space measure including the flux factor, A F Born is the Born amplitude for the production of the colour-singlet system, and µ and ν denote the renormalisation and rapidity regularisation scales, respectively. All quantities on the second line of Eq. (2.1) depend explicitly on the jet radius R. The index F indicates the flavour configuration of the initial state, i.e. either qq (F = q) or gg (F = g), and for simplicity we will drop it when referring to the ingredients of the factorisation theorem (2.1). In Eq. (2.1), the hard function H describes the dynamics at large virtuality scales µ ∼ Q, and therefore contains purely virtual contributions. It is defined as the square of the matching coefficient between QCD and SCET, i.e., (2. 2) The two beam functions B n and Bn describe collinear dynamics of radiation along the lightcone directions n µ andn µ , which correspond to the beam directions. These are defined by matrix elements of collinear operators in SCET. Finally, the soft function S describes the dynamics of soft radiation emitted off the two initial-state partons. The resummation of the logarithms ln Q/p veto T is achieved by solving evolution equations for each of the functions in the factorisation theorem (2.1) between their canonical scales and the two scales µ and ν. Specifically, the hard function satisfies the Renormalisation Group Equation (RGE) related to that of the matching coefficient C(Q; µ) (see e.g. Refs. [32,33] and references therein) where Γ F cusp and γ F H are the cusp and hard anomalous dimensions of the quark (F = q) or gluon (F = g) form factor in the MS scheme. The boundary condition of the evolution is set at the canonical scale µ = Q. Unlike the hard function, the beam functions B i also depend on the rapidity regularisation scale ν, and satisfy a system of evolution equations (see e.g. Refs. [3,4]). The first is the RGE where γ F B is the collinear anomalous dimension, and the second is the rapidity evolution equation where γ F ν denotes the observable-dependent rapidity anomalous dimension. The boundary condition for the joint {µ, ν} evolution is set at the canonical scales µ = p veto T and ν = Q. Finally, the soft function S satisfies the system of evolution equations with canonical scales µ = ν = p veto T . In the framework of the rapidity renormalisation group, the dependence on the rapidity anomalous dimension cancels between the soft and beam functions, while the soft (γ F S ) and collinear (γ F B ) anomalous dimensions are related to the hard anomalous dimension γ F H by the RG invariance of the physical cross section, that is The resummation of the jet-vetoed cross section at N k LL will require the cusp anomalous dimension Γ F cusp up to k + 1 loops, and the quantities γ F H , γ F S , γ F B , γ F ν up to k loops. Finally, the boundary conditions of the RGEs introduced in this section are needed up to k − 1 loops. In this article, we focus on the calculation of such boundary conditions for the soft function up to the two-loop order, which, together with the two-loop beam functions and the three-loop rapidity anomalous dimension γ F ν , are the missing ingredients to formulate the N 3 LL computation.
Before we carry on with the computation of the soft function, a remark is in order. The validity of the factorisation theorem in Eq. (2.1) at arbitrary logarithmic order has been the subject of debate in the SCET literature. Specifically, while its validity up to NLL is uncontroversial, going to higher logarithmic orders requires some extra care. The authors of Refs. [4,34] argue that, starting from NNLL, the existence of soft-collinear mixing terms may break the factorisation of the jet-veto cross section into soft and beam functions, while the authors of Ref. [3] disagree with this statement and argue that the cancellation of such terms requires the multipole expansion of the phase-space constraint in the computation. In this article, we do not discuss this matter any further as it is not directly relevant for the results presented here. We will present a discussion of such soft-collinear mixing terms in a future publication.

The soft function
In this section we focus on the soft function, whose operatorial definition reads: and analogously for Yn. The colour charge operator T a is understood to be in the adjoint (fundamental) representation of SU (3) for gluon (quark) initiated processes. The trace in Eq. (3.1) is performed over colour indices and the normalisation factor d F denotes the dimension of the representation, i.e. d F = N 2 C − 1 for gluon-initiated reactions and d F = N C for quark-initiated ones. The operator M(p veto T , R 2 ) acts on a given state of X s soft particles |X s by applying a veto p veto T on the final-state jets of jet radius R, obtained with a generalised k T family of clustering algorithms.
To proceed we introduce a complete set of soft states in Eq. (3.1) and recast it as where M(p veto T , R 2 ) is the phase-space constraint, which schematically reads is the generic clustering condition on the X s soft particles in the final state. This is defined algorithmically for the generalised-k T family of jet algorithms with distance measures with p = −1 for the anti-k T [35], p = 0 for the Cambridge-Aachen [36,37], and p = 1 for the k T [38] algorithm. Here k i⊥ denotes the transverse momentum of particle i with respect to the beam direction, and ∆η ij and ∆φ ij are the relative rapidity and azimuthal angle between particles i and j, respectively. The clustering condition Θ cluster (R 2 ) specifies how particles are sequentially clustered according to the above distance measures. Its expression will be given below when discussing the calculation of the soft function. The computation can be decomposed into the sum of the soft function for a reference observable and a correction factor. A similar approach for the calculation of two-loop soft functions has been exploited also in Refs. [18,39]. Here, for the reference observable we take the transverse momentum of the colour singlet system, which we denote by S ⊥ (p veto T , µ, ν) (defined in Ref. [40] for the regularisation scheme adopted here). The correction factor ∆S(p veto T , R 2 ; µ, ν) accounts for the effect of the clustering algorithm. The soft function S ⊥ (p veto T , µ, ν) is well understood and known up to O(α 3 s ) [40,41]. We recompute it up to O(α 2 s ) in this article to validate our approach. The jet-veto soft function is then: This equation defines the remainder term ∆S(p veto T , R 2 ; µ, ν), whose evaluation at the twoloop level is the main result of this article.
The remainder function ∆S depends on the jet algorithm and only contributes for two or more real emissions. This means that the difference ∆S starts at O(α 2 s ) and at this accuracy is purely determined by double-real diagrams. The decomposition in Eq. (3.6), together with the definition (3.4) allows us to determine the measurement function for the remainder term ∆S as (3.7)

Regularisation of rapidity divergences
We now discuss the regularisation procedure used to deal with the rapidity divergences that feature in SCET II problems. We adopt the exponential regulator initially proposed in Ref. [40], which is defined by modifying the phase-space integration measure for each real emission, such that where ν denotes the rapidity regularisation scale discussed in section 2. On the other hand, the integration measure for virtual corrections remains unchanged. The role of the regulator is to suppress the integrand as the light-cone components become large, hence regulating the phase-space region affected by rapidity divergences. The regularised soft function is obtained by taking the Laurent expansion around ν → +∞ after the phase-space integrals over k i have been evaluated. In the expansion, one neglects terms of O(ν −2 ). The exponential regulator introduced above has several useful features. In particular, it preserves non-abelian exponentiation, and allows for the resummation of rapidity logarithms by means of evolution equations given in section 2, whose solution is independent of the path chosen in the (µ, ν) plane. In the following sections we present the results for the one-loop and two-loop soft functions, and discuss how the above regulator is dealt with in the explicit computation.

The renormalised one-loop soft function
We define the perturbative expansion of the soft function as The leading order result is S (0) (p veto T , R 2 ; µ, ν) = S implying that In fact, the above considerations trivially hold true for all single-emission contributions to all orders in α s . By explicit computation, and performing the renormalisation in the MS scheme, one finds the result (3.12)

The renormalised two-loop soft function
At two loops, we start by computing S ⊥ (p T , µ, ν) using an analytic implementation of sector decomposition and the Mathematica packages HypExp [42] and PolylogTools [43]. Our results agree with the findings of Refs. [41,44], which provides a first consistency check on our computational approach. For completeness, we include them in the ancillary file SoftFunctionPerp.m. At this order we also get the first non-zero contribution to ∆S(p T , R 2 ; µ, ν), coming from double-real diagrams describing the emission of either two soft gluons or a soft quark-antiquark pair. The corresponding double-real squared amplitudes were obtained in Refs. [45][46][47] (see also Refs. [48][49][50] for calculations in the context of SCET). At two-loops, the µ dependence in ∆S(p T , R 2 ; µ, ν) is entirely encoded in the strong coupling constant α s (µ), while the ν dependence is induced by the exponential regulator of which we will eventually consider the expansion around ν → +∞.
It is convenient to split the calculation of ∆S (2) (p T , R 2 ; µ, ν) into the sum of a correlated and an uncorrelated term (see also e.g. Refs. [18,51]). We write which follows from a similar decomposition of the squared amplitudes, (3.14) The uncorrelated term A uncor. (k 1 , k 2 ) encodes the limit of the double-soft squared amplitude |A(k 1 , k 2 )| 2 for relative rapidity η ≡ ∆η 12 = η 1 − η 2 → ∞. This term is proportional to the colour factor C 2 R , where C R = C F or C R = C A for quark or gluon initiated processes. Conversely, the correlated term encodes the remaining contributions to |A(k 1 , k 2 )| 2 , which vanish in the limit where the two emissions are widely separated in rapidity. This term can be further separated into two colour factors, C R C A and C R T R n F : The utility of the decomposition in Eq. (3.13) lies in the fact that the two terms on the right-hand side require different treatments of the rapidity regulator. We have performed two independent computations for ∆S (2) (p T , R 2 ; µ, ν), one analytical and one numerical, which we discuss in more detail below. Before doing so, however, we note that ∆S (2) (p T , R 2 ; µ, ν) is finite in all infrared and collinear limits of the two-particle phase space (because the r.h.s. of Eq. (3.7) vanishes in these limits), and therefore the calculation can be directly performed in d = 4 space-time dimensions. In particular, this implies that ∆S (2) (p T , R 2 ; µ, ν) does not depend on µ.

Analytic computation of the two-loop soft function
We start by discussing the analytic calculation of ∆S (2) (p T , R 2 ; µ, ν). The momenta of the two real particles (either gluons or quarks) are denoted by k i , and we adopt the following parametrisation for the phase space in the r.h.s. of Eq. (3.8): in terms of the (pseudo-)rapidities η i , the azimuthal angles φ i , and the transverse momenta k i⊥ ≡ | k i⊥ |. We then perform a change of variable in the parametrisation of k 2 , (3.17) in order to express its kinematics relatively to the ones of k 1 . With this change of variable, the measurement function (3.7) takes the simple form where we used the explicit form of Θ cluster (R 2 ) in the variables defined above, namely the relation With the above parametrisation, the correlated and uncorrelated contributions to the squared amplitude factorise as (see e.g. Refs. [1,51]) where the function D cor./uncor. (ζ, η, φ) is regular in the limit ζ → 0. In particular one has The calculation of ∆S (2) (p T , R 2 ; µ, ν) will then involve the evaluation of phase-space integrals of the type with ∆M(p veto T , R 2 ) being defined in Eq. (3.18), and where we have already performed the integration over φ 1 . We now discuss separately the correlated and uncorrelated contributions.
The correlated correction. We start with the correlated corrections. The function D cor. (ζ, η, φ) vanishes in the limits η → ±∞ by definition of this contribution, and therefore the only source of rapidity divergences in Eq. (3.23) are the limits η 1 → ±∞. We can then integrate over η 1 and expand the result around ν → ∞ and, according to our regularisation procedure, retain only the leading term. We write We then decompose the Θ(η 2 + φ 2 − R 2 ) function in ∆M(p veto T , R 2 ) given in Eq. (3.18) as (3.26) The integral I A associated with part A is simple to evaluate by direct integration using the package PolyLogTools [43]. We obtained the exact R 2 dependence and subsequently expanded the final result in a Laurent series about R 2 = 0.
To compute the integral corresponding to part B, we consider its derivative with respect to the (squared) jet radius R 2 . Equation (3.25) leads to the differential equation We compute the right-hand side of the above equation as a Laurent series in R 2 , up to and including R 8 terms. As a boundary condition, we take the limit of I B for R = 0, which is allowed by the fact that the entire collinear divergence in Eq. (3.26) is encoded in part A. Some care is needed in applying the method of expansion by regions [52] when taking this limit. Starting from O(R 2 ) more than one region contributes which requires the introduction of additional regulators. The calculation of the O(1) term, however, only receives contributions from a single region and the integral is well defined. We apply the above strategy to both colour structures contributing to ∆S cor. (p T , R 2 ; µ, ν). The resulting expressions are provided in electronic form in the ancillary files SoftFunctionFermion.m and SoftFunctionNonAbelian.m for C R n F T R and C R C A , respectively.
The uncorrelated correction. We now consider the uncorrelated case. Thanks to Eq. (3.22), the integrand takes a simpler form than in the correlated case: Despite its apparent simplicity, this integral has a more complicated structure of rapidity divergences, which now originate both from the η 1 → ±∞ and η → ±∞ limits. To handle this situation, we start by taking a Laplace transform of the exponential regulator with respect to the variable e −γ E /ν and denote its conjugate by τ . Introducing the more convenient variables the Laplace transform gives (including the Jacobian corresponding to the above change of variable) We then perform an expansion of the exponential regulator in distributions in both x and w in Laplace space, where the exponential regulator becomes the rational function in Eq. (3.30). Finally, we take the leading term in the limit of τ → ∞ (corresponding to ν → ∞) and take the inverse Laplace transform of the result. This procedure yields where we also evaluated the x integral using the fact that the remaining factors in the integrand of Eq. (3.28) are x independent. For the remaining integrals, we follow the same procedure as for the calculation of the correlated term ∆S cor. (p T , R 2 ; µ, ν). We report the resulting expression in the ancillary file SoftFunctionAbelian.m.

Numerical computation of the two-loop soft function
In this section we briefly outline the procedure used for the numerical evaluation of the quantity ∆S (2) , retaining the full dependence on the jet radius. This independent calculation will provide a non-trivial check of the analytic calculation discussed in the previous section. Numerical calculations were performed to a precision at the permille level.
The correlated correction. We express the integral for ∆S (2) in terms of the variables , (3.32) in addition to φ, as defined in Eq. (3.17). For the correlated contribution, the integrand is strongly peaked around η = 0, meaning that no rapidity regulator is needed for this integration. The integrand is however flat in η t , leading to a rapidity divergence when integrating over this variable. The structure of the integral over η t , including the exponential regulator factor, is exactly that of Eq. (3.24), and we use this equation to perform the integral analytically (in the ν → ∞ limit). The structure of the resulting integrand in the variable K 2 T is very simple, containing only terms of the form ln n K 2 T /ν 2 /K 2 T , and we can also perform the integration over K 2 T analytically. Finally, we are left with (finite) integrations over the η, φ, and z variables. These are performed numerically using the GlobalAdaptive NIntegrate routine from Mathematica.
The uncorrelated correction. Here we re-express the constraint in the measurement function for ∆S (2) as (3.33) For the second term on the right-hand side, there is guaranteed to be no rapidity divergence in the integration over η by definition, and we may use exactly the same techniques as were used for the correlated correction. For the first term, we use the integration variables η 1 , η 2 , K 2 T , z and φ. The integrand does not depend on either η 1 or η 2 , such that we experience rapidity divergences in the integrations over both of these variables. After inserting the exponential regulator, we may perform the integrations over both η 1 and η 2 analytically by making use of the result in Eq. (3.24) (where here we only need the result for the case where ζ = 0). As for the correlated case, the structure of the integral in K 2 T is simple and may also be performed analytically. This then leaves us with the z and φ integrations, which are performed numerically via the GlobalAdaptive NIntegrate routine from Mathematica.

Results and dependence on the jet radius
We now discuss some consistency checks of our results and the validity of the small-R expansion for phenomenologically relevant values of the jet radius.
We start by extracting the rapidity anomalous dimensions, which can be compared with the results of Refs. [2][3][4]. Specifically, starting from the results for ∆S (2) (p T , R 2 ; µ, ν) computed in the previous section we can obtain the difference between the rapidity anomalous dimension entering the resummation of p veto T and that of p T . At O(α 2 s ) this amounts to computing:

∂∆S
(2) −131 + 12 π 2 + 132 ln(2) ln(R) + 8 27 108ζ 3 − 805 + 33 π 2 + 396 ln 2 (2) + 420 ln (2) + 1429 + 3600 π 2 + 12480 ln(2) 2700 where C R is the quadratic Casimir in the representation of the initial-state partons (either quarks or gluons). Equations (4.1) and (4.2) are obtained from the correlated corrections, while Eq. (4.3) is obtained from the uncorrelated corrections. For the C 2 R colour factor we present the full R 2 dependence, since the series terminates at order R 4 for the anomalous  dimension if R < π (see the discussion in the appendix of Ref. [1] for R > π). For the other two colour factors, higher-order terms in R 2 are not given here for simplicity but can be obtained from the results in the ancillary files. These results reproduce those given in Refs. [2][3][4]. A numerical calculation of the same soft function defined in Eq. (3.1) has been previously presented in Ref. [26] using the SoftServe code [53]. This computation uses a different rapidity regularisation procedure to ours. As a consequence, while we find agreement at the level of the anomalous dimensions, the boundary conditions to Eqs. (2.6) are scheme dependent and thus differ from the result of Ref. [26]. We next assess the validity of the small-R expansion used in the calculation of the soft function, and specifically whether this expansion is sufficiently accurate for typical values of the jet radius R ∈ [0.2, 0.8]. We first compare the analytic expansion in R 2 to the numerical calculation of the soft function, whose R 2 dependence is exact. Fig. 1 displays the R dependence of the two results for ν = p T , which agree well in the considered range of R. As a second check, we consider the quantity ∆S (2) (p T , R 2 ; µ, ν) truncated at different orders in R 2 for the three different colour structures. More precisely, we define the relative difference of the expansions at sixth and eighth order in R (for ν = p T ), and plot the quantity for the three different colour factors, C ∈ {C R n F T R , C R C A , C 2 R } in Fig. 2. The plot shows an excellent convergence of the small-R expansion all the way up to R = 1.0, with residual corrections being at the sub-permille level. As a benchmark, in Table 1 we present the values of δ C (R 2 ) for R = 0.4 and R = 0.8, showing that O(R 8 ) terms are indeed numerically negligible. The two figures 1-2 indicate that the truncation error associated with our small-R expansion is well under control for phenomenologically relevant values of the jet radius. We note nevertheless that the same strategy that was used here to obtain the results up to O(R 8 ) allows us to compute higher-order terms should they be required.

Conclusions
In this article, we presented the first complete analytic calculation of the two-loop soft function which enters the factorisation theorem for the leading-jet transverse momentum resummation, or equivalently the jet-vetoed cross section for the production of any colour singlet system. We carry out two independent calculations: an analytic expansion for small values of the jet radius R up to and including O(R 8 ) terms; and a numerical calculation for selected values of R where the full dependence on R is retained. The two calculations ingredients which are currently unknown, one needs the two-loop beam functions as well as the three-loop rapidity anomalous dimension. Moreover, going beyond NNLL requires also a careful formulation of the factorisation theorem in SCET. We will address the above points in forthcoming publications.