Threshold factorization of the Drell-Yan process at next-to-leading power

We present a factorization theorem valid near the kinematic threshold z=Q2/ŝ→1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ z={Q}^2/\hat{s}\to 1 $$\end{document} of the partonic Drell-Yan process qq¯→γ∗+X\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ q\overline{q}\to {\gamma}^{\ast }+X $$\end{document} for general subleading powers in the (1 − z) expansion. We then consider the specific case of next-to-leading power. We discuss the emergence of collinear functions, which are a key ingredient to factorization starting at next-to-leading power. We calculate the relevant collinear functions at Oαs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}_s\right) $$\end{document} by employing an operator matching equation and we compare our results to the expansion-by- regions computation up to the next-to-next-to-leading order, finding agreement. Factorization holds only before the dimensional regulator is removed, due to a divergent convolution when the collinear and soft functions are first expanded around d = 4 before the convolution is performed. This demonstrates an issue for threshold resummation beyond the leading-logarithmic accuracy at next-to-leading power.


Introduction
The study of soft emission in the threshold regime z = Q 2 /ŝ → 1 of the Drell-Yan (DY) process A + B → γ * (Q) + X has a long history. The all-order summation of the leading-power (LP) logarithms in (1−z) was pioneered in [1,2] and was later studied using soft-collinear effective theory (SCET) methods [3][4][5]. Currently LP threshold logarithms can be resummed up to next-to-next-to-next-to-leading logarithmic accuracy [5,6]. In comparison, the structure of factorization and resummation at the next-to-leading power (NLP), that is, the next order in the expansion in (1 − z), is not very well understood.
The DY process, given it is the simplest hadron-hadron collision process, has also been the target of several new calculations at subleading power. In this direction explicit computations of partonic cross sections at NLP up to next-to-next-to-leading order (NNLO) and partly beyond were performed in the coupling expansion by employing the expansion-by-regions method [7,8] and diagrammatic factorization techniques [9][10][11][12][13]. The leading logarithmic (LL) resummation of the Drell-Yan processes qq → γ * + X and gg → H + X was first achieved using SCET methods [14,15] and soon after in the diagrammatic framework [16]. Besides the threshold regime, the analysis of subleading power corrections for DY and single Higgs production has been investigated at fixedorder for resolution variables such as N-jettiness [17][18][19][20][21][22] and the q T of the lepton pair or the Higgs boson [23,24]. The resummation of NLP LLs for an event shape can be found in [25].
The resummation of NLP leading logarithms [14,15] relies on a factorization formula that was anticipated in these papers, and is also a prerequisite for taking the non-trivial step beyond LLs. In the present work we fill the theoretical gaps and provide the details of the derivation of the factorization formula beyond LP for qq → γ * + X. The factorization formula, which achieves the separation of scales through operator definitions of the relevant functions, and its check against the existing NNLO NLP results from the expansion-by-regions approach, is the first main result of this paper. Nevertheless, it must be regarded as a formal result, because it applies to bare regularized quantities. As will be explained, when one attempts to renormalize these quantities by subtracting the divergent parts, the convolution of the various factors becomes itself divergent. This complicates the resummation of NLP logarithms beyond the LL accuracy with SCET renormalization group methods.
The second main result of this paper, already extensively used in [14,15], is the identification of collinear functions or radiative jet functions at the amplitude level in the factorization formula at NLP. We discuss their origin, why they do not appear in the well-known LP factorization formula, and provide their precise operator definition in SCET. We also calculate the collinear functions at O(α s ), which illustrates the concept at the practical level and is required for the above mentioned NNLO comparison.
The concept of a jet function radiating a soft gluon was originally introduced in [26] by way of extending the Low-Burnett-Kroll formula in QED. It has been extensively discussed in the diagrammatic factorization approach [12,13,27] for the production of a colourless final state in hadronic collisions. While these functions are closely related to the collinear functions above, since they describe the same physics, they are yet different. The collinear functions in SCET are defined at the operator level as the matrix elements of collinear fields. They are single scale objects, excluding soft contributions, and therefore appear suitable for the formulation of the NLP factorization formula.
The paper is structured as follows. In Sec. 2 we discuss the emergence of the collinear functions and provide their definition through an operator matching equation. The factorization formula valid at general subleading powers is derived in Sec. 3. We then specialize this formula to NLP and identify the relevant soft and collinear functions which appear at this order in the power expansion. We present one of the main results of our paper in Sec. 4 where we extract the collinear functions at O(α s ) through a matching calculation. We compare the result that we obtain for the factorized cross section, where we employ the newly computed collinear functions, to the expansion-by-regions results up to NNLO in fixed-order perturbation theory [12,13] in Sec. 5. In particular, we find agreement with the collinear-soft NNLO contribution, if the convolution of collinear and soft function is performed in d dimensions. In Sec. 6 we demonstrate the appearance of a divergent convolution when we expand the collinear and soft function in d − 4 before performing the convolution between the two. We conclude in Sec. 7. The NLP SCET Lagrangian and supplementary results for the NLP one-loop soft emission amplitude are provided in Appendices A and B, respectively.

Threshold dynamics and collinear functions
The object of our investigation is the partonic DY process qq → γ * [→ ¯ ] + X in the kinematic region z = Q 2 /ŝ → 1, whereŝ = x a x b s is the partonic centre-of-mass energy squared, x a , x b are the momentum fractions of the partons inside the incoming hadrons and Q 2 is the invariant mass squared of the lepton pair. The factorization of the Drell-Yan process near threshold at NLP will be conducted within the positionspace formulation [28,29] of SCET [30,31]. Four-momenta will often be decomposed using light-like vectors n µ + and n µ − satisfying n + · n − = 2 and n 2 − = n 2 + = 0, according to At the partonic level, the threshold configuration constitutes a SCET I problem, hence to capture the dynamics, collinear, anticollinear, and soft fields are required. The scaling of the corresponding momenta, written in component notation (n + p, n − p, p ⊥ ), is Q(1, λ 2 , λ), Q(λ 2 , 1, λ), and Q(λ 2 , λ 2 , λ 2 ), respectively, where Q is the hard scale of the process and λ is the small power-counting parameter given by λ = √ 1 − z . We note that thresholdcollinear modes cannot be radiated into the final state X, since there is not enough energy available in threshold kinematics.
At the hadronic level, (anti)collinear-PDF modes with transverse momentum scaling p ⊥ ∼ Λ, where Λ denotes the strong interaction scale, exist in addition to the above threshold-collinear modes. These can be radiated into the hadronic final state. The q p l A0 Figure 1: Example of a threshold-collinear loop attached to the external PDF-collinear line together with a LP (A0) current. This diagram, and the ones with more loops attaching to the (anti)collinear leg, yield scaleless integrals and therefore vanish in dimensional regularization.
ordinary parton distribution functions are defined in terms of these modes. Concretely, the c-PDF modes have momentum scaling (Q, Λ 2 /Q, Λ), whereas it is (Λ 2 /Q, Q, Λ) for thec-PDF modes. We assume that the scale Λ of the strong interaction is parametrically much smaller than the threshold-collinear scale, Λ Qλ = Q(1 − z) 1/2 . We consider power corrections in λ, but we always work at leading power in Λ/Q.
In the derivation of the LP factorization theorem the threshold-collinear fields are usually ignored, since they can be trivially integrated out. This can be traced to the soft-collinear decoupling transformation [31], which removes completely the soft-collinear interactions from the LP Lagrangian. As is well known, the LP partonic cross section is then factorized at threshold into the convolution of a hard function, which is the square of a hard matching coefficient, and a soft function, which is a vacuum matrix element of soft Wilson lines [32]. At subleading power, soft-collinear interactions remain after the decoupling transformation, resulting in time-ordered product operators [33]. Threshold-collinear loops no longer vanish, and the threshold-collinear fields must now be matched to c-PDF collinear fields. The non-trivial matching coefficients constitute the amplitude-level collinear functions. In the following we will make these qualitative statements more precise.

Leading power and decoupling
We begin the discussion by considering the purely threshold-collinear 1 loop corrections to the DY process, as depicted in Figure 1. The LP current, denoted J A0,A0 following the notation of [34], consists of a collinear quark field in the n µ − direction and an anticollinear antiquark field in the n µ + direction. The first important observation is that on-shell such loops yield scaleless integrals, which vanish in dimensional regularization.
In order to obtain non-vanishing corrections, the introduction of an additional scale is necessary, for example, through the injection of a soft momentum. This is possible in threshold kinematics since the final state is composed of soft radiation. To this end, we consider the LP SCET Lagrangian written in terms of standard SCET fields, where the quark part is written explicitly as it will serve as an example. In this form, soft-collinear interactions are present at LP since The n − component of the soft gluon field is unsuppressed with respect to the corresponding component of the collinear field, resulting in the well-known eikonal form of the soft-collinear interaction. This means diagrams of type shown in Figure 2 exist. The external soft line provides a scale to the collinear loop, and indeed, individually, such diagrams are non-vanishing. Following the labeling in Figure 2, k, p, and l are soft, collinear, and anticollinear momenta, respectively. One can then form the collinear invariant (n − k)(n + p) ∼ λ 2 , resulting in dimensionally regulated results proportional to [µ 2 /((n − k)(n + p))] . It therefore appears that there should be collinear functions already at LP.
However, at LP, the decoupling transformation [31] can be applied, where the soft Wilson line is defined as

6)
2 Similar definitions apply to the collinear gluon field, and to anticollinear fields with n + ↔ n − .
The partonic cross sectionσ ab (z) factorizes into a hard function, originating from squaring the hard matching coefficient C A0,A0 (t,t ) in (2.8), and a soft function: The leading power DY soft function is given by [32]

Emergence of collinear functions
The analysis becomes more involved when subleading-power effects are studied. The framework employed here for the power-suppressed corrections in SCET was developed in [33][34][35][36]. It makes use of collinear gauge-invariant building blocks, which consist of collinear quark and gluon fields in a particular collinear direction, and non-local operators with insertions of terms from the power-suppressed SCET Lagrangian to systematically include subleading-power contributions in perturbative calculations. In what follows, we use this general framework to derive power corrections to the LP factorization formula for DY production at threshold. We find that the new physical ingredients, the collinear functions, arise from soft-collinear interactions present in the power-suppressed Lagrangian. These technically appear as a consequence of Lagrangian insertions in timeordered product operators.
As an illustrative example, we consider the insertion of the NLP soft-collinear interaction Lagrangian from (A.1). The decoupling transformation has already been performed (and the superscript (0) on the collinear gauge-invariant quark field χ c dropped), and the B ± field is a soft building block formed by a soft covariant derivative and soft Wilson lines (we also define the soft quark building block for completeness) 14) In contrast to LP, the decoupling transformation does not remove completely the soft-collinear interactions. In fact, the insertions of Lagrangian terms appear in nonlocal operators with an integral over the position of the insertion, where the field χ c (tn + ) arises from the LP J A0,A0 current. See 2ξ into a collinear quark line.
B ± (z − ) field on the other hand has dependence only on z µ − = (n + z) n µ − 2 due to multipole expansion, but this dependence links the collinear and soft fields and leads to a collinear invariant for collinear loop integrals. Concretely, consider the DY matrix element with an insertion of the above Lagrangian, Compared to the LP expression (2.8), there are extra collinear fields in the c-PDF matrix element and there is a convolution in z − between the soft and collinear matrix elements. It is precisely the presence of this extra convolution, injecting momentum with a soft scaling into the collinear matrix element, which induces a scale and leads to the emergence of collinear functions. The soft matrix element in the last line now contains an explicit gauge field insertion in addition to the Wilson lines, and will form a part of the generalized soft function. The anticollinear matrix element is the same as before, and will form part of a parton distribution function (PDF) upon squaring. We now focus on the collinear matrix element, which appears in the second line. Due to threshold kinematics the threshold-collinear modes are forbidden from entering the final state. At leading power in the Λ/Q expansion, the threshold-collinear fields in the collinear matrix element must be integrated out and matched to c-PDF mode operators consisting of a single quark (or gluon) field, which after squaring the amplitude will lead to the standard PDFs. A prototype for this matching step is the equation (refined later) where L (2) c refers to only the collinear pieces of the Lagrangian insertion. The perturbative matching coefficientJ(t, u; z − ) is the collinear function. It contains the collinear physics at the amplitude level. We stress once more that it appears first in powersuppressed corrections to the DY process. The above equation provides an operator definition of the concept of the "radiative jet amplitude" [12,13,26]. The matching should be performed in the presence of soft structures which, acting as projectors, define independent collinear functions. We give formal definitions in Sec. 2.3 below.
For the above example, we calculate the tree-level contribution to the collinear terms in the second line of (2.17). To this purpose it is convenient to introduce the momentumspace operator where we have contracted two of the collinear fields to form the collinear quark propagator, performed the z-derivatives and the integral over t, and used for the incoming quark. The z-integral can next be performed, yielding delta functions which remove the remaining integral over the momentum k. Then we find We have underbraced the matching coefficient which defines the tree-level collinear function. The appearance of collinear functions beyond LP is generic and constitutes a key concept in NLP investigations. In Sec. 4 we will calculate the one-loop corrections to these functions.

Collinear matching: formal definitions
The general collinear matching equation, suppressing indices, is given by for the collinear quark and gluon field in i-th direction, respectively. Furthermore, is a soft operator and the sum over i runs over a basis of soft structures, Here [µ, ν] denotes antisymmetrization µν − νµ, and the ellipses indicate all possible independent soft structures after utilizing the equation of motion Eq. (2.23) is a formal all-order and all-power matching equation, and it will be used extensively in the following sections. A graphical illustration is given in Figure 4.

Factorization near threshold
We now turn to the formal derivation of the factorization formula beyond leading power. We recall that the SCET derivation of factorization at LP [5] involves matching the Figure 4: A momentum-space pictorial representation of the matching equation (2.23). The oval labelled J is a collinear function. The ω i variables are conjugate to the respective positions of the insertions of subleading-power Lagrangians. Many threshold-collinear fields may join to the (possibly power-suppressed) SCET currents of the A, B, C... type [34], but there is only a single c-PDF field at leading twist in the Λ/Q expansion.
coupling to the virtual photon to the LP SCET current, where (prior to use of decoupling transformation [31]) The matching coefficient is related to the corresponding momentum-space coefficient by The fields, denoted by χ c (and A µ c⊥ further on), are the collinear-gauge-invariant collinear quark (and gluon) fields out of which building blocks of general N -jet operators are formed [34]. To derive the factorization formula valid at subleading powers, the matching equation (3.1) must be modified to include higher orders in the (1 − z) expansion. In general, this is accomplished through inclusion of all possible combinations of powersuppressed currents and subleading Lagrangian insertions. We obtain such general factorization formula in the following, before specializing to the case of NLP in later sections where we provide explicit results for the objects appearing in the factorization formula.

Factorization at general subleading powers
Omitting the index structure for clarity, the general, all power, hard matching of the vector current is given bȳ The sizes of the sets {dt k }, {t k } (and the barred sets for the anticollinear direction) for each term on the right-hand side of the matching equation depend on the type of current present in that term. Inclusion of all contributions is accounted for by the sum over indices m 1 and m 2 , which label the basis of SCET operators (and their corresponding short-distance matching coefficients C m 1 m 2 ) depending on the content of its building blocks using the formalism and notation developed in [34,35], for example m 1,2 = A0 in the LP case (3.1). Explicitly, the DY process consists of a collinear and an anticollinear direction both of which can contain sources of power suppression, hence the SCET currents are built as follows (3.5) As mentioned above, J m 2 c ( {t k }) is constructed using collinear-gauge-invariant collinear building blocks given in (2.24). In this general construction, the letter A, B, C etc. used to label the operator denotes the number of fields in a particular collinear direction, and the number 0, 1, 2 etc. denotes the overall power of λ of the current with respect to the LP, which is labelled 0. Hence, the A-type current consists of one field in one direction and derivatives of that field, the B-type current contains two fields and their derivatives, and so on. Γ m 1 ,m 2 ρ in (3.5) stands for the appropriate spinor and Lorentz structure of the operator. For instance, in (3.2) Γ A0,A0 In addition, there exist time-ordered products of currents with subleading terms in the SCET Lagrangian. These are denoted by T n, for example YM are the power-suppressed terms in the SCET Lagrangian [29]. As discussed in Sec. 2, the decoupling transformation does not remove the soft-collinear interactions in the subleading SCET Lagrangian and the injection of soft momentum into collinear loops is necessary to form non-vanishing collinear functions. For this reason the time-ordered product terms are crucial ingredients of the factorization of the DY process at NLP. To yield a non-zero subleading power amplitude at least one leg must have such time-ordered product. The other leg can then contribute to power suppression through power-suppressed currents of A, B, C etc. type or another operator containing a time-ordered product. Starting from O (λ 3 ), in addition to collinear fields, the current operator can contain purely soft building blocks [34], denoted by J s (0) here.
As discussed above, only the (c)c-PDF modes can be radiated into the final state. Eq. (3.5), however, contains threshold (anti)collinear modes, hence a second collinear matching onto a (c)c-PDF field must be performed using (2.23). The first line of (2.23) corresponds to the time-ordered product of the collinear part J m 2 c ({t k }) of a general SCET operator (3.5) with a subleading-power Lagrangian, hence applying (2.23) to the collinear and anticollinear sectors, we obtain the DY matrix element of (3.4) in the form In this equation, the index k (k) counts the number of building block fields in the collinear (anticollinear) direction within each current, and we sum over all currents. The index j (j ) refers to the number of Lagrangian insertions in the collinear (anticollinear) sector, where we also sum over all possibilities. As discussed in Sec. 2 at least one Lagrangian insertion is necessary to yield a non-vanishing subleading-power amplitude. Finally, thẽ J i (Jī) are the (anti)collinear functions and s i ({z j− }) (sī ( {zj + })) are made up of explicit B + , q + (B − , q − ) field products and their derivatives, as indicated in (2.25). 3 The further derivation of the general factorization formula follows closely the steps presented in [14] for the derivation of the NLP leading logarithmic resummation. Suppressing the m 1,2 labels, the hard matching coefficients and c-PDF fields are Fourier-transformed using and respectively. For the collinear functions we define (z j− = n + z j /2) The set {ω j } is a set of variables with soft scaling conjugate to { z j− }, and Einstein's summation convention is implied in the exponents. Equations analogous to (3.9) and (3.10) are used for the anticollinear direction. Implementing these in (3.7) we arrive at generalized version of Eq. (3.16) of [14]: For brevity, we define the coefficient function that contains both, the hard and collinear matching functions at the amplitude level. The next step in the derivation of the factorization formula is to square the amplitude, which gives the hadronic tensor W µρ defined below. Combined with the transverse lepton tensor belonging to the final-state lepton pair, for which the phase-space integrals are computed in d dimensions, we obtain an expression for the cross section (3.14) At this point we transform the c-PDF fields back to coordinate space and use the standard definition for the PDF. After performing the integrations over n + p a , n − p b and some further manipulations, we extract the convolution with the PDFs from the hadronic DY spectrum (2.10), and obtain the expression for the qq-induced partonic cross section near threshold including power corrections in (1−z) in the most general form. In the last line we introduced the generalized multi-local soft function, This concludes the derivation of the general formula for the DY cross section near threshold including power corrections. Note that these results were stated in Eqs. (2.1) and (2.2) of [14] without details, which are given here.

Factorization at NLP
We next focus on the next-to-leading power effects where certain simplifications in the general formula (3.16) can be made. We first note that since the ω variables are connected to the soft emissions from collinear functions, and therefore come from insertions of subleading-power Lagrangians in a time-ordered product, at NLP their total number is highly constrained. On the one hand, there must be at least one ω present due to the fact that at least one time-ordered product operator must appear in the SCET amplitude in order to provide a threshold-collinear scale and not lead to a trivial null result, as explained earlier in the text. On the other hand, the total power suppression at NLP is O(λ 2 ), which means that there can be at most two separate ω variables which correspond to two L (1) insertions, each contributing O(λ) suppression. The constraint on the number of subleading power interactions also limits the number of soft structures s i from the set (2.25), required at NLP. In the position-space SCET framework, soft fields in the current operators appear only from O(λ 3 ) [34]. Hence, at NLP, the soft part J s (0) is not present, and the soft structures come only from single insertions of the O(λ 2 ) SCET Lagrangian, L amplitude also cannot exist in the qq channel. 4 In consequence, at cross section level at NLP, the O(λ 2 ) suppression must be generated in the amplitude which then interferes with the LP amplitude according to (3.14), yielding the O(λ 2 ) suppressed cross section. This still leaves the possibility of O(λ 2 ) suppression to be generated by the J T 2 (t) operator formed by a L (1) insertion and a subleading current of A1 or B1-type. Due to chirality and helicity conservation in QCD, the possible currents are (3.19) and corresponding ones with power suppression in the anticollinear direction. The important detail to note is that both currents are each proportional to n ±ρ . However, the power-suppressed amplitude in which these currents could appear, is interfered with the LP amplitude, which is proportional to γ ⊥ρ , as can be seen in (3.2). Contraction of these two structures makes such contribution vanish at the cross section level to all orders in perturbation theory. This means that at NLP the sum over indices in m 1,2 in the formula derived in Sec. 3.1 contains only the A0-type current, along with time-ordered products of the LP current with Lagrangian insertions. Hence, only the hard matching coefficient C A0,A0 of the LP current appears in the NLP factorization formula. These considerations lead to the conclusion that the soft structures relevant at NLP are in fact the terms already explicitly presented in (2.25) (dropping the ellipsis) after the use of the equation of motion to eliminate the redundant n + B + structure.
The simplifications outlined above make it possible to write down a NLP version of the general subleading-power factorization formula (3.16) in a more compact way.
Namely, up to NLP (3.16) simplifies tô +c-terms , (3.20) The set notation, with {ω j } = {ω 1 , ω 2 }, is only required for terms i = 4, 5 where the soft structures consist of insertions of fields at different positions, as can be seen in the explicit expressions below. All other terms require only a single ω variable, aside from the LP position-space soft function In (3.20) the terms with power suppression placed on the anticollinear leg, both in the amplitude and its conjugate, are indicated by "c-terms" and not written explicitly, since eventually they contribute a factor of 2 to the power-suppressed terms in the above formula.
As explained above, the general structure Γ ρ defined in (3.5) is simply γ ρ ⊥ at NLP, since only the J A0,A0 current needs to be used in time-ordered products with Lagrangian insertions in the matching to the DY current. Furthermore, the anticollinear func- The index i, which is summed over in (3.20), stands in place of all indices-Dirac, Lorentz, and colour-required by each term depending on the specific soft structure appearing in the collinear matching ( Eq. (3.20) still contains the unexpanded final-state phase-space integral over the lepton-pair momentum q. This means that in addition to the dynamical power corrections to the amplitude, there is a kinematic power correction from the phase-space integration over the LP amplitude, which will be discussed in more detail below.
Next, we would like to draw attention to the collinear functions themselves. Since, as noted above, at NLP only the LP current J A0,A0 is needed in time-ordered products with Lagrangian insertions, the set {ψ c (t k n + )} in the general collinear matching equation (2.23) consists of a single quark (or antiquark) collinear field, and the set We also use momentumspace collinear functions as defined in (3.10), hence the collinear matching equation at NLP is either 23) or the one with two L (1) insertions, in which case the collinear function (soft function) has two arguments {ω 1 , ω 2 } ({z 1− , z 2− }) and the corresponding integrations must be added. The indices µ and d carried by s represent the collective Lorentz and colour indices appropriate for the given soft structure. For each independent soft structure s i there exists a corresponding collinear function J i as shown on the right-hand side of the above equation. Thus far we have focused on the derivation of the factorization formula for the bare partonic cross sectionσ. Indeed,σ still contains collinear singularities, which are usually subtracted by PDF renormalization. Therefore, care has to be taken when dealing with this d-dimensional quantity. For instance, the spin trace which appears at leading power gives a factor Tr In order to compare with with literature we find it convenient to consider the quantity ∆(z) defined through with the factor (1 − ) divided out compared to [37]. We next simplify the factorization formula in (3.20) further by discussing separately the kinematic and dynamical NLP correction.

NLP kinematic correction ∆ kin
NLP (z) In the partonic centre-of-mass frame of the DY process, where x a p A + x b p B = 0, the three-momentum of the DY boson has to be balanced by the soft radiation, q = − p Xs . The soft radiation energy is expanded in powers of λ as follows: where the first term has a further expansion in (1 − z) , Starting with the LP soft function term in (3.20), contributions to the NLP cross section come from expanding the kinematic factors. Focusing on this LP soft function term and recalling the simplification of the D coefficients for this case noted after (3.22), we start from In the above equation a number of kinematic corrections can be identified. The first is due to power suppression provided by second term in the exponent. The second, originates in the expansion of Ω * itself. The expansion of the 1/z factor gives the third kinematic correction, and a fourth kinematic correction comes from expansion of the argument of the hard function H(ŝ). After expanding out these terms, the integral over q can be performed, yielding a delta function, which sets x = 0 in the soft function. We write the four corrections in order as where S DY (Ω, x ) is the LP soft function defined in (2.12), but with argument x 0 generalized to non-zero x. The full NLP kinematic correction, ∆ kin NLP (z), is given by the sum of these four terms. In Sec. 5 we present the result of evaluating these expressions up to NNLO.

Dynamical NLP power correction ∆ dyn
NLP (z) Next we consider the contribution to the NLP cross section due to insertions of subleadingpower Lagrangians and LP kinematics. Thus we keep only the first term in the expansions (3.25), (3.26). The d d−1 q integral then gives a delta function for the spatial part of x, hence in the soft functions we can immediately set x = 0.
As opposed to the kinematic correction, the collinear functions appearing here are non-trivial. Note that the collinear functions will carry the same indices as the corresponding soft function. On top of the indices connecting to the soft function, the collinear functions carry two Dirac and two colour indices, γβ and f b, from the threshold-collinear and c-PDF fields in the matching equation (3.23). It is understood that the collinear functions, J i , in (3.20) carry indices as prescribed by (3.23). For instance, the first soft structure in the set (2.25) has one B + field and therefore carries an adjoint index A connecting to the collinear function. This means that J 1 carries one additional adjoint index corresponding to the colour generator. Explicitly, J 1 (n + p, x a n + p A ; ω) in (3.20) stands for J A 1;γβ,f b (n + p, x a n + p A ; ω). In order to simplify the ∆ dyn NLP (z) part of the factorization formula (3.20) further, we decompose the collinear functions into all possible colour and spinor structures. Continuing with the example from above, this particular collinear function must be proportional to T A f b since this is the only structure which carries one adjoint, A, and two fundamental, f b, colour indices. At this point we can define a scalar collinear function multiplied by T A f b and move the colour factor into the soft function where it forms part of the trace over the colour indices. In a similar way, the colour factors in other collinear functions can be absorbed into their corresponding soft functions. The dynamical NLP part of (3.20) can then be simplified to {dω j } J i,γβ (n + p, x a n + p A ; {ω j }) S i (Ω; {ω j }) + h.c. , (3.32) where here Ω = Q(1 − z). As in (3.20), the double-valued set {ω j } = {ω 1 , ω 2 }, is only required for terms i = 4, 5. For i = 5, in addition to the Dirac indices βγ written explicitly, J i and S i contain further indices, see the definition of S 5 below, which are contracted among them. As mentioned above, a factor of 2 in this formula comes from thec-terms. Furthermore, one of the D coefficients always reduces to the LP expression, since at NLP there is no O(λ) amplitude in the qq-channel, as discussed above. We point out again the main difference to the LP factorization formula, namely the presence of the convolution of a jet function with multi-local, generalized soft functions. 6 We define the multi-local, generalized soft functions in momentum space as the Fourier transforms The position-space soft functions appearing at NLP are given by There exists in principle another soft function, with the soft structure given by the second term in (2.25). This soft function is required to obtain the NLP one-soft-gluon emission amplitude, see Appendix B, but does not contribute to the DY cross section at any order in perturbation theory. This is because the soft functions are vacuum matrix elements of Wilson lines and soft field insertions, hence, the only structure which can carry the Lorentz indices of the anti-symmetric structure is the epsilon tensor. However, this is excluded in QCD by parity conservation. Therefore, only the five soft functions given in (3.34) to (3.38) and their corresponding collinear functions appear in the factorization formula.
The above all-order formulation of NLP threshold factorization and the operator definition of the appearing jet and soft functions is one of the main results of this paper.

Expansion up to NNLO
In Sec. 5 we will check the NLP factorization formula by comparing to existing fixedorder O(α 2 s ) results in the literature and to own expansion-by-region calculations. To prepare this discussion we consider here the terms that arise in the NNLO expansion of (3.32). Each of the objects in the formula, the hard matching coefficient C A0,A0 (n + p, n −p ), the collinear functions J i (n + p, x a n + p A ; {ω j }), the soft functions S i (x; {ω j }), has a perturbative expansion in the strong coupling. Since at NLP the generalized soft functions contain explicit soft field insertions, as opposed to simply being composed of Wilson lines as at LP, the lowest order at which they can contribute is α 1 s . The hard and the collinear functions can have tree-level contributions. This means that in order to reproduce NLO results, only one combination is needed, tree-level hard and collinear functions and a NLO soft function. Then, to reproduce the NNLO fixed order results, there are three contributions: (1) Tree-level hard function together with one-loop collinear and soft functions, (2) one-loop hard function, tree-level collinear and one-loop soft function, and finally, (3) the soft functions at O(α 2 s ). Before proceeding, it is important to note that since the kinematic set-up allows only for soft radiation, the large component of the incoming PDF-collinear momentum must be identical to the sum of the large components of the outgoing threshold momenta of the collinear function. Since for the A0 current there is only one outgoing momentum, the collinear functions relevant to NLP will be be proportional to δ(n + p − x a n + p A ). However, due to the presence of n − z in the soft-collinear interactions, which translates into a n + p derivative in momentum space, the momentum-space collinear functions can also contain derivatives of the momentum-conserving delta function. This occurs for J 1 (n + p, x a n + p A ; ω). Since it is also diagonal in the Dirac indices we write this collinear function in terms of two scalar components as follows: J 1;γβ (n + p, x a n + p A ; ω) = δ γβ J 1,1 (x a n + p A ; ω) δ(n + p − x a n + p A ) + J 1,2 (x a n + p A ; ω) ∂ ∂(n + p) δ(n + p − x a n + p A ) . (3.40) In the factorization formula, we can integrate by parts the derivative such that it acts on the amplitude hard-scattering coefficient. Hence the derivative collinear-function term contributes only when the hard matching coefficient is momentum-dependent, which happens only from the one-loop order on for C A0,A0 . Once the derivative on the coefficient function is taken, one can perform the remaining d(n + p) integral using the extracted delta function.
The only soft function from the list in (3.34)-(3.38) which begins at lowest, next-toleading order is S 1 . The others contain at least two insertions of subleading soft fields, which implies that the leading contribution to the cross section is NNLO. Therefore, expanded up to NLO, we have ∆ dyn (1) 1,1 (x a n + p A ; ω) S Moving on to NNLO accuracy, the three contributions discussed above take the following expressions: • Collinear: one-loop collinear and NLO soft functions 1,1 (x a n + p A ; ω) S In ∆ dyn (2) NLP−soft (z) the derivative terms in the collinear functions do not contribute, since the hard function is taken at tree level.
All of the above formulas can be simplified by using tree-level values for the relevant objects. In particular since H (0) (Q 2 ) = 1, we have   for the collinear term. Next, the hard contribution in (3.44) can be simplified using treelevel values for the collinear functions in (4.15) and (4.16). Care has to be taken when dealing with this expression, since it refers to d-dimensional regularized objects. The one-loop d-dimensional hard matching coefficient depends on Q 2 = x a x b n + p A n − p B only through an overall factor (−Q 2 /µ 2 ) − . Performing the derivative therefore gives back the hard matching coefficient multiplied by a factor of − /Q. Together with the hermitian conjugate term in (3.44), we obtain − /Q × (C * A0(0) C A0(1) + C * A0(1) C A0(0) ) = − H (1) /Q from the derivative term. Then we arrive at The tree-level collinear functions can also be utilized to simplify the soft term, but we do not present it here.

Calculation of collinear functions
In this section we present the computation of the collinear functions to one-loop accuracy. The presence of these functions at NLP is one of the main results of this paper, and we will need the one-loop calculation in the subsequent section to verify the NLP factorization formula to NNLO. The collinear functions are defined through the non-perturbative operator matching equation (3.23). The left-hand side includes the threshold-collinear fields originating from time-ordered products of the LP current with subleading-power Lagrangian terms. We introduce the abbreviatioñ for the left-hand side of (3.23), and define its Fourier transform by T γf (n + q) = dt e i(n + q) tT γf (t) .

(4.2)
The momentum-space matching equation reads For soft structures with two soft gluon emissions the generalization explained below (3.23) applies. We recall that since the collinear scale Q 2 (1 − z) Λ 2 by assumption, the collinear function is a perturbatively calculable short-distance coefficient in the matching of (4.1) and (4.3). We can therefore extract the collinear functions J i by taking an appropriate matrix element between partonic states. For example, in case of collinear functions with a single external soft gluon, the simplest choice is the matrix element g(k)|...|q(p) with a soft gluon and PDF-collinear quark. We then compute both sides of the matching equation with the LP collinear Lagrangian with soft fields decoupled, in which case the soft fields on both sides act only as external fields. Hence, the soft matrix element g(k)|s i;µ,d (z − )|0 takes its tree-level expression (since only soft loops could contribute). The same is true for 0|χ PDF c,βb (un + )|q(p) , because loop corrections are scaleless in this case.

Collinear functions at O(α 0 s )
For the qq-induced DY process, only the insertions of the quark-gluon subleading SCET Lagrangian but not the Yang-Mills terms contribute at tree level to the collinear functions. Indeed, at least one collinear gluon loop would be needed to which a L YM insertion could be attached via a triple-gluon interaction.
We use momentum-space Feynman rules for the soft-collinear interactions vertices from the power-suppressed SCET Lagrangian given in Appendix A of [35] to perform the computation. The collinear-quark soft-gluon interaction vertex is given by and The momentum k + , which appears in the argument of the delta function above, is defined as k µ + = (n − k) n µ + 2 . The three terms in the O(λ 2 ) vertex (4.4) correspond directly to the three terms in the power-suppressed SCET Lagrangian given in Eq. (28) of [29]. This has been rewritten in terms of gauge-invariant building blocks in (A.1) such that the first term in (4.4) corresponds to L (2) 1ξ , the second term to L (2) 2ξ , and the third term to L (2) 4ξ . In Appendix A we also provide the soft-quark and Yang-Mills SCET Lagrangian in this notation, which are needed for the one-loop calculation and for soft structures with soft quarks.
An important feature of the NLP Feynman rules is that they contain derivatives of momentum-conservation delta functions at the subleading-power vertices. This is due to the appearance of explicit position-space arguments, x µ , in the SCET Lagrangian terms owing to multipole expansion [28]. These derivatives must first be integrated by parts to act on the rest of the amplitude before imposing momentum conservation.

Single soft gluon structures
Inspection of the subleading-power SCET Lagrangian (A.1) shows that only the two soft gluon structures can have tree-level single-gluon matrix elements at O(λ 2 ). Hence the sum over i in (4.3) reduces to i = 1, 6. Explicitly, (4.3) turns into where K, q refer to the colour of the external state and the superscript 1g reminds us that we consider the collinear functions for single soft gluon emission. The c-PDF collinear matrix element on the right-hand side equals 0|χ PDF c,βb (un + )|q(p) q = δ bq Z q,PDF u c,β (p) e −i(n + p) u , (4.8) where Z q,PDF is the on-shell wave renormalization factor of the c-PDF field. The soft matrix elements are found to give Inserting these results into (4.7), we obtain This is the final expression for the right-hand side of the matching equation (4.3) for single soft gluon structures for the chosen partonic state. We note that this expression is exact to all orders in perturbation theory, since, as mentioned above, there are no loop corrections to the above matrix elements.
We next turn our attention to the computation of the left-hand side of the matching equation (4.3). The relevant terms in L (2) are L (2) 1ξ,2ξ,4ξ , which give rise to the NLP soft-gluon vertex (4.4). A straightforward tree-level calculation gives where Z q,c | tree = 1 is the tree-level value of the on-shell wave function renormalization factor of the quark field in the effective theory including the threshold-collinear mode.

Double soft parton structures
We now consider the collinear functions multiplying soft structures with at least two soft fields. In the graphical representation of Figure 4, these correspond to diagrams with one external quark to the left and right, and two external soft gluons or a soft quark-antiquark pair attaching to J. Specifically, we consider the single insertions of the L ξq , see Appendix A for the definition of these terms.
A closer inspection of the subleading-power SCET Lagrangian (A.1) shows that after the soft fields are stripped off, the remaining collinear parts of L (2) 3ξ and L (2) 5ξ are identical to those of L (2) 2ξ and L (2) 4ξ , respectively. This means that the collinear functions are the same, that is J 3,1 is equal to J 1,1 , and J 2 to J 6 .
Therefore, we only need to consider the double soft parton emission matrix elements for the double L (1) insertions. The collinear matching equation for double Lagrangian insertions is The partonic matrix elements to be calculated here are g(k 1 )g(k 2 )|...|q(p) for the soft structure with emitted gluons and q(k 1 ) k 1q (k 2 ) k 2 |...|q(p) for the structure with emitted soft quarks. The right-hand side of the matching equation is then obtained as The left-hand side is calculated as for the single soft gluon case, where for T 2g γf (T 2q γf ) we use L   Let us emphasize that the previous two equations result from one-soft-particle irreducible tree diagrams formed by contracting the collinear fields in the double Lagrangian insertion. When the equation-of-motion relation (2.26) is used to eliminate the n + B + soft field, the two soft parton collinear functions receive another contribution from the timeordered product (4.1) and the two soft field terms (gluon-gluon and quark-antiquark) in the equation-of-motion relation, which correspond to one-soft-particle reducible tree diagrams, in which the n + B + gluon from the L (2) 1ξ insertions splits into a gluon or quarkantiquark pair through a LP interaction. The additional contribution can be obtained from the one-gluon matrix element (4.12) before eliminating the n + B + structure through the on-shell and transversality relation (4.25) below.

Collinear functions at O(α s )
In this section we focus on demonstrating the consistency of the concept of collinear functions by calculating J 1 and J 6 at the one-loop level. J 1 is also the only collinear function which is needed at the one-loop order to verify the NLP factorization formula at NNLO accuracy, see (3.46). We do not calculate the loop correction to the collinear functions of the two soft-parton structures, since it is a next-to-next-to-leading order (NNNLO) effect.
The right-hand side of the matching equation has already been obtained in (4.11), which is valid to all orders in α s . The on-shell wave function renormalization factor should now be evaluated with one-loop accuracy. However, when dimensional regu- Figure 6: One of the diagrams contributing to the one-loop collinear functions. Through calculation of this diagram using Feynman rules from [35] we can obtain the J 1 and J 6 collinear functions, corresponding to insertions of L (2) 1ξ and L (2) 2ξ , and L (2) 4ξ , respectively. larization is used for ultraviolet and infrared divergences, Z q,PDF = 1 to all orders, because the loops are scaleless. 7 The coupling renormalization is also the same on both sides of the matching equation, and drops out at the one-loop order.
We therefore focus on the calculation of g(k) K |T 1g γf (n + q)|q(p) q on the left-hand side of (4.3), which requires the calculation of the Feynman diagrams with one collinear loop and a single soft emission, generated by insertions of the power-suppressed Lagrangian. The relevant SCET diagrams are shown in Figure 5. The circled vertex denotes the subleading-power Lagrangian insertion, while all other vertices are LP interactions.

Detailed computation
We illustrate the computation by considering as an example the top-left diagram in Figure 5, which we draw again with momentum labels in Figure 6. All necessary Feynman rules were provided in Appendix A of [35] and (4.4). Applying them to the diagram under consideration leads to After substituting the expression for S σδ (−k, p − p 2 , p 1 ) from (4.4) and performing an integration by parts of the derivative with respect to p 1 contained in S, the derivative acts on the integrand including the delta function in the second line. At this point, the momentum conservation delta function at the subleading-power interaction vertex can be imposed by performing the integral over p 1 . This identifies p µ 1 = p µ − p µ 2 − k µ + and results in connected via the equation-of-motion identity (2.26). We can use the transversality and on-shell conditions k · * = 0 and k 2 = 0, respectively, for the emitted gluon, which has not yet been done in obtaining (4.24). The relation k · * = 0 can be written in light-cone components as at which point we see that indeed we can express the first soft structure in the curly bracket of (4.24) in terms of the second, This is expected as we know that the insertions of L (2) 1ξ and L (2) 2ξ contribute to the same collinear function J 1 , since the soft structures are connected via (2.26).

Amplitude calculation results
The calculation of all diagrams in Figure 5 gives the following result, after using the on-shell and transversality relations: These results constitute the left-hand side of the matching equation, that is, the extension of (4.12) to one-loop accuracy. Remarkably, we find that the one-loop correction to the derivative delta-function term cancels exactly when all diagrams are added, which explains the absence of such term in the above equation.

Fixed-order results
There exists a number of NLP results for the DY process at NLO and NNLO in the strong coupling [7,[11][12][13]40], obtained from direct expansions of the QCD diagrams. In this section we verify the correctness of the NLP factorization formula by comparing to these results and own results from the expansion-by-regions method [41].

NLO
Expanding the NLP factorization formula to NLO, one finds only one dynamical contribution to the cross section, since the soft function begins at O(α s ). The one-loop soft function is given by 9

S
(1) Using this result in (3.42) and performing the convolution integral over ω gives ∆ dyn (1) where we set µ = Q.
In addition, at NLO we need to take into account the kinematic corrections ∆ K1 NLP (Ω), ∆ K2 NLP (Ω), and ∆ K3 NLP (Ω) in (3.28), (3.29), and (3.30). For the latter two, we can use (5.1) and H (0) (Q 2 ) = 1, since no derivatives with respect to the coordinate x needs to be taken. To compute ∆ K1 NLP (Ω) we use the result for the one-loop soft function with full x dependence from [14,42]. Upon summing the three kinematic corrections we obtain ∆ kin (1) Results for the NLO NLP contribution to DY production have been presented in [11] within a diagrammatic approach, in which power-suppressed soft radiation is described it terms of generalized next-to-soft Wilson lines. Our result (5.2) agrees with the corresponding expression Eq. (6.17) of [11]. The kinematic correction (5.3) is provided in Eq. (6.13) of [11] as a correction to the LP matrix element. Agreement can be easily checked. After summing (5.2) and (5.3), and applying the subtractions that arise from PDF renormalization, we also find agreement with the NLO NLP result reported in Eq. (B.29) of [40].

NNLO
In Sec. 3.3 the three possible dynamical NLP contributions to the cross section at NNLO have been discussed. These are collinear, hard, and soft contributions presented in (3.46), (3.47) and (3.45), respectively. In this section we explicitly compute and check the first two of these. The soft contribution requires a full NLP NNLO soft function computation, which is beyond the scope of this work. However, we present the one-virtual, one-real soft contribution to the cross section here. Also, in Appendix B we present complete results for the one-loop power-suppressed amplitude with one real soft emission, including the soft loop contribution. The latter forms part of the virtual-real contribution to the NNLO soft function. The missing contribution comes from double real soft emission, which we leave for future work.

Collinear contribution
This contribution comes from the one-loop collinear functions combined with the NLO soft function and tree-level hard function, see (3.46). We recall that the delta-function derivative term in the collinear function, spelled out in (3.40), vanishes after partial integration, since the hard function at tree level is a constant. The one-loop collinear function that is required is then given by (4.28) with colour generator and Dirac-index Kronecker-symbol removed.
For the purpose of deriving the NNLO fixed-order result, we keep must use the ddimensional expression of the collinear function and perform the convolution with the d-dimensional soft function. Then expanding in and setting Ω = Q(1 − z), µ = Q yields Notice that, as expected, there are no leading logarithms O(α 2 s ln 3 (1−z)) in the collinear contribution, since the highest power of the logarithm in the finite terms in the second line is NLL accuracy, ln 2 (1 − z).
Results describing virtual collinear radiation at one loop with emission of a soft gluon have been derived in [7] within the expansion-by-regions approach [41], and in [12,13] within a diagrammatic approach, in which the effect of collinear loops is described in terms of a "radiative jet function". The C 2 F term in (5.4) is in agreement with the corresponding contribution in Eqs. (13), (14) of [7] and Eq. (4.22) of [12], where the abelian contribution only is considered. 10 The C A C F term in our result (5.4) is not provided separately in literature, but only in sum with the hard and soft contribution, that we consider in the following.

Hard contribution
Next we check the contribution composed of the one-loop hard function, the tree-level collinear functions, and the one-loop soft function. In contrast to the collinear contribution, here the collinear function with the derivative contributes, since the hard matching coefficient is momentum-dependent beyond tree level.
The relevant formula is now (3.47), which already made use of the expressions for the collinear functions at tree level. The one-loop soft function was given in (5.1). The d-dimensional hard matching coefficient at the one-loop order can be found in Eq. (2.23) of [43], where Q 2 = (n + p)(n −p ). Taking care of the imaginary part, the one-loop hard function Performing the ω-integration in (3.47), setting µ = Q, and expanding in leads to where ζ(3) is a Riemann zeta value. In contrast to the NLP collinear contribution, LLs appear in this expression. Resummation of the hard function is therefore necessary in order to sum LLs to all orders in α s , as was done in [14].
In the literature, the hard one-loop times one real soft gluon result has been considered before within the expansion-by-regions method. The expression for the abelian C 2 F term has been given in Eq. (12) of [7], and agrees with (5.7). 11 Within the diagrammatic approach [12,13], the hard contribution (5.7) arises from dressing the non-radiative amplitude by a one real soft gluon, according to the LBK theorem. For a discussion of the LBK theorem in the present approach, see [44].

Soft contribution
The soft contribution provided here is the one-real, one-virtual piece of the full NNLO soft function as mentioned in the introduction of Sec 5.2. In (3.45) one can see that there are contributions to the NNLO soft function from different soft structures. However, as detailed in Appendix B.3, only one soft structure, S 1 , and corresponding tree-level collinear function actually contribute to this piece. Hence the simplified factorization formula is (Ω, ω).
The result for one-real, one-virtual contribution to the two-loop soft function reads Using (4.15) for the tree-level collinear function in (5.8), integrating over ω and expanding in yields In the literature the non-abelian C F C A term of the one-real, one-virtual contribution has been provided as a sum of the collinear, soft and kinematic correction (see Eq. (4.6) of [13]), thus (5.10) cannot be compared directly. We performed an independent calculation of the full one-real, one-virtual correction within the expansion-by-regions method, and (5.10) agrees with the soft region, as it should be.

Kinematic contribution
The kinematic correction from the sum of terms (3.28)-(3.31) can also be obtained at NNLO by using the NNLO soft function with full x-dependence presented in [42]. We find ∆ kin (2) We note that there are no LLs due to kinematic corrections. The kinematic contribution has been calculated previously within the expansionby-regions or the diagrammatic approach as the NLP phase-space corrections to the LP matrix element, but the expression corresponding to (5.11) has not been provided explicitly. (It is part of Eq. (4.6) and (5.2) of [13], but it cannot be separated from the NLP matrix element.) We thus compare (5.11) with an own independent expansion-byregions calculation, in which we take the matrix element at leading power (both one-real, one-virtual and two-real diagrams), and integrate it against the NLP phase space, finding agreement.

Ill-defined convolution
One of the primary uses of factorization formulas in SCET is to perform resummation using renormalization group equations. Soft-collinear factorization often involves convolutions C ⊗ F of hard functions with collinear factors, for example, in deep-inelastic scattering or in convolutions with PDFs for any hadronic scattering cross section, or J ⊗S of jet with soft functions, for example in the description of radiation from final-state jets. Resummation relies on defining renormalized factors by subtracting their poles in dimensional regularization and deriving a renormalization group equation for the renormalized function, which usually also has a convolution form. Large logarithms are then summed by evolving one function to the characteristic scale of the other. Finally, the convolution of the two factors is done.
This procedure evidently requires that the final convolution integral of the renormalized factors is well defined. As we discuss now, this important requirement is not satisfied by the NLP factorization formula for the DY process.
The issue is most clearly exposed when we focus on the functional form of the objects appearing in the one-loop collinear times one-loop soft NNLO term in factorization formula given in (3.46). The one-loop collinear function J It is evident that the integral is well defined when keeping the exact dependence in the integrand, as was done in the previous section in order to obtain and reproduce the fixedorder NNLO NLP results. However, as explained above, for resummation we would like to treat the parts originating in the collinear function, (n + p ω) − , and the soft function pieces, ω −1− (Ω − ω) − , independently. That is, we wish to expand each in and define renormalized functions. However, it is clear that there is a problem when this procedure is attempted in (6.1). Concretely, one encounters a divergent integral, or dω δ(ω) ln(ω) and other ill-defined integrals after introducing the standard plus distribution for the 1/ω 1+ factor. 12 In order to make the issue even more explicit, we take the -expanded collinear function given in (4.29) and also expand the one-loop soft function (5.1) in , The convolution of this expression with (4.29) according to (3.46) (at µ = Q as in the section above) gives where only the pole terms in are shown. There are two issues with this result. First, one of the terms with 1/ pole is ill-defined as we encounter the integral dω δ(ω) ln(ω). Second, the coefficient of the C 2 F / 2 , C F C A / pole terms which are not divergent have changed with respect to the correct result from (5.4) obtained from expanding in after performing the convolution in d dimensions.
It is clear from the above that it is not possible to obtain the NLP logarithms of (1 − z) correctly from the standard renormalization procedure and four-dimensional convolutions. The leading logarithms in the qq (gg) channels in DY (Higgs) production summed in [14,15] form an exception, since they require only tree-level collinear functions and since the loop corrections to the collinear functions do not contribute leading logarithms. The ill-defined convolution, however, hampers the extension of resummation to NLL. The convolution itself requires subtraction, and contributes to the logarithms, which can therefore not be obtained from the separate renormalization group equations for the renormalized collinear and soft functions. Nevertheless, the NLP formula derived in this paper factorizes the different momentum scales of the DY process consistently at the level of regularized matrix elements of the soft and collinear operators, and therefore can be justifiably called a factorization formula. It may be hoped that it provides the starting point for understanding how to renormalize d-dimensional convolutions, which appear to be a generic feature of NLP factorization. 13

Summary
In this work, we derived for the first time a factorization formula for DY production near threshold in the qq-channel at general powers in the (1 − z) expansion. We then focused on the next-to-leading power, which entails several simplifications. The main result is the NLP factorization formula (3.32), which generalizes the LL-accurate formula in [14].
As one of the new key ingredients of the subleading-power factorization formula, we identify and discuss the emergence of collinear functions at the amplitude level. While the related concept of a "radiative jet function" [26] has been known to be relevant to power corrections at the DY threshold from diagrammatic studies [12,13,27], the benefit of the present SCET treatment is an operator definition, which renders the function gauge-invariant by construction. More precisely, see (2.23), the collinear functions are the perturbative matching coefficients, when threshold-collinear fields are matched to c-PDF fields in the presence of external soft structures that describe the emission of one (or several) soft gluons. Due to the strict scale separation and systematic power expansion, the collinear functions are single-scale objects. They are extracted from partonic matrix elements since the threshold-collinear scale is assumed to be much larger than the scale of strong interactions, Q(1 − z) 1/2 Λ. The tree-level collinear function required for LL resummation in threshold DY (and Higgs) production has already been used in [14,15]. In this work, we computed the oneloop O(α s ) corrections (4.28) and (4.31) to the collinear functions, which can contribute to the DY cross section at NNLO, and to the one-loop one-gluon emission amplitude. These results confirm explicitly the observation made in [14,15] that the DY collinear function cannot contain LLs. The one-loop calculation demonstrates the validity of the definition of these NLP objects and allows us to verify the correctness of the factorization formula at NNLO by comparing its expansion in powers of α s with existing results obtained at this order with the expansion-by-regions method.
However, our investigation also highlights that factorization at NLP is not yet understood at a similar level as at LP. The factorization formula separates the scales relevant to the DY threshold in the form of well-defined, dimensionally regulated collinear and soft functions, which have to be convoluted in the soft momentum variables ω i . The O(α 2 s ) calculation makes explicit what can already be seen from general scaling arguments that the convolutions exist only for the d-dimensional functions. When the expansion in is performed before the convolution, the latter is ill-defined and leads to a divergence. This implies that the formula is not yet in a form suitable for the resummation of large threshold logarithms beyond the LLs through the renormalization group equations for renormalized hard, collinear and soft functions. Nevertheless, it may be hoped that it provides the starting point for understanding how to renormalize d-dimensional convolutions, which would open the path to NLL resummations beyond LP.

A.2 YM subleading SCET Lagrangian
The subleading-power gluon self-interaction terms of the soft-collinear Yang-Mills Lagrangian [29] expressed in terms of the collinear and soft gauge-invariant fields are given by

B One-loop single soft real emission amplitude
In the main body of the text we focused on the factorization formula at the crosssection level. As a by-product of the computation of the collinear functions, which are amplitude-level objects, we also obtained the power-suppressed one-loop one-soft emission DY amplitude, which we summarize here. The results below, computed directly in SCET, were shown to agree with in-house results obtained by applying the expansionby-regions method to the same quantity. We consider the following operator, which is the right-hand side of (3.4) without the soft current J s : as in (3.5). The variables appearing in this expression are defined in Sec. 3.1, and the sum is performed over the different power-suppressed currents in the N -jet SCET operator matched to the QCD current. Below we focus solely on the case in which the power suppression is in the collinear sector, thereby setting m 1 = A0, and allow for structures which give power-suppression up to O(λ 2 ) (NLP). Specifically, we consider the time-ordered product of J m 1 ,m 2 ρ with subleading-power Lagrangian insertions between an emitted soft gluon g(k) K |, and an incoming collinear quark and anticollinear antiquark, |q(p)q(l) . This defines the amplitude that we calculate at the one-loop order. Concretely, we consider only the time-ordered products of the collinear operator part J m 2 c in (B.2) with subleading-power soft-collinear (not: soft-anticollinear) Lagrangian insertions. The complete result for the amplitude is obtained by adding to the contributions given below the corresponding ones with n + and n − interchanged.
In the following sections we present the different contributions to this object. Partial results obtained when the virtual loop is collinear (soft) carry a subscript c (s), M c (M s ). The NLO contributions from the one-loop hard matching coefficient are marked with h, M h . Moreover, we further split the results according to the polarization of the off-shell DY photon γ * produced by the vector current, that is, we separate the amplitude into the terms proportional to γ ⊥ρ , n +ρ , and n −ρ . Notice that the γ ⊥ρ structure appears due to the LP current in (3.2), while n ±ρ terms arise from the power-suppressed A1 and B1 currents in (3.18) and (3.19), respectively.

B.1 Collinear loop: γ ⊥ρ
We begin with the results for the set of diagrams in which the virtual loop has collinear momentum scaling and the virtual photon created by the vector current has a transverse ρ index. In (B.3) this means taking the LP current, and index m spans over time-ordered product insertions of the L (2) Lagrangian. The equations below are in fact related to the results presented in (4.27) and come from calculation of the diagrams in Figure 7.
We separate the resulting expression into the amplitude with colour factor C F and C A . The former receives contributions from the diagrams in the top line of Figure 7, the latter from those in the bottom line and the non-abelian part of the last two diagrams in the top line. We find duce the corresponding virtual-real contribution from the expansion-by-regions method. Hence only non-abelian contributions arise here and we find Details on the vanishing of numerous other a priori possible diagrams are provided in Figures 10 and 11. Note that the latter figure also includes diagrams that represent insertions of both, the collinear (on the upper leg) and anticollinear (on the lower leg) subleading soft-collinear interactions, when a = b = 1. However, as all these terms vanish, there is a unique separation of contributions from collinear Lagrangian insertions and from anticollinear Lagrangian insertions. In (B.8) have we have given the a = 2, b = 0 contribution from the last diagram in Figure 11, while the a = 0, b = 2 anticollinear one is obtained by exchanging n + ↔ n − . We further note that the absence of a contribution of the second diagram in Figure 11, containing a power-suppressed two-soft gluon vertex, implies the statement made in Sec. 5.2.3 that only the single soft-gluon structures with their corresponding soft functions S 1 , S 6 contribute at NNLO, of which only S 1 is relevant at cross-section level as explained in the main text.

B.4 Soft loop: n ρ +
The relevant diagram is again the topology of Figure 9. However, since one power of λ is used up by the power-suppressed current, at the soft-collinear vertex we now insert the L (1) term from the SCET Lagrangian.
The J A0,B1 current cannot give a contribution here since it produces a collinear gluon, that cannot be contracted to form a soft loop.
The diagrams shown in Figures 10 and 11 are also present here. The only change is that the LP hard current is replaced by J A0,A1 and the sum of a + b (+ c) = 1 only. Once  Figure 10: Diagrams with one soft emitted gluon and one soft loop. Since all the diagrams here include the LP J A0,A0 current, the O(λ 2 ) power suppression must be provided by Lagrangian insertions. This means using all possible insertions such that a + b (+ c) = 2 at the indicated vertices. Out of the 20 possibilities, many vanish immediately due to contractions which yield n 2 ± = n ± · γ ⊥ = 0 or propagators which give zero due to the vanishing external transverse momentum. The remaining integrals, where the integrand does not immediately vanish, are either scaleless or vanish by Cauchy's theorem, because all propagator poles lie in one half of the complex momentum plane.
again only the last diagram in Figure 11 does not vanish, and we find M n + ρ K s,C A =vc(l) n ρ + ig s α s (4π) −(n − k)(n + k) µ 2 There is no term proportional to n ρ − .
B.6 Hard loop: n ρ + This contribution comes from the one-loop correction to the matching coefficient C A0,A1 of the J A0,A1 current together with an insertion of the O(λ) piece of quark SCET Lagrangian. C A0,A1 is related to C A0,A0 by reparametrization invariance [45]. With the definition (3.18) the relation reads C A0,A1 = −1/(n + p) C A0,A0 . We then find There is no term proportional to n ρ − .