On next-to-leading power threshold corrections in Drell-Yan production at N$^3$LO

The cross-section for Drell-Yan production of a vector boson has been previously calculated at next-to-next-to-leading order, supplemented by enhanced logarithmic terms associated with the threshold region. In this paper, we calculate a large set of enhanced terms associated with the colour structure $C_F^3$ at N$^3$LO, for the double real emission contribution in the quark-antiquark channel, as an expansion around the threshold region up to and including the first subleading power. We perform our calculation using the method of regions, which systematically characterises all contributions according to whether the virtual gluon is (next-to) soft, collinear or hard in nature. Our results will prove useful for developing general formalisms for classifying next-to-leading power (NLP) threshold effects. They are also interesting in their own right, given that they constitute a previously unknown contribution to the Drell-Yan cross-section at ${\cal O}(\alpha_s^3)$.


Introduction
The ongoing Large Hadron Collider programme, together with related experimental facilities, necessitates the calculation of scattering processes in perturbative quantum field theory to ever greater precision. The state of the art in such calculations typically evolves on two complementary fronts. Firstly, there is the calculation of specific processes at fixed order in perturbation theory (including both QCD and electroweak corrections). Secondly, there is the inclusion of successive infinite towers of kinematically enhanced contributions, and the matching of these so-called resummed predictions with fixed order results. The state of the art for most processes of interest is next-to-leading order (NLO) in perturbation theory, supplemented by next-to-next-to-leading logarithmic (NNLL) resummed contributions. A few processes are known beyond this order, and in this paper we focus on inclusive quantities in the production of heavy particles, which depend on a single ratio ξ of kinematic scales, such that ξ → 0 at threshold. Examples include the Drell-Yan production of a vector boson, which is currently known to NNLO [1][2][3][4][5][6][7][8], and the closely related process of Higgs boson production via gluon-gluon fusion, which has been calculated up to an impressive N 3 LO [9][10][11][12][13][14][15][16] in the large top mass limit. The differential cross-section in QCD for these and other single-scale quantities assumes the generic form nm log m ξ ξ + + c (δ) n δ(ξ) + c (0) nm log m ξ + . . . , (1) where K ew collects electroweak coupling and normalisation factors, α s = g 2 s /(4π) is the strong coupling, and n 0 denotes the power of the strong coupling in the Born interaction. Commencing at NLO, each order in α s is accompanied by a series of divergent contributions as the threshold variable tends to zero, associated with QCD radiation that is soft and / or collinear with the hard particles in the underlying scattering process. The first set of terms in the square bracket in eq. (1) constitutes the leading power (LP) in the threshold variable ξ, which mixes with the second set of terms, that originates also from purely virtual corrections. The third set of terms is next-to-leading power (NLP) in a systematic expansion in ξ, and formally divergent as ξ → 0, albeit integrably so. Finally, the ellipses in eq. (1) denotes higher power corrections in ξ which vanish at threshold.
The practical significance of threshold contributions is well-known, and a variety of approaches exist for resumming LP terms to all orders in perturbation theory [17][18][19][20][21][22][23] in order to obtain meaningful comparisons of theory with data. In recent years, the NLP terms in eq. (1) have also received a great deal of attention, for a number of reasons. Firstly, they can dominate the theoretical uncertainty in the threshold region once the first few powers of LP logarithms have been resummed (see e.g. [24], and [25] for a more recent discussion). Secondly, the origin and general structure of NLP terms -including whether or not they share similar universality properties with their LP counterparts -is an interesting problem of quantum field theory in its own right. Thirdly, the classification of NLP contributions in cross-sections is closely related [26] to the study of so-called next-to-soft theorems, which have been explored in both a gauge theory [27] and gravitational context [28][29][30] due to their intriguing relation with asymptotic symmetries.
Whether or not a general resummation prescription exists for NLP terms is still an open question, that has been explored using an assortment of methods , some of them building upon the earlier work of refs. [56][57][58]. In order to further develop and test such formalisms, it is crucial to have detailed theoretical data -namely, explicit results for threshold logarithms up to NLP power in specific processes. Furthermore, it is extremely useful to classify separately contributions to each individual NLP term that come from real or virtual radiation that is soft and / or collinear (or hard, in the case of multiple emissions). Drell-Yan production offers a particularly clean testing ground in this regard, given that all threshold logarithms associated with purely real radiation are manifestly (nextto-) soft in origin (see e.g. [22]). Virtual gluons, however, can indeed be collinear with one of the incoming parton legs, as well as hard or soft, thus leading to a nontrivial structure of threshold logarithms. A convenient way to classify each individual contribution is to carry out the integration over virtual momenta using the method of regions [59][60][61], which explicitly separates out the modes of the loop momentum into non-overlapping soft, collinear or hard configurations. This method was heavily used in the calculation of the total cross-section for Higgs boson production at N 3 LO [14,15], and was also used in ref. [34] to reanalyse the 1-real, 1-virtual contribution to the NNLO Drell-Yan cross-section, first calculated in refs. [3][4][5][6][7][8], to obtain the contribution associated with each separate virtual region. This data proved essential when deriving a factorisation formula for next-to-soft effects [35,36], which generalises the well-known soft-collinear factorisation formula at LP (see e.g. ref. [62]), and which may pave the way for a NLP resummation formalism (see refs. [45-51, 54, 55] for an alternative approach based on effective field theory).
Reference [34] focused specifically on abelian-like contributions to the qq initial state, which in QCD are associated with the colour structure C n F at O(α n s ). At any given order, such terms are amongst the most complicated in terms of the number of different NLP effects that underly their structure. Furthermore, the development of factorisation formulae and / or resummation prescriptions for threshold corrections can be made systematically simpler by beginning with the abelian-like theory (as in refs. [34][35][36]58]), before generalising to the non-abelian case. We will thus restrict ourselves to abelian-like contributions in this paper, but our aim is to extend the classification of threshold contributions, up to NLP in the threshold variable, to diagrams involving one virtual gluon and two real emissions. As in ref. [34], the presence of the virtual gluon means that there are non-trivial regions to analyse. Furthermore, the results will have a direct bearing on how to generalise the factorisation formula of refs. [35,36] to include the effects of more than one gluon emission, which is clearly a necessary component for resummation. Although this is our main motivation, it should be stressed that the results of this paper constitute part of the Drell-Yan crosssection at N 3 LO, which is not yet known, although leading power threshold terms have been previously evaluated in refs. [63][64][65].
The structure of our paper is as follows. In section 2, we review necessary facts regarding Drell-Yan production, and outline the various steps used in our calculation. In section 3, we present results for the abelian-like contribution to the Drell-Yan K factor, before discussing their structure. We conclude in section 4. Some technical details are contained in the appendices.
2 Outline of the calculation

