Quark and gluon two-loop beam functions for leading-jet $p_T$ and slicing at NNLO

We compute the complete set of two-loop beam functions for the transverse momentum distribution of the leading jet produced in association with an arbitrary colour-singlet system. Our results constitute the last missing ingredient for the calculation of the jet-vetoed cross section at small veto scales at the next-to-next-to-leading order, as well as an important ingredient for its resummation to next-to-next-to-next-to-leading logarithmic order. Our calculation is performed in the soft-collinear effective theory framework with a suitable regularisation of the rapidity divergences occurring in the phase-space integrals. We discuss the occurrence of soft-collinear mixing terms that might violate the factorisation theorem, and demonstrate that they vanish at two loops in the exponential rapidity regularisation scheme when performing a multipole expansion of the measurement function. As in our recent computation of the two-loop soft function, we present the results as a Laurent expansion in the jet radius $R$. We provide analytic expressions for all flavour channels in $x$ space with the exception of a set of $R$-independent non-logarithmic terms that are given as numerical grids. We also perform a fully numerical calculation with exact $R$ dependence, and find that it agrees with our analytic expansion at the permyriad level or better. Our calculation allows us to define a next-to-next-to-leading order slicing method using the leading-jet $p_T$ as a slicing variable. As a check of our results, we carry out a calculation of the Higgs and $Z$ boson total production cross sections at the next-to-next-to-leading order in QCD.


Introduction
Precision measurements at hadron colliders often rely on jet vetoes to reduce the impact of background due to QCD radiation. This is done by rejecting events containing jets with a transverse momentum exceeding some cutoff value p veto T that is often much smaller than the large momentum transfer of the hard scattering process Q. This strategy finds common applications in the field of Higgs physics, as well as in a number of electro-weak and QCD measurements at the Large Hadron Collider (LHC), and has therefore motivated a large number of studies [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. For instance, the state-of-the-art prediction for the jet-vetoed Higgs production cross section involves the resummation of the large logarithms ln(p veto T /Q) up to next-to-next-to-leading logarithmic (NNLL) order [3][4][5][6] matched to fixed order calculations up to next-to-next-to-next-to-leading order (N 3 LO) [7]. The above NNLL calculations also include a numerical extraction of the O(α 2 s ) non-logarithmic terms relative to the Born (often referred to as NNLL ) from fixed-order calculations. In this article, we will complete the direct calculation of such constant terms in view of pushing the resummation accuracy to N 3 LL, which is demanded by the outstanding experimental precision foreseen at the LHC in the coming years.
Specifically, we consider the calculation of the two-loop beam functions entering the factorisation and resummation of the jet-vetoed cross section. Beam functions describe the dynamics of radiation collinear to the beam direction in high-energy hadron collisions. Together with our recent calculation of the two-loop soft function in Ref. [16], the results presented here constitute the last missing ingredient for the calculation of the jet-vetoed cross section at small value of the veto scale (i.e., up to power corrections in p veto T /Q) at the next-to-next-to-leading order (NNLO), and an important step towards the resummation of the leading-jet transverse momentum distribution (and the related jet-vetoed cross section) at N 3 LL order. Moreover, our results allow us to construct a slicing method for NNLO calculations in QCD using the leading-jet transverse momentum as a slicing variable.
We work in the framework of soft-collinear effective field theory (SCET) [17][18][19][20][21]. More specifically, the jet-veto cross section belongs to the class of SCET II problems, which are affected by the so-called factorisation (or collinear) anomaly [22,23], 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. Here we use the exponential regularisation scheme [24], consistently with our recent calculation of the two-loop soft function [16].
The validity of SCET factorisation for this observable at arbitrary logarithmic order has been the subject of debate in the literature [5,6,25]. Indeed, particular attention must be paid to the presence of soft-collinear mixing terms that might violate the factorisation theorem for this observable. While some groups [6,25] suggest that the SCET factorisation would be broken already at the next-to-next-to-leading logarithmic (NNLL) order, the authors of Ref. [5] present an argument as to why such terms should cancel if one performs a consistent multipole expansion of the measurement functions. Following the observation of Ref. [5], we explicitly show here that within the exponential regularisation scheme such a multipole expansion of the jet-clustering measurement function leads to the cancellation of the factorisation breaking terms, in the regime in which the jet radius R is treated as an O(1) parameter. 1 This explicitly confirms the validity of the factorisation theorem at NNLO, and constitutes an important step towards the resummation at the N 3 LL order.
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 we discuss the definition of the beam functions in the presence of a rapidity regulator. Section 3 contains a discussion of our analytic computation of the two-loop beam functions as a small-R expansion, as well as the details related to the zero-bin subtraction and the cancellation of the soft-collinear mixing terms. Section 4 reviews our numerical calculation of the beam functions that retains the full-R dependence, and compares the results obtained by the analytic and numerical computations, finding good agreement between the two. Finally, in section 5 we construct a phase-space slicing scheme based on the leading-jet transverse momentum for fully differential NNLO calculations of the production of colour-singlet systems. We test the scheme, and our results, by calculating the NNLO total cross section for Higgs and Z boson production at the LHC. Finally, our conclusions are given in section 6.

Factorisation of leading-jet transverse momentum in SCET
We begin by recalling the factorisation theorem for the jet-vetoed cross-section. We consider the production of an arbitrary colour-singlet system of total invariant mass Q in proton-proton collisions. The cross section, differential in the system's kinematics dΦ Born and with a veto on the transverse momentum of the leading jet p jet T < p veto T , factorises in the limit p veto Born is the Born amplitude for the production of the colour-singlet system, and µ and ν denote the renormalisation and rapidity scales, respectively. 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 from now on when referring to the ingredients of the factorisation theorem (2.1). The hard function H describes the dynamics at the hard scale, i.e. with virtuality µ ∼ Q. This scale is integrated-out in the SCET construction. Therefore, the hard function contains purely virtual contributions, and it is defined as the squared matching coefficient of the leading-power two-jet SCET current, i.e. (2. 2) The soft function S describes the dynamics of soft radiation off the initial-state partons and is defined as a matrix element of soft Wilson lines. Lastly, the beam functions are denoted by B n and Bn. Their two-loop calculation is the main focus of this paper. They are defined by matrix elements of collinear fields in SCET and describe the (anti-)collinear dynamics of radiation along the light-cone directions n µ andn µ of the beams. Within the SCET formalism, the resummation of the logarithms ln p veto T /Q appearing in Eq. (2.1) is achieved by evolving each of the functions in the factorisation theorem from their canonical scales to two common µ, ν scales. The hard matching coefficient obeys the renormalisation group equation (RGE) where Γ F cusp and γ F H are the cusp and hard anomalous dimensions of the quark (F = q) or gluon (F = g) form factors, renormalised in the MS scheme. The hard function's canonical scale is the hard scale µ = Q. The beam functions B i depend, in addition to the renormalisation scale µ, on the rapidity regularisation scale ν. They satisfy a system of coupled evolution equations (see e.g. Refs. [5,6]). The RGE is given by 4) and the rapidity evolution equation reads where γ F ν denotes the observable-dependent rapidity anomalous dimension. The boundary condition for the {µ, ν} evolution is set at the canonical scales µ = p veto T and ν = Q. Finally, the evolution equations for the soft function read with canonical scales µ = ν = p veto T . The dependence on the rapidity anomalous dimension cancels between the evolution of the soft and beam functions in the framework of the rapidity renormalisation group [26,27]. The soft (γ F S ) and collinear (γ F B ) anomalous dimensions are related to the hard anomalous dimension γ F H by the invariance of the physical cross section under a change of the renormalisation scale, that is The resummation of the jet-vetoed cross section at N k LL requires the cusp anomalous dimension Γ F cusp up to k + 1 loops, and the anomalous dimensions γ F H , γ F S , γ F B , γ F ν up to k loops. The boundary conditions (non-logarithmic terms) of the evolution equations need to be known up to k − 1 loops. Achieving N 3 LL accuracy for the jet-vetoed cross section requires the knowledge of the non-logarithmic terms in C(Q; µ) at two loops, which is given by the QCD on-shell form-factor [28] and has been known for a long time [29,30]. The two loop computation of the S function was presented in our earlier article [16]. Below, we focus on the evaluation of the two-loop beam functions. While the two-loop anomalous dimensions are known, the non-logarithmic terms are computed here for the first time. Moreover, since the anomalous dimensions featuring in Eq. (2.7) are known up to three loops [31][32][33], with the results presented in this work, the only missing ingredient for the N 3 LL computation is the three-loop rapidity anomalous dimension γ F ν .

