NNLO zero-jettiness beam and soft functions to higher orders in the dimensional-regularization parameter $\epsilon$

We present the calculation of the next-to-next-to-leading order (NNLO) zero-jettiness beam and soft functions, up to the second order in the expansion in the dimensional regularization parameter $\epsilon$. These higher order terms are needed for the computation of the next-to-next-to-next-to-leading order (N$^3$LO) zero-jettiness soft and beam functions. As a byproduct, we confirm the $O(\epsilon^0)$ results for NNLO beam and soft functions available in the literature.


Introduction
To find signals of physics beyond the Standard Model, many interesting processes at the LHC are being studied with ever increasing precision. An important part of these efforts is the development of methods that enable N 3 LO QCD calculations, at least for the simplest processes where colorsinglet final states are produced. In the absence of fully-developed N 3 LO subtractions schemes, a promising approach is the slicing method [6][7][8][9] that has seen a recent resurgence in the context of LHC physics.
Any slicing method is based on the idea that one can split the phase space for a process of interest into partially-resolved and fully-unresolved parts. The fully-unresolved contribution originates from virtual, real-soft and real-collinear emissions. Conversely, the resolved one requires a final state that contains at least one additional QCD jet in comparison to the lowest order final state and, for this reason, it must be computed through lower order in the perturbative expansion in QCD than the unresolved one.
Phase-space separation into fully-unresolved and resolved parts can be accomplished using different kinematic variables. The two most popular ones are p ⊥ and N -jettiness variables that have been used recently in many NNLO QCD computations [10][11][12][13][14][15][16][17][18][19][20][21]. In this paper we will deal with the so-called zero-jettiness variable that can be used to perform a slicing computation of N 3 LO QCD corrections to the production of a colorless final state V (H, W , Z, γ * , W W , ZZ, γγ, etc.) in hadron collisions. This variable reads [22] τ = m min i∈{1,2} where p i are the four-momenta of incoming partons, k m are the momenta of final state QCD partons and Q i are the so-called hardness variables. In the limit of small τ , the cross section factorizes [23] into a product of hard H, beam B and soft S functions All quantities that appear in Eq. (1.2) are known through NNLO QCD. Moreover, the hard function H is known through N 3 LO QCD for single vector boson and Higgs boson production [24,25] and, recently, the three-loop quark-to-quark matching coefficient, needed to relate the beam function to parton distribution functions, was computed in the generalized large-N c approximation in Ref. [26]. 1 The computation reported in Ref. [26] required the knowledge of certain NNLO beam functions through second order in the dimensional regularization parameter . These functions were calculated in Ref. [28] and the results of that computation were used in Ref. [26]. The goal of this paper is twofold. First, we aim to extend the calculation reported in Ref. [28] and to compute all NNLO QCD matching coefficients through the second order in , as required for the calculation of matching coefficients through N 3 LO QCD. Second, we will compute the NNLO QCD soft function through the second order in , as required for the calculation of the N 3 LO QCD soft function. We note that NNLO QCD zero-jettiness beam functions were computed in Refs. [1][2][3] through zeroth order in , whereas the NNLO soft function was originally calculated in Refs. [4,5].
To extend the calculation of beam and soft functions to higher orders in , we use methods that may be of interest in their own right. Indeed, we employ collinear and soft limits of QCD amplitudes [29,30], reverse unitarity [31] and integration-by-parts identities [32] to show that computation of soft and all NNLO beam functions for zero-jettiness can be significantly simplified. In the case of the soft function, we rewrite step functions that arise from the definition of the zero-jettiness variable as integrals of delta functions over auxillary parameters before applying reverse unitarity. We note that these methods allow one to express any NNLO zero-jettiness beam function through just twelve and the NNLO soft function through just nine simple (phasespace or loop) integrals. In case of the soft function, integrations over auxillary parameters turn out to be remarkably simple.
The remainder of the paper is organized as follows. In Section 2 we describe the computation of the partonic beam functions through O( 2 ) starting from collinear limits of scattering amplitudes and explain how the master integrals are calculated. We discuss the calculation of the bare soft function through O( 2 ) in Section 3. We conclude in Section 4. Finally, we note that results for the NNLO bare soft function and beam function matching coefficients are collected in an ancillary file provided with this submission.

