Leading-logarithmic threshold resummation of the Drell-Yan process at next-to-leading power

We resum the leading logarithms $\alpha_s^n \ln^{2 n-1}(1-z)$, $n=1,2,\ldots$ near the kinematic threshold $z=Q^2/\hat{s}\to 1$ of the Drell-Yan process at next-to-leading power in the expansion in $(1-z)$. The derivation of this result employs soft-collinear effective theory in position space and the anomalous dimensions of subleading-power soft functions, which are computed. Expansion of the resummed result leads to the leading logarithms at fixed loop order, in agreement with exact results at NLO and NNLO and predictions from the physical evolution kernel at N$^3$LO and N$^4$LO, and to new results at the five-loop order and beyond.


Introduction
The weak-coupling expansion of QCD high-energy scattering fails near kinematic thresholds due to the restricted phase space for real emission. The logarithmic enhancements in the kinematic variable that characterizes the threshold must be resummed to all orders in the coupling expansion to arrive at a reliable approximation. This has been studied first [1,2] and in greatest detail for the simplest such situation, the production of a single uncoloured particle DY (Drell-Yan process) in the collision of two hadrons, A(p A )B(p B ) → DY(Q)+X, where X denotes an unobserved QCD final state. The DY process has always provided the first physically very relevant case on which to push the accuracy of resummation to the next level, or explore new approaches to resummation [3].
The DY spectrum dσ DY /dQ 2 is given by the convolution of parton distributions in the incoming hadrons with partonic short-distance cross sectionsσ ab in partonic channels ab. The parton scattering cross sections can be regarded as functions of z = Q 2 /ŝ, wherê s = x a x b s is the partonic center-of-mass (cms) energy squared, and In this expression the series with coefficients c n , c nm encompass the leading power (LP) singular terms, and, more specifically, the terms c 0 and c n(2n−1) constitute the leading logarithms (LL). The terms multiplied by d nm are suppressed by one power of (1 − z) and are referred to as next-to-leading power (NLP). The NLP LL series is given by the highest power NLP logarithms with coefficients d n(2n−1) for n = 1, 2, . . .. Existing approaches to soft gluon resummation of the DY threshold apply only to the LP terms. The key result is the factorization of the partonic cross section into the product of a hard function and the DY soft function [4] S DY (Ω) = dx 0 4π expressed in terms of Wilson lines, as defined below. Both functions depend on a renormalization scale µ. This dependence is important to perform the resummation via a renormalization group equation, but will not be indicated explicitly unless necessary. In principle it is possible to sum arbitrary subleading logarithms at LP by computing the hard and soft function and the evolution equation to sufficiently high order. Presently, LP logarithms can be summed to the next-to-next-to-next-to-leading logarithmic order [3,5].
In contrast, much less is understood at NLP. The structure of NLP logarithms has recently received increased interest with explicit calculations at fixed order n = 1, 2 using the method of region approach [6][7][8] and diagrammatic factorization techniques [9][10][11]. However, an all-order resummation has not yet been performed for NLP threshold logarithms in the Drell-Yan process.
In the present work we accomplish this task for the leading logarithmic terms at NLP for the quark-antiquark production channel of the DY process within the framework of soft-collinear effective theory (SCET). This framework has the distinct advantage of providing precise operator definitions of the factors appearing in the factorization formula extended to NLP, which converts the resummation problem into a renormalization problem of SCET operators and soft functions. In section 2 we discuss the factorization formula for the DY process at NLP. While the form of the result is rather intuitive, this section draws heavily on a companion publication [12], where this formula is derived in detail and verified at the one-loop order. The core result of NLP LL resummation is contained in section 3, which identifies the sources of NLP LLs and derives the hard, soft and collinear functions needed for resummation, as well as the renormalization-group equation and its LL solution. In section 4 we expand the resummed result in α s , which provides both, a check of the resummed result by comparison with existing fixed-order expressions, and the JHEP03(2019)043 so far unknown logarithmic terms at higher order. The technical details of the one-loop anomalous dimension calculation for the soft functions and the derivation of an evolution equation for the NLP partonic cross section are given in two appendices.
In the following we summarize some general results on the DY threshold at NLP from [12]. We first note from (1.2) that the LP factorization formula does not contain a collinear function with threshold momentum scaling for two reasons: 1) Due to threshold kinematics generic collinear modes cannot be radiated into the hadronic final state. 2) Virtual collinear loops are scaleless, because a threshold-collinear scale can only be formed by attaching a soft momentum to the collinear loop. However, after the soft-gluon decoupling transformation [14] of the collinear fields, there are no soft-collinear interactions in the leading-power SCET Lagrangian.
The factorization of the DY threshold at NLP is obtained by matching the coupling to the DY particle to higher order in the (1 − z) expansion and by including subleading interactions from the SCET Lagrangian as perturbations. The decoupling of collinear and soft fields in the LP Lagrangian guarantees that the amplitude factorizes into a convolution of hard, collinear and soft amplitudes. However, reason 2) no longer holds at NLP, since the NLP SCET Lagrangian contains soft-collinear interactions. The time-ordered products with the hard vertex inject soft momentum into the collinear loops, resulting in a nonvanishing collinear function at the amplitude level.
The DY spectrum for the production of a lepton-antilepton pair with invariant mass Q through a virtual photon can be written as (leaving out quark electric charge factors) (2.1) 1 The above classification of regions refers to the 'standard' treatment of factorization near threshold, which neglects momentum regions which lead to scaleless integrals in dimensional regularization. If one aims at the correct identification of ultraviolet and infrared singularities, the situation is more subtle, as discussed in [17]. The whole effect of the additional collinear-soft region introduced in this reference is to convert the IR singularities in the soft function of the standard treatment into UV singularities, ultimately leading to the same factorization formula at leading power. If this holds at LP, it must also hold at NLP, since the additional region interpolates between the soft scale and the scale of the parton distributions, which are the same leading-twist parton distributions at NLP in the threshold expansion.