The beam functions
The quark and gluon beam functions are defined as matrix elements of non-local collinear operators between the proton states |P (p) carrying large momentum p. Specifically, we have [2,5,6] where the collinear-gauge invariant collinear building blocks are defined in terms of the fields Here ξ(x) is the collinear quark field and D µ ⊥ is the covariant ⊥ collinear derivative. Gauge invariance is achieved by introducing the collinear Wilson line W (x) defined as (2.10) The operator M(p veto T , R 2 ) acts on a given state of X c collinear particles |X c by applying a veto on the final-state jets of radius R such that p jet i T < p veto with the phase space constraint M(p veto T , R 2 ) being Here max{p jet i T } is the transverse momentum of the hardest jet, where jets are defined in the E recombination scheme [34]. The constraint Θ cluster (R 2 ) is the generic clustering condition on the X c collinear final state particles, defined for a k T -class of jet algorithms [35] with jet distance measures Specific choices of the parameter p correspond to the anti-k T [35] algorithm (p = −1), the Cambridge-Aachen [36,37] algorithm (p = 0), and the k T [38] algorithm (p = 1). The results obtained in this article are valid for any of these choices. In Eq. (2.13), k ⊥i is 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 particles are clustered sequentially with respect to the above distance measure, as specified by the clustering condition Θ cluster (R 2 ). Since the definition (2.8) involves matrix elements between proton states, the beam functions are in general non-perturbative objects. However, for p veto T Λ QCD , it is possible to perturbatively match them onto the standard parton distribution functions (PDFs) f F/P (x, µ) as follows [2,5,6], where the standard proton PDFs are defined as for the quark and the gluon, respectively. The perturbative matching kernels I F F (z, Q, p veto T , R 2 ; µ, ν) are short-distance Wilson coefficients, whose computation is the focus of this article. At the LO, the collinear partons do not radiate and we have I F F (z, Q, p veto T , R 2 ; µ, ν) = δ F F δ(1 − z). The computation of radiative corrections to this relation up to the two-loop order will be the subject of Section 3.
The formal definition of the matching in Eq. (2.14) requires splitting the generic collinear fields into the perturbative collinear modes with the virtuality of the order of p veto T and the PDF-collinear modes with virtuality Λ QCD . The non-perturbative PDFs are then defined in terms of the matrix elements of the PDF-collinear modes only, while the perturbative collinear modes are integrated out from the theory. This distinction can often be ignored in practice at the leading power in SCET. Hence, in the rest of this article, we will refer to both types of collinear fields as simply collinear modes.
Even though Eq. (2.14) is written as an identity relating specific matrix elements, it represents an operatorial identity, which does not depend on the specific choice of the external states. Thus, to compute the matching kernels we can replace the external proton states by perturbative partonic states. With this replacement, the bare partonic PDF for finding a parton i inside parton j becomes which is valid to all orders in perturbation theory and the partonic matrix elements in Eq. (2.8) are then directly equal to the bare perturbative marching kernels. For this reason, in what follows, we will refer to the I F F coefficients interchangeably as matching coefficients or beam functions.