Drell-Yan production
In this section, we review some necessary facts about the Drell-Yan process, and the method of regions, that will be needed for what follows. Throughout, we focus on the quark-antiquark Drell-Yan production of a colour singlet vector boson, corresponding to the LO process For our purposes, we may take V to be an off-shell photon, and let e q denote the electromagnetic charge of the incoming quark. We further define the variable where Q 2 is the virtuality of the vector boson, and s = (p +p) 2 the squared centre of mass energy. At leading order, z = 1, such that the cross-section may be written where and N c is the number of colours. At higher orders, one has 0 ≤ z ≤ 1, such that the upper limit corresponds to threshold production. We may then define the K factor where the right-hand side contains the differential cross-section at O(α n s ). The complete K factor for Drell-Yan production, including all partonic channels and full z dependence, has been previously calculated up to NNLO (n = 2) [2][3][4][5][6][7][8], and leading power threshold contributions at N 3 LO have been evaluated in refs. [63][64][65]. At any given order, one must include the effects of additional radiation, that may be real or virtual. Reference [34] reanalysed the 1-real, 1-virtual contribution to K (2) (for the qq channel), up to the first subleading power in a threshold expansion about z = 1. In this limit, the K factor assumes a form similar to eq. (1), containing plus distributions and logarithms of the threshold variable ξ = 1 − z. As discussed in the introduction, ref. [34] focused on all contributions up to next-to-leading power (NLP) in ξ, that are proportional to the colour factor C 2 F , where C F is the quadratic Casimir in the fundamental representation. Such contributions are similar to those one would obtain in an Abelian theory, upon replacing g s C F with the relevant electromagnetic charge of the quark, and the aim of ref. [34] was to classify the precise origin of all such contributions, according to whether the virtual gluon is hard, soft or collinear with one of the incoming (anti-)quarks. Here, we carry out a similar analysis for the case of one virtual gluon, and two real emissions. This contributes to the N 3 LO factor K (3) (z), and the virtual gluon has a number of non-trivial momentum regions that give rise to NLP terms.
The amplitude we consider is shown schematically in figure 1, and corresponds to the process at one-loop order. Labelling this by A 2r,1v , its contribution to the differential cross-section occurs through interference with the pure two real emission amplitude A 2r : where the prefactors originate from colour / spin averaging and the Lorentz-invariant flux factor, we work in d = 4 − 2 spacetime dimensions throughout, and dΦ (3) is the differential phase space for the 3-body final state. There are 48 distinct Feynman diagrams that contribute to the abelian-like one-loop amplitude (where we define abelian-like diagrams to be those that contribute to the C 3 F colour structure in the cross section, thereby also excluding diagrams with a fermion loop). We have generated all such diagrams using QGRAF [66], and subsequently used Reduze [67,68] (version 2) to construct the interference term appearing in eq. (8). At this stage, one must carry out the integration over the loop momentum k appearing in eq. (8) and figure 1. To this end, we also use Reduze to reduce the one-loop integration to a set of scalar master integrals, using integration by parts identities. These integrals may themselves be represented as scalar Feynman diagrams with topologies of increasing complexity. The box and pentagon master diagrams are shown in figure 2, where the simpler bubbles and triangles are omitted for brevity.
As stated above, the aim of our paper is to classify the structure of the K factor up to NLP terms in the threshold expansion. We must then consider each master integral, and elucidate its corresponding contribution to threshold behaviour, according to whether the loop momentum is hard, soft or (anti-)collinear to one of the incoming partons. Here we follow the standard approach of the method of regions [59][60][61], which we describe more fully in the following section.