JHEP03(2019)043
The factorization formula forσ ab (z) for z → 1 covering power-suppressed terms in (1 −z) iŝ The formula holds for each partonic channel ab, and we dropped the parton indices ab for notational convenience. We also used the same symbol i for four separate index sets i j and a sum over various terms of this form is implicit. This will be made more precise in the following section for the terms relevant to LL resummation. S(x; ω i ,ω i , ω i ,ω i ) denotes a generalized multilocal soft function, which depends on the soft momenta radiated from the collinear and anticollinear functions. The coefficient function D(−ŝ; ω i ,ω i ) is defined at the amplitude level, and summarizes convolutions of the hard functions with the initial-state collinear and anticollinear functions at the amplitude level: Here C(n + p i , n −pi ) is the momentum-space coefficient function of a general two-jet operator as defined in [18]. Beyond LP, its Fourier transform can depend on several momentum fractions as indicated by the unspecified generic index i. Note that colour and Lorentz indices on the hard, collinear and soft functions are suppressed for the purpose of this generic discussion. The LP factorization formula is recovered from the general formula as the special case when there are no collinear functions. In this case the index set i is empty, the convolutions over the various ω i variables in (2.2) are absent and D(−ŝ; ω i ,ω i ) → C A0 (−ŝ) ≡ C A0 (x a n + p A , x b n − p B ), where the latter denotes the matching coefficient of the LP SCET current,ψ Here χ c and A µ c⊥ (below) denote the collinear-gauge-invariant collinear quark and gluon fields that form the building blocks of general N -jet operators as defined in [18]. Eq. (2.2) then turns into an intermediate result derived in [3], which after further kinematic simplifications valid at LP leads to the LP factorization formula (1.2).
We provide some further details on the definition of the generalized soft and the collinear functions. We define the soft Wilson lines

JHEP03(2019)043
in the fundamental representation, and perform the standard soft-decoupling transformation from (anti)collinear fields with Y + (Y − ). In terms of the decoupled collinear fields, which will be used below, the current reads We then introduce two sets of soft gluon and soft quark building blocks in terms of which the NLP soft-collinear SCET quark-gluon interaction Lagrangian [16] can be written as The soft-collinear interactions in the SCET Yang-Mills Lagrangian can be rewritten in a similar way. We recall that the soft fields B + µ , q + are evaluated at the multipoleexpanded position. A corresponding expression holds for the anticollinear sector. The generalized soft functions are the vacuum matrix elements of Wilson lines from the DY current similar to (1.3) and in addition collections of B µ ± (z ∓ ) insertions, where the argument z µ ∓ = (n ± z)n µ ∓ /2 arises from the integration over d 4 z in the time-ordered product and the multipole expansion of soft fields in soft-collinear interactions. Similar multilocal soft functions have appeared as "subleading shape functions" in early applications of SCET to power-suppressed effects in semi-inclusive heavy-quark decay [19][20][21], and recently in the resummation of NLP LLs of the thrust distribution in Higgs decay to two gluons [22].
A novel feature of the NLP factorization formula for the DY process is the appearance of collinear functions at the amplitude level [12]. They are defined as the matching coefficients of a product of generic collinear fields to collinear-PDF fields in the presence of soft field operators. Equivalent definitions hold for the anticollinear sector. To this end, we construct a basis of gauge-covariant soft field operators

JHEP03(2019)043
from all the possible insertions of the subleading Lagrangian. We then collect the remaining collinear fields, and define the collinear functions J i as the matching coefficients in equations, such as where ab refer to colour and αβ to Dirac indices, and d (µ) represents a collective colour (Lorentz) index from the soft operator. In the DY process the c-PDF modes can be radiated into the final state without violating threshold kinematics, since their transverse momentum of order Λ is negligible, hence the amplitude involves the matrix element X c,PDF | . . . |A(p A ) of the above equation. After squaring the amplitude and summing over the unobserved X c,PDF particles along the beam direction, the matrix element X c,PDF |χ PDF c,βb (un + )|A(p A ) can be associated with the parton distribution function f a/A (x a ). The matching coefficient J i;αβ,µ,abd (u, t; ω) is a perturbative short-distance coefficient dominated by virtualities Q 2 (1 − z) Λ 2 . Since (2.13) is an operator equation, the short-distance coefficient can be extracted by computing the partonic matrix element 0| . . . |q(x a p A ) . These considerations generalize to all collinear matrix elements that appear upon working out the time-ordered products with the subleading SCET Lagrangian. At leading twist in the Λ/Q expansion, but to all orders in the (1 − z) expansion, only c-PDF operators with a single quark or a single gluon building block on the right-hand side of (2.13) are needed. Eq. (2.13) provides the SCET analogue and operator definition of the concept of 'radiative jet functions' in the diagrammatic approach [9,10,23]. At NLP, the first subleading power in the (1 − z) expansion, the hard DY vertex should in principle be matched to SCET current operators with up to three collinear building blocks. The general basis of these operators and their renormalization is discussed in [18,24]. In the classification of these papers the operators in question correspond to A0, A1, B1 and A2, B2, C2-type operators at orders 1, λ, λ 2 , respectively. However, a non-scaleless collinear function can arise only in conjunction with a time-ordered product, from the subleading SCET Lagrangian. Since the subleading Lagrangian insertions start at O(λ), the O(λ 2 ) suppressed A2, B2, C2-type can be dropped from the beginning. Further simplifications can be made when one is interested only in the leading logarithms, as described in the next section.

Relevant terms
For the following discussion, we adopt a frame where the transverse momenta of the colliding partons vanish, p µ a = x a √ sn µ − /2, p µ b = x b √ sn µ + /2, and write the NLP factorization formula in the schematic formσ (3.1)