Analytic computation of the quark and gluon beam functions
In this section, we discuss the computation of the renormalised matching coefficients I F F (x, Q, p veto T , R 2 ; µ, ν), obtained after the renormalisation of the collinear PDFs and of the remaining UV divergences. The perturbative expansion of the matching kernels in powers of the strong coupling constant α s is defined as To efficiently compute the matching coefficients I F F (x, Q, p veto T , R 2 ; µ, ν), we decompose each beam function into the sum of the beam function for a specifically chosen reference observable and a remainder term. The reference observable is chosen so that it has the same single-emission limit. With this choice, the two-loop matching coefficients of the reference observable have the same divergences in the dimensional regularisation parameter around d = 4 as those for leading-jet p T , and the remainder term can then be computed directly in four dimensions. As in our calculation of the soft functions [16], we take the transverse momentum of the colour singlet system as the reference observable. The associated beam functions are denoted by I ⊥ F F (x, Q, p veto T ; µ, ν). These are known up to O(α 3 s ) [39][40][41][42][43][44] and we consider the renormalised O(α 2 s ) result of Refs. [41,42] as a reference in our calculation as these are also computed within the exponential rapidity regularisation scheme. The remainder term ∆I F F (x, Q, p veto T , R 2 ; µ, ν) accounts for the effects of the jet clustering algorithm. The perturbative matching coefficients are then expressed as The functions I ⊥ F F (x, Q, p veto T ; µ, ν) are obtained from the beam functions of the transversemomentum of the colour singlet system [41,42], and we include the one-loop and two-loop contributions in our ancillary files. 2 We stress here that the decomposition (3.2) is simply a convenient way of organising the calculation, and the ingredient I ⊥ F F has different physical properties from the actual beam functions entering transverse momentum resummation. A first difference is that the latter are defined in impact parameter space, since they are sensitive to the vectorial nature of transverse momentum factorisation. A second, related difference concerns the gluon beam functions, which in transverse momentum resummation receive a correction from different Lorentz structures including a linearly polarised contribution (see e.g. [39,42,45,46]). This leads to peculiar azimuthal correlations between radiation collinear to the two initial state (beam) legs. The above effect is absent in the gluon beam functions defined in Eq. (2.8), which are already integrated over the azimuthal angle of the emitted radiation. A simple physical explanation for this observation is that, unlike in the transverse momentum case, the jet algorithms considered in the factorisation theorem (2.1) will never cluster together emissions collinear to opposite incoming legs, therefore leaving no phase space for the azimuthal correlations to occur.
The term ∆I F F (x, Q, p veto T , R 2 ; µ, ν) defined in (3.2) contributes only when two or more real emissions are present and consequently ∆I F F starts at O(α 2 s ). At this order it can be computed directly in d = 4 space-time dimensions and only real emission diagrams contribute, with the measurement function To set-up the calculation, we need to take into account that the phase-space integrals, as in a typical SCET II problem, exhibit rapidity divergences which require additional regularisation. A consistent computation of the soft and beam functions requires that the same regulator is used in both calculations. We therefore adopt the exponential regulator defined in Ref. [24], which we also used in our recent calculation of the two-loop soft function [16]. This regularisation procedure is defined by altering the phase-space integration measure for real emissions, such that where ν is the rapidity regularisation scale discussed in section 2. The regularised beam functions are obtained by performing a Laurent expansion about ν → +∞ and neglecting terms of O(ν −1 ). In the rest of this section, we outline the main technical aspects of the calculation of the matching coefficients and present our results.
2 To be precise, we first perform the inverse Fourier transform of the beam functions obtained in Refs. [41,42], which are defined in impact parameter space. Then we integrate them up to p veto T to obtain

The renormalised one-loop beam functions
At O(α s ), the jet algorithm does not play a role and the correction term in Eq. (3.2) vanishes After the renormalisation of the collinear PDFs and of the remaining UV divergences, the one-loop result reads (see e.g. [41,42]) where the space-like splitting functions are