The method of regions
In the method of regions, singular parts of integrals in perturbative amplitudes are partitioned, according to physical criteria on the loop momenta. In the case of the threshold expansion considered in this paper, it is possible to separate completely the singular behaviour into non-overlapping regions, whose individual contributions reconstruct the full integral (itself expanded about the threshold limit) when summed. As an example, consider the diagram (B 1 ) of figure 2, where we have associated the loop momentum k with a particular internal line. One may expand this momentum in a Sudakov decomposition where we have defined dimensionless lightlike vectors in the directions of the incoming particles, as well as the vector k ⊥ transverse to the beam direction i.e. such that k ⊥ · n − = k ⊥ · n + = 0.
Denoting the Sudakov components of the loop momentum via k µ = (k + , k ⊥ , k − ), we may define the various regions by different scaling behaviours of these components. That is, one may introduce a book-keeping parameter λ ∼ √ 1 − z, such that the regions we need to consider are given by momenta of the form where the terms collinear and anti-collinear denote collinearity with respect to p andp respectively. In any given (scalar) master integral, the denominators can be systematically expanded in λ in each region, keeping the first subleading power where necessary to achieve NLP order in the final expression for the K factor. The integral in each region can then be carried out, and the results from all regions added together to reproduce, in principle, the threshold expansion of the full integral. Note that these are not the only possible scalings: in principle, it is also possible to consider momenta scaling as and so on. It is possible, however, to show that the only regions relevant for the threshold expansion are the hard, (anti-)collinear and soft regions defined by the scalings of eqs. (12). All other regions give scaleless integrals, which vanish in dimensional regularisation, such that we may discard them in the following. By definition, the incoming momenta are (anti-)collinear: while the gluon momenta are soft, i.e.
There is an interesting subtlety in the above procedure, if one wants to be sure of having characterised all possible regions of a given master integral. Before the region expansion, a given master integral possesses a symmetry under shifts of the loop momentum, such that one may associate the loop momentum k with an arbitrary internal line of the master diagram. However, the decomposition of k into regions breaks Lorentz invariance, leading to a violation of the shift symmetry. It may then be the case that particular choices of k are such that one cannot unambiguously identify all possible regions. To illustrate this point, let us consider diagram (B1) of figure 2, which we redraw in figure 3 so as to label the internal lines in what follows. In this particular case, a naïve choice of loop momentum will indeed lead to an important region being missed. Furthermore, this is a problem that arises for the first time at N 3 LO, due to requiring the presence of a virtual gluon, and two real emissions. Given that this problem does not seem to have been spelled out in related calculations in the literature (e.g. [15]), we believe it is instructive to examine this example in detail here. We consider the expansion in regions of the box integral represented in figure 3. The integral is defined as where D i represents the propagator associated with line i in figure 3, and we have introduced the convenient notation where d = 4 − 2 , and µ MS = µ e −γ E /2 (4π) 1/2 . Choosing the loop momentum k to correspond to line a seems natural, because in this way the regions are directly associated with having a hard, collinear or soft "gluon" exchange in the loop, which should be easily interpreted in the context of an effective field theory containing soft and collinear gluons. We can then define the denominators and expand the loop momentum k in regions using the Sudakov decomposition of eq. (9). One obtains (writing a · b ≡ ab in places so as to compactify expressions), Table 1: Scaling associated with the terms in the propagators D b and D d , as defined in eq. (18), where we set s ∼ 1. Leading terms in each region are highlighted in grey. Table 2: Scaling associated with the terms in the propagators D c , as defined in eq. (18). Leading terms in each region are highlighted in grey.
The scaling in λ of the various terms in the different regions is provided in tables 1 and 2.
In the following we keep only the leading terms for each propagator, thus getting the leading power contribution to the box integral. The hard region turns out to give Following the same criterion, a naive expansion in the collinear region, assuming the scaling assigned in table 1 and 2 gives, to leading power, Note that the hard region gives a subleading power contribution compared to the collinear region. Within a consistent expansion to leading power the hard region is thus zero, even if it is not scaleless. Furthermore, is it possible to show that integration in the anti-collinear and soft regions give scaleless results: Table 3: Scaling associated with the terms in the propagators D a and D d , as defined in eq. (23).
Thus, the leading power contribution to the integral in eq. (15) seems to be given by the collinear region in eq. (20). This conclusion is erroneous, however, as an important contribution has been missed, where the latter can be revealed easily by shifting the loop momentum to k = k + p. As discussed above, shift symmetry is broken by the region expansion, such that shifting the loop momentum can lead to inequivalent regions in general. With the new choice of loop momentum, the propagators read so that applying the Sudakov decomposition of eq. (9) gives The scaling of the various component in the different regions is provided in tables 3 and 4. Notice that we label the new regions with a prime, to distinguish them from the regions considered with the previous parameterization. It is easy to check that the new hard, collinear and anti-collinear regions still give the same result as the old corresponding regions: Table 4: Scaling associated with the terms in the propagators D c , as defined in eq. (23).
The new soft region, however, is not scaleless, and gives a new contribution which was not present in the old parameterization: In order to reconcile these results, note that the problem with the original choice of loop momentum is that the external scales are not well separated: both the "collinear" scale √ s n + (k 1 + k 2 ) ∼ λ 2 and the "soft" scale 2k 1 · k 2 ∼ λ 4 appear in the same propagator D c .
This causes problems in the collinear region because, even if the leading power terms in D c scales as λ 2 (see table 4), the loop integration is still over the full domain. There is therefore a region of the integration domain in which k ∼ −p, so that one has i.e. the leading power terms ∼ λ 2 cancel, causing the subleading power term 2k 1 · k 2 ∼ λ 4 to become leading. Considering this term small in the expansion of the propagator thus leads to the wrong analytic structure of the integral in this limit. The consequence is that the propagator D c cannot be expanded in the collinear region when parametrizing the loop momentum as in eq. (17). Rather, one needs to consider a more general collinear region "C", in which the propagator D c is kept unexpanded: Evaluating I C exactly and expanding at threshold after integration, indeed one finds that it contains both the contribution from the collinear and the soft region associated with the alternative loop momentum choice of eq. (22): An independent check can be performed with the program Asy [60,69], which provides a geometrical method to reveal the regions contributing to a given integral. Using the program with the integral in eq. (15) reveals the existence of three non-scaleless regions, which correspond to the hard, collinear and soft regions found within the second parameterization of the loop momentum in eq. (22). The same program can be used to verify that we have captured all regions in every other diagram. Some readers may be wondering why the hard region exhibits infrared singularities in the above results, which can be a common point of confusion in the method of regions. The approach we have taken above is to perform all required momentum scalings, and to set to zero any integrals which remain scaleless in dimensionless regularisation. In the soft region, expansion of the propagators changes the ultraviolet scaling behaviour of the integral, and thus introduces (spurious) ultraviolet divergences, whose effect is to cancel infrared divergences associated with exchange of multiple gluons between the incoming (anti)-quark legs, i.e. associated with the scale s. One can instead choose to isolate these UV divergences and absorb them into the hard function, and the effect of this procedure is to transfer poles in from the hard to the soft region. Given that this has no bearing on the final result for the K factor (which is a sum of all regions), we do not do this here. However, it should be remembered throughout that poles appearing in the hard region are indeed of soft origin.
Despite the above cancellation between UV and IR divergences, there remains the abovementioned nonzero contribution to the soft region, which is particularly interesting in that it is new at N 3 LO in perturbation theory. To see this, note that we need a virtual gluon in order to discuss decomposition of the loop momentum. Furthermore, the new soft region involves the momentum scale k 1 · k 2 , which can only be formed if there are at least two soft gluons in the final state. Detailed scrutiny of the region expansion applied to each of our Feynman diagrams reveals that the sole contribution to the soft region stems from physical configurations similar to those of figure 4. In the example shown, the incoming collinear quark turns into a soft quark by emitting a collinear gluon, where the soft quark then emits two soft gluons. As is well-known, soft quarks are subleading (in the momentum expansion) relative to soft gluons. Thus, we expect the soft region to contribute (if at all) at NLP level only. Furthermore, the somewhat complicated structure of soft and collinear emissions, together with the fact that this region occurs for the first time at N 3 LO, suggests that it will be suppressed by a number of powers of , so as to give subleading logs in the final result for the K factor. We will see in what follows that both of these expectations are borne out. It is also worth mentioning that a similar soft region was seen in the N 3 LO Higgs boson computation of ref. [15], where it was found to indeed be nonzero. We expect an essentially identical contribution to appear within the framework of soft collinear effective theory (SCET).
In summary, careful application of the method of regions to the process of figure 1 reveals the presence of hard, (anti-)collinear and soft regions. The latter crucially relies on the presence of a virtual gluon (giving rise to a loop momentum expandable in regions), as well as two real emissions, to provide the nonzero scale k 1 ·k 2 associated with the soft region. After expanding all propagators in each region, all integrals over the loop momentum k can be carried out analytically. Given that such integrals at one-loop order are quite standard in the literature, we do not report intermediate results here. Results for the squared matrix element in each region can be found in the following section. In order to cross-check our results, all steps of this calculation (e.g. diagram calculation, reduction to master integrals, expansion in regions, loop integration) have been carried out twice, in two completely independent implementations, and with full agreement.