JHEP03(2019)043
Every factor includingσ on the left-hand side depends on the renormalization and factorization scale µ, but the scale dependence cancels after convolution with the parton densities within a given accuracy, if all factors are computed consistently. Let us assume that µ is chosen at the collinear scale O(Qλ). Since each factor depends only on a single scale, with this choice, the collinear functions J do not contain large logarithms at any order in the expansion in the strong coupling α s . The hard (C) and soft (S) functions exhibit large logarithms of (1 − z), which we aim to resum at LL. The NLP LL series is given by the terms α n s ln 2n−1 (1 − z). The soft functions start at O(α s ) or higher, since they involve at least one soft gluon radiated into the final state. A NLP LL can be generated at one-loop only if the one-loop soft function contains α s ln(1−z) and if the product [C ⊗ J ⊗J ] 2 has an O(α 0 s ) term. We can therefore drop all terms for which the hard and collinear functions do not have tree-level contributions.
Since at least one power of λ comes with the necessary time-ordered product, at NLP O(λ 2 ), we can allow for at most one further factor λ from the hard matching of the DY current. Hence, only one or none of the two factors of C in (3.1) can be the coefficient function of an O(λ) suppressed current. The available structures arē The collinear gauge field in the B1-type currents must be contracted within the amplitude, which leads to collinear functions that start at the one-loop order O(α s ). We can therefore neglect these operators. The A1-operators could produce NLP LL terms together with a soft function generated by a single L ξ , whereas the leading-power SCET soft gluon interaction does not involve the transverse direction. The time-ordered product must then vanish in the adopted reference frame, since there is no external vector available to contract an odd number of transverse indices. We conclude that the NLP LL series arises entirely from the time-ordered products of subleading power soft-collinear Lagrangian interactions with the leading-power SCET current (2.4). Moreover, its coefficient function can be taken in the tree-level approximation at the hard scale.

Quark-antiquark channel
The possible time-ordered products of the A0-operator at NLP can be inferred from (2.11) as well as the terms in L  [16]. We now argue that many of them do not contribute LLs.
The single insertion of L YM vanishes, since the associated collinear function depends on a single transverse index, but an external transverse vector is not available in the chosen reference frame. As expected, there is no O( √ 1 − z) power correction. The same argument is also valid for the case of two insertions of L

JHEP03(2019)043
We turn now to the single insertions of the O(λ 2 ) interactions in (2.11). The first step, as discussed before (2.13), is to find a minimal basis of soft fields operators. This is done systematically by considering Lagrangian insertions with one soft gluon, two soft gluons, and so on. For LL resummation only contributions from operators with a single soft gluon are necessary. This is because the soft function associated with the insertion containing two soft gluons starts at O(α 2 s ). Such soft functions could contribute to the NLP LL series only if the O(α 2 s ) contribution has a ln 3 (1 − z) term, which could arise, if there is logarithmically enhanced one-loop mixing with a one-gluon soft function. No such logarithmic enhancements of the cusp anomalous dimension type are known for off-diagonal operator mixing and we assume the absence of such logarithmically enhanced mixing in the following. This argument also excludes the possibility of LL terms arising from a double insertion of L (1) ξ and the double insertion of L (1) ξq into the amplitude as the corresponding soft functions also start at O(α 2 s ). The square of the single insertion of latter interaction contributes only to the quark-gluon channel.
We restrict then to the O(λ 2 ) Lagrangian insertion with one soft gluon field. Inspecting (2.11) and the YM Lagrangian we note that the possible soft gluon structures are n + B + , ∂ ν B + ⊥µ for the collinear direction. The soft structure n + B + , contained for example in L 1ξ , can in fact be related to ∂ µ ⊥ B + ⊥µ using the equation of motion for the soft field, namely . The first, symmetric part must be proportional to g µν ⊥ , because the soft function is a vacuum matrix element of Wilson lines and the soft gluon field insertions. g µν ⊥ is then the only symmetric structure which can carry two transverse Lorentz indices. The only remaining soft structure is given by the second, antisymmetric combination i∂ ⊥ν B + ⊥µ − i∂ ⊥µ B + ⊥ν . This combination does not contribute, because by the above argument the NLP soft function with this soft field insertion must be proportional to the epsilon tensor, which is excluded by parity conservation of QCD.
We are thus left with only one collinear function proportional to a single soft building block, given by We note at this point that L 1ξ contains n − z, which in momentum space can be converted into a derivative ∂ n + p on C A0 (n + p, n −p ). But at tree-level the hard coefficient is momentum-independent, and the derivative evaluates to zero. Moreover, the insertions of the YM Lagranian can only start contributing through one loop. As discussed in section 3.1, only the tree level collinear function contribution is necessary for LL resummation, when µ is chosen at the collinear scale. Hence, the only Lagrangian insertion contributing

JHEP03(2019)043
at LL accuracy is given by L (2) 2ξ , since it is the only piece which contributes at tree level to the collinear function in (3.4).
We therefore find that for LL resummation at NLP in the quark-antiquark channel only the single time-ordered product needs to be considered when the renormalization scale µ is chosen at the collinear scale.
To NLP LL accuracy the matching equation (2.4) is then simply extended tō We shall consider explicitly the insertion on the incoming collinear quark line. There is the corresponding term, where the insertion is placed on the anticollinear antiquark line. The NLP factorization formula (2.2) for the qq channel is obtained after extracting the collinear function from the above equation, followed by squaring the amplitude. One also needs to keep track of a kinematic correction, which arises from evaluating the d 3 q, d 3 x integrals in (2.2) with NLP accuracy.

(Anti)quark-gluon channel
For gq → γ * (→¯ ) +q to occur near threshold, the PDF-collinear gluon needs to be converted into a collinear quark carrying almost all its momentum by emission of a soft antiquark. This process vanishes at LP, but can be realized at NLP by inserting L (1) ξq once at the amplitude level. The matching equation for the (anti)quark gluon channel is The cross section follows from the interference of this amplitude with itself. A non-vanishing interference requires that the two L (1) ξq insertions are either both on the collinear quark line, or both on the anticollinear antiquark line. No kinematic corrections need to be considered in this channel, since there is no LP amplitude.

Gluon-gluon channel
The gluon-gluon channel at threshold does not contain NLP leading logarithms for the production of a lepton-antilepton pair through a virtual photon that couples only to quarks. As shown above, NLP LLs could come only from the A0 quark-antiquark SCET current, but to turn the external PDF-(anti)collinear gluons into (anti)collinear (anti)quarks, would require at least four insertions of the L suppression relative to the LP qq partonic channel. We do not discuss here the case of Higgs production for which this channel is relevant for the leading NLP logarithms.
In the following we focus on the quark-antiquark production channel and defer the further discussion of the quark-gluon channel to future work.