The renormalised two-loop beam functions
In this section we outline the procedure used in our computation of the two-loop beam functions. In particular, we do not discuss further the calculation of the already known 2), whose expression, after the renormalisation of the collinear PDFs and of the UV divergences, can be found in Refs. [41,42]. The following discussion refers to the correction ∆I F F for a generic flavour channel, although specific channels may present a simpler structure, and some of the contributions discussed below may vanish in their calculation.
The momenta of the two real particles (either gluons or quarks) are denoted by k i , with i = 1, 2, and we adopt the following parametrisation for the phase space on the r.h.s. of Eq. (3.4): in terms of the (pseudo-)rapidities η i , the azimuthal angles φ i , and the magnitudes of transverse momenta k i⊥ ≡ | k i⊥ |. We then perform a change of variables in the parametrisation of k 2 , in order to express its kinematics relative to that of k 1 . With this change of variables, the measurement function (3.3) takes the simple form where we used the explicit form of Θ cluster (R 2 ) in the variables defined above, namely the relation followed by the identity We now consider the squared amplitudes for the radiation of two collinear partons in a generic flavour channel |A F F | 2 , which have been derived in Refs. [47,48]. Without loss of generality, they can be expressed as is the contribution that survives in the limit in which the two emissions k 1 , k 2 are infinitely far in rapidity, while the remaining part A correlated F F encodes configurations in which the two emissions are close in rapidity (see also e.g. Refs. [49,50]). The above decomposition is useful in that the two contributions give rise to integrals with a different structure of rapidity divergences, and as such they require slightly different treatments. In the parametrisation (3.15) each contribution to the squared amplitude factorises as with x denoting the longitudinal momentum fraction. We parametrise the light-cone components in the delta function by means of Eq. (3.15) and the Sudakov parametrisation The large momentum component considered in the argument of the delta function of Eq. (3.21) (either + or −) depends on whether we consider the beam functions along then µ = (1, 0, 0, 1) or n µ = (1, 0, 0, −1), respectively. Without loss of generality, we will consider n µ as the hard-collinear direction, but the same considerations apply to the calculation of the beam functions alongn µ . In the following, we will discuss separately the treatment of the exponential regulator for the correlated and uncorrelated contributions to the beam function.
Rapidity regularisation for the correlated correction. The integrand of the correlated contribution vanishes by construction in the limit of large rapidity separation. Therefore, the only rapidity divergence in the integrals (3.21) arises when x = 1, namely when the rapidity of both emissions is unconstrained. We can handle the exponential regulator by integrating over η 1 using the delta function, and then expanding in distributions of (1 − x) as follows As usual, we kept only the leading power terms in the limit ν → ∞.
Rapidity regularisation for the uncorrelated correction. The integrand of the uncorrelated contribution does not vanish asymptotically in the regime of large rapidity separation. Therefore, it features two types of rapidity divergences, which emerge in the limits η → ±∞ and x = 1. Eq. (3.23) must be modified accordingly to deal with this more complicated structure, and should now involve the divergence in η as well. We proceed by expanding in distributions also in the variable w ≡ e η as follows 3 The integral over dw can only be performed after inserting the squared amplitude and the measurement function.
Laurent expansion in the jet radius. To proceed, in each of the contributions listed above, we consider the differential equation derived by taking the derivative of the integrals (3.21) in R. Since only the measurement function depends on R, this amounts to the replacement where R 2 > φ 2 . The resulting integral can be evaluated as a Laurent expansion in the jet radius R, that we obtain analytically in Mathematica with the help of the package PolyLogTools [51].
To calculate the boundary condition, we decompose the Θ(η 2 + φ 2 − R 2 ) function in ∆M(p veto T , R 2 ) given in Eq. (3.16) as The contribution stemming from part A contains the collinear singularity proportional to ln(R), while that arising from part B is regular in the R → 0 limit. The collinear singularity in part A does not directly allow us to take the boundary condition at R = 0. We then take an expansion by regions around R = R 0 1 and neglect terms of O(R 2 0 ). All boundary conditions are calculated analytically with the exception of the O(1) constant terms arising from part B of Eq. (3.26) for the correlated corrections, which are obtained numerically as a grid in the x variable. This is the only non-analytic ingredient in our calculation. We use the resulting boundary conditions to solve the differential equation in R, and afterwards we take the limit R 0 → 0. This procedure allows us to obtain the Laurent expansion to any order in R. In this article we present results up to and including O(R 8 ) terms, which are sufficient to reach very high precision in the numerical evaluation of the beam functions. Higher order terms in R could be in principle included in our expansion. 4