Phase space integration
Applying the methods of the previous section, one obtains the interference term appearing in the integrand of eq. (8), expanded in regions and integrated over the loop momentum. The results are compact enough to report here, and it is first convenient to define the invariants as well as the combination consisting of the 1-loop double real contribution contracted with the conjugate tree-level result, integrated over the loop momentum. Results for the hard region at (next-to) leading power are then as follows: and the various functions {f X i } are defined in appendix A. Likewise, the squared matrix element in the collinear region turns out to be The anti-collinear region can be straightforwardly obtained through the exchange p ↔p. Finally, there is the soft region, which yields To compute the contribution of eqs. (31)-(34) to the differential cross-section or K factor, we must integrate over the Lorentz-invariant three-body phase space associated with the final state, as stated in eq. (8). One is free to choose a particular momentum frame for the phase space integration. Furthermore, given that each separate term in eqs. (31)- (34) is Lorentz invariant, we are free to choose different frames for different types of contribution, according to convenience. For the hard and collinear regions, expanding the right-hand side of eqs. (31)-(33) before substituting into eq. (8) reveals a series of terms, all containing the master integral where δ ∈ {0, 1}. For these values of δ, it is possible to obtain a result for this integral as an expansion the threshold variable (1 − z) for any value of the spacetime dimension d, by decomposing each real gluon momentum k i in a Sudakov decomposition, similar to eq. (9). We spell out this derivation in appendix B, and here present the results which are sufficient to integrate eqs. (31) and (33) to NLP order in (1 − z). Here we have defined as well as the total solid angle in (d − 2) spatial dimensions For the soft region, we rely on the symmetry of eq. (34) under the (combined) exchange of p ↔p and k 1 ↔ k 2 to reduce the number of distinct terms that need to be integrated. There remain two types of terms: (i) those involving the hypergeometric function 2 F 1 (1, 1; 1 − ; t 2 /(t 2 + t 3 )); (ii) those without the hypergeometric. Terms of the latter form are similar to those that occur in the double real emission contribution to the NNLO Drell-Yan crosssection [3,4] (see also ref. [33] for a recent derivation in the present notation). To integrate them, one may apply straightforward algebraic identities such as (and similarly for {u i }) to create a series of terms of the form of eq. (35), with at most one α i and at most one β i non-zero. Furthermore, δ will have a fractional power that depends on , due to the presence of the factor s − 12 in eq. (34). As described in refs. [3,4,33], this integral can be carried out exactly in the centre of mass frame of the two final state gluons. We review this derivation in appendix C.
The most difficult phase space integrals occur in terms of type (i) above, namely those in the soft region involving a hypergeometric function. All such terms involve the master integral and we note that similar integrals have been carried out for Higgs boson production in refs. [12,70], whose methods prove very useful for the present study. We proceed as follows.
We first apply identities similar to eq. (39) to put all terms in the form where at most one α i and at most one β i is nonzero, finding in all cases that α 2 = 0. As we explain in appendix C, for integrals involving only (α 1 , β 1 ) potentially nonzero, one may use the centre of mass frame of the outgoing gluons to derive the analytic result (valid for arbitrary d) where the ellipsis denotes subleading powers of (1 − z). This expression can be easily expanded in using the HypExp package for the hypergeometric function [71,72]. All necessary values of the parameters {α i , β i , γ i , δ} are collected in appendix C, together with results for each integral, where for convenience we define For integrals involving (α 1 , β 2 ) non-zero, we were not able to find any comparable closed form expression. However, they can be evaluated using Mellin-Barnes techniques, and the "energies and angles" phase space parametrisation described in refs. [12,70]. We describe this method in appendix C, but note here that in order to apply it to integrals involving negative powers of γ 1 and / or γ 2 , one must reexpress them in terms of other integrals, some involving more than two nonzero values of (α 1 , α 2 , β 1 , β 2 ). Results are collected in appendix C, again using the notation of eq. (42). All aspects of the phase space integration, including the calculation of all relevant master integrals, have been carried out twice and completely independently, with full agreement.

