Two-Loop Beam and Soft Functions for Rapidity-Dependent Jet Vetoes

Jet vetoes play an important role in many analyses at the LHC. Traditionally, jet vetoes have been imposed using a restriction on the transverse momentum pTj of jets. Alternatively, one can also consider jet observables for which pTj is weighted by a smooth function of the jet rapidity yj that vanishes as |yj | → ∞. Such observables are useful as they provide a natural way to impose a tight veto on central jets but a looser one at forward rapidities. We consider two such rapidity-dependent jet veto observables, TBj and TCj , and compute the required beam and dijet soft functions for the jet-vetoed colorsinglet production cross section at two loops. At this order, clustering effects from the jet algorithm become important. The dominant contributions are computed fully analytically while corrections that are subleading in the limit of small jet radii are expressed in terms of finite numerical integrals. Our results enable the full NNLL′ resummation and are an important step towards N3LL resummation for cross sections with a TBj or TCj jet veto.


Introduction
Jet vetoes find frequent application at the LHC, e.g. in Higgs property measurements as well as in searches for physics beyond the Standard Model. They are used to cut away backgrounds, and more generally to classify the data into exclusive categories, or 'jet bins', based on the number of jets in the final state, in order to increase the signal sensitivity. When a generic jet veto observable T is constrained to be much smaller than the hard scale of the process Q, large logarithms of T /Q appear in the perturbative expansion of the jet-vetoed cross section, and should be resummed to obtain precise predictions. The default jet variable by which jets are currently classified and vetoed is the transverse momentum p T j of a jet. However, there are some drawbacks to using this variable. In harsh pile-up conditions, it can be difficult to identify (and veto) low-p T jets in the forward region (beyond |η| 2.5) when a large part or all of the jet lies in a detector region where no tracking information is available. One way to get around this problem is by introducing a hard cut on the jet (pseudo)rapidity and only consider jets at |η j | < η cut . However, such a cut also changes the logarithmic structure [1], and none of the extant resummations of p T j take account of such a rapidity cut. Another way one might try to avoid the problem is to raise the cut on p T j , but then one loses the potential benefits of a tight jet veto (such as its utility to identify the initial state of heavy resonances [2]).
An alternative approach is to consider a generalized jet-veto variable [3] T f j = p T j f (y j ) (1.1) that is smoothly dependent on the jet rapidity y j , where the function f (y) is decreasing for increasing |y| such that limiting T f j to small values tightly constrains central jets, but only loosely constrains the forward ones. In addition to the above practical considerations, given the importance of jet binning it is highly desirable to have several different options for performing jet vetoes, both experimentally and theoretically. This avoids having to rely exclusively on p T j , and provides important complementary information on the pattern of additional jets produced, for example in Higgs production [4].
In this paper, we consider the two representative jet observables 1 (1.2) We have introduced light-cone coordinates, where an arbitrary four-vector q µ is decomposed as q µ = q − n µ /2+q +nµ /2+q µ ⊥ with n µ andn µ being light-like vectors (n 2 =n 2 = 0, n·n = 2) along the beam directions. T Bj gives the plus (minus) momentum of the jet j if the jet lies in the right (left) hemisphere with p −(+) j > p . That is, it has the same rapidity weighting as the global beam thrust hadronic event shape [5]. T Cj has the same rapidity weighting as the C-parameter event-shape for e + e − → dijet processes. It becomes equal to T Bj at forward rapidities and approaches p T j /2 at central rapidities. The T Cj spectrum has been measured in Higgs production by ATLAS [4].
The purpose of the present paper is to provide the full two-loop corrections to the T Bj -dependent and T Cj -dependent beam and soft functions that are required to bring the resummation for T Bj and T Cj to the NNLL level. These corrections also provide the necessary fixed-order boundary conditions for the N 3 LL resummation, with the remaining missing ingredients being the three-loop clustering correction to the noncusp and the fourloop correction to the cusp anomalous dimensions.
We compute, partly analytically and partly numerically, the beam and dijet soft functions for both rapidity-dependent jet veto observables in eq. (1.2), for all possible color and parton channels. To be precise, for the beam functions we compute the two-loop perturbative matching coefficients between the beam functions and the standard parton distribution functions (PDFs). This is the first explicit calculation of the full set of two-loop singular matching corrections for a jet-algorithm dependent jet veto. This includes the complete set of corrections arising from the clustering of two independent emissions. (In ref. [19], the full two-loop soft function for p T j was calculated, while the two-loop beam functions required for the NNLL resummation was extracted numerically from fixed-order codes.) Our results for the two-loop beam functions can also be used in computations of N -jet cross sections with a veto on further jets being imposed via T Bj or T Cj . This paper is structured as follows: In section 2, we define precisely the soft and beam functions for T Bj and T Cj vetoes, and give an outline of their calculation at two loops. Additional technical details are relegated to the appendices. In section 3, we give the obtained results for the two-loop beam and soft functions, and we conclude in section 4.