Zero-bin subtraction and cancellation of soft-collinear mixing at two loops
The structure of SCET reproduces that of an expansion by regions [52] of the relevant integrals occurring in the (real and virtual) radiative corrections to the observable under study. This method requires a full expansion of the integrals, in such a way that any expansion of each region within the scaling corresponding to a different region leads to scaleless integrals. 5 In practical applications of SCET, the presence of scales related to either additional regulators or the observable itself might render the overlap contributions non-vanishing. This can be overcome by subtracting by hand such overlapping regions to avoid doublecounting using the so-called zero-bin subtraction procedure [53]. In this sub-section we consider this procedure in the context of the factorisation theorem in Eq. (2.1), more precisely in its application to the matching coefficients I F F (x, Q, p veto T , R 2 ; µ, ν). Starting from the definition of the renormalised beam function given in Eq. (3.2), we observe that I ⊥ F F (x, Q, p veto T , R 2 ; µ, ν), which we extract from Refs. [41,42], already underwent zero-bin subtraction and thus does not contain any overlap between soft and collinear modes. It is thus sufficient to discuss the new contribution ∆I F F (x, Q, p veto T , R 2 ; µ, ν) computed in this article, which still contains contamination from soft modes due to the presence of additional scales such as the jet radius R and the rapidity regulator ν.
The starting point of the zero-bin procedure is the subtraction from ∆I F F of its own expansions when either one or both partons become soft. To be more precise, we introduce the operator SC which acts on ∆I F F by taking the expansion in the region in which emission k 1 is soft (k µ 1 ∼ (k + ∼ λ, k − ∼ λ, k ⊥ ∼ λ)) and k 2 is collinear (k µ 2 ∼ (λ 2 , 1, λ)).
This expansion affects both the squared amplitudes and the phase-space constraint in the ∆I F F integrals. Similarly, we introduce the operators CS and SS which, when acting on ∆I F F , perform the expansion of the beam function in the limit in which k 2 is soft and k 1 is collinear or soft, respectively. For the problem under consideration, we observe that the CS and SS operations commute, that is (CS)(SS) = (SS)(CS) (and similarly for SC). At two loops, we can then define the zero-bin subtracted beam functions as where the terms (1 − SS) are responsible for subtracting the soft-soft limit of the softcollinear subtraction. For a generic channel F F , some of the terms in Eq. (3.27) might vanish at leading power in the counting parameter λ.
In this procedure, a crucial role is played by the terms CS(1 − SS)∆I F F and SC(1 − SS)∆I F F , which describe an interplay between the soft and collinear modes. The only overlap between soft and collinear modes predicted by the SCET factorisation theorem (2.1) at a given perturbative order amounts to products of terms arising from the lower-order soft and beam functions. The absence of any other type of overlap is a necessary requirement for the observable under consideration to factorise and therefore to be resummable in SCET.
In the case at hand, the terms CS(1 − SS)∆I F F and SC(1 − SS)∆I F F in Eq. (3.27) are responsible for subtracting the overlap between soft and collinear regions. Performing the multipole expansion of the phase-space constraint is necessary to demonstrate that these terms have the expected form and thus that the factorisation theorem in Eq. (2.1) is indeed correct. This must be explicitly verified in the presence of the exponential regularisation scheme, since it introduces an additional scale that prevents integrals that would otherwise be scaleless from vanishing.
To be concrete, let us focus on the term CS(1 − SS)∆I F F on the r.h.s. of Eq. (3.27). To compute this term, we start by acting with the operator CS(1−SS) on the measurement function in Eq. (3.16), namely on We then rewrite The first term on the right-hand side of the above equation leads to a factorising contribution, in that the contribution associated with the theta func- (3.28) will cancel exactly against I ⊥ F F (x, Q, p veto T , R 2 ; µ, ν) when considering the full beam function I F F (x, Q, p veto T , R 2 ; µ, ν) as defined in Eq. (3.2). One is then left with the term (3.30) The trivial action of the CS(1 − SS) operator on this term amounts to simply replacing the transverse momenta k i⊥ with those of the collinear and soft particles. The resulting phase-space integral reduces to the product of the one loop beam and soft functions, in line with what is predicted by the factorisation theorem (2.1). Instead, the second term on the right-hand side of Eq. (3.29) seemingly leads to a term that is not captured by the factorisation theorem, featuring the measurement function To proceed, we act with the CS(1 − SS) operator on the above measurement function. The action of CS corresponds to taking the limit k µ 1 ∼ (λ 2 , 1, λ) and k µ 2 ∼ (λ, λ, λ). Following Ref. [5], one must then expand the clustering condition in Eq. (3.31). Noticing that |η| 1, this amounts to The phase-space constraints on the r.h.s. of the above equation lead to vanishing integrals as in this region |η| 1. We conclude that the terms CS(1−SS)∆I F F and SC(1 −SS)∆I F F vanish for the mixing configuration arising from the measurement function (3.31), in line with the prediction of the factorisation theorem (2.1). This result demonstrates that this factorisation is formally preserved at the two-loop order in the exponential regularisation scheme.
We then perform an explicit calculation of the remaining non-vanishing integrals entering the definition of the zero-bin subtraction (3.27). The calculation is performed using the same approach discussed in the previous sub-section, where now all the boundary conditions for the R 2 differential equation are evaluated fully analytically. The only technical difference with the calculation discussed in the previous section is the treatment of the exponential regulator, which is now modified by the fact that either one (for the softcollinear zero bin) or none (for the double-soft zero bin) of the momenta k 1 , k 2 appear in the longitudinal δ function in the integrals corresponding to Eq. (3.21). The analogues of Eqs. (3.23) and (3.24) for these contributions are given in Appendix A.
An alternative approach to the zero-bin subtraction procedure, adopted in Refs. [6,25], is to evaluate Eq. (3.27) using an alternative set of operators CS, SC and SS = SS, where the bar indicates that they do not act on the clustering condition Θ(η 2 + φ 2 − R 2 ) present in the measurement function defining ∆I F F , which is then left unexpanded. This leads to the following alternative definition of the zero-bin subtracted beam function ∆I Within this approach, the mixing terms originating from the integrals given by CS(1 − SS)∆I F F with the measurement function (3.31) do not vanish any longer. As such, to reproduce the QCD result, one has to add them back by hand to the factorisation theorem (2.1). We denote these terms by ∆I mix F F . Refs. [6,25] claim that such mixing terms constitute a violation of the SCET factorisation theorem (2.1) already at the NNLO (and NNLL) level, making it impossible to resum the jet-veto cross section at N 3 LL in SCET. Ref. [6] proposes to add these terms back by hand to Eq. (2.1) in order to achieve NNLL accuracy, but no fix is proposed beyond this order. One can however note that the soft-collinear mixing terms at the two-loop order have the same logarithmic structure as the zero-bin subtracted beam function ∆I F F which allows one, at this perturbative order, to absorb them into a re-definition of the subtracted two-loop beam functions as We verified by explicit calculation that As such, we conclude that the factorisation theorem (2.1) is preserved at NNLO and the mixing terms in Eq. (3.34) vanish upon performing the multipole expansion discussed above. The procedure leading to Eq. (3.34) can be used as a way to compute ∆I subtracted F F without performing the multipole expansion of the measurement function, but the apparent presence of mixing terms does not constitute a breaking of the SCET factorisation theorem at this perturbative order. For the interested reader, together with the results for the twoloop subtracted beam functions, we also provide in the ancillary files the expressions for the mixing terms of Eq. (3.34), ∆I mix F F , obtained without performing the multipole expansion discussed above.
We note that our findings are also consistent with the QCD formulation of the resummation of the jet-vetoed cross section [1,3]. Eq. (3.31) predicts a clustering between soft k µ ∼ (λ, λ, λ) and collinear k µ ∼ (λ 2 , 1, λ) modes which is absent in the QCD formulation due to the nature of the jet algorithms belonging to the generalised k t family. By construction, these do not cluster together partons that fly at very different rapidities, as it is the case for a soft and a collinear parton which feature a large rapidity separation |η| ∼ | ln 1/λ 2 | 1. The only possible clustering between a collinear and a soft parton in QCD is when the latter is also collinear, i.e. k µ ∼ (λ 2 , 1, λ) albeit with a small energy, which is entirely accounted for in the definition of the beam functions. Therefore, the absence of mixing terms in the SCET formulation is consistent with the QCD expectation.