Results
We now have all the necessary ingredients for assembling the abelian-like terms (∼ C 3 F ) in the 2-real, 1-virtual contribution to the K factor of eq. (6), in the qq channel up to NLP order 1 . We will present separate results for the hard, (anti-)collinear and soft regions. For the hard region, one has (in the normalisation of eq. (6)) where (given the focus of our study) we report only enhanced (non-constant) terms in the finite part, and we have made the conventional choice for the dimensional regularisation scale in the MS scheme. NLP terms will be sensitive to this choice, given that the K factor contains the dimensional combination Note that we have identified everywhere, i.e. we have neglected the delta function contribution that arises from rewriting LP terms in terms of plus distributions. The delta function terms mix with virtual corrections not included here, and thus are not worth reporting. For the collinear region, we find where the anti-collinear region gives an identical contribution. Finally, we have the soft region, whose contribution is The total result for the (unrenormalised) K factor up to NLP order in the threshold expansion can be obtained from the above results through the combination Equations (43,45,46) constitute the main results of this paper. As discussed above, our main motivation for presenting them is as a prerequisite for formulating and testing general prescriptions for classifying (and potentially resumming) NLP threshold corrections in arbitrary processes. Although a full study in this regard is beyond the scope of this paper, it is worthwhile to make a few remarks regarding the implications of our results. Following a detailed analysis of the 1-real, 1-virtual K factor in the qq channel [34], refs. [35] considered a general amplitude with N hard particles, which is then dressed by an extra gluon emission. A process-independent factorisation formula was presented, building on the earlier work of ref. [58], which captured all abelian-like contributions to the amplitude up to NLP order in the threshold expansion, in the absence of final state jets. This formula was generalised to include all fully non-abelian contributions in ref. [36], and extends the well-known soft-collinear factorisation formula for LP threshold effects (see e.g. [62]). It includes a number of universal functions describing soft and collinear behaviour, whose operator definitions involve (generalised) Wilson lines [29,32]. A new function occuring at NLP level is the so-called jet emission function, first introduced in ref. [58]. As its name suggests, it describes the dressing of a jet function (collecting virtual collinear effects) with an additional radiative gluon. A fully non-abelian operator definition for this quantity has been proposed for (anti-)quark jets in ref. [36], and calculated at one-loop order. A similar calculation is in progress for gluons, which would have immediate applications in e.g. Higgs boson production via gluon-gluon fusion.
In processes containing two or more additional gluons, an open question is whether the functions appearing in the one-emission case are sufficient to capture all physics up to NLP order in the threshold expansion, or whether new functions should appear. For example, one may consider generalising the jet emission function to a family of quantities representing the dressing of a nonradiative jet with arbitrary numbers of additional gluons. For resummation of NLP effects to be possible, it should ideally be the case that these higher multiplicity jet emission functions are related by an iterative property to those with lower numbers of emissions (for a preliminary discussion in a purely abelian context, see ref. [58]). Or, this may be possible only up to a given subleading logarithmic order.
At NNLO in Drell-Yan production, it was already noticed that, perhaps unsurprisingly, the (anti-)collinear region in the method of regions maps straightforwardly to the contribution of the jet emission functions associated with the incoming (anti-)quark legs in the factorisation approach. Furthermore, this contribution started only at next-to-leading-logarithmic (NLL) order, and at next-to-leading power (NLP) in the threshold variable. In the present calculation, we also see that the (anti-)collinear regions start only at O( −4 ) rather than O( −5 ). Thus, again we find that collinear effects are NLP, and give only subleading (NLL) threshold logarithms. Indeed, the only source of leading LP or NLP effects is the hard region, as can be clearly seen in eq. (43). This observation will certainly be a useful guide when examining the extent to which (multiple) jet emission functions are relevant at higher orders in perturbation theory. Furthermore, there is much existing evidence (most notably in ref. [38]) that the highest power of the NLP log exponentiates in Drell-Yan. The observation that collinear effects do not affect this log at N 3 LO provides a significant hint regarding how to formally prove this property.
The new soft region at this order depends crucially on the presence of two gluons, and so would seem to be a correction to factorisation formulae of the type presented in refs. [35,36], in that it cannot be composed iteratively from lower-order information. However, it is worthwhile to note that the soft region itself is heavily suppressed in the expansion, so that it only contributes logarithmic terms at N 4 LL level. If this behaviour persists at higher orders, such a region is unlikely to trouble realistic efforts to resum NLP effects, but it should of course be fully understood, as it will be present in the exact Drell-Yan K factor at higher orders.
Further insights into the iterative structure of our results can be obtained by examining the squared matrix elements before integration over the final state phase space, but after the integration over the loop momentum of the virtual gluon. In the case of the hard region (eq. (31)), we find that the coefficient f H 1 matches the similar function found in the oneloop quark form factor, such that the leading power term agrees with what one obtains from applying the well-known eikonal Feynman rules to the non-radiative one-loop Drell-Yan process. At NLP, we noted that the second coefficient f H 2 already appears in the 1-real, 1-virtual contributions at NNLO. Thus, there is strong evidence that the hard region can indeed be understood using the existing tools of refs. [32,33,35,36]. In the collinear region, we find that the function f C 1 in the first line of eq. (33) occurs already at NNLO, such that this contribution factorises into a one-loop jet emission on the quark leg, dressed by a tree-level emission from the anti-quark (and vice versa for the anti-collinear region). The remaining collinear contributions, involving the additional coefficients f C 2,3 , lack such a straightforward interpretation, leaving open the possibility that one must consider a separate jet emission function for pairs of gluons. Finally, as discussed already above, the soft region is not expected to be iteratively obtainable from lower order information.

Conclusion
In this paper, we examine abelian-like contributions to Drell-Yan production in the qq channel at N 3 LO, namely those with the colour structure ∝ C 3 F . We have classified all logarithmically enhanced contributions near threshold when one gluon is virtual, and the other two real, up to next-to-leading power (NLP) in the threshold variable (1 − z). Our motivation is to work towards a systematic classification of NLP threshold effects, building on e.g. the factorisation formulae of refs. [35,36] (see refs. [45][46][47][48][49][50][51] for similar work within the context of effective field theory). To this end, we present results for the unrenormalised K factor, using the method of regions [59][60][61] to separate contributions according to whether the virtual gluon is hard, soft or collinear with one of the incoming particles. Our hope is that this provides a great deal of useful information for elucidating the general structure of NLP effects, similar to how previous methods of region analyses at NNLO [34] directly informed the construction of factorisation formulae valid to subleading order in the threshold expansion.
There are a number of noteworthy features in our result. Firstly, there is a nonzero soft region that appears for the first time at N 3 LO, and which we find persists upon integration over the final state phase space. The presence of such a contribution requires at least one virtual gluon and two real gluons, and thus does not appear to be iteratively relatable to lower order information. A similar region was found to be nonzero in the recent (and closely related) calculation of Higgs boson production via gluon-gluon fusion [70], whose methods prove very useful for the present analysis. The overall contribution of this region to the Drell-Yan K factor is highly subleading, in that it contributes with a single pole in the dimensional regularisation parameter at O(α 3 s ), corresponding to a N 4 LL NLP logarithm in the finite part of the K factor. It would be interesting to see what effect such a region has at higher orders in perturbation theory, and indeed whether it has a straightforward counterpart in SCET.
Unlike the hard region, the collinear region does not contribute to the leading NLP logarithm, suggesting that collinear effects are not relevant to the potential resummation of the highest power of NLP logs to all orders in perturbation theory. Both the hard and collinear regions in our analysis show signs of an iterative structure, whereby parts of the results can be obtained from lower order information. These observations will prove highly useful in generalising factorisation formulae for NLP effects to higher orders in perturbation theory.
There are a number of directions for further work. Immediately related to the present study would be the calculation of threshold contributions in the triple real emission contributions to Drell-Yan production at N 3 LO, or in the double-virtual, single real channel. Furthermore, one can generalise the calculation to include all possible colour structures, involving fully non-abelian corrections. Finally, the implications of our results for developing a fully systematic classification of NLP threshold effects in arbitrary scattering processes will be the subject of much further study.

A Coefficients entering the matrix element
In this appendix, we collect results for the various coefficients appearing in eqs. (31)- (34). Starting with the hard region, we have The coefficients for the (anti-)collinear regions are For the soft region, we have

B Phase space integrals in the hard and (anti-)collinear regions
In this appendix, we spell out the derivation of eq.(36), using the Sudakov decomposition of eqs. (9)(10)(11). Furthermore, we define the quantities k i+ = n − · k i and k i− = n + · k i , using a slightly different convention to the Sudakov decomposition of the loop momentum in section 2.2, so as to make factors of 2 more convenient in the following. The 3-body phase space in d dimensions is given by where and θ is the Heaviside function We may carry out the integral over the photon momentum q using the delta function in eq. (51), obtaining where in the second line we have used eq.
where we can Taylor expand the exponential in k 1 · k 2 , given that higher order terms will be suppressed by powers of 1 − z: Putting things together, the phase space becomes where we have transformedω = iω. We can now use this result to carry out the integral of eq. (35) for the two special cases of δ ∈ {0, 1}. For δ = 0, we may note that the integrand of eq.(35) has no transverse momentum dependence, such that the linear term k 1⊥ · k 2⊥ in eq. (57) leads to an odd integrand, and can be neglected. Using polar coordinates for the k i⊥ integrals, one may use the onshell delta functions to eliminate the integral over |k i⊥ |, such that eq. (35) becomes After a variable changek i± =ω √ s k i± , we may recognize the inverse Laplace transform The integrals overk i± will be of the form: for which the variable transformatioñ Substituting these results, we obtain eq. (36) as required.
The integral of eq. (35) with δ = 1 appears only at NLP level, such that we may entirely neglect the term k 1 · k 2 in eq. (56), as it will lead to terms suppressed by further powers of (1 − z). Carrying out similar steps to the δ = 0 case, one again finds eq. (36).

C Phase space integrals in the soft region
In this appendix, we describe various integrals (of increasing complexity) that occur when integrating the squared matrix element in the soft region (eq. (41)) over the final state phase space.

C.1 Integrands with no hypergeometric function
First, we need integrals of the form of eq. (35), in which at most one parameter {α i } and at most one parameter {β i }. The Sudakov decomposition of appendix (B) turns out not to be helpful here, due to the fractional power of δ. Instead, one may simplify the calculation by working in the centre of mass frame of the two outgoing gluons [3,4,33]. In this frame, one writes wheret and λ is the Källen function λ(a, b, c) = a 2 + b 2 + c 2 − 2ab − 2ac − 2bc. The Mandelstam invariantst andũ can in turn be expressed as functions of the photon energy fraction z = Q 2 /s and of two further variables 0 < x < 1 and 0 < y < 1, such that where (1 − z) is the threshold variable. The 3-body phase space in d dimensions now takes the form In terms of the above definitions, one finds where These quantities satisfy the relation such that upon defining and using the above definitions, the angular integral may be carried out using the result [3] (first derived in [73]) At this stage, one must carry out the integrals over the variables x and y appearing in eq. (64). These can all be carried out in terms of beta functions, or using the identity

C.2 Integrands with a hypergeometric function
Next, we must consider phase space integrals such as those of eq. (40), where the integrand contains a hypergeometric function. As is the case for the similar integrals in refs. [12,70], we have not found it possible to obtain a useful closed form analytic result for arbitrary values of the parameters. However, for a certain subclass of the parameters, we can indeed find such a result, valid for any d. Let us present this case first.
If α 2 and β 2 are both zero, eq. (40) reduces to where In the centre of mass frame of the two outgoing gluons (see section C.1), this becomes Next, one can use the Mellin-Barnes representation for the hypergeometric function so that eq. (73) becomes The angular integrals can be carried out using eq. (69), to get At this point, we may expand the integrand in (1 − z), taking the leading power only. After some work, we end up with The y integral can be carried out immediately in terms of Gamma functions. The x integral would give a 3 F 2 , but then the remaining Mellin-Barnes integral could be cumbersome. Instead, we can introduce a second Mellin-Barnes representation, after which the x integral can be carried out in terms of Gamma functions, yielding We must now carry out the double Mellin-Barnes integral. However, this can be done straightforwardly, by recognising the w 2 integral as i∞ −i∞ where we have used Gauss' identity At this stage we are left with Using eq. (74) we can recognise the w 1 integral as Putting everything together, we obtain the result of eq. (41).

C.3 General parameter values
As stated above, for other necessary values of the parameters, we are not able to find a closed form solution for the integral of eq. (40), valid for any spacetime dimension d. Instead, we may settle for an expansion in the dimensional regularisation parameter . To this end, it is useful to use an alternative phase space parametrisation, as discussed in refs. [12,70]. We first write eq. (40) as where a = − and which differs from eq. (72) in having arbitrary powers of all two-particle invariants. Reference [12] starts by scaling momenta according to 2 so that eq. (84) becomes × (p 1 · p 3 + p 1 · p 4 ) −γ 1 (p 2 · p 3 + p 2 · p 4 ) −γ 2 2 F 1 1, 1; a + 1; where The integral in the second line is now dimensionless. Furthermore, if one wants the leading behaviour in (1 − z), then this has already been extracted, so that one can set z = 1 in the integral itself. In practice this is done by using a particular parametrisation for the rescaled momenta, and a particular expression for the soft phase space. The momenta are parametrised in the lab frame, which immediately implies Furthermore, we can choose to write p 3 and p 4 in terms of a d-velocity β i : (88) 2 Our notation {p i } coincides with the notation used in ref. [12] after the rescaling has taken place.
Note that, despite appearances, E i is dimensionless due to the rescaling introduced above. Upon substituting eqs. (87) and (88) where following ref. [12] we have defined where the current notation s 12 should not be confused with the scale s 12 = 2k 1 · k 2 used in the main text. At this point one may introduce the Mellin-Barnes representation (see e.g. ref. [74]) as well as the identity The phase space integral now has the form of multiple products of two-particle invariants, thus is of the same form as the integrals considered in refs. [12,70]. The invariants can be rewritten using the parametrisation of eqs. (87) and (88): Furthermore, the leading behaviour of the phase space measure as z → 1 is given by 3 (see e.g. ref. [12]) where dΩ (d−1) i is the differential solid angle associated with particle i. Using eqs. (94) and (95) in eq. (93), one may carry out the E i integrals using yielding 2 2C−δ+5−4d π 3−2d Γ(a + 1) Γ 2 (a)Γ(γ 2 )Γ(2d − C + 2δ − 4) Next, we must carry out the angular integrals. Given that each d-velocity β 3 and β 4 occurs thrice rather than twice, we can no longer use eq. (69). Unfortunately, there is no known closed form for the angular integral involving three angular quantities. There is, however, a triple Mellin-Barnes form [75] (see also eq. (5.17) of ref. [12]) in d = 4 − 2 dimensions: × Γ(λ 1 + z 4 + z 5 )Γ(λ 2 + z 4 + z 6 )Γ(λ 3 + z 5 + z 6 )Γ(1 − λ 1 − λ 2 − λ 3 − − z 4 − z 5 − z 6 ) Upon using this result, the remaining integral over the angular variables of particle 4 can be carried out using eq. (69), which it is more convenient to write as Our general phase space integral now has the form of a six-fold Mellin-Barnes integral, which applies if γ 1 and γ 2 are both non-zero. If either of them is zero, we do not need to apply eq. (92) for the relevant combination of invariants, and thus we will obtain a lower order Mellin-Barnes integral from the outset. Our strategy for carrying out an integral for general (α 1 , α 2 , β 1 , β 2 , γ 1 , γ 2 , δ, a) is as follows: Integrals on the right-hand side that only involve powers of p · k 1 and / orp · k 1 can be carried out using the analytic result of eq. (41). Remaining integrals can be carried out using the Mellin-Barnes approach outlined in this section. Note, however, that for the second term in the last line of eq. (101), it is straightforward to derive a closed form, valid for any d.