Collinear functions
The generic collinear modes with virtuality Q 2 λ 2 are integrated out at the amplitude level by matching collinear field operators to PDF-collinear fields. In this section we define the collinear functions relevant to NLP LL resummation and calculate the tree-level coefficient functions.
So far we have shown that the only time-ordered product with insertion of L 2ξ contributes to LLs when µ is chosen at the collinear scale. We collect the collinear fields from the time-ordered product (3.5) in the operator where α is a Dirac index, and Latin indices (except for c) are colour indices. The derivative prefactor (in − ∂ z ) 2 is conventional and conforms with the definition of the soft field product in (3.18) below as well as being consistent with equation of motion in equation (3.3). We find it convenient to define the analogue of the matching equation (2.13) in momentum space with the soft field already stripped off. Introducing and the Fourier-transform of the PDF-collinear quark field, the matching equation takes the form J 2ξ;α,ade (n + p, ω) = d(n + p ) J 2ξ;αβ,abde (n + p, n + p ; ω)χ PDF c,βb (n + p ). (3.12) The collinear function J 2ξ;αβ,abde (n + p, n + p ; ω) is defined as the matching coefficient in this operator equation. The matching coefficient is governed by the large scale Qλ Λ. We determine it by calculating the 0| . . . |q(p a ) matrix element of the above equation. At tree level we find J 2ξ;αβ,abde (n + p, n + p ; ω) ≡ J 2ξ;αβ,abde (n + p; ω)δ(n + p − n + p ) For time-ordered products originating from A-type operators, which contain only a single collinear building block by definition, collinear momentum conservation implies that one can always extract a delta function, or derivatives of the delta function from the collinear function, as done above. We can use the colour-Fierz relation with SU(N c ) generators T A in the fundamental representation to write

JHEP03(2019)043
which is also the most general colour decomposition. Since from the definition the indices de are contracted with the matrix soft gluon field B ± ⊥µ;de , the first term will never contribute. General arguments also imply that this collinear function is diagonal in the Dirac indices, hence to all orders in perturbation theory only a single scalar jet function defined through J 2ξ (n + p; ω) arises from the time-ordered product (3.5). We now calculate the matrix element of (3.6) after the (second) matching step that makes the collinear function explicit. The integrals over t,t can be performed by introducing the momentum-space hard function. After some manipulations, we obtain The above expression includes the LP term and the LP matching on the anticollinear antiquark leg. The new NLP contribution appears in the second-and third-to-last line. The corresponding power-suppressed contribution from the insertion on the anticollinear antiquark leg is thec-term not written explicitly. Making use of the delta function in the collinear factor and its colour and spin decomposition, we can simplify (3.16) to Eq. (3.16) as it stands is still valid to all orders in perturbation theory. In the tree approximation for the hard and collinear functions, which is sufficient to resum the NLP LLs, we can further set C A0 (n + p a , −n − p b ) → 1 at the hard scale and J into the (anti)quark parton distribution after squaring the above amplitude and summing over the hadronic final state X, which on the right-hand side has been split into its soft and PDF-(anti)collinear components.

Soft functions
The soft functions are defined after squaring the amplitude and summing over the soft final state X s . We introduce the soft operator 18) and the Fourier transform of its (colour-traced) vacuum matrix element To see how this generalized soft function appears, we recall some standard steps in the derivation of the factorization formula for the DY process. We first express the phase-space delta function (2π) 4 δ(p A + p B − q − p Xs − p X c,PDF − p Xc ,PDF ) in terms of the space-time integral of the exponential, then use the translation operator to absorb the hadronic final state momenta into translations of the field arguments, then square the amplitude. At this point, the sum over the PDF-(anti)collinear state can be performed, and the matrix element of the PDF-(anti)collinear fields expressed in terms of the parton distributions, (3.20) A similar standard definition applies to the anti-quark distribution in hadron B. 2 Applying these steps to (3.17), performing the integrations over n + p a , n − p b and stripping off the convolution with the parton distribution functions, we arrive at 2ξ (x a n + p A ; ω) S 2ξ (x, ω) +c-term . respect to x 0 of S 2ξ (x, ω) | x=0 . The factor of two in (3.21) arises from the two identical (see appendix A.2) NLP terms in the square of the amplitude. We introduced the hard function which is the same for the LP and NLP term, and indicated the dependence on the hard renormalization scale for later convenience. While the collinear function J 2ξ (x a n + p A ; ω) is non-zero at tree level, the soft function starts only at O (α s ). The straightforward one-loop calculation gives where Eq. (3.23) exhibits a divergence despite being the lowest-order contribution. The divergence can be interpreted as mixing into the soft function S x 0 , defined as In position space, this is simply the LP soft function with an extra factor −2i/x 0 , which leads to O(λ 2 ) power suppression and the presence of the θ(Ω) factor in the tree-level expression S x 0 (Ω) = θ(Ω), required to cancel the divergence. At first sight, such soft functions might appear peculiar. However, similar objects with collinear fields were required in the renormalization of subleading gluon jet functions [25] and the above is the positionspace and Drell-Yan process equivalent of the "θ-soft functions" introduced in [22] for the subleading-power thrust distribution. We renormalize the soft functions by writing where ω i , ω j denote (possibly empty) sets of continuous variables that parameterize the non-locality of the soft function beyond the dependence on Ω. In general, the number of arguments ω i can be different than the number of ω j , and the integration is over all ω j that the bare soft function depends on. If it depends only on Ω, the integration dω j can be omitted. Explicitly, the 2ξ soft function satisfies

JHEP03(2019)043
The mixing term in the second line of (3.27) subtracts the divergent part of the first term on the right-hand side, resulting in a finite, renormalized soft function at O(α s ). The complete one-loop anomalous dimension matrix for the above soft functions required for LL resummation at NLP is derived in appendix A.2.