Results and convergence of the small-R expansion
The two-loop zero-bin subtracted beam functions are included as Mathematica-readable files in the ancillary files accompanying this article. For each channel, we decompose the result into the different colour structures contributing at two loops. The final corrections to the (zero-bin subtracted) matching coefficients are obtained from their own colour decompositions as follows (we drop here the superscript subtracted used in the previous sub-section to simplify the notation): (3.36) The full beam functions are then obtained with Eq. (3.2). We stress that in our conventions the ∆I matching coefficients already contain a factor of n F , while the ∆I We now discuss some consistency checks on our results and on the validity of the small-R expansion for phenomenologically relevant values of the jet radius. As a first check, we verified that the dependence of the beam functions on ln ν matches the prediction from the evolution equation (2.5). As a second check, to assess the validity of our expansion in R we considered the quantity ∆I (2) F F (before performing the zero-bin subtraction discussed in Sec. 3.3) truncated at different orders in R 2 . More precisely, we defined the relative difference of the expansions at sixth and eighth order in R, and plotted the quantity for each different flavour channel (we set ν = p 2 T /p − = p 2 T /Q to remove the rapidity logarithms in ∆I (2) F F ). We find that in all cases δ F F (R 2 ) is vanishingly small up to R = 1, where the convergence of the R 2 expansion is not necessarily guaranteed. The convergence of the series is drastically improved at smaller values of R which are relevant for collider phenomenology. As an example, we plot in Fig. 1 the correlated corrections at R = 1. The plot shows an excellent convergence of the small-R expansion with residual corrections well below the permille level.

Numerical computation of the quark and gluon beam functions
We now discuss a numerical evaluation of the quark and gluon beam functions, which provides a crucial consistency check of the analytic results discussed in the previous section. In this section, we first outline the steps followed in the numerical calculation and then present a comparison between our numerical and analytic results.