Calculation of the beam function
In this section we describe the calculation of the bare partonic beam function. We split the discussion into two parts. In Section 2.1 we explain the general set up and relate the calculation of the beam functions to collinear limits of QCD amplitudes. We also use reverse unitarity to express bare beam functions through master integrals. In Section 2.2 we describe the calculation of these master integrals. We present some results in Section 2.3.

General setup
It was pointed out in Ref. [33] that a bare partonic beam function B b ij , that describes the transition of a parton j to a parton i, can be obtained by integrating spin-and color-averaged collinear splitting functions P j→i * {m} over an unresolved m-particle phase-space The phase-space measure is defined as follows where {m} is the set of collinearly-radiated partons. In Eq. (2.2) we denote the momentum of the incoming parton j as p, its complementary light-cone momentum asp and the momenta of final state partons as k m . Furthermore, t is the so-called transverse virtuality of the off-shell parton i, z · p is its longitudinal momentum and s = 2p ·p. It was explained in Ref. [29] how splitting functions P j→i * for all parton-to-parton transitions can be calculated. This requires the use of a physical (axial) gauge for gluons and projection operators that decouple collinear emissions from hard matrix elements. These projection operators act on matrix elements M j→i * {m} describing the process of a parton j splitting into on-shell partons {m} and an off-shell parton i * . Following Ref. [29], we write  6) where N m are symmetry and averaging factors. To compute all beam functions it is sufficient to consider i's and j's from the following set (i, j) ∈ {(q l , q m ), (q l , g), (q l ,q m ), (g, g), (g, q m )} [1,2], where the indices l and m denote quark flavours. We note that a flavour-preserving transition in B b q l qm is obtained by setting l = m. Similar to regular splitting functions, all other beam functions can be obtained from the above set. Examples of diagrams that are required for the calculation of beam functions are shown 2 in Fig. 1.
The bare partonic beam functions B b ij Eq. (2.6) can now be calculated as standard phase-space and loop integrals with the projection operator P as a special Feynman rule. To facilitate this computation, we apply reverse unitarity [31] and rewrite delta functions in Eq. (2.2) as differences of two "propagators" with opposite signs in the i0 prescription, mapping phase-space integrals in Eq. (2.6) onto loop integrals. We then use integration-by-parts (IBP) identities [32] to express the beam function through master integrals. The IBP reduction is performed using FIRE [35].
We find that all five beam functions can be expressed through just 12 master integrals. They include nine double-real master integrals , , , (2) , , (2.7) (2) , , 2 We use FeynGame [34] to draw Feynman diagrams. (2) , and three real-virtual master integrals , (2.8) where for a given integrand f we write [f ] (2) = dPS (2) f, We note that phase-space measures dPS (2,1) are defined in Eq. (2.2). We also note that the gluon and quark beam function share the same set of master integrals. We describe the calculation of these master integrals in the next section.

Master integrals
The master integrals shown in Eqs. (2.7) and (2.8) are sufficiently simple to be evaluated directly.
To illustrate the computation, we discuss three representative examples. All other master integrals can be calculated along similar lines. We begin with the master integral (2.10)  We start by rescaling the momenta p,p, k 1 and k 2 in such a way that the dependencies of the integrals on t and s factor out. To this end, we write 3 (2.11) and obtain I 6 (s, t, z) = t d−5 s −1 I 6 (1, 1, z) . (2.12) To simplify the notation we drop tildes over momenta and turn to the calculation of the following integral into the integrand and change the order of integration. We find We first compute the function F in Eq. (2.15) in the rest frame of the time-like vector Q. In that frame Q = (Q 0 , 0, 0, 0) and F becomes (2.16) We parameterize the two light-like momenta as p = p 0 (1, n p ),p =p 0 (1, np), 4 and introduce spherical coordinates for k 1 . We obtain where we introduced the notation p 1 = (1, n p ), p 2 = (1, − np) and k n = (1, n k ). The angular integral in Eq. (2.17) was discussed in Refs. [36,37]. The result reads 3 We note that for real-virtual master integrals we also rescale the loop momentum l →l √ t. 4 Note that in the rest frame of Q, p andp are not in a back-to-back configuration. where is the d-dimensional solid angle, 2 F 1 is the Gauss hypergeometric function and ρ 12 = (1 − n p 1 · n p 2 ). Finally, we rewrite ρ 12 in a Lorentz-invariant way .

(2.21)
To integrate over Q we employ the Sudakov decomposition In Eq. (2.22) we used the fact that Q 2 > 0, Q · p > 0 and Q ·p > 0 to constrain integrations over α and β. After eliminating the delta functions δ (β − 1/z) and δ (α − (1 − z)) by integrating over α and β, we obtain integrate over u and find where 3 F 2 is the generalized hypergeometric function [38]. The expansion of the hypergeometric function in is easily obtained using the program HypExp [39,40]. Our next example is the master integral into the integrand and write the integral as We compute the integral Eq. (2.27) in the rest frame of the vector Q. To this end, we parameterize the phase space as shown in Eq. (2.16), integrate over k 1 to remove the delta function, introduce spherical coordinates for k 2 and integrate over the absolute value of k 2 to remove the remaining delta function. We obtain the angular integral where we introduced the notation p 1 = (1, 1 λ np), with λ = 1/(Q 0p0 ) − 1, p 2 = (1, n p ) and k n = (1, n k ). The angular integration in Eq. (2.28) was discussed in Ref. [37]; the result reads (2.29) In Eq. (2.29) F 1 is the Appell hypergeometric function (see e.g. Ref. [41]). Writing Eq. (2.29) in a Lorentz-invariant way, we obtain .
We substitute Eq. (2.30) into Eq. (2.26), introduce the Sudakov decomposition Q µ = αp µ + βp µ + Q µ ⊥ and integrate over Q . We substitute Q 2 ⊥ = l(1 − z)/z and find To perform the l-integration we use the integral representation of the Appell function [38] We find We would like to expand the integrand in a Laurent series in and compute the integral order by order in this expansion. This can be done if the integrand remains integrable at = 0. It is easy to see that this is not the case; while the integral over l in Eq. (2.33) converges if we Taylor expand around = 0, the integral over u diverges at u = 1.
We remove the divergence by performing an end-point subtraction at u = 1, splitting the integral into two pieces. To write the result, we define two functions and write Eq. (2.33) as (2.36) The u = 1 singularity in the first term on the right hand side of Eq. (2.36) is now regulated, while the last term in Eq. (2.36) can be easily integrated over u. We find (2.37) All remaining integrands in Eq. (2.37) can now expanded to the required order in and integrated using the HyperInt package [42]. The final result reads where H( m w , z) are harmonic polylogarithms (HPLs) [43]. Finally, we consider the real-virtual master integral I 10 . It reads (2.39) We perform the l-integration first. To this end, we combine the propagators 1/l 2 and 1/l ·p . We write and obtain the standard loop integral over l (2.41) The integration is now straightforward and we obtain (2.42) The remaining integration over the on-shell momentum k is performed by introducing the Sudakov decomposition k µ = αp µ + βp µ + k µ ⊥ . We find This concludes the discussion of the evaluation of the master integrals. All manipulations with hypergeometric functions that appear in master integrals, including their expansions in , are performed with the help of the HypExp package [40]. We describe some results for the beam functions in the next section.

Results
We are now in a position to present the bare partonic beam functions By performing the renormalization procedure and matching onto partonic distribution functions, as discussed in Refs. [1,2,23,26], we also obtain the matching coefficients I q i q j , I q i g , I q iqj , I gg and I gq i . To present the results, we write the beam functions and the matching coefficients as a series in the renormalized MS coupling constant Since the expressions for the bare partonic beam functions B b ij and the matching coefficients I ij through O 2 are lengthy, we only discuss some features of the most complicated coefficient I gg ; complete expressions for all other matching coefficients are given in an ancillary file provided with this submission. We write the matching coefficient in the following form where we define the plus distribution (2.47) For brevity, we only show the coefficient C −1 as well as the function F δ,h (z) in pure gluodynamics (n f = 0). For the coefficient C −1 we find

Calculation of the soft function
In this section we describe the calculation of the bare zero-jettiness soft function S at NNLO in QCD. We begin by discussing the general setup in Section 3.1, relating the calculation of the soft function to soft limits of QCD amplitudes for color singlet production. We re-write step functions, that originate from the zero-jettiness measure, as integrals of delta functions over auxillary parameters. We then use reverse unitarity to express the soft function through master integrals. In Section 3.2 we describe the calculation of master integrals as functions of the auxillary parameters and explain in Section 3.3 how the remaining integrations over auxillary parameters can be performed.

General setup
The zero-jettiness bare soft function can be calculated by considering soft limits of scattering amplitudes for colour singlet production. These soft limits, described by eikonal functions, were calculated through NNLO QCD in Refs. [29,30]. We extract them from that reference and integrate the obtained expression over the m-particle unresolved phase space dPS (m) S including the m-particle zero-jettiness measure M m for the set of radiated partons {m} with momenta k m .
We begin by writing the bare soft function as a series in the bare strong coupling constant where we defined The lower order results read where C a = C F (C A ) if the incoming particles are quarks(gluons), respectively. At NNLO we need to consider the following contributions to the soft function where the functions ξ (2) g,qq,gg denote various eikonal functions and for m = 1, 2 we introduced the short-hand notation The first term in Eq. g has been calculated to arbitrary order in in Ref. [4]. It reads and we thus focus on the double-real emission pieces. The zero-jettiness measure for two real partons reads [4] where the momenta p andp are again two complementary light-like vectors and we set p ·p = 1/2. We refer to different sets of delta functions and step functions in Eq. (3.7) as "configurations". Since the integrands in Eq. (3.4) are invariant under exchange of p andp, it is sufficient to only consider two configurations M 2 = 2 M A + 2 M B , which we refer to as A and B. Hence, we write For color-singlet production, the quantities ξ , (3.14) with p 1 = p, p 2 =p. We note that S (2) qq and S (2) gg were obtained in Refs. [4,5] by directly integrating T ij and ξ ij over the relevant phase space. We will discuss an alternative to this approach, that is in line with the beam function calculation discussed in Section 2. We hope that this approach can be extended to enable an N 3 LO calculation of the zero-jettiness soft function.
To this end, we would like to employ reverse unitarity and IBP technology to simplify calculation of the soft function. To do so, we map step functions on to delta functions, using the following identity which holds for a, b ∈ [0, ∞). Since k 1,2 · p, k 1,2 ·p ∈ [0, ∞), Eq. (3.15) is applicable. We therefore rewrite Eqs. (3.8) and (3.9) as follows To illustrate this point, we discuss the computation of S (2) qq in detail; the computation of S (2) gg is analogous. According to our earlier discussion, contributions to the soft functions due to an emission of a qq pair read qq,B . (3.18) We note that we have split Eq. (3.18) into two contributions, stemming from configurations A and B. They read We proceed by writing all delta functions in Eqs. (3.16) and (3.17) as linear combinations of the corresponding "propagators" and performing partial fractioning. We find that in configuration A all integrals can be mapped onto two integral families (3.20) where for a given integrand f we write We perform the IBP reduction using FIRE [35] and obtain the following master integrals I qq,3 00 = 1 (1) , , . (3.23) For configuration B, we obtain two integral families (3.24) , (3.25) that are mapped on the following master integrals I qq,2 00 = 1 (2) , , I qq,2 01 = (k 1 · k 2 ) −1 (2) , , (3.26) . In Eqs. (3.24) -(3.26) we used We describe the calculation of the master integrals in the next section.

Master integrals
The master integrals shown in Eqs. (3.23) and (3.26) can be evaluated directly. When describing this calculation below, we will always assume that z 1 > z 2 since all contributions to the soft function are symmetric with respect to z 1 ↔ z 2 permutation.
To illustrate the simplicity of the computation, we discuss the calculation of the most complicated master integral. We provide explicit solutions to all other master integrals in Appendix A. We consider the master integral The computation proceeds as follows. We begin by performing the Sudakov decomposition of the two light-like momenta k 1,2 k 1,2 = α 1,2 p + β 1,2p + k 1,2⊥ . (3.29) The integration measure dPS (2) S is then written as Note that integrations over α and β extend from zero to infinity with constraints imposed by δ-functions. We write I qq, 3 11 = The angular integrations in Eq. (3.31) were discussed in Ref. [4]. The result reads Because of the delta functions in Eq. (3.31), we need Eq. (3.32) for (3.33) The hypergeometric function can be simplified using the following identity 2 F 1 1, which is valid for |z| < 1. Since we work in the region where z 2 < z 1 , we can immediately use Eq. (3.34) to simplify Eq. (3.33). We obtain Remarkably, the hypergeometric function in Eq. (3.35) is independent of the parameters α i and β i , allowing for a straightforward integration. We substitute Eq. (3.35) back into Eq. (3.31), integrate over α 1 , α 2 , β 2 and change the integration variable β 1 → β 1 = β 1 /τ . We find Upon integrating over β 1 , we obtain the following result for the most complicated of the nine master integrals needed to describe the NNLO soft function (3.37) A complete list of master integrals can be found in Appendix A. We note that for the gluon emission contribution S (2) gg no further master integrals are required. This concludes our discussion of the evaluation of master integrals. We discuss the remaining integration over the auxillary parameters z 1,2 in the next section.

Integration over auxillary parameters
We express the double-real contributions in terms of master integrals and write S (2) qq,A as follows . (3.38) It appears that upon substituting solutions for the master integrals Eqs. (A.1) -(A.4) into Eq. (3.38), we will have to perform non-trivial integrations over z 1 and z 2 . However, after changing variables z 2 = t z 1 , the z 1 integration factors out. The remaining t integration seems to include terms that are proportional to (1 − t) −4 . However, upon taking the limit t → 1 we find that the most singular term actually scales like (1 − t) −1−2 and, therefore, can be easily subtracted.
We perform an endpoint subtraction at t = 1, expand the integrand in a Laurent series in and compute the integral order by order in with the help of HyperInt [42]. The final result reads (3.39) The physical meaning of the auxillary variables z 1 and z 2 can be understood by considering the Sudakov decomposition of k 1,2 . The singularity at z 1 = z 2 describes the limit were the quark and the anti-quark become collinear to each other, while the z 1 = 0, singularity describes the kinematic configuration in which the gluon, that emits the qq pair, becomes collinear to the light-like directions p µ . The τ → 0 limit controls the double-soft divergence. Next, we discuss the contribution S anti-quark are emitted into different hemispheres. Thus both z 1 = 0 and z 1 = z 2 (collinear) singularities should be absent. We therefore expect that we can simply expand the integrand in Eq. (3.40)   The calculation of S (2) gg can be performed in the same way. While the gluon emission amplitudes include an additional singular configuration compared to the qq case, the z 1,2 singularity structure remains unchanged. The additional "single-soft" divergence, which is absent in qq emission, is accounted for by an additional factor −1 that originates from the IBP reduction, and thus the complexity of the z 1,2 integrations remains unchanged. We present our results for the soft function in the next section.

Results
We now present our final result for the bare soft function S (2) through O 2 at NNLO QCD. To this end, we write A + C a T F n f S

Conclusion
We computed all NNLO zero-jettiness beam functions and the soft function expanded through O( 2 ) using soft and collinear limits of QCD amplitudes, reverse unitarity and IBP relations.
Our results provide one of the building blocks for calculating the N 3 LO soft function and beam function matching coefficients; some results for the beam functions described here have already been used in Ref. [26]. While the N 3 LO QCD computations of beam functions [26,27] are the first steps towards implementing zero-jettiness slicing to describe color-singlet production in hadron collisions, a significant amount of work remains to be done. Indeed, in addition to going beyond the large-N c approximation other matching coefficients I q i ,g , I q i ,q j , I g,g and I g,q i have to be calculated. Furthermore, the N 3 LO zero-jettiness soft function is currently unknown. Since the computation of the soft function is complicated by step functions in the phase-space measure, it is important to understand how to connect it to modern computational methods that involve IBP reductions and differential equations. The method discussed in this paper is a first step in that direction.