Kinematic corrections
Kinematic corrections arise from the evaluation of the second line of (2.2) with NLP accuracy. In the partonic center-of-mass frame x a p A + x b p B = 0, the three-momentum p Xs of the soft hadronic final state is balanced by the lepton-pair, q + p Xs = 0. This implies the counting q ∼ λ 2 , q 0 = √ŝ + O(λ 2 ). The energy of the soft hadronic final state is expanded as with We then find for the second line of (2.2) the approximation valid to NLP. It is understood that x = 0 is set after the spatial derivatives are taken. The result is general and holds for any soft function. However, at NLP we need to apply the kinematic correction only to the LP term. Eq. (3.21) therefore simplifies tô 2ξ (x a n + p A ; ω) S 2ξ (Q (1 − z), ω) +c-term . (3.33) Since S 0 (x) = 1 at tree level, the derivative soft function starts at O(α s ). For dimensional reasons x 0 ∂ 2 x S 0 (x) | x=0 ∝ 1/x 0 . Hence, like the soft function S 2ξ , the derivative soft function times x 0 mixes into S x 0 and produces leading logarithms at NLP. Further details on the renormalization of S 0 (x) and ∂ 2 x S 0 (x) | x=0 are provided in appendix A.1. The kinematic corrections contained in the first line of (3.33) can be made more explicit. First, additional corrections arise from expanding the hard matching coefficient , but these terms do not contribute to the NLP LL series. Indeed, H (Q 2 ) starts with O(α s ln(1 − z)), but the tree-level LP soft function is δ(1 − z), which sets this term to zero. Any further term arising from the product of (1 − z)H (Q 2 ) with the LP soft function is explicitly at most of NLL accuracy. Hence we can set H(ŝ) → H(Q 2 ) in (3.33). Another implicit kinematic correction arises from the (1 − z) 2 term in the definition (3.31) of Ω * . Defining the sum of all three kinematic corrections is finite, and hence there is no leading logarithm. The diagonal renormalization of all three kinematic soft functions involves the same cusp anomalous dimension, since they descend from the same S 0 (x). The general structure of the renormalization group equation then implies that the cancellation of kinematic NLP leading logarithms for ∆ qq (z) (but notσ qq (z)) holds to all orders in perturbation theory, see appendix A.1 for further details.

Resummation
We are now in the position to sum the leading logarithms at NLP to all orders in the α s expansion. The logarithms arise from the ratio of the scales involved in the process. We shall sum the logarithms by evolving the hard function from the hard scale µ h ∼ Q and the soft functions from the soft scale µ s ∼ Q(1 − z) to a common scale µ c ∼ Q √ 1 − z using the renormalization group equations (RGEs) for the hard and soft functions. Choosing µ c to be of order of the collinear scale, we do not need the RGE of the collinear function. We shall consider the expansion of ∆(z) =σ(z)/z, since, as discussed above, the kinematic corrections cancel for this quantity, which simplifies the discussion.

JHEP03(2019)043
The hard matching function H(Q 2 , µ) obeys the RGE which follows from the anomalous dimension of the LP SCET A0 operator. 3 Here and α s denotes the MS QCD coupling at the scale µ. The general solution to (3.43) reads where [3] S(ν, µ) = − αs(µ) To LL accuracy, a Γ and a γ can be set to zero, and S(ν, µ) evaluated with the one-loop approximation to the cusp anomalous dimension and the beta function To evolve the soft function S 2ξ we have to solve the coupled system of RGEs derived in (A.47) of appendix A.2,

JHEP03(2019)043
Here µ does not have to be chosen of order of the soft scale, in which case the solution sums the leading large logarithms ln(µ/µ s ) to all orders. With the LL evolved hard and soft function at hand, we can proceed with the evaluation of the NLP partonic cross section (3.33). We write this equation for ∆(z) =σ(z)/z, in which case we can drop the kinematic correction at LL accuracy, hence (3.33) simplifies to 2ξ (x a n + p A ; ω, µ c ) S 2ξ (Q (1 − z), ω, µ c ) +c-term . 2ξ (x a n + p A ; ω, µ) = −2/(x a n + p A ) and using x a n + p A = Q + O(Q(1 − z)) in the NLP term, we find At this point thec-term takes an identical form to the second term in the curly bracket and the two can be added. In effect, the four insertions of the L 2ξ Lagrangian on the four external legs of the squared amplitude all contribute the same amount to the leadinglogarithmic next-to-leading power correction. We refer to appendix A.2 for the details of the argument.
With the explicit LL solutions (3.45), (3.51) for H(Q 2 , µ) and S 2ξ (Ω, ω, µ), respectively, inserted, the above equation reads where we also used H(Q 2 , µ h ) = 1 + O(α s ). ∆ LL LP (z) represents the LL-resummed leadingpower partonic cross section, in the present formalism given in [3]. We can set µ h = Q, µ s = Q(1−z) and µ c = Q √ 1 − z, since the precise choice is irrelevant for the leading logarithms. Note, however, that (3.54) is not of the most general form, since it implies that the factorization scale µ is set to µ c = Q √ 1 − z in the parton distributions. In order to translate the result to arbitrary µ we use the evolution equation for the partonic cross section, where . For the qq channel, no mixing terms in the splitting kernels need to be considered up to O(λ 2 ), because P ab starts at O(λ 2 ) for a = b, and the contributions fromσ ab for ab = qq yield additional power-suppression, starting at O(λ 2 ) for the qg channel. For the quantity ∆(z) =σ qq (z)/z this implies where we used also P qq = Pqq.
In appendix B we derive from this equation the evolution equation for the NLP part ∆ NLP (z, µ) in the expansion ∆(z, µ) = ∆ LP (z, µ) + ∆ NLP (z, µ) + . . . , and find with the inhomogeneous term given by Since we are interested in LL accuracy only, all terms in the square bracket of (3.59) except for Γ cusp (α s ) ln(1 − z) can be dropped, as well as the second term on the right-hand side of (3.60), that arises from the kinematic correction. The solution of (3.59) under this approximation is given by We use the LL approximation to this solution to evolve the NLP term in (3.54) from µ c to an arbitrary µ. To calculate the inhomogeneous term to LL accuracy, we use the expression for the LL resummed LP partonic cross section, where η = 2a Γ (µ s , µ) and the expression must be understood as a distribution [3]. Recalling that the expansion of S LL yields JHEP03(2019)043 double logarithms at every order, while that of η produces only single logarithms, we obtain from (3.60) To leading-logarithmic accuracy which also implies S LL (µ h , µ c ) − S LL (µ s , µ c ) = 0 in (3.54). Hence the two terms on the right-hand side of (3.61) are given by where we usedŜ(z, µ , µ) +Ŝ(z, µ c , µ ) =Ŝ(z, µ c , µ). Adding the two terms we find which is identical to the NLP term in (3.54) except that µ in the above expression is not restricted to the collinear scale, but can take any value. The scale of the parton luminosity that multiplies the above expression is now manifestly independent of z, and the logarithms of (1−z) are generated by setting µ s ∼ Q(1−z). The fact that the form of (3.54) and (3.68) are identical implies that the collinear function cannot contain leading logarithms when evaluated at a scale µ different from µ c . The simple result (3.68) for the summed NLP LLs is the main result of this paper. Let us briefly summarize the main steps of the derivation. Starting from the NLP factorization formula (3.1), we obtained the explicit LL-accurate form (3.33), which contains NLP terms from a kinematic correction and the convolution of a radiative jet function with a generalized soft function. We then noted that the kinematic correction cancels for the quantity ∆(z) and obtained (3.53). The LL resummed NLP correction (3.54) follows from inserting the LL resummed generalized soft function (3.51) and the well-known hard function. The final step to arrive at (3.68) consisted in evolving (3.54) from the collinear scale µ c to an arbitrary µ, which turned out to preserve the form of the resummed result.

JHEP03(2019)043 4 Fixed-order expansion: predictions and checks
We expand the resummed NLP cross section in the strong coupling. This serves both, as a check of the equation by comparing the logarithms at fixed order with the known NLO and NNLO partonic cross sections and provides new results beyond these orders.
To generate the fixed-order logarithms from (3.68), we expand the ratios of the running strong coupling into a series in α s (µ) and logarithms. To LL accuracy, we may approximate The scale of α s on the right-hand sides of these equations is not specified, since its precise value is a NLL effect. The NLP term in (3.68), with µ taken to be the free renormalization and factorization scale as discussed above, then reduces to Interestingly, for the special choice µ = µ c , the two exponentials cancel to LL accuracy and the NLP LL series becomes trivial, In this case, the z-dependence of the hadronic cross section is hidden in the scaledependent parton densities, evaluated at the collinear scale. For arbitrary µ, we use (4.2) with µ h = Q and µ s = Q(1 − z), and define L µ = ln(µ/Q). Then where (log) 11 stands for some combination of the two logarithms to the 11th power. The first two lines can be compared to known exact results [26]. 4 In particular, the NLO result in the first line agrees with eq. (B.29) in [26], and the NNLO contribution

JHEP03(2019)043
in the second line with eq. (B.31) in [26], up to subleading terms in the expansion in logarithms. The N 3 LO term confirms the conjecture [27,28] that the leading logarithm at this order can be simply obtained from including the NLP term in the Altarelli-Parisi splitting kernels in the standard LP resummation formalism. The N 3 LO and N 4 LO terms for L µ = 0 have been given in eqs. (B.4) and (B.5) of [29] based on the observation that the "physical evolution kernels", which express the scale dependence of a given quantity in terms of the quantity itself, exhibits only single logarithms as z → 1 to the respective order in the α s expansion. Our direct derivation from the all-order resummed partonic cross section agrees with these results as well. The N 5 LO term is a new result and the expansion to any order can be obtained without effort from (4.2).

Summary
Historically, soft gluon resummation was first studied for the threshold of the DY process [1,2] and then extended to increasingly higher logarithmic accuracy at leading power in the expansion in the threshold variable (1 − z), based on the factorization into a hard and soft function. In this work, we considered the next-to-leading power (NLP) in (1 − z) and provided an all-order resummation of the leading NLP logarithms of the form α n s ln 2n−1 (1− z), n = 1, 2, . . . near the kinematic threshold z = Q 2 /ŝ → 1.
This result is made possible by the systematic treatment of the threshold in the framework of position-space soft-collinear effective theory, which implies a generalized factorization formula including collinear functions at the amplitude level beyond leading power. We sketched and explained the factorization formula, which will be further elaborated in [12], and derived a simple expression for the leading logarithmic terms at NLP in the quark-antiquark production channel. Our main results are the resummed partonic cross sections (3.68) and (4.2), which follow from the solution to the renormalization group equations for the hard function and generalized subleading-power soft functions. The final result is stunningly simple. When the parton distributions are evaluated at a factorization scale of order Q √ 1 − z, the LL series terminates at O(α s ). For general factorization scale, re-expansion in α s shows agreement with known exact results at NLO and NNLO and with predictions from the physical evolution kernel at N 3 LO and N 4 LO, and leads to new results at the five-loop order and beyond.
We are aware of two other resummations at subleading power, both also to LL accuracy. The thrust distribution in Higgs decay to two gluons for (1−T ) → 1 was recently considered, also in the SCET framework [22]. There are differences of conceptual nature in the collinear physics, as thrust is infrared-safe, while the parton distributions must be factorized for the DY process, which requires the introduction of the PDF collinear modes and the matching of the collinear functions at the amplitude level. Nevertheless, at the LL level, the NLP resummation has a very similar structure and the result a similar level of simplicity. In particular, the mixing pattern of the required subleading soft functions is similar, and the "theta soft-function" (here the 1/x 0 soft function) appears in both cases. A number of quark-mass suppressed form factors has been resummed with LL accuracy in [30,31] through a diagrammatic all-order analysis. Here the NLP suppression is provided by the JHEP03(2019)043 quark mass m Q rather than the kinematics itself. It would be interesting to reproduce and generalize this result in the SCET framework.
All NLP resummations are presently restricted to the leading logarithmic accuracy. We expect extensions of these results to NLL to reveal the full complexity of NLP resummation. Indication of this is already provided by the form of the non-cusp terms of the one-loop anomalous dimension kernels of subleading-power N -jet operators [18,24], which enter at this order. In addition, the renormalization of subleading-power soft and collinear functions must be better understood to determine the single-pole terms needed for NLL resummation. A related question is whether the convolution integrals in the ω i variables of collinear and soft functions exhibit singularities at ω i → 0, requiring extra renormalization.

A Anomalous dimensions of the soft functions
In this appendix, we consider the renormalization of the soft functions that contribute to the renormalization group equations for the NLP DY cross section in the qq channel, and calculate their one-loop anomalous dimensions. The general form of the RGE reads with the anomalous dimension matrix defined in terms of the Z-factor (3.26) as The inverse Z-factor is defined such that where δ (ω − ω ) = j δ ω j − ω j . For further convenience, we introduce the perturbative expansion of the renormalization constants where Z (n) AB ∝ α n s .

A.1 Kinematic soft functions
The basic object from which all kinematic soft functions can be derived is the LP positionspace soft function generalized to arbitrary position x, The leading-power soft function is recovered for x = 0, while a Taylor expansion in x yields power-suppressed contributions. Lorentz invariance and invariance under reparametrization of the collinear basis vectors, n − → an − , n + → a −1 n + , for any a > 0, implies that S(x) is a function of two dimensionless variables. The soft function is an example of a closed Wilson loop with two cusps, light-like segments and no crossing point, whose anomalous dimension arises only from the cusp points (see, e.g., [4,[32][33][34]). In particular, they are renormalized multiplicatively in position space. At the one-loop order in dimensional regularization with d = 4 − 2 , the bare soft function must have a simple dependencẽ on the position variable with some function f ( , u). The explicit evaluation gives (see also [35]), where we defined The renormalized position-space soft function satisfies the renormalization group equation with Γ cusp (α s ) defined in (3.44), and at the one-loop order, as follows from (A.8).

JHEP03(2019)043
S(x) is analytic around the point x 0 = (x 0 , 0, 0, 0) (corresponding to x 2 /(n + xn − x) = 1). We can therefore obtain the RGEs for the kinematic soft functions by taking the appropriate derivatives of (A.10) and their one-loop expressions from (A.8). Let us therefore define the following vector of NLP position-space soft functions where (assuming the transverse plane to be the x 1 -x 2 -plane) where now The mixing of S 3 (x 0 ), S K2 (x 0 ), and S K3 (x 0 ) into S x 0 (x 0 ) arises from the expansion and from x 0 -derivatives on L 0 , and is given by the cusp anomalous dimension to all orders. The diagonal entries of the anomalous dimension matrix are identical for all functions and equal to the one of the LP DY soft function.
We are now in the position to justify the statement made in the main text that in the sum of the kinematic corrections to ∆(z) =σ(z)/z the leading NLP logarithms cancel to all orders in perturbation theory. From (A.18) we obtain, by summing,

JHEP03(2019)043
that is, the mixing into the 1/x 0 soft function vanishes. Since the kinematic soft functions all start at the one-loop order, the evolution can only produce terms at most of the order of α s (α s ln 2 (1 − z)) n , which correspond to next-to-leading logarithms. This conclusion does not hold forσ(z) defined as in (2.1) itself, for which the relevant kinematic soft function is the sum S K1+K2 (x 0 ). For this reason, and for further use in section 3.5, it is instructive to solve the RGE system For the present case of interest S NLP (x 0 ) = S K1+K2 (x 0 ) and the off-diagonal anomalous dimension is given by ∆ = Γ cusp . Since the RGE is local in position space, for fixed x 0 it is of the same form as (3.45), except for the off-diagonal term. The general solution can be written as with S(ν, µ), a Γ (ν, µ) as given in (3.46), (3.47). a γ W , a ∆ are defined analogously to a γ with the obvious replacement of γ by the respective anomalous dimension function. The 1/x 0 soft function does not appear directly in the NLP DY cross section, but only through mixing. Also, the initial condition S NLP (x 0 , µ s ) in (A.24) is O(α s ), hence to LL accuracy, we may approximate with S x 0 ,tree (x 0 , µ s ) = −2i/(x 0 − iε) and The function S LL (ν, µ) was defined in (3.49), and ∆(α s ) = ∆ (0) α s /π + O(α 2 s ). In this approximation the momentum space solution for S LL NLP (Ω, µ) is simply given by substituting S x 0 ,tree (x 0 , µ s ) → θ(Ω).

A.2 S 2ξ (Ω, ω) soft function
Here we consider the renormalization of the soft function S 2ξ (Ω, ω) defined in (3.19), which arises from the L (2) 2ξ Lagrangian insertion. As discussed in section 3.3, this soft function mixes into the 1/x 0 soft function S x 0 (Ω). We therefore determine the one-loop anomalous dimension matrix Γ AB for A, B = x 0 , 2ξ, more precisely the terms required JHEP03(2019)043 for LL resummation, namely the cusp-logarithmic terms in the diagonal entries, and the off-diagonal entries.
The diagonal entry Γ x 0 x 0 has been shown above to be identical to the one of the leading power soft function in position space. Written as momentum-space kernel corresponding to the definition (A.1), and keeping only the cusp-logarithm part, we have where µ s is a soft scale of order Ω. Its precise value is not needed, since it only affects the non-logarithmic term. The off-diagonal term Γ 2ξ x 0 can be inferred from (3.29) to be The 1/x 0 soft function cannot mix into the non-local 2ξ soft function, hence Γ x 0 2ξ = 0. We use two different methods to extract the diagonal part Γ 2ξ 2ξ . First, we consider the renormalization of the soft operator instead of the soft function which allows us to extract Γ 2ξ 2ξ from a one-loop calculation of the matrix element g|S 2ξ |0 involving an external gluon. Second, we compute the soft function (A.31), involving the vacuum matrix element 0|S 2ξ |0 , at the two-loop level.
First method. We generalize (3.26) to the corresponding operators, The operator (S 2ξ ) ab carries open colour indices, and in general the renormalization factor (Z 2ξ,2ξ ) ab,cd has a matrix structure with respect to these colour indices. The renormalization factor of the soft function is obtained by projecting on the colour singlet part, For the leading 1/ 2 pole we find that (Z 2ξ 2ξ ) ab,cd ≡ δ ac δ bd Z 2ξ 2ξ + O( −1 ) is diagonal in colour indices, such that Z 2ξ 2ξ = Z 2ξ 2ξ + O( −1 ). Since for the purpose of LL resummation we are interested in the leading pole only, we do not discriminate between Z 2ξ 2ξ and Z 2ξ 2ξ in the following. Figure 1. Tree-level diagrams for the one-gluon matrix element of the soft operator S 2ξ . The part to the left (right) of the cut corresponds to the time-ordered (anti-time-ordered) part of the diagram, and lines labeled by n ± with in (out)-going arrow correspond to soft Wilson lines Y ∓ (Y † ∓ ). The filled square and the two solid lines connected to it stand for the soft covariant derivative and the Wilson lines contained in

JHEP03(2019)043
In order to extract Z 2ξ 2ξ we consider a matrix element of S 2ξ with a non-vanishing tree-level contribution. We find it convenient to consider a single gluon with momentum p and colour index A in the external state, g A (p)|S 2ξ |0 . The gluon can be emitted from both of the fields inside the time-and anti-time ordered part. For a general polarization vector , the tree-level matrix element is (see figure 1) The components of the polarization vector are related by For the computation of the one-loop matrix element, we find it convenient to choose a polarization vector satisfying n − = 0, such that the tree-level matrix element is proportional to p ⊥ · * ⊥ (this implies that only the left diagram in figure 1 contributes). The one-loop computation yields p ⊥ · ⊥ as well as n + terms. The latter can be eliminated in terms of the former using The computation is done in dimensional regularization for both UV and IR singularities of the on-shell matrix element. As mentioned above, it will be sufficient to extract the leading 1/ 2 pole. The relevant diagrams are shown in figures 2 and 3. They can be divided into "real" (internal gluon connecting time-and anti-time-ordered parts of the operator, see figure 2) and "virtual" contributions (internal gluon connecting two fields within the time-or antitime-ordered parts, respectively, see figure 3). In both cases, additional diagrams that vanish due to our assumption n − = 0 are not shown. In particular, this implies that the external gluon cannot be directly attached to the Wilson lines contained in B µ + . In addition, virtual diagrams with an internal gluon line in the anti-time-ordered part are not shown. They yield scaleless integrals, analogous to the leading-power soft function. First consider diagrams a − c in figure 2. In these diagrams, the ⊥ gluon from the B µ + insertion is connected directly to the external state, while the loop is formed by the gluons connecting different Wilson lines. We find

JHEP03(2019)043
Diagrams d) and e) can be shown not to lead to a 1/ 2 pole. Next we compute real diagrams that involve a triple-gluon vertex shown in figure 2 f − i. We find that the 1/ 2 pole cancels among these diagrams. There exist four additional diagrams, that are

JHEP03(2019)043
and therefore the diagonal part of the anomalous dimension matrix is found to be Γ 2ξ 2ξ Ω, ω; Ω , ω = 4 α s C F π ln µ µ s δ Ω − Ω δ ω − ω , (A.40) again omitting non-logarithmic, µ-independent contributions. This method of computation does not allow us to determine the exact value of µ s or the non-logarithmic terms associated with the single 1/ poles, because we find that the single poles depend on the momentum components of the external state. This points to the need for a better understanding of the renormalization properties of generalized soft functions that arise from time-ordered products at subleading power. The issue does not appear to be relevant to LL accuracy (the double pole). Nevertheless, to check the above calculation, we determine the logarithmic term of Γ 2ξ 2ξ by another method in the following.
Second method. We derive Γ 2ξ 2ξ in an alternative way from the two-loop calculation of the soft function S 2ξ = 1 Nc Tr 0|S 2ξ |0 . In the following, S The renormalization condition (3.26) for the two-loop matrix element, restricted to A = 2ξ, x 0 , reads where we omit convolutions with respect to Ω and ω j for brevity. We use this equation to extract the double pole part of Z 2ξ,2ξ . For that purpose it is sufficient to focus on the leading pole 1/ 3 term of the equation.
The counterterm Z 2ξ x 0 has already been determined in (3.29) by requiring that the one-loop matrix element (3.23) is finite. We further use that the renormalization of the S x 0 soft function is the same as for the leading-power soft function. Thus we have S (1) x 0 = −Z (1) x 0 x 0 S (0) where Z 2ξ x 0 and Z 2ξ 2ξ are unknown at this point. Now we use that the relevant entries of Γ AA have explicit, linear dependence on ln µ, while the off-diagonal terms depend on µ only implicitly through α s . Under these assumptions we can solve (A.2) perturbatively for Z AB and find a relation that constrains the highest poles of the renormalization factor,  We could also determine the 1/ 2 pole of Z AB , but it is not required for LL resummation. The above relation implies that the double pole part of the renormalization constant Z 2ξ to cross-check the previous method and we obtained a result that agrees with (A.39).
Combining the results we find the following RGE equation for the generalized soft functions at LL accuracy, where µ s = O(Q(1 − z)).
As discussed in the main text, the L 2ξ insertion appears four times. In addition to the insertion on the incoming quark leg, discussed explicitly in the main text, it may also occur along the n + direction ("c-term"), or within the complex conjugated amplitude contributing to the DY cross section, in which case the Lagrangian insertion occurs within the anti-time-ordered part ("T -term"). Finally, there is a contribution from L (2) 2ξ along the n + direction within the conjugated amplitude ("cT -term").
These contributions can be treated in a way analogous to the one discussed in the main text. Thec-term involves the soft operator The corresponding soft function is Sc 2ξ (Ω, ω) = 1 Nc Tr 0|Sc 2ξ (Ω, ω) |0 . Its value at the oneloop order is identical to (3.23) up to a sign, such that Γc 2ξ x 0 = −Γ 2ξ x 0 . To determine the corresponding diagonal part of the anomalous dimension, we proceed as in the first method outlined above. However, we find it convenient to use a polarization vector for the external gluon with n + = 0 and n − = 0. The relevant Feynman diagrams can then be obtained from figures 2 and 3 by interchanging n − ↔ n + as well as the arrows, and replacing z − → z + . We refer to the labels of the corresponding "flipped" diagrams below, and describe the changes. Due to the different operator ordering, one needs to replace C F → C F − C A /2 in a) and b). Diagram c) flips its sign, since the virtual gluon is attached to Y † + (x 0 ) instead of Y − (x 0 ). These changes compensate in the sum of a − c, which gives