Method for the numerical computation
All numerical integrations discussed below are performed using the GlobalAdaptive NIntegrate routine from Mathematica, to an accuracy at the permyriad level or better. This will allow for precise numerical tests of the analytic calculation. The correlated correction. For this part of the beam function, the integrand can only have a rapidity divergence as η t ≡ 1 2 (η 1 + η 2 ) → −∞ at finite η ≡ η 1 − η 2 , corresponding to x → 1. Note that the restriction k − 1 + k − 2 < p − forbids us from approaching η t → ∞ and encountering a rapidity divergence in this limit. We find it convenient to choose integration variables that remain finite as η t → −∞ (at finite k 2 1⊥ , k 2 2⊥ , η). In this way, x controls the approach to the rapidity divergence. Explicitly, we choose along with η and the azimuthal separation between the two emitted partons, φ. Using these variables, the rapidity divergences manifest themselves as a factor of 1/(1 − x) in the squared amplitude. This is regulated by inserting the exponential regulator factor (3.4), where in the exponent we may drop the k − ≡ k − 1 + k − 2 as we do not encounter any rapidity divergences associated with k − 1 , k − 2 → ∞. We take the limit ν → +∞ in the regulator and make use of the distributional expansion given in Eq. (3.5) of [41]: The structure of the integrand in T is simple, containing only terms of the form log n (T )/T , and integration over this variable may be done analytically. This just leaves the integrations over Z, η, and φ, which are performed numerically at fixed values of x and R.
The uncorrelated correction. After symmetrisation of the integrand in partons k 1 and k 2 , we may choose to integrate only over η < 0 and then multiply the final result by 2. With this restriction, for the uncorrelated contribution to the squared amplitude we encounter rapidity divergences as η 1 → −∞, as well as when η t → −∞. For our integration variables we should choose two variables that control the approaches to these two limits, and then other variables that remain finite in these limits. We choose to use: as well as x and φ. 6 Then, the limit η 1 → −∞ corresponds to Z → 0, whilst η t → −∞ corresponds to x → 1 and the rapidity divergences manifest themselves as a factor 1/[Z(1 − x)] in the integrand. We insert the exponential regulator (again, we can drop the k + in the exponent), and make use of the distributional expansion given in Eq. (3.30) of Ref. [41]:

