The NNLO gluon beam function for jet-veto resummation

We compute the gluon beam function for jet-veto resummation to next-to-next-to-leading order (NNLO) in the strong-coupling expansion. Our calculation is based on an automated framework that was previously used for the computation of the respective quark beam function, and which we significantly extended for the present calculation. In particular, the perturbative matching kernels are directly calculated in momentum space, without the need to perform an additional Mellin transform. We present results for both gluon and quark-initiated processes, which we cross-checked with an independent semi-analytical method that exploits the similarity of the beam functions to the more familiar case of transverse-momentum resummation. Our computation is relevant for jet-veto resummations at NNLL$'$ accuracy.


Introduction
Many experimental analyses at the Large Hadron Collider (LHC) apply jet vetoes to enhance the signal sensitivity.In these analyses, the jet veto is usually imposed as a cut on the transverse momenta p veto T of the reconstructed jets.If this scale is much smaller than the typical hard scale Q of the process, the QCD radiation gets constrained to be either soft or collinear to the beam directions, which induces large Sudakov-type corrections of the form α n S ln 2n (p veto T /Q).These logarithms challenge the convergence of the perturbative expansion, and they need to be resummed to all orders to obtain reliable predictions.
The starting point of the SCET analyses is a factorisation theorem, which at leading power in p veto T /Q ≪ 1 can be written in the schematic form where the sum runs over all partonic channels, Q is the invariant mass and Y the rapidity of the colourless final state with x 1,2 = (Q/ √ s) e ±Y .Here the purely virtual corrections to the Born process are contained in the hard functions H ij , whereas the beam functions B i/h and the soft functions S ij describe the effects from collinear and soft radiation, respectively.The factorisation theorem furthermore assumes that the jet radius R satisfies λ ≪ R ≪ ln(1/λ) with λ = p veto T /Q [2].While the hard functions H ij (Q, µ) are process-dependent, the soft functions S ij (p veto T , µ) trivially depend on the process via the colour representations of the partons i, j.They are currently known to next-to-next-to-leading order (NNLO) accuracy [6,27,28].The beam functions B i/h (x, p veto T , µ), on the other hand, are universal and depend on the parton i within the hadron h that carries the longitudinal momentum fraction x after emitting collinear radiation that passes the jet-veto constraint.They are by themselves non-perturbative objects, but as long as the jet-veto scale satisfies p veto T ≫ Λ QCD , they can be matched onto the standard parton distribution functions as we will review below.The corresponding matching kernels can then be computed perturbatively, and they were determined to NLO in [2,9].A subset of the current authors extended this calculation to NNLO for quarkinitiated processes [29], whereas the full set of matching kernels was determined afterwards in [30].These results were subsequently used to perform NNLL ′ resummations in [15,16]. 1he purpose of this work is twofold.First, we will complete the NNLO calculation initiated in [29] by computing the matching kernels for gluon-initiated processes.This serves, on the one hand, as a cross-check of the calculation in [30], which was based on a different technique and, in particular, a different prescription to regularise rapidity divergences.We will also study the endpoint behaviour of the matching kernels in some detail since these are often divergent and therefore hard to control numerically.In addition we for the first time present results for the gluon beam function in Mellin space.Second, we stress that the method we use for the calculation of the jet-veto kernels is generic, and it can be applied to a much broader class of observables.It actually follows a strategy that was successfully applied for the calculation of soft functions within the SoftSERVE approach [27,[31][32][33], which relies on a universal parametrisation of the observable-dependent measurement function in Laplace space.In our previous study [29], we performed a Mellin transform to resolve the distributions associated with the momentum fraction x, whereas we extract these distributions directly in momentum space in the current work.We view this extension as a major improvement of our automated approach, and we anticipate many further results of this novel framework in the future.We, in fact, already recalculated the matching kernels for transverse-momentum resummation and the event-shape variable N-jettiness in this setup [34][35][36].
To validate our results, we performed an independent semi-analytical calculation that is closer to the method that was used in [30].Specifically, it exploits the similarity of the jetveto beam functions to the ones relevant for transverse-momentum resummation, which are known analytically to the considered NNLO and beyond [37][38][39][40][41].In the difference between the two beam functions, many contributions drop out and the calculation can moreover directly be performed in four space-time dimensions.While similar in spirit, our method is not equivalent to the one that was used in [30] as we will explain below.
The remainder of the paper develops as follows: We first introduce the theoretical foundation of our framework in Sec. 2. In Sec. 3 we provide the computational details of both our automated numerical approach and the semi-analytical method that was used for cross-checks.In Sec. 4 we present our results for the gluon beam function as well as a comparison to the results from [30].We conclude in Sec. 5, collect the relevant anomalous dimensions and splitting functions in App.A, and provide details on the reference observable that was used in the semi-analytical approach in App.B. Finally, in App.C we present our novel results for the quark beam function in momentum space.

Theoretical framework
To set up the notation, we introduce two light-like vectors n µ and nµ that obey n 2 = n2 = 0 and n • n = 2. Any four-vector k µ can then be decomposed according to its projections In this notation, the gluon beam function for jet-veto resummation is defined as where is the n-collinear gluon field operator, W n a collinear Wilson line, and the sum over X represents the phase space of the final-state partons with momenta {k i }.At the tree level, when there is no collinear emission, this is simply the vacuum state, and at NNLO it denotes the phase space of up to two massless partons.The hadronic state with collinear momentum P µ = P − n µ /2 is furthermore indicated by |h(P )⟩, whereas the jet-veto constraint is imposed by the measurement function M(p veto T ; {k i }), which we will specify in Sec. 3 below.As a matrix element of hadronic states, the beam function is not directly accessible in perturbation theory.But as long as the scale that constrains the collinear emissions satisfies p veto T ≫ Λ QCD , it can be matched onto the usual parton distribution functions f i/h (x, µ) according to which holds at leading power in Λ QCD /p veto T .Our goal consists in computing the matching kernels I i←k (z, p veto T , µ) for gluon-initiated processes (i = g) to NNLO in QCD.However, since our previous results for the quark channels (i = q) were provided only in Mellin space [29], we will for completeness also determine those kernels directly in momentum space in this work.The matching calculation can, in fact, be most easily performed using partonic on-shell states, since the parton distribution functions then evaluate to f i/j (x, µ) = δ ij δ(1 − x) to all orders in perturbation theory when dimensional regularisation is used to regularise both ultraviolet (UV) and infrared (IR) divergences.The partonic calculation therefore directly yields the matching kernels in this case.
As stated in the introduction, the jet veto is usually imposed on the transverse momenta of the reconstructed jets.It is well known, however, that transverse-momentum dependent observables belong to a special class of observables that suffer from rapidity divergences that are not captured by the dimensional regulator ϵ = (4 − d)/2.This follows from the fact that the relevant soft and collinear modes have the same virtuality in the effective theory, which is also known as SCET-2.To regularise these rapidity divergences, we use the analytic phase-space regulator from [42], which is introduced for each emission with momentum k µ i via where α is the rapidity regulator and ν the associated rapidity scale.The same prescription is also implemented for the soft integrals in SoftSERVE, but for collinear emissions with k − i ≫ k + i it can be further simplified at leading power.The very fact that the rapidity regulator in (2.3) respects the n-n symmetry of the process, ensures that the beam functions for collinear and anti-collinear emissions are equivalent to all orders in perturbation theory.
The matching kernels defined in (2.2) thus depend on the scheme that is used to regularise the rapidity divergences and, in particular, the rapidity scale ν.To obtain a result that is scheme independent, we follow the collinear-anomaly approach [24,25], which states that the product of the soft and collinear functions can be refactorised as where the matching kernels I i←j (z, p veto T , µ) on the right-hand side of this equation are independent of the rapidity scale ν.The dependence on the hard scale Q is furthermore resummed through the collinear-anomaly exponent F gg (p veto T , µ), which satisfies the RG equation with the cusp anomalous dimension in the adjoint representation Γ A cusp (α s ).Up to two-loop order, its solution is given by where L = ln µ/p veto T , and Γ A i and β i are the expansion coefficients of the cusp anomalous dimension and the β-function, respectively, that are given explicitly to the required order in App. A. As the anomaly exponent describes the physics from the overlap of the soft and collinear regions, the non-logarithmic terms d A i in (2.6) are constrained by Casimir scaling, and they are therefore related to the corresponding quantities for the quark channels, which were given in our previous work [29].We also note that the anomaly exponent renormalises additively, F bare gg = F gg + Z F gg , and the corresponding MS counterterm Z F gg obeys a similar RG equation as in (2.5), whose two-loop solution reads The remaining ingredients in the collinear-anomaly relation (2.4) are the refactorised matching kernels I i←j (z, p veto T , µ), which are independent of the rapidity regularisation scheme.Their scale dependence is determined by the RG equation where Γ R i cusp (α s ) is the cusp anomalous dimension in the representation of the parton i, γ i (α s ) is the collinear anomalous dimension for quarks (i = q) or gluons (i = g), and P k←j (z, α s ) are the DGLAP splitting functions.The renormalisation of the refactorised matching kernels is, in fact, slightly more complicated, and it is convenient to introduce two types of counterterms that subtract the UV divergences of the beam function (Z B i ) and the IR divergences associated with the matching onto the parton distribution functions (Z f k←j ) for them [29].Specifically, we write where the UV counterterm obeys the RG equation (2.10) Up to two loops, its solution reads and the specific values for the coefficients of the anomalous dimensions are also summarised in App. A. The IR counterterm, on the other hand, satisfies the RG equation whose two-loop solution is given by where the coefficients of the DGLAP splitting functions P (i) k←j (z) are defined in App. A. We furthermore introduced a short-hand notation for the convolutions that appear frequently in the calculation.Note that this notation implies a sum over intermediate partonic channels with index l.For completeness, we give the explicit expressions for these convolutions also in App. A. Finally, we state the solution of the RG equation (2.8) for the renormalised matching kernels.Up to NNLO, it takes the form As the logarithmic terms in this expression are already known, we will focus on the nonlogarithmic terms I (m) i←j (z) in later sections.

Computational aspects
As stated in the introduction, we followed two different approaches to compute the renormalised matching kernels for jet-veto resummation.The first approach is based on an automated numerical framework that can be used for a generic class of observables, while the second approach is a semi-analytical method that exploits the similarity of the jet-veto beam functions to the ones relevant for transverse-momentum resummation.We will now describe each of these methods in turn.

Numerical approach
In the first approach, we start from the definition of the beam function (2.1) with the measurement function replaced by the ansatz where τ is a Laplace variable of dimension 1/mass, and the function ω({k i }) specifies a generic observable.Following [29,35], we find it convenient to work in Laplace space, since both the calculation of the bare beam function and its renormalisation simplify in this case.We note that for the specific jet-veto beam function defined in (2.1), the measurement function is not of the form (3.1), but it is instead given by M(p veto T ; We then follow the procedure described in [27] to convert this into the form (3.1) by taking a Laplace transform, and we extract the results in the original p veto T space by inverting the Laplace transformation at the very end of the calculation.
To factorise the implicit divergences of the collinear matrix elements, we use universal phase-space parametrisations as in [29,35].Particularly, for a single emission with momentum k µ , we use the magnitude of its transverse momentum relative to the beam axis, which is set by P µ ∝ n µ , and we allow for a non-trivial azimuthal dependence of the observable, where θ k is the angle between ⃗ k ⊥ and a reference vector ⃗ v ⊥ that could be imposed by the observable.The remaining momentum components are then fixed by the on-shell condition and the explicit delta function in (2.1), yielding k + = k2 T /k − and k − = zP − with z = 1 − z.This leads to the following measurement function for a single emission, For two emissions, we proceed similarly and parametrise the momenta k µ and l µ of the massless final-state partons in the form where again and we in addition now have three non-trivial angles in the transverse plane, with θ k and θ l referring to the external vector ⃗ v ⊥ as before, and θ kl being the angle between ⃗ k ⊥ and ⃗ l ⊥ .We then rewrite these angles in terms of variables t k , t l , and t kl that are defined on the unit interval, similar to (3.2).
In terms of these variables, the two-emission measurement function for jet-veto resummation takes the form where R is the jet radius and ∆ = ln 2 a + arccos 2 (1 − 2t kl ) is the distance measure of the jet algorithm.The twofold structure in (3.5) arises as follows: If the distance between the two emissions is smaller than R, the emissions are clustered together and the veto is imposed on the recombined pseudo-particle.If, on the other hand, this distance is larger than R, the veto constrains both of the reconstructed jets.Similar to our previous work on the quark beam function [29], we consider the general class of k T -type jet algorithms in this work, for which the expression in (3.5) turns out to be independent of the specific clustering prescription (k T , Cambridge/Aachen, anti-k T , . . . ) [2].
Having defined the phase-space parametrisations and the structure of the measurement function, the partonic beam functions can now be computed.The required collinear matrix elements are related to the spin-averaged d-dimensional time-like splitting functions [43].Specifically, we used the expressions from [44][45][46] for the real-virtual contributions, and the ones from [47,48] for the double real-emission part, and we applied standard crossing relations to convert these to time-like kinematics [36].The main task then consists in factorising the implicit phase-space divergences, which allows one to perform the expansion in the two regulators α and ϵ directly on the integrand level.While this can easily be achieved in the one-emission case, this is only true in certain cases for the double real emissions in the parametrisation (3.4).We therefore employed a number of further techniques that involve sector-decomposition steps and non-linear transformations to factorise all divergences for this contribution (details can be found in [36]).
In comparison to our previous work on the quark beam function [29], our novel framework is capable of computing the beam functions directly in momentum space.In particular, within the parameterisations (3.2) and (3.4) it is possible to factorise the divergence associated with the longitudinal momentum fraction z.We may then introduce the corresponding distributions via After expanding the integrand first in the rapidity regulator α and subsequently in the dimensional regulator ϵ, the coefficients in this double expansion can be integrated numerically.The final expressions for the matching kernels then consist of distributions with numerical coefficients and a non-trivial z-dependent 'grid' contribution that we sample for different values of z.To perform these steps, we have implemented our formalism in the publicly available program pySecDec [49], and we use the Vegas routine of the Cuba library [50] for the numerical integrations.By following this procedure, one obtains the regulator-dependent matching kernels on the left-hand side of the refactorisation condition (2.4).In the last step, we then combine these expressions with the corresponding NNLO soft function that is provided by the SoftSERVE distribution to extract the I (m) i←j (z) part of the refactorised and renormalised matching kernels that are defined in (2.15).

Semi-analytical approach
In the second approach, we use the fact that the beam functions for jet-veto resummation are related to the ones for transverse-momentum resummation.As the latter are known analytically at the considered NNLO, one focuses on the difference where the first term refers to the jet-veto matching kernels that appear on the left-hand side of the refactorisation condition (2.4), and the second term indicates the corresponding kernels for a so-called reference observable.Following [5], we choose for the latter the integrated q T -spectrum, and we provide details on the derivation of the corresponding refactorised matching kernels At NLO the difference in (3.7) vanishes by construction, and at NNLO it only receives contributions from diagrams with double real emissions.The measurement function capturing these contributions can be expressed as where the angular separation ∆ between the two emissions is quantified as before.
As has been further argued in [5], there is a price to pay when working with the difference (3.7), i.e. one has to explicitly compute mixing terms that originate from different sectors of the effective theory.Notably, these contributions are independent of the jet radius R and they emerge solely within the uncorrelated contribution to the amplitude.To manage these contributions effectively, we decompose the double-emission amplitude in the collinear region A (2) c into its uncorrelated and correlated components, wherein the uncorrelated term is constructed from the product of one-emission collinear and soft amplitudes, Therefore only specific colour structures receive mixing contributions.The starting point for the evaluation of (3.7) is the representation where the rapidity regulator α has been factored in.As the divergences in the dimensional regulator ϵ are absent for the measurement (3.9), one can immediately set d = 4 in the calculation.To evaluate this contribution we employ the following phase-space parametrisation [5], In the computation of the correlated term, we then utilize integration techniques as delineated in [5].Initially, we determine the asymptotic behaviour in the limit R → 0 by expanding the integrands around this limit.Moreover, the integral includes non-singular terms, which can be represented as a power series in R 2 .To proceed, we first subtract the leading singularities as ∆ϕ and ∆y approach zero, and we evaluate the coefficients of the power series in R 2 numerically.The evaluation of the uncorrelated contribution necessitates consideration of the mixing terms as discussed above.For this part, we rewrite the theta function θ(∆ − R) in (3.9) as 1 − θ(R − ∆), which bifurcates the analysis into R-independent and R-dependent terms.The R-dependent terms are free of mixing contributions, and their computation parallels the treatment of the correlated part.Conversely, for the R-independent terms one starts from the measurement function and the mixed collinear/soft contribution can in this case be calculated from where we stress that the rapidity regulator is retained in its original form (2.3) in the soft region.Due to the n-n symmetry of the process, there is no need to evaluate the mixing between anti-collinear and soft contributions explicitly, and for the collinear/anti-collinear mixing one finds a scaleless integral that vanishes in the applied regularisation prescription.Finally, for the double soft emissions one has where A s denotes the soft double-emission amplitude, whose explicit expression can be found e.g. in [51].
Focusing on concreteness on the diagonal gluon channel, the various contributions can be put together and the final correction to the refactorised matching kernel in (2.15) can be obtained from ∆I (2)  g←g (z) = ∆I (2)  g←g (z, p veto T , µ, ν) + ∆I (2,cs)  g←g (z, p veto T , µ, ν) where the last term subtracts the contribution from the two-loop anomaly exponent, which we decompose as ∆d A 2 (R) = ∆d which agrees with previous calculations [5] 3 .The extraction is less involved for the offdiagonal channels, and in each case, we can write the final expression for the refactorised matching kernels in the form where the coefficient of the plus-distribution is given by the anomaly exponent in the representation R i of the parton i, and the coefficient of the delta function has a similar expansion in terms of colour factors and R-dependence.Particularly, we obtain for the diagonal gluon channel and the corresponding expressions for the diagonal quark channel can be found in App. C. Finally, we obtain purely numerical results for the non-distributional terms H i←j (z, R) in (3.19) that depend on both the momentum fraction z and the jet radius R.
We conclude this section by noting that our semi-analytical approach is similar in spirit to the method that was used in [30], but it differs in several technical aspects.First and foremost, our calculation uses a different rapidity regulator, which is purely analytic, and therefore all zero-bin subtractions are scaleless and vanish in our setup.Moreover, we use slightly different techniques for computing the uncorrelated-emission contribution, as well as a different reference observable than [30] as described above.

Results
In this section, we present our results for the non-logarithmic contributions to the renormalised matching kernels I as a function of the jet radius R. The dots show the result of our numerical approach, and the lines represent the expression in (4.2).Right: The same for the coefficient of the delta function where the lines refer to (4.5).
we can directly compute these kernels in momentum space, but for completeness, we also present the respective Mellin-space results that were obtained with the method described in [29].While we focus here on the gluon channels with i = g, our novel momentum-space results for the quark kernels (i = q) are collected in App. C. As QCD is invariant under charge conjugation, the anti-quark kernels (i = q) directly follow from these expressions.

Momentum-space kernels
At NLO, the algorithmic nature of the jet veto does not show up yet, and the matching kernels can be obtained analytically.Specifically, they read [2] whereas the one-loop anomaly coefficient d A 1 vanishes at this order.At NNLO, we first verify if our setup reproduces the known results for the two-loop anomaly coefficient d A 2 and the non-cusp anomalous dimension γ g 1 .As the former depends on the jet radius, we evaluate it numerically for three different values of R ∈ {0.2, 0.5, 0.8} using the method described in Sec.3.1.In the left panel of Fig. 1, we compare these numbers that are indicated by the dots against the semi-analytical expression where the first two terms originate from the reference observable (3.8), and the latter two terms were given in (3.18) above.The plot shows that the two calculations are in excellent agreement within the uncertainties that are too small to be visible on the scale of the plot.We also remark that the anomaly exponent exhibits Casimir scaling, and it is therefore related to the anomaly coefficient in the fundamental representation d F 2 (R) that is relevant for Drell-Yan production.On the other hand, the non-cusp anomalous dimension is independent of the jet radius R and it is constrained by RG consistency.Writing where the numbers in square brackets show the known analytic results.The agreement provides another strong check of our computation.
Coming to the renormalised matching kernels I (2) i←j (z), we first consider their distribution-valued component.To this end, we decompose the kernels in the form where the coefficient of the plus-distribution is given by the two-loop anomaly exponent in the representation R i of the parton i, and it thus carries no new information.The coefficient of the delta function, on the other hand, is for the gluon channel given by with explicit expressions for the last two terms given in (3.20).Our numerical results for these coefficients are compared to these semi-analytical expressions in the right panel of Fig. 1, which again shows perfect agreement between the two calculations.We finally turn to the grid contributions, which we decompose further in terms of their colour structure, where we kept the R-dependence of the individual kernels on the right-hand-side of this equation implicit.We evaluated these kernels numerically for three different values of the jet radius R ∈ {0.2, 0.5, 0.8} and 125 values of the momentum fraction z using the method described in Sec.3.1.The results are displayed in Fig. 2, and they are also provided in electronic form in the file that accompanies the present article.Specifically, one notices that the off-diagonal kernel I (2,C F T F n f ) g←q (z) does not show any dependence on the jet radius, since it only receives real-virtual contributions, whereas all the other kernels exhibit a dependence on R, which is more pronounced for small values of the momentum fraction z.Similar to the plots in Fig. 1, the numerical uncertainties are not visible on the scale of the plots, with the only exception being the low z-region of the off-diagonal I    uncertainties are generically at the sub-percent level, we observe accidental cancellations in this case, which make them somewhat more pronounced for this particular kernel.We next compare these numbers against the output of our semi-analytical approach.For this purpose, we define the following ratios involving the grid contributions of the two approaches, where the subscripts 'Numerical' and 'Semi-analytical' refer to the two setups of our calculation described in Sec.3.1 and Sec.3.2, respectively, and the index X labels one of Lower: The same for the ratio of the grid contributions between our numerical approach and the results of [30] as defined in (4.8).
the colour structures in the notation of (4.6).We present a template comparison for one diagonal and one off-diagonal kernel in the upper panels of Fig. 3.In general, we find a very good agreement between the two calculations, which supports our previous claim that the uncertainties of our numerical approach are at the sub-percent level.This can clearly be seen in the bulk of the distributions in the first panel of Fig. 3.In these plots one also notices that the relative uncertainties are enhanced in those regions where the central values of the kernels are close to zero.As argued above, the situation is special for the off-diagonal channel that is shown in the right panel of Fig. 3 because of accidental cancellations, and in this case one finds relative uncertainties of a few percent even in regions where the central values of the kernels are not small.We remark, however, that these uncertainties are likely to be overestimated, since the shifts of the central values are in all plots much smaller than what is indicated by the error estimate.
As stated in the introduction, the jet-veto matching kernels were previously determined to NNLO in [30].Although our Mellin-space results for quark-initiated processes were already available at that time [29], a comparison between the two calculations has not yet been presented.Here we perform such a comparison.Note that this requires to recombine Figure 4: Endpoint regions z → (upper row) and z → 1 (lower row) for the same matching kernels as in the lower panels of Fig. 3, where our are shown in colour and the ones from [30] in gray.The smallest z values (upper row) and the largest z values (lower row) of the grids provided in [30] are indicated by the dashed vertical lines.
the scheme-dependent kernels from [30] with the corresponding soft function [28] in order to extract the refactorised matching kernels according to (2.4).Using (2.15) and (4.4) one then determines the grid contributions I (2,Grid) i←j (z, R), which we use for the comparison.Specifically, we define the ratios Ref. [30] , which are displayed for two template kernels in the lower panels of Fig. 3.As is evident from these plots, we find a very good agreement between the two calculations for all values of the momentum fraction z (the pattern around the zeroes of the matching kernels being similar to the one in the upper panels).In order to verify if this agreement persists in the endpoint regions z → 0 and z → 1, where most of the kernels are enhanced, we show the same matching kernels in those regions in Fig. 4.While it becomes hard to distinguish our results (in colour) from the ones of [30] (in gray) in these plots, we find (i) that the agreement extends nicely into both endpoint regions, and (ii) that our numbers are stable even for extreme values of z, which allows one to sample these regions to very high accuracy. -300 -60

Mellin-space kernels
In a previous study [29], a subset of the current authors computed the jet-veto matching kernels for quark-initiated processes at NNLO directly in Mellin space.These kernels are related to the refactorised matching kernels defined in (2.4) via and one similarly denotes the Mellin transform of the non-logarithmic contribution in (4.4) by I (2) i←j (N ).While our results from the previous section can in principle be used for this conversion, the discretised grids would introduce systematic uncertainties in addition to the statistical errors from the Monte-Carlo integration.We therefore follow the strategy outlined in [29] here, and compute the Mellin-space kernels directly within our numerical approach by including an additional integration over the variable z according to (4.9).At NLO this integration can be performed analytically, and one obtains At NNLO, we then sample the Mellin-space kernels for ten values of the Mellin parameter N ∈ {2, 4, 6, 8, 10, 12, 14, 16, 18, 20} as in [29].Writing the Mellin-space kernels in the form we show our results for the individual coefficients in this decomposition in Fig. 5.We note that the numerical uncertainties are again at the sub-percent level in this case, and they are therefore not visible on the scale of the plots.The Mellin-space kernels are particularly useful when the resummation is performed in Mellin space, as was done e.g. for joint threshold and transverse-momentum resummation in [52].

Conclusion
We computed the beam-function matching kernels for jet-veto resummation to NNLO in QCD.Building on our earlier study of the respective quark beam function [29], we provide the full set of matching kernels in this work.To this end, we applied two different computational methods, and we found perfect agreement between the two approaches.Our final results are displayed in Fig. 2 and Fig. 7, and they are also provided in electronic form as supplementary material to this paper.The jet-veto matching kernels were previously determined to NNLO in [30].We performed a detailed comparison to these results, paying special attention to the endpoint regions z → 0 and z → 1, and found very good agreement for all kernels.This represents the first confirmation of the results in [30].In addition, we for the first time computed the NNLO matching kernels for the gluon channels in Mellin space.
Whereas our semi-analytical method is specific to the jet-veto observable, our numerical approach is generic, and it can therefore be applied to a much broader class of observables.In contrast to our previous study [29], in which we determined the matching kernels in Mellin space to avoid distribution-valued expressions, we have extended our formalism in this work such that it becomes capable of calculating the matching kernels directly in momentum space.We view this extension as a major improvement of our framework, and we anticipate many further applications in the future.In the long term, we plan to provide a public code for the computation of NNLO beam functions in the spirit of the SoftSERVE distribution.

A Anomalous dimensions and splitting functions
We define the coefficients in the perturbative expansion of the anomalous dimensions appearing in the RG equations (2.5) and (2.8) via and likewise for the QCD β-function The relevant coefficients of the cusp anomalous dimension are given by where i = F refers to the fundamental and i = A to the adjoint representation.Up to two-loop order, the collinear quark and gluon anomalous dimensions read [53,54], whereas we only need the one-loop coefficient of the QCD β-function, β 0 = 11 3 C A − 4 3 T F n f , in this work.
The two-loop solutions of the RG equations given in (2.13) and (2.15) involve convolutions of one-loop splitting functions.In the notation introduced in (2.14), which implies a sum over partons with index l, one has g←l ⊗ P (0) where the factor n f in the last line arises from a sum over intermediate massless quarks.
Similarly, the convolutions of the one-loop splitting functions with the one-loop jet-veto matching kernels read g←l ⊗ P where the factor n f in the last term again arises from a sum over all massless quark flavours.

B Details on the reference observable
Within the approach described in Sec.3.2, the jet-veto matching kernels are computed starting from the difference (3.7) with respect to a reference observable, for which we use the integrated q T -spectrum according to (3.8).One then needs to add the matching kernels for this reference observable to the final expressions for the refactorised matching kernels given in (3.19).
The matching kernels for the reference observable can easily be derived from available results for transverse-momentum resummation.As the latter are usually given in positionspace, one needs to Fourier-transform these expressions and to integrate them up to the p veto T cut.Specifically, logarithms in the transverse-momentum framework then translate into p veto T logarithms using the relation where J 1 (x) is a Bessel function.For n = 0, 1, 2 there is thus a one-to-one correspondence between the two types of logarithms, whereas for n = 3, 4 they differ by certain ζ 3 terms.This implies that the refactorised matching kernels for the reference observable are not precisely equal to the ones from transverse-momentum resummation for the diagonal channels.Particularly, we find where i←j (z) are the matching kernels for transverse-momentum resummation, for which we use the explicit expressions provided in [38] (setting the RG logarithms therein to zero).

C Quark beam function in momentum space
In [29] we presented the first computation of the quark beam function for jet-veto resummation in Mellin space.Here we complement these results by calculating the matching kernels directly in momentum space using the two approaches described in Sec. 3. At NLO the matching kernels can be obtained analytically.They read [9] whereas we decompose the NNLO kernels in terms of distributions and a grid contribution according to (4.4).Our results for the coefficients of the distributions are shown in Fig. 6, where the numbers indicated by the dots have been obtained within our numerical approach from Sec.   for the coefficient of the delta function.We recall that both of these coefficients are constrained by Casimir scaling, and they are therefore related to the expressions in (3.18) and (3.20).The agreement between our numerical and semi-analytical results provides a check of our calculation.We finally turn to the grid contributions, for which we use the same colour decomposition as in [29], In total there are thus seven independent kernels in this case, which implicitly depend on the jet radius R. Our results for these kernels are shown in Fig. 7, and they are also contained in the accompanying electronic file.We note that these numbers, which were obtained with our numerical setup, have sub-percent uncertainties, and they show a similar level of agreement in comparison with both our semi-analytical approach and the results from [30] as for the gluon channels.

( 2 ,
C F C A ) g←q (z) kernel.While these

Figure 2 :
Figure 2: Grid contributions to the NNLO matching kernels defined in (4.4) and (4.6) for three different values of the jet radius R.

8 Figure 3 :
Figure3: Upper: Ratio of the grid contributions between our numerical and our semi-analytical results as defined in (4.7).We show two sample matching kernels for three values of the jet radius R. Lower: The same for the ratio of the grid contributions between our numerical approach and the results of[30] as defined in (4.8).

Figure 5 :
Figure 5: NNLO matching kernels for three different values of the jet radius R in Mellin space.

Figure 6 :2 = d F 2 2 F + d C A 2 C F C A + d n f 2 C 2 F
Figure 6: Left: Two-loop anomaly exponent dF 2 = d F 2 2 F + d C A 2 C F C A + d n f2 C F T F n f as a function of the jet radius R. The dots show the result of our numerical approach, and the lines represent the expression in (C.2).Right: The same for the coefficient of the delta functionC F δ (R) = C C F δ C 2 F + C C A δ C F C A + C n f δ C F T F n f, where the lines refer to (C.3).

Figure 7 :2 2 +
Figure 7: Grid contributions to the NNLO matching kernels defined in (C.6) for three different values of the jet radius R.