Calculation
The measurement function for a jet veto using a generic variable T f j is given by θ(T f j < T cut ) . (2.1) This constraints all jets to have T f j < T cut and thus vetoes any jets with T f j > T cut . The jets J(R) in eq. (2.2) are identified with a specific jet algorithm with jet radius R. At NNLO, we can have at most two emissions, and the precise form of the jet algorithm is irrelevant. Our results are valid for any jet algorithm that clusters the emissions together if they are closer in ∆R 2 = ∆φ 2 + ∆y 2 than R 2 , where ∆y is the separation in rapidity, ∆φ is the separation in azimuthal angle and R is the jet radius. In experimental analyses at the LHC typically R = 0.4 or 0.5. For two real emissions in the final state with momenta k 1 and k 2 , eq. (2.2) reduces to That is, for ∆R < R the two emissions are clustered into a jet and T f j is computed from the sum p j = k 1 + k 2 . For ∆R > R, each emissions forms its own jet and T f 1 and T f 2 are computed with p j = k 1 and p j = k 2 , respectively. The measurement function M jet f (T cut ) is inserted into the usual SCET operator matrix elements defining the beam and soft functions. In this case, the jets J(R) are obtained purely from the collinear or soft radiation within each sector. For details we refer to refs. [1,3]. The practical implementation of such a measurement in beam and soft function calculations is discussed below.
The factorization of the cross section with a T Bj or T Cj veto in refs. [1,3] is strictly speaking valid only to lowest order in an expansion in R. The possibility of clustering independent soft and collinear emissions into the same jet breaks the soft-collinear factorization of the measurement function with the corresponding corrections starting at O(R 2 ) [1]. Since this only affects the measurement itself but not the soft-collinear factorization of the amplitudes and SCET Lagrangian, these soft-collinear clustering corrections can be computed in the effective theory and are included in our results. We separate out the corresponding contributions in the two-loop beam and soft functions that are associated with the clustering of independent emissions. They are denoted with the subscript 'indep' and together with the corrections from soft-collinear clustering reproduce the two-loop clustering behaviour of independent emissions in full QCD (in the singular limit). In section 3, we give two prescriptions as to how this collection of terms can be treated in the NNLL resummation.
In addition, at O(R 2 ) (potentially) factorization breaking effects due to Glauber interactions can play a role [28]. At the perturbative level they first appear at O(α 4 s ) and are not discussed further here.
Our calculation of the beam and soft functions is organized in an expansion in R as well. We will give terms in this expansion up to orders high enough for all practical purposes. This expansion only involves even powers of R (up to few exceptional terms at lower orders). We find the R 2 expansion to converge very quickly, suggesting that the relevant expansion parameter is (R/R 0 ) 2 with R 0 2. (Similar observations have been made recently also in other contexts involving small-R expansions, see e.g. [29,30].) As pointed out in ref. [1], in the small-R limit one should also consider resumming the corresponding logarithms ln R appearing in the jet-vetoed cross section. The dominant contribution beyond O(α 2 s ) was obtained in ref. [31]. Their resummation at the LL was obtained in refs. [26,32], and very recently methods have been developed [33,34] that allow to systematically carry out their resummation to higher orders.
To perform the computation of the jet-algorithm dependent soft and beam functions we follow the same strategy used in refs. [1,19] by computing the difference to a reference soft or beam function defined with a global (jet-algorithm independent) measurement, with the reference functions having been computed elsewhere in the literature. A key property of the global reference measurements we use is that they coincide with the jet-dependent measurements for the case of one real emission. Then when we compute the differences, we only need to consider the double-real emission amplitudes. Since the measurement functions are different for T Bj and T Cj , we must perform a separate computation of the soft function for the two observables. Since T Cj is equal to T Bj at forward rapidities, the beam function is in fact the same for both observables and there is only computation to be performed.

Soft function
We discuss explicitly the calculation of the quark soft function, i.e., where the two partons initiating the hard process are quarks. The results for the two-loop gluon soft function are obtained in the usual way by replacing C F → C A due to Casimir scaling. Following ref. [1], the calculation is done in a different way for the 'uncorrelated' C 2 F and the 'correlated' C F C A and C F T F n f color channels. The notion of correlated and uncorrelated soft emissions is closely connected to the notion of webs in the context of non-Abelian exponentiation of eikonal matrix elements [35,36]. Without the potentially exponentiation breaking measurement operator the soft amplitude factorizes into sums of products of webs. For example in our two-loop calculation the C 2 F part of the double-real soft emission amplitude can be written as the product of two identical one-loop amplitudes (webs) proportional to C F . We therefore regard the two emissions in the C 2 F channel as independent or uncorrelated. In contrast, the C F C A and C F T F n f parts of the total soft amplitude are nonfactorizable two-loop webs, and we refer to the corresponding emissions as correlated. The method of calculation for the correlated and uncorrelated channels is described in the following.

C 2
F piece For the uncorrelated C 2 F part of the bare soft function, we conveniently compute the two-loop correction relative to the expectation from non-Abelian exponentiation for the unclustered case, i.e. half of the one-loop soft function squared. This difference is directly the two-loop clustering correction for independent emissions, so we denote it by ∆S bare (2) f,indep (T cut , R). The total C 2 F contribution to the bare soft function is then The clustering correction ∆S f,indep starts at order R 2 , as discussed earlier, i.e., non-Abelian exponentiation for the jet-veto soft function works up to terms of order R 2 . The first term in eq. (2.3) corresponds to using a reference measurement θ(T f 1 < T cut ) θ(T f 2 < T cut ), which separately restricts each emissions irrespective of their separation. The difference to eq. (2.2) then gives the measurement function ∆M f,indep that corresponds to ∆S (2) f,indep , For the T Bj veto, we split the measurement eq. (2.4) into two pieces: The jet and total rapidities are defined by (2.8) For the T Cj veto, we split eq. (2.4) as follows: In both cases, the contribution to ∆S associated with the first part of the measurement, ∆M f,indep,1 , starts to contribute at order R 2 and contains a 1/ divergent piece (using dimensional regularization with d = 4 − 2 ), and thus produces a single logarithm of T cut at O(R 2 ). The second part is found to start at order R 4 and is finite in four dimensions. 2 Both contributions can be obtained analytically order by order in the R 2 expansion. This calculation is performed in appendix A, and we give the results up to O(R 8 ) in section 3.2. The jet veto soft function is renormalized multiplicatively [1,3]. To convert the bare results of the above computation to the renormalized S (2) , we expand the relation S bare = Z × S to second order, For the case of the uncorrelated soft function contribution we can directly use eq. (2.3) together with the one-loop relation S bare(1) = Z (1) + S (1) and find This shows that the renormalized S (2) for the C 2 F channel is equal to the expectation from non-Abelian exponentiation for the unclustered jet veto, 1 2 (S (1) ) 2 , plus the finite part ∆S (2) indep of the clustering correction.
The correlated C F C A and C F T F n f contributions to the soft function are split into three pieces, G,f (T cut ) + ∆S (2) base (T cut , R) + ∆S (2) rest,f (T cut , R) . (2.15) Here, S G,f is defined for a (known) global reference measurement. For T Bj , S G,B (T cut ) is given by the cumulant of the single-differential thrust soft function. It has been computed in refs. [37][38][39], and the explicit T B -differential expression can be found e.g. in ref. [40].
S G,C is the cumulant of the C-parameter soft function, which is defined e.g. in eq. (28) of ref. [41] and has been obtained numerically in refs. [41,42]. Since the reference measurements are chosen to coincide with the jet measurements for a single real emission, S (2) G,f already contains the correct real-virtual contributions. The measurement function ∆M f corresponding to the total difference ∆S G,f is thus given by where the first line is the full jet measurement for two emissions in eq. (2.2) and the second line subtracts the reference measurement M G,f = θ(T f 1 + T f 2 < T cut ).
The divergences of S f and S G,f differ by a 1/ term and are the same for T Bj and T Cj . The second quantity ∆S (2) base in eq. (2.15) is designed to capture this divergence and is the same for both T Bj and T Cj . The remaining piece ∆S (2) rest,f is then a finite correction. The measurement function for ∆S (2) base is defined as The global reference measurement essentially amounts to always clustering the emissions, and ∆M base corrects this to a constraint on the individual emissions when they are further apart than R. The reference measurement thus already captures most of the singularities. The remaining 1/ divergence of ∆S (2) f is associated with the limit of large jet rapidity |y j | → ∞. The correlated amplitude has support only in a finite range of ∆y = y 1 − y 2 around zero. The limits |y j | → ∞, |y t | → ∞ and y 1 , y 2 → ±∞ simultaneously are therefore effectively equivalent. In particular, in the latter limit ∆M f in eq. (2.16) becomes equal to ∆M base in eq. (2.17). Thus the divergence of ∆S (2) f is equal to the one of ∆S (2) base .
As we will see below, the 1/ divergence in ∆S (2) base can be isolated analytically, which makes ∆S (2) base much easier to compute than the full difference ∆S (2) f . In this sense, ∆S (2) base performs a very similar role as a subtraction term in a conventional fixed-order calculation. In particular, the remainder ∆S (2) rest,f = ∆S (2) f −∆S (2) base can be computed numerically setting d = 4 from the start.
The soft function S (2) f contains terms proportional to ln R and ln 2 R. These logarithms are entirely contained in ∆S (2) base , such that ∆S (2) rest,f is finite as R → 0, see appendix B. We compute these logarithms analytically as detailed below, while all remaining contributions are computed numerically. Some details regarding the setup for the numerical calculations are given in appendix B.
In the remainder of this section we describe explicitly the computation of ∆S (2) base . We first write it in terms of the momenta k 1 and k 2 of the two emitted gluons, corr (k 1 , k 2 ) is the amplitude obtained by adding all C F C A and C F T F n f terms of the double-real diagrams as given explicitly in appendix B of Ref. [39], ∆M base is the measurement function in eq. (2.17), and we have included the usual MS factor withμ = µ exp{ 1 2 [γ E − ln(4π)]}. We can rewrite the phase-space integral in terms of light-cone components in d = 4−2 dimensions and arrive at To incorporate the constraint on ∆R we write k ± 1 and k ± 2 in terms of the variables ∆y, y t , ∆φ, T T , and z, as in Appendix C of ref. [1], (2.20) The measurement function in eq. (2.17) can be expressed in these variables as .
As mentioned previously, the integration of A corr with ∆M base gives terms proportional to ln n R. We calculate these terms analytically by expanding the full integrand of the ∆y, ∆φ, z, T T integration (including the amplitude and phase space factors) for small ∆R ∼ ∆y ∼ ∆φ, keeping only the lowest order terms of order ∆R −2 . The integral may then be performed analytically by using the relations To compute the remainder of ∆S (2) base , we subtract the expanded integrand from the full one and integrate the result. The integrand has a sufficiently simple dependence on y t and T T that we can perform the integrations over these variables first analytically. The integral over y t simply gives a 1/(4 ) factor. The result can then be expanded in up to order 0 and integrated numerically over the remaining variables ∆y, ∆φ, z order-by-order in .
The ln R and ln 2 R terms in the soft and beam functions all arise from eqs. (2.22) and (2.23). They are remainders of the collinear divergence between particles 1 and 2, which is regulated by R. For a p T j veto there is only a single ln R. The ln 2 R terms we find for the T f j veto are related to the different treatment of rapidity divergences compared to the p T j case. In the latter an additional rapidity regulator is required, whereas in our T f j veto calculation rapidity-type divergences are effectively regulated (along with the usual IR and UV divergences) by dimensional regularization. The ln 2 R terms only arise from the O( ) terms in eqs. (2.22) and (2.23) and therefore cannot appear in the anomalous dimensions. For the same reason, their coefficient in the fixed-order terms is the same as the coefficient of the ln R terms in the anomalous dimensions. Consequently, they cancel between the beam and soft functions at fixed order due to RG consistency.
To convert the bare result into a renormalized one, we again use eq. (2.12). For the correlated channels, the cross term Z (1) S (1) however vanishes, implying that S (2) is just the finite part of S bare (2) here. In particular, we can obtain the C F C A and C F T F n f contributions to the renormalized jet-dependent S (2) f exactly as shown in eq. (2.15), namely by adding to the corresponding reference S (2) G,f the finite part of ∆S (2) base as well as ∆S (2) rest,f .

Beam Function
For the beam function, the calculation is done in the same way for all color and partonic channels. The jet-veto beam functions B i can be computed as a convolution of perturbative matching coefficients I ij with the standard parton distribution functions (PDFs) f j as [5,43,44] To obtain the matching coefficients at two loops, we compute the two-loop partonic beam functions B ij , see refs. [5,13,14,44] for details and related definitions in terms of SCET operator matrix elements. The partonic jet-dependent B ij are calculated via the difference ∆B ij between them and the reference beam functions as where t cut ≡ xp − T cut with p − being the large light-cone momentum of the incoming parton. The reference functions B G,ij (t cut , x) are defined using the global reference measurement θ(T < T cut ), where T is the total plus momentum of all real emissions. They are given by the cumulant of the virtuality-dependent beam functions calculated in refs. [13,14], The measurement function corresponding to ∆B ij is given by the difference of eq. (2.2) and the global reference measurement, In the second step we used that in the collinear sector we always have T i = k + i and T j = k + 1 + k + 2 for both T Bj and T Cj . In addition, we have the usual measurement δfunction that fixes the large light-cone momentum fraction of the parton that enters the hard process to x.
Although collinear matrix elements, like the beam functions, do not exponentiate, we nevertheless introduce a notion of 'correlated' and 'uncorrelated' emissions also for their calculation. We note, however, that this distinction is purely technical and there is no one-to-one correspondence to the different color factors. It is inspired by the behavior of the amplitude in the limit where (at least) one emission becomes soft.
We start by taking the double-real emission amplitudes A for the partonic beam function for each parton and color channel, previously calculated in refs. [13,14]. These amplitudes are gauge invariant and have been calculated in both Feynman and axial (light-cone) gauge as a cross check. We again denote the momenta of the emitted partons as k 1 and k 2 , and (without loss of generality) we have symmetrized the amplitudes such that they are symmetric under interchanging k 1 ↔ k 2 . For each color channel we separate A into two pieces, where A A is defined by and A B = A − A A is the remainder. Labeling the parton entering the hard process with i and the incoming parton with j, the term A A always has the form The color factors are defined as and theP (0) ij (x) are the unregularized one-loop splitting functionŝ Note that A A is nonzero only for certain parton and color channels. As can be seen from eq. (2.29), A A is essentially obtained by taking the soft limit for one of the partons. In this limit the amplitude, since it contains a sum over all diagrams, falls apart into a product of a one-loop soft amplitude A S(1) i and a one-loop collinear amplitude A B(1) ij by Ward identity arguments. As should be clear from its structure, the term A A corresponds to the uncorrelated emission part of the amplitude. This term is the part of the amplitude that when expressed in terms of the ∆y, ∆φ, z, T T , y t variables in eq. (2.20), and when multiplied by the appropriate Jacobian becomes flat in ∆y, ∆φ as well as in y t in four dimensions. It is also the only part of the amplitude that remains in the zero-bin subtraction terms (see below).
The integration of the term A B yields only a 1/ divergence, and is essentially identical in character to the integral of the soft correlated amplitude. In particular there are no divergences associated with the z, ∆y, ∆φ or T T integrations, so the same integration variables and techniques as for the correlated soft contributions can be used to integrate this piece. Integrating the A A piece yields a much deeper divergence structure, but this piece can be integrated in a straightforward way analytically.

'Correlated' piece
To integrate A B , the y t integration can be performed using the δ-function for the large minus momentum, after which the expressions for k ± 1 and k ± 2 in terms of the remaining integration variables z, ∆y, ∆φ, T T are The integrand contains an explicit factor of (1 − x) −1−2 which we pull out and expand in terms of distributions: The remaining part of the integrand is then simply expanded in and integrated over the remaining variables. There are no additional divergences associated with these integrals, and they yield a function that is finite in the limit x → 1. Note that for the channels with i = j, this function in fact vanishes for x → 1, such that for these channels we can replace the right hand side of eq. (2.42) by 1/(1 − x). This means that there is only a 1/ piece for the i = j channels proportional to δ(1 − x), as one may anticipate from the fact that the anomalous dimension for the beam function is diagonal in flavour. The integration of A B gives terms proportional to ln R and ln 2 R for some color structures. We note in passing that in axial gauge the ln n R terms arise from only one diagram topology, namely the 'bubble insertion' graph shown in figure 2f ) and 2o) of ref. [13] for the quark beam function and in figure 1i) and 1o) of ref. [14] for the gluon beam function. This topology is also sketched on the left-hand side of figure 1. As for the soft function, we calculate the ln n R terms analytically by expanding the integrand for small ∆R, and performing the integration using eqs. (2.22) and (2.23). Compact expressions for the required small ∆R limits of the beam function amplitudes are given in appendix D.
The remaining contributions from A B can then be obtained by subtracting the expanded integrand from the full integrand and integrating this numerically, as was done for the correlated soft function contributions. However, here this must be done for various points in both the x and R directions, requiring to generate a 2D grid of points. We did compute such a grid -however, since such results are not so straightforward to use or present, we also follow an alternative approach, the results of which are given in section 3 below. In this alternative approach we explicitly compute the leading R-dependent terms proportional to R 2 of the non-δ(1 − x) parts of the beam function analytically. The R 0 terms in the non-δ(1−x) piece are computed numerically, along with the full R dependence of the δ(1 − x) piece. For these pieces we can then give simple one-dimensional functions (in either x or R) obtained by interpolating the numerically generated points. The results from this approach are sufficiently accurate for smaller R values. More quantitative statements and a comparison of these approximated results to the exact numerical results are given in section 3.
To compute the terms proportional to R 2 , we expand the integrand up to order ∆R 0 . We then take the θ(∆R > R) in the measurement eq. (2.27) and divide it into two pieces, which yields the integral Here dΦ denotes the phase-space integration over all of the relevant integration variables, in particular over ∆φ and ∆y. The I −2 term denotes the part of the integrand of order ∆R −2 , whilst I 0 is the part of O(∆R 0 ). It is clear that the integrations involving the '1-term' in the right factor will give a R-independent result, which we can drop if we are only interested in the R 2 terms. The integration of the remaining O(∆R 2 ) terms with θ(∆R < R) gives contributions of O(R 4 ), so we can drop these too. The integration of the I −2 with the θ(∆R < R) generates the above mentioned ln n R plus some O(R 0 ) terms. The R 2 contributions we want to isolate thus exclusively come from the integral 3 − dΦ I 0 θ(∆R < R) . (2.45) We performed this integral analytically for the non-δ(1 − x) parts of each beam function. 4 The non-δ(1 − x) piece ∝ R 0 we obtained by numerically computing the integral in eq. (2.44) at R = 0.2, and then extrapolating to R = 0 using the leading R-dependence just obtained, which works very well for the small R values between 0 and 0.2. We do not directly compute the integral at R = 0 to avoid numerical instabilities.

'Uncorrelated' piece
Next, to integrate the uncorrelated A A part of the amplitude, our method is as follows. We first write again θ(∆R > R) = 1 − θ(∆R < R) as in eq. (2.43). The integration of A A with the '1-term' can be straightforwardly done analytically without a change of variables. This yields terms constant in R that can be as divergent as 1/ 2 . Note that the '1-term' here corrects the reference measurement to the fully unclustered case. The θ(∆R < R) term then is the analogue to the clustering correction for the soft uncorrelated emissions in eq. (2.4).
The integration of A A with θ(∆R < R) can be done analytically order-by-order in R after changing variables to ∆y, ∆φ, z, T T , and like in the soft case yields terms starting at O(R 2 ). For the qq C 2 F and gg C 2 A channels only, this piece yields a divergent R 2 / term that is exactly the one anticipated in eq. (50) of ref. [1], but with the opposite sign.
For this piece, the zero-bin subtractions are nonvanishing. Any of the zero-bin limits, k 1 soft, k 2 soft, or both soft, picks out only the term A A from the full amplitude A. The 3 This technique can easily be extended to obtain the R 2n corrections for arbitrary n, and we have indeed obtained R 4 expressions for the diagonal channels as discussed below. 4 Of course the analogous analytical calculation can be done for the O(R 2 ) pieces in the δ(1 − x) terms.
For these we can however perform a 1-parameter fit for the full R dependence with x = 1, which is accurate enough that analytical results of the R 2 contributions are not needed.
zero-bin computation yields terms of O(R 2 ). For the qq C 2 F and gg C 2 A channels the total zero bin contribution gives a divergent R 2 / term, which is exactly twice the naive result from the θ(∆R < R)×A A term above. The full R 2 / contribution, given by subtracting the zero-bin contribution from the naive result, then exactly reproduces and confirms eq. (50) of ref. [1]. Note that the zero-bin contribution for the '1-term' is scaleless and vanishes, i.e., the zero-bin for the fully unclustered contribution is scaleless just like in the inclusive case. We discuss the full computation of the θ(∆R < R) × A A term and the corresponding zero-bin terms in more detail in appendix C.

Renormalization
The result of the above calculation yields the difference between the bare jet-veto and reference beam functions, ∆B bare . To convert this to a difference between renormalized beam functions, we must expand the relation B bare = Z ⊗ B for the two beam functions to two loops. Even though Z (1) and B (1) are the same for both (B (1) G = B (1) ), the two beam functions renormalize differently, i.e., the nature of the convolution between the counter terms Z and the renormalized beam functions B is different. In particular, for the jet-veto beam function this convolution is just the simple product. Taking the difference between the renormalized beam functions according to eq. (2.25) we thus obtain where we define The second term in square brackets in eq. (2.46) accounts for the different renormalization of the global reference beam function. The result for ∆Z (2) i can be converted to the corresponding clustering correction for the two-loop anomalous dimensions, ∆γ (2) Bi . The relation needed to achieve this can be obtained by expanding Z ⊗ γ B = −dZ/d ln µ to two-loop order for the two beam functions and taking the difference, . (2.48) The result we find for ∆γ (2) B is precisely as was predicted in ref. [1]. Namely, it is (−1/2) times that of the soft function, plus an additional term related to soft-collinear mixing as required by RG consistency. This is an important check of our calculation.
Finally, the difference in two-loop renormalized partonic beam functions ∆B (2) ij can be converted into a difference in two-loop matching coefficients ∆I (2) ij by expanding eq. (2.24) for a partonic incoming state to two loops for the jet-veto and reference beam functions, and taking the difference. Since ∆I (0) ij and ∆I (1) ij are both zero, this immediately yields ∆I (2) ij (t cut , x, R, µ) = ∆B (2) ij (t cut , x, R, µ) . (2.49) The contributions to ∆I (2) ij from the integration of A A with the unclustered '1 term', as well as the terms in square brackets in eq. (2.46), are of O(R 0 ). As mentioned above, their purpose is to correct the running of the beam function to be multiplicative rather than involving convolutions. We therefore present these pieces together with the integration of the correlated amplitude A B below. On the other hand, the integration of A A with θ(∆R < R), together with corresponding the zero-bin terms, are of O(R 2 ). They are in fact the clustering corrections for independent emissions (analogous to ∆S (2) indep ), and so we separate these pieces off, and denote them using the subscript 'indep' in the results below.

Results
Here we present the two-loop results for the T f j -dependent beam and soft functions, together with the finite part of the two-loop soft-collinear mixing terms discussed in ref. [1]. For completeness, we also present the corresponding anomalous dimensions of these functions. In certain places in the results we give functions of R obtained by fitting numerically generated points. These points were generated using the GlobalAdaptive NIntegrate routine from Mathematica, and have a numerical uncertainty at the sub-per-mille level. We checked this by using the Suave routine from the Cuba library [45] to perform the same integrations and confirmed that the results agreed at the sub-per-mille level with the GlobalAdaptive method once a sufficiently large number of integration points were used. The number of data points we used for the fits was 30 in each case, spanning the range R = 0.05 − 1.5. We checked that the fit functions reproduce each point at the sub-per-mille level, and that there are no visible interpolation artifacts (such as 'polynomial wiggle') in these fits. In section 3.1.2 we give functions of x that have also been fitted. The points for these fits have also been generated using Mathematica NIntegrate with sub-per-mille precision (and checked using Suave). The fitting procedure for these functions is slightly more involved and described at the end of section 3.1.2.
Our final results for the renormalized two-loop beam and soft functions to be used in the NNLL resummation are given by 5 where the one-loop coefficients do not yet depend on R and are simply given by the cumulants of the corresponding differential functions. The two-loop contributions are written as

+ ∆S
(2) f (T cut , R, µ) + ∆S (2) f (T cut , R, µ) , I G,ij (t cut , x, µ) + ∆I (2) ij (t cut , x, R, µ) + ∆I (2) ij (t cut , x, R, µ) . Here, I G,ij is the cumulant of the virtuality-dependent beam function matching coefficients, which can straight-forwardly be obtained from the results of refs. [13,14]. For f = B the function S G,f is the cumulant of the renormalized thrust soft function and for f = C it is the cumulant of the renormalized C-parameter soft function.
The results for the corrections due to the jet clustering are split into two parts. The corrections ∆S (2) f and ∆I (2) ij are not associated with the clustering of independent emissions. For simplicity we call them 'correlated' pieces. They are given in section 3.1. Finally, the ∆S (2) f,indep and ∆I (2) ij are the corrections from the clustering of independent emissions, which start at O(R 2 ). They are given in section 3.2, together with the associated two-loop softcollinear mixing terms. There, we also give two alternative prescriptions for incorporating these terms in the resummation at NNLL .

We find ∆S
(2) with the anomalous dimension correction and where C i = C F for the quark soft function and C i = C A for the gluon soft function. The functions F are fitted from numerically generated points using the form and the fitted coefficients a i are given in table 1. The functions f are different for T Bj and T Cj . They are fitted from numerically generated points using the form and the fitted coefficients c i are given in table 2.
R-dependence of the δ(1 − x) terms in ∆I gg and ∆I qqV . They are fitted using the form in eq. (3.8), and the fitted coefficients are given in table 3. The function h(x) gives the x-dependence of the 0 non-δ(1 − x) piece in ∆I at the point R = 0.2. Where appropriate we subtracted out the analytically-calculated small-R result and/or the endpoint x = 1 behaviour. For this function, a more complicated form is used. In fact, it is fitted separately in two regions: 0.0009 < x < 0.6 and 0.6 < x < 1, where x = 0.0009 is the smallest x value for which we generate data. In the high-x region 0.6 < x < 1, we use a form that is equal to the first few terms of a Taylor expansion around x = 1. In the low-x region 0.0009 < x < 0.6, we use a form equal to the first few terms of a Taylor expansion around x = 0, plus terms with powers of ln x up to the third power. For the gg and gq channels only, a final term proportional to 1/x is also included. It can be shown that only these channels have a piece in h(x) proportional to 1/x by expanding the respective integrands for small x. Summarizing, we use the following form to fit h(x): This functional form is satisfactory for every h(x) function in eqs. (3.16) - (3.12). The only exception is h qqS , which falls to zero so steeply near x = 1 that it cannot be described well by the above high-x fit function. For this function only we therefore use a high-x fit function given by n=3,9d n (1 − x) n . Including these higher powers of (1 − x) ameliorated the fit, and in this case including the lower powers of (1 − x) up to the second power is not necessary for a good fit. We have performed the fits for the coefficients d n andd n separately for each h(x), each using 30 data points in the relevant x range. We checked that the resulting fit reproduces all the data points at the sub-per-mille level, and confirmed by eye that it does not contain interpolation artifacts. We also checked by eye that the low-x fit transitions smoothly onto the high-x fit. The fit coefficients for all h(x) functions are presented in appendix E.
As mentioned in section 2.2 and indicated above, our results for ∆I ij only contain terms up to O(R 2 ) for the non-δ(1 − x) pieces. For the reader interested in the more precise behaviour in R of these pieces, we have also obtained full 2D grids by numerically integrating the beam function amplitude minus its small ∆R expansion, spanning 0.0009 < x < 1 and 0.2 < R < 1.2. These results can be provided upon request.
For most parton and color channels, the level of error of the O(R 2 )-accurate expressions above compared to the full numerical results is a few permille for R values up to around 0.6 − 0.8, and a few percent at the largest R values ∼ 1.0 − 1.2. The exceptions to this are the qqV C 2 F and gg C 2 A channels, where the errors hit the percent level already for R = 0.6, and are much bigger for larger R values. 6 For these two channels we have also computed the R 4 corrections, and their inclusion improves the agreement with the full numerical result to a similar level as the other channels. The expressions for these R 4 terms may also be provided upon request.

Independent emission pieces
The corrections to the two-loop soft and beam functions associated with the clustering of independent emissions, together with the soft-collinear mixing term, are given by U B (R) and U C (R) are the terms associated with the measurements ∆M B,indep,2 and ∆M C,indep,2 in eqs. (2.7) and (2.11), respectively, cf. appendix A. Further terms in the R expansion may be easily computed, but have a negligible impact for R < 1.5. From the above expansions it seems clear that the expansion parameter is R/2 or even smaller. As mentioned already, the factorization of the jet-veto measurement, does not hold at O(R 2 ) due to the soft-collinear mixing contributions. As a result, the scale dependence from the independent emission terms only cancels, when all three contributions in eqs. There are two possibilities as to how one can treat these O(R 2 ) terms at NNLL . In the first, we note that at the two-loop order we work here, the soft-collinear term has the same type of logarithms as the beam function. This allows us, at this order at least, to absorb the soft-collinear mixing terms into the beam functions, such that in eqs. (3.2) and (3.3) we have ∆I (2) ij (t cut , x, µ) = ∆I (2) ij,indep (t cut , x, R, µ) + SC (2) ij (t cut , x, R, µ) , ∆S (2) f (T cut , R, µ) = ∆S (2) f,indep (T cut , R, µ) , (3.28) and the µ-dependence cancels between ∆I (2) ij and ∆S (2) f . In practice, this scheme is effectively equivalent to the one employed in ref. [18] for a p T j veto.
Let us define the anomalous dimensions of the beam and soft functions (for either T Bj or T Cj ) as where γ G,S is the cumulant of the anomalous dimension for thrust (and C-parameter), and γ G,B is the cumulant of the anomalous dimension for the virtuality-dependent beam function. Then in the scheme defined by eq. (3.28), the jet clustering corrections to the (noncusp) anomalous dimensions are given by where the last C 2 i R 2 term is from the independent emission contributions. Alternatively, we can simply exclude eqs.
f (T cut ) = 0 , (3.31) and in addition the last C 2 i R 2 term is removed from eq. (3.30). Instead, we combine all terms into a common independent emissions contribution and add the following correction term to the 0-jet cross section: (3.32) This equation is somewhat schematic -the ⊗ symbols represent the Mellin convolution in light-cone momentum fractions along with the appropriate flavour sums, and all pieces are evaluated at a common scale µ. In the resummed cross section this term can then be multiplied with the overall evolution factor of the O(R 0 ) cross section. This corresponds to what is done in refs. [17,19] in the context of a p T j veto.

Conclusion
In this paper, we have computed the two-loop beam and dijet soft functions for the rapiditydependent jet-veto observables T Bj and T Cj . Each function is computed as the difference from the known results for the corresponding global jet-independent observable. The beam and soft functions have been computed for the complete set of color and parton channels.
The dominant parts of the functions in the small-R limit are computed fully analyticallye.g. the ln n R terms and the R 2 terms in the non-δ(1 − x) part of the beam functions. The remaining parts are determined in terms of finite numerical integrals, for which we provide simple interpolation functions.
Our results enable the resummation of color-singlet production cross sections with T Bj or T Cj jet vetoes at full NNLL order, and are a necessary ingredient for the N 3 LL resummation (with the missing ingredient being the corrections to the three-loop anomalous dimension). Our computed beam functions can also be used in the resummation of N -jet cross sections for which a central jet veto is imposed using T Bj or T Cj . Such measurements will provide important information on the production pattern of additional jets in hard processes at the LHC.

Acknowledgments
This work was supported the DFG Emmy-Noether Grant No. TA 867/1-1. J.G. acknowledges financial support from the European Community under the FP7 Ideas program QWORK (contract 320389). M.S. is supported by GFK and by the PRISMA cluster of excellence at Mainz university. The authors thank the Erwin Schrödinger Institute program "Challenges and Concepts for Field Theory and Applications in the Era of the LHC", and J.G. thanks JGU Mainz, for hospitality while portions of this work were completed. The figure in this paper was produced using JaxoDraw [46].

A.1 T Bj veto measurement
We first compute the contribution to the soft function corresponding to ∆M B,indep,1 in eq. (2.6). We have We express this integral in terms of light-cone momentum components, perform the transverse integrals using the delta functions, and change variables from k + 1 , k + 2 , k − 1 , k − 2 to k + 1 , k + 2 , y t , ∆y, with y t and ∆y defined in eq. (2.20). This gives For the calculation of ∆S B,indep,2 we use the variables in eq. (2.20) and express the measurement as After performing the T T , ∆φ, y t integrations we expand in small ∆y ∼ R and finally integrate over z and ∆y. The result is given in section 3.2.

A.2 T Cj veto measurement
For T C , the starting expression is identical to eq. (A.1), albeit without the 2θ(y t > 0) piece, and with k + 1 and k + 2 replaced by T C1 and T C2 in the theta functions of the second line. We change variables in this case from k + 1 , k + 2 , k − 1 , k − 2 to T C1 , T C2 , y t , ∆y, with The measurement function for ∆S C,indep,2 in eq. (2.11) can be written as (0 < X < 1) with X = T T T Cj = 2e 2(∆y+yt) (2z − 1) sinh(∆y) sinh(2y t ) + cosh(∆y) cosh(2y t ) + 1 e ∆y+2yt + e 2∆y (1 − z ) + z e 2yt (e 2∆y − 1)z + 1 + e ∆y , (A. 8) and We can now directly carry out the intergrations over T T and ∆φ. After expanding in small ∆y ∼ R we then integrate over y t , z , and ∆y to arrive at the result given in section 3.2.

B Soft function correction: ∆S rest
For completeness we document in this appendix how we organize the numerical computation of ∆S (2) rest,f in eq. (2.15). In particular we show the split of the corresponding measurement function ∆M rest,f , we have chosen in order to conveniently divide the calculation into separate parts. We note however that other decompositions (or no split at all) might work just as well in practice.
For the T Bj veto we write where ∆M rest,B1 = 4θ(∆R < R) θ(y j > 0) θ(y 1 > 0) θ(y 2 < 0) The numbers (4 or 2) in front of the θ functions indicate that we have exploited phase-space symmetries and summed over equivalent contributions from different regions. We note that ∆M rest,B3 amounts to the difference between the (inclusive) cumulant beam thrust and cumulant double-hemisphere soft function measurements. Correspondingly, ∆S is a finite R-independent constant. A simple dimensional analysis also shows that it is independent of T cut . The piece ∆S (2) rest,B1 is of O(R) and ∆S (2) rest,B2 contains O(R 0 ) terms. For the T Cj veto we write Here we have defined A = k + 1 /T T = z (1+e −∆y−2yt ) and B = k + 2 /T T = (1−z )(1+e ∆y−2yt ). The corresponding part ∆S (2) rest,C2 turns out to be of O(R 2 ), whereas ∆S (2) rest,C1 is of O(R 0 ). The form of the different parts of the measurements in ∆M rest,f already explains why ∆S (2) rest,f does not contain terms proportional to ln n R. For the measurements ∆M rest,B1 and ∆M rest,C2 this is obvious, because they are proportional to θ(∆R < R). A term in the integrand proportional to 1/∆R 2 , which could potentially generate such a logarithm, must be absent, because otherwise it would cause a divergence for ∆R → 0. The measurements ∆M rest,B2 and ∆M rest,C1 on the other hand effectively vanish linearly in the small ∆y ∼ ∆R limit. To see this for the latter measurement one has to take into account that in four dimensions the associated integrand is proportional to 1/T T . We are therefore free to rescale T T → T T /(A + B) in the second term of eq. (B.4). As the two-loop soft function integrands contain at most 1/∆R 2 poles, these measurements thus prohibit logarithmic terms in R.
The integrations for ∆S

C Beam function correction: independent emission clustering part
We first calculate the O(R 2 ) contributions associated with the θ(∆R < R) × A A term in eq. (2.43). We have Let us first consider the evaluation of the 1/ piece of this expression. For this purpose many of the terms in eq. (C.3) can be set to 1, and we obtain where again t cut ≡ xp − T cut . Plugging in the values for C i andP ij (x) from eqs. (2.32) and (2.33), one observes that the divergent part is zero when i = j, and for the case of i = j corresponds to half of eq. (50) in ref. [1], but with the opposite sign, as stated in the main text. Now let us move to consider the zero-bin subtractions [47]. There are several zero bins one must consider: k 1 alone soft, k 2 alone soft, and both k 1 and k 2 soft simultaneously. The first two we must subtract, and the final one we add. Our procedure for computing the zero bins coincides with that of refs. [1,48]. Namely, we take the appropriate soft limit of the full amplitude A, we take the limit of the minus momentum delta function and leave the rest of the measurement function in eq. (2.27) untouched. As mentioned in the main text, taking the soft limit of the amplitude A in fact picks out only the term A A for any of the zero bins.
We now make the decomposition in eq. (2.43). However, for any of the zero bins, the '1-term' contains scaleless integrals and is therefore zero. For example, the '1-term' for the k 1 soft zero bin is proportional to Next we consider the θ(∆R < R) term. For the k 1 soft zero bin it equals It is clear to see that the 1/ piece of this is identical to eq. (C.4). The k − 2 soft gives an identical contribution. There remains the contribution in which both k − 1 and k − 2 go soft simultaneously. Here the minus delta function just becomes δ[(1 − x)p − ] and factors out of the expression, and the remainder looks very similar to a soft function calculation, so we can use the same variables (∆y, ∆φ, y t , z) as we use in that case. The integration for y t then looks like 7 Thus, the net effect of the zero bin in the 1/ terms is to invert the sign of eq. (C.4), to the sign predicted in ref. [1]. We can also compute the full expression for the integral of the θ(∆R < R) × A A and zero-bin terms, including the finite pieces. Here we give the analytic expressions for the θ(∆R < R) × A A term: and the sum of the three zero bin terms, which must be subtracted from the result of the naive beam function computation.

D Beam function amplitudes in the small-R limit
It is possible to write compact formulae for the small-R limit of the beam function amplitudes A in terms of one-loop splitting functions. This is expected since at small R the 7 Note that in the computation of ∆S (2) base and ∆S (2) B,indep,1 , see e.g. eq. (A.3), we obtain the same integral, except accompanied by an explicit theta function constraining yt > 0. Thus in the latter case we obtain 1/(4 ) rather than zero. amplitude factorizes into two splitting processes. The first is a splitting i → j + m, where i is the initial state parton in A, and j is the physical parton that goes into the operator with momentum fraction x. The second is the splitting m → k + l, where k and l are the partons that go into the final state and generate the ∆R. m is an intermediate parton whose nature can be determined by quark number conservation from the other partons. As discussed in the main text, there is only one diagram that contributes in the small-R limitnamely, the one in figure 1 that already has the topology of a diagonal i → j + (m → k + l) process. Figure 1 also illustrates the two splitting processes that this graph decomposes into at small R.
There are two distinct cases, corresponding to whether the intermediate parton m is a gluon or a quark. We must treat the two cases differently because if one has a quark in the initial state, the plane of the splitting is not correlated to whether the quark has positive or negative helicity. However, when one has a gluon, there is a well-known correlation between the plane of splitting and the gluon polarization.
The formula for the case when m = q (orq) is i→jq (x, ) originates from the first graph on the right hand side of figure 1, whilst P (0) q→kl (z, ) originates in the second graph. The final factor in eq. (D.1) arises from converting the transverse momentum denominators (etc.) in the splitting amplitudes relative to the initial-state parton direction (for the first amplitude on the right hand side of figure 1), or relative to the parton m (for the second amplitude), to our variables. In this formula one must always sum over all possibilities for (k, l) -i.e. (q, g) and (g, q). To get the full expression for A R in d dimensions, one must use the d-dimensional splitting functions given, e.g. in ref. [49] The result when m is an intermediate gluon is a little more complex and given by g(in)→kl (z, )∆y 2 +P The splitting functionsP i→jg(in/out) (x, ) are splitting functions depending on whether the final state gluon is polarized in or out of the splitting plane. These functions may be computed from results in ref. [50]. The two nonzero cases are given bŷ Here, the second outgoing parton is polarized, while the remaining spin/polarization indices are summed over or averaged as appropriate, and x is the momentum fraction of the first outgoing parton relative to the incoming quark. SimilarlyP where x is again the momentum fraction of the first outgoing parton relative to the incoming gluon. In the case where the intermediate gluon decays into quarks one must remember to sum over the possibilities (j, k) = (q,q) and (j, k) = (q, q). Say, without loss of generality, that the initial parton travels along the z direction, and gluon m is emitted somewhere in the x − z plane. Then the first term corresponds to the case in which the gluon polarization lies in the x − z plane (call this plane P 1 ). If the gluon splits in plane P 1 , resulting in jk having a separation in ∆y, then the gluon polarization and splitting planes are coincident, so we should weight the ∆y 2 factor with the g(in) → kl splitting function. On the other hand, if the gluon splitting is in the plane which contains the m direction and the y direction (call this plane P 2 ), the gluon polarization and splitting planes are perpendicular. Here, partons jk gain a separation in ∆φ, so we weight the ∆φ 2 factor with the g(out) → kl splitting function. The second term in eq. (D.4) corresponds to the case where gluon m has polarization in the plane P 2 . Here we must swap the weighting factors multiplied by ∆φ and ∆y around, for obvious reasons. The final term is actually where the gluon polarization lies in one of the 'extra' −2 dimensions. Here, the gluon polarization is outside both splitting planes. So both splitting function weightings are for the 'out' polarization.

E Fit function coefficients for beam function
In these tables, the notation used for the numbers is such that a −b means a × 10 −b -e.g. 4.197602 −3 means 4.197602 × 10 −3 . The functional form for the fit is given in eq. (3.21).
Coefficients for the low-x fitting functions: x 3 x 4 x 5 gg CACA