(4.4)
We perform the integration over χ 1 analytically, and the integrations over the χ 2 , Z and φ variables numerically (for terms containing a δ(Z), we perform the trivial Z integration analytically).
The soft-collinear zero bins. We use the approach of Refs. [6,25] as a way to compute the zero-bin subtraction without performing the multipole expansion of the measurement function, see Sec. 3.3 and in particular Eq. (3.35). Let us, without loss of generality, take parton k 1 to be soft. Then k − 1 is no longer restricted by the delta function on the minus light-cone momentum, and we may have rapidity divergences for η 1 → ±∞ as well as for η 2 → −∞. We handle this calculation by re-expressing the clustering constraint in the measurement as: (4.5) In the second term on the right hand side, the two partons are restricted to be close together in rapidity, such that we only have rapidity divergences corresponding to η t → −∞ (or x → 1). The same strategy may then be used for this term as was used for the correlated corrections. Note that there is no collinear divergence here associated with η, φ → 0, due to the form of the squared amplitude for the soft-collinear zero bin (which coincides with the squared amplitude for the uncorrelated correction). For the first term on the right-hand side of Eq. (4.5), we choose to use the same variables we used in Ref. [16]: [ΔIFF'] Num.
[ppm]  along with η 1 , φ and x. The approach to η 2 → −∞ is then controlled by x, whilst η 1 directly controls the approach to η 1 → ±∞. We introduce the exponential regulator, and split it in a straightforward way into two factors depending on k 1 and k 2 respectively; we drop the k + 2 in the exponent as before, but now may no longer drop k + 1 . The integrand does not depend on η 1 (except in the regulator factor), and we may perform the integral over η 1 using Eq. (3.24) from Ref. [16]. We utilise Eq. (4.4) for the rapidity divergence corresponding to x → 1. The integration over K 2 T is performed analytically, and the z and φ integrals are done numerically.
The soft-soft zero bins. This calculation coincides exactly with that performed in Ref. [16] (up to a prefactor of δ(1 − x) that appears here), and we use the results presented in that paper for this contribution.

Comparison with the analytic results
As a further assessment of the quality of our analytic small-R expansion, we compare the numerical calculations (which have exact R dependence) with the analytic results obtained in the previous section at different values of R. Due to the many flavour channels, we choose to show here only the worst-case scenario, namely the comparison between the two calculations for the most complicated contributions, corresponding to the correlated part of the squared amplitudes in Eq. (3.19). Fig. 2 shows the outcome of this comparison at R = 1, and we can see that the difference between the two computations is at the level of parts per million, which is the level of accuracy of the numerical calculation. This demonstrates that the R expansion converges extremely well up to R = 1.

Leading-jet p T slicing at NNLO
The computation of the two-loop beam functions for the leading-jet p T constitutes the last missing ingredient to construct a non-local subtraction scheme for colour singlet production at NNLO based on p jet T . In analogy with non-local subtraction schemes such as q T -subtraction [54], jettiness subtraction [55,56], and k ness T subtraction [57], we can formulate an NNLO slicing fully differential in the Born phase space for the production of a colour singlet F , as (dσ ≡ dσ dΦ Born ) where the first term on the right hand side coincides with the non-logarithmic terms of the jet-veto cross section Eq. (2.1), the last term is its derivative with respect to p veto T expanded through O(α 2 s ) relative to the Born, and the second term is the NLO cross-section for the production of the colour singlet in association with a jet. The above formula formally reduces to the NNLO result in the limit p jet T,cut → 0. However, since the second and the third terms are both divergent logarithmically in this limit, Eq. (5.1) can be computed numerically only by choosing a finite value of p jet T,cut > 0. This introduces a slicing error O((p jet T,cut /Q) m ), where m is an integer value to be determined by studying the p jet T,cut → 0 behaviour of the non-singular contribution contained within brackets in Eq. (5.1).
The comparison of the NNLO results obtained using Eq. (5.1) to the known NNLO cross sections provides a very robust check of the correctness of the results presented in this work. We perform this test by considering on-shell Z and H production, which allows us to independently check the quark and gluon beam functions, respectively. We use MCFM 9.1 [58] to compute the NLO result for Z + j [59] and H + j [60][61][62] production, while we use the implementation of the jet-veto resummation [2,3,5,6] in the RadISH code [14,63,64] to compute the factorised expression (2.1) and its expansion up to NNLO. We compare our results with the analytic NNLO cross section for Z [65,66] and H [61,67,68] production which we computed using the n3loxs code. 7 For our numerical checks, we consider proton-proton collisions at a centre-of-mass energy of 13 TeV and R = 0.4. We adopt the LUXqed plus PDF4LHC15 nnlo 100 parton distribution functions [69] through the LHAPDF interface [70]. We choose factorisation and resummation scales equal to µ R = µ F = m Z , m H for Z and H production, respectively, with m Z = 91.1876 and m H = 125 GeV. In Fig. 3 we study the dependence of the NNLO correction on p jet T,cut /Q for Z and H production for different partonic channels by normalising it to the analytic result. We compare the results obtained using p jet T -slicing (in orange) with those obtained using q T -slicing (in blue) to assess the performance of the two methods. For Z production we are able to lower the value of p jet T,cut down to 0.1 GeV, whereas we stop at p jet T,cut = 0.5 GeV for Higgs production as the fixed order H + j calculation becomes slightly unstable in some channels below this  value. 8 We observe that in all the channels the results obtained using leading-jet p T slicing converge to the exact cross section in the p jet T,cut → 0 limit, thus providing a powerful check of the validity of our computations. By comparing the results obtained with p jet T -slicing to those obtained using q T -slicing we notice that the convergence towards the analytic result is comparable between the two methods, with q T -slicing converging slightly faster in most cases for R = 0.4. Smaller values of the jet radius R appear to improve the convergence of the p jet T subtraction, possibly due to the reduced size of the subleading power corrections. Further investigations on the size of subleading power corrections deserve dedicated studies.

Conclusions
In this article, we presented the first calculation of the complete set of two-loop beam functions relevant for the leading-jet transverse momentum resummation in colour singlet production. The results were obtained using two independent methods: a semi-analytical expansion for small jet-radius R up to and including terms of O(R 8 ), and a fully numerical evaluation for several fixed values of R. The small-R expansion is analytical with the only exception being a set of R-independent regular terms. The numerical calculation retains the complete R dependence and shows perfect agreement with the analytical expansion in the range R ∈ [0, 1] which is relevant for collider phenomenology. We further checked our computation by performing an NNLO calculation of the total cross section for Higgs and Z boson production using a slicing subtraction scheme based on the leading-jet p T . Our 8 We thank A. Huss for providing results calculated with the NNLOJET code [71] at p jet T,cut = 0.1 GeV for Higgs production, which we used as an independent cross-check.
calculation reproduces known analytic predictions for the NNLO total cross section in all flavour channels, thus validating our results.
When describing the technical aspects of the calculation, we discussed in detail the complications related to zero-bin subtraction and soft-collinear mixing. In particular, we explicitly showed that if one performs a multipole expansion of the measurement functions there exist no mixed soft-collinear contributions which break the SCET factorisation theorem at NNLO. This observation is non-trivial in the presence of the exponential rapidity regulator in that it adds a new scale to the problem, which leads to the presence of non-vanishing integrals that would otherwise be scaleless.
Our complete results are provided in Mathematica-readable files attached to the arXiv version of this article. Together with our earlier analytic results for the leading-jet p T soft function [16], they constitute an important step towards the N 3 LL resummation of this observable, with the only missing ingredient being the three-loop rapidity anomalous dimension.

Acknowledgments
We are grateful to Thomas Becher for helpful discussions on the cancellation of soft-collinear mixing terms in SCET factorisation. We would also like to thank Julien Baglio, Claude Duhr and Bernhard Mistlberger for providing us with a preliminary version of their computer code n3loxs used in our checks of the total cross section, and Alexander Huss for kindly providing a cross check of the differential distributions with the NNLOJet code. The work of JRG is supported by the Royal Society through Grant URF\R1\201500. LR has received funding from the Swiss National Science Foundation (SNF) under contract PZ00P2 201878. RS is supported by the United States Department of Energy under Grant Contract DE-SC0012704.
Note added In the final stages of the preparation of this article, Ref. [72] appeared with a numerical calculation of the beam functions in the quark channel. These results are obtained with a different rapidity regulator and computed for a discrete set of real points in the Mellin variable N conjugate to the longitudinal momentum x. For this reason, it is not immediately clear how to compare the results of Ref. [72] with the ones presented here.

A Expansion of the exponential regulator in zero-bin integrals
In this appendix we provide the ingredients to calculate the integrals contributing to the zero-bin subtraction discussed in Sec. 3.3. Specifically, we provide the analogues of Eqs. (3.23), (3.24) needed for the calculation of the correlated and uncorrelated contributions, respectively.
Soft-collinear zero-bin. We consider the limit in which one of the two partons is soft (say k 2 ) and the second is collinear. The exponential regulator in the correlated corrections can be expanded as: Similarly, we can use the following formula to deal with the uncorrelated contribution (see footnote 3): Analogous expressions hold for the case in which k 1 is soft.
Double-soft zero-bin. In the limit in which both partons are soft, the exponential regulator in the correlated corrections can be expanded as: