Threshold factorization of the Drell-Yan quark-gluon channel and two-loop soft function at next-to-leading power

We present a factorization theorem of the partonic Drell-Yan off-diagonal processes gq¯qg→γ∗+X\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ g\overline{q}(qg)\to {\gamma}^{\ast }+X $$\end{document} in the kinematic threshold regime z=Q2/ŝ→1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ z={Q}^2/\hat{s}\to 1 $$\end{document} at general subleading powers in the (1 − z) expansion. Focusing on the first order of the expansion (next-to-leading power accuracy with respect to the leading power qq¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ q\overline{q} $$\end{document} channel), we validate the bare factorization formula up to Oαs2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}_s^2\right) $$\end{document}. This is achieved by carrying out an explicit calculation of the generalized soft function in d-dimensions using the reduction to master integrals and the differential equations method. The collinear function is a universal object which we compute from an operator matching equation at one-loop level. Next, we integrate the soft and collinear functions over the convolution variables and remove the remaining initial state collinear singularities through PDF renormalization. The resulting expression agrees with the known cross section in the literature.


Introduction
It is well documented that the perturbative expansion of QCD fails near the kinematic threshold, as the phase space for real emission is restricted to contain only low-energy (soft) radiation.Considering the Drell-Yan (DY) process A(p A ) + B(p B ) → γ * (Q)[→ ℓ l ] + X, with X being the unobserved QCD final state, the threshold regime of the partonic cross section is characterised by the limit z ≡ Q 2 /ŝ → 1, where Q 2 represents the invariant mass of the final state lepton pair and ŝ the partonic centre-of-mass energy squared.In this region, physical observables are expressed as a power expansion in (1 − z) → 0 and feature large logarithmic corrections in this variable.Reliable results can only be obtained by resumming these logarithms to all orders in perturbation theory.This was first achieved for the leading-power (LP) contribution more than thirty years ago using diagrammatic techniques in [1,2] and later using soft-collinear effective theory (SCET) methods in [3][4][5].
Nowadays the resummation of LP threshold logarithms is understood in DY up to next-to-next-to-next-to-leading logarithmic (N 3 LL) accuracy [6][7][8][9].However, comparatively less is known about terms which are suppressed by a power of (1 − z), the so-called next-to-leading power (NLP) logarithms, despite the fact that studies of amplitudes with next-to-soft emissions have quite a long history [10][11][12].Focusing on the diagonal (q q) channel of the Drell-Yan process, calculations of partonic cross sections at NLP up to next-to-next-to-leading order (NNLO) in the strong coupling expansion, and partly beyond, have been carried out using the expansion-by-regions method [13,14] and diagrammatic factorization techniques [15][16][17][18][19]. Investigations of the NLP terms for the DY process within the effective field theory (EFT) framework were initiated in [20] and [21].In the former, the NLP logarithms were resummed to leading logarithmic (LL) accuracy, and in the latter, the complete subleading power factorization theorem was derived.In the present work, we complete these investigations by deriving and validating up to NNLO the NLP factorization theorem for the off-diagonal g q (qg) channel of the DY production process at threshold, which constitutes the last missing piece up to NLP accuracy.
The development of this framework is timely as plenty of attention has recently been given to NLP studies of the analogous off-diagonal channels in deep-inelastic scattering (DIS) at large Bjorken-x and "gluon thrust" in hadronic e + e − annihilation, their corresponding diagonal channels, and Higgs production in gluon fusion at threshold [19][20][21][22][23][24][25][26][27][28][29][30][31][32].Progress beyond LP has also been achieved for variables such as N-jettiness [33][34][35][36][37][38] and the q T of the lepton pair or the Higgs boson [39][40][41].However, the extension of the standard LP factorization theorems to NLP has not proved straightforward.The main stumbling block being the ubiquitous appearance of endpoint divergences in the convolution integrals that connect the hard, (anti-) collinear and soft functions in NLP factorization theorems [21].This conceptual issue has been investigated in a number of contexts: for instance, conjecture regarding the form of the leading double logarithms related to soft quark emission has been presented in [26], and refactorization ideas were developed for DIS [29] and for Higgs decay to two photons (and gluons) through bottom-quark loops [42][43][44][45].Combination of standard SCET factorization and endpoint factorization was put forward in [32] to arrive at a factorization formula for the "gluon thrust" valid in d = 4, such that standard systematically improvable Renormalization Group (RG) methods can be applied to perform the resummation of large logarithms.Endpoint divergences have also recently been explored in QED [46] and B physics [47][48][49].
Concerning the study of the all-order resummation for the g q channel of DY, results at LL accuracy have recently been obtained using diagrammatic techniques in [31] and d-dimensional consistency relations in [50].In general, the development of a resummation framework at higher orders will require a genuine prescription for the treatment of the endpoint convolution divergences in the presence of hadronic initial states.In this respect, it is crucial to derive a consistent factorization theorem valid in d = 4.As a first step toward this goal, however, one needs to develop a consistent factorization theorem at the level of the bare cross section, and calculate all of the contributing functions to the relevant perturbative order in dimensional regularization.While this task has been completed for the DY q q channel in [21,51], the g q channel of DY has not been systematically address yet.The purpose of this paper is to fill this gap: in particular, starting from time-ordered power suppressed operators in SCET, we derive the factorization of the partonic cross section in terms of short-distance coefficients times the convolution of collinear and soft matrix elements, valid at any subleading power.We then focus on the NLP contribution, and evaluate the collinear and soft functions appearing in the factorization theorem respectively at O(α s ) and O(α 2 s ), which is the required accuracy to reproduce the known inclusive DY cross section result at NNLO.This provides a strong sanity check of our framework and allows us to reach the same accuracy as for the diagonal q q channel in [21,51].
The collinear function has a purely virtual origin and we evaluate it at one-loop order.Its appearance is the consequence of the insertion of the power suppressed interaction which couples a soft quark with collinear gluons and a collinear quark.This function is practically equivalent (with just some minor differences in its definition) to the radiative collinear function computed to O(α 2 s ) in [52].The calculation of the collinear function which appears in the factorization formula for the g q channel turns out to be much simpler in comparison to the evaluation of the collinear functions which appear in [21] since the derivative terms acting on the momentum conserving delta function are not present in this case.The soft function depends on the total energy of the soft emissions as well as two additional convolution variables.We evaluate the soft contributions at NNLO accuracy in exact d-dimensions by employing the reduction to master integrals and the differential equations methods that were also employed in [51].Other soft functions which appear in the h → γγ decay NLP factorization formula were computed in [53,54] at O(α s ).
Working with d-dimensional expressions allows us to carry out the convolution integrals at fixed-order accuracy.This is the first step needed for the application of the refactorization ideas developed in [32,44].We also study asymptotic limits of the obtained functions which is useful for any such program.However, additional insights concerning the operatorial structure of the incoming hadronic states are required to carry out this program, which we leave to a future investigation.
The paper is organised as follows: in section 2 we derive the general subleading power factorization formula for the partonic off-diagonal channels of the DY process and we specialize it to NLP.Then, in sections 3 and 4 we calculate the required NLP collinear and generalized soft functions up to O(α s ) and O(α 2 s ), respectively.As a last step in section 4, we consider asymptotic limits of the soft function.The obtained accuracy of the soft and collinear functions enables the validation of the factorization formula derived in section 2 to NNLO, which is carried out in section 5. We summarize and discuss possible future developments in section 6.

Factorization near threshold
We consider the Drell-Yan invariant mass distribution ŝ are respectively the hadronic and partonic threshold variables, with s = (p A +p B ) 2 and ŝ = (p a +p b ) 2 = x a x b s the hadronic and partonic centre of mass energy; Q 2 ≡ q 2 is the invariant mass of the final state photon; last, Λ is the confinement scale of QCD.In the present analysis, we do not consider power corrections in Λ/Q.
Near threshold the partonic cross section has the following power expansion where the power-counting parameter is defined as λ = √ 1 − z ≪ 1; the leading power term in the equation above is O(λ −2 ) and the next-to-leading power term is O(λ 0 ).Only the diagonal production channel ab = q q contributes at LP.At NLP more channels open up: in addition to the NLP correction to the q q channel, which has been studied at length in [21,51], there are also contributions from the gluon-antiquark g q and the quark-gluon qg channels, which are the topic of the present work.In this section we derive their factorization structure near threshold, valid at general subleading powers.Given that the two channels g q and qg give identical contributions, in what follows we focus for simplicity on the gluon-antiquark g q channel.
We derive the factorization theorem for the bare cross section in d = 4−2ϵ dimensions: in this case, dimensional regularization regulates the endpoint divergences arising in convolution integrals, hence calculations at fixed orders are well defined.In order to facilitate comparison with literature, we will consider an equivalent form of Eq. (2.1), given by where the parton luminosity function L ab (y) is defined as and ∆ ab (z) is related to the partonic cross section σab (z) in Eq. (2.1) by In order to obtain the factorization theorem for ∆ g q(z) let us start from the standard expression for the hadronic d-dimensional differential cross section.After integration over the phase-space for the final state leptons one has where W µρ is the hadronic tensor and in turn J ρ = q e q ψq γ ρ ψ q represents the electromagnetic quark current.In Eq. (2.7), = p X is the total momentum of the final state radiation X.In order to not obscure the derivation of the factorization theorem, in what follows we work with a single quark flavour and set the electromagnetic charge e q = 1.
The partonic kinematic threshold corresponds to the case where almost all of the energy in the interaction between the two incoming partons is carried away by the intermediate off-shell boson, γ * .The relevant modes consist thus of PDF-collinear and PDF-anticollinear modes, threshold-collinear and threshold-anticollinear modes, and soft modes.Momenta are decomposed along two light-like directions n µ ± , defined by The scaling of a given momentum l µ is thus indicated by the scaling of the components (n + l, l ⊥ , n − l).The PDF-collinear and anticollinear modes scale respectively as p µ represents the confinement scale of QCD, while threshold collinear, anticollinear, and soft modes scales respectively as and p s ∼ Q(λ 2 , λ 2 , λ 2 ).We work in position-space SCET [55,56] and describe each mode by its own set of fields.We have therefore PDF-and threshold-collinear quarks and gluons, and soft quarks and gluons.These fields are expressed in terms of gaugeinvariant building blocks, after the soft-collinear decoupling transformation has been applied.For instance, we express collinear fields modes in terms of the gauge-invariant building blocks where ] represent the original non-decoupled collinear fields.In turn, the covariant derivative is defined as In what follows we will always work with decoupled fields, and drop the index (0).In Eq. (2.8) the soft Wilson lines are defined as while the collinear Wilson lines reads In both of the above equations P is the path ordering operator.Notice also that in eq.(2.8) the soft Wilson lines are evaluated at z µ − ≡ (n + z)n µ − /2, in accordance with the fact that soft fields are multipole expanded when multiplying collinear fields [56].Concerning soft fields, we write them in terms of the gauge invariant building blocks where the soft covariant derivative is iD µ s = i∂ µ + g s A µ s .The power-suppressed softcollinear Lagrangian YM is then expressed in terms of these fields, as listed in appendix A of [21].We refer to [55,56] for a thorough derivation of the SCET Lagrangian and SCET fields definition.
Given that the final state is forced to only contain soft radiation, the hard matching to SCET fields can be performed at amplitude level, since there are no contributions to the hadronic tensor where the currents at positions 0 and x are connected by hard partons.For the process g q → γ * + X to take place in the z → 1 limit, the incoming c-PDF gluon must be converted to a threshold collinear quark through the emission of a soft antiquark.The threshold-collinear quark retains almost all of the momentum of the incoming c-PDF gluon.The soft-quark interaction with collinear fields is inherently a subleading power effect.This can be deduced from the fact that soft quarks at leading power only appear in the soft term, qs i / D s q s , in the SCET Lagrangian; the soft quarks first appear in soft-collinear interaction terms at O(λ) in L (2.12) Therefore, unlike in the case of the q q-channel discussed in [21], the contribution to the Drell-Yan cross-section from the g q-channel begins at NLP (i.e., at O(λ 2 )), through a time-ordered product insertion of L ξq into the amplitude and its complex conjugate, as represented in Fig. 1.
Figure 1: Schematic representation of the lowest-power contribution to the gluonantiquark channel in Drell-Yan.It involves the conversion of a (PDF) collinear gluon into a (threshold) collinear quark by emission of a soft anti-quark.The emission proceeds first at NLP by means of a time-ordered product, involving the leading power SCET current J A0 times the insertions of the power suppressed Lagrangian L (1) ξq , indicated by the index "1" in the picture.The matching of threshold modes onto PDF modes is trivial for the anti-collinear sector, as there is no power suppression on this leg [21].

Factorization at general subleading powers
The formal factorization formula for the g q-channel at general subleading can be obtained following closely the derivation for the q q-channel given in [21], therefore we provide here a rather concise discussion and draw attention to the differences between the two cases.Our discussion below follows the derivation given in [57].In general, one needs to take into account the contribution of time-ordered products involving the (multiple) insertions of YM Lagrangian terms and operators arising from the expansion of the Drell-Yan current, with the constraint that the time-ordered products must involve a collinear gluon and at least one soft quark emission.Omitting the index structure for clarity, the general, all power, hard matching of the vector current is given by where the symbols {dt k }, {t k } (and the barred ones for the anticollinear direction) indicate sets of convolution variables associated to the number of fields in each collinear direction.The indices m 1 and m 2 label the basis of SCET operators (and their corresponding short-distance matching coefficients C m1 m 2 ), following the formalism and notation developed in [58][59][60] 1 .The function J s involves only soft fields and starts at O(λ 3 ) [58].Given that the DY process involves a collinear and an anticollinear direction, the currents J m 1 ,m 2 ρ explicitly read where for instance the terms J m 2 c ( {t k }) are constructed using collinear-gauge-invariant collinear building blocks given in Eq. (2.8).In this general construction, the indices m 1 , m 2 consist of letters A, B, C etc., labelling the number of fields in a particular collinear direction, and numbers 0, 1, 2 etc., denoting the overall power of λ of the current with respect to the LP, which is labelled 0. The function Γ m 1 ,m 2 ρ in Eq. (2.14) gives the appropriate spinor and Lorentz structure of the operator.For instance, the structure of the LP operator reads Γ A0,A0 ρ = γ ⊥ρ ; at O(λ) one has Γ A0,A1 ρ = n +ρ , etc. Hard modes are integrated out in the matching between QCD and SCET, and appear as short-distance coefficients C m 1 ,m 2 in Eq. (2.13).In the next step, one integrates out threshold-collinear and anticollinear modes.As discussed in [20,21], this gives rise to a matching equation, which involves on the left-hand-side the SCET Lagrangian with threshold-collinear and soft fields, and on the right-hand side PDF-collinear and soft fields, times a short-distance coefficient, defined as "collinear function", which arises as a consequence of integrating out threshold-collinear modes.The collinear function is the analogue of the "radiative jet" discussed in [12,18,19], and we refer to section 2 of [21] for a detailed discussion.The specific form of the matching equation is dictated by gauge invariance, momentum and Fermion number conservation.For the g q channel, the general collinear matching equation reads where {d d z j } = m j=0 d d z j and {dz j− } = m j=0 dn + z j 2 .{z j− } denotes the set of m positions at which the soft building block insertions are located, and T represents the time-ordering operator.On the left-hand side, L (l) (z j ) is a set of m O(λ l )-suppressed Lagrangian insertions and {ψ c (t k n + )} denotes a set of n collinear fields, ψ c = χ c or A µ c⊥ , each dependent on one variable from the n-sized set {t k }, originating from the collinear part of the Drell-Yan current operators.On the right-hand side, the function Gi represents the NLP collinear function resulting from the matching of threshold collinear fields onto the initial state PDF-collinear gluon A PDF c and is analogous to the NLP collinear function Ji in [20,21], representing the matching onto an initial state PDF-collinear quark, which appears in the factorization of the q q channel.Last, as in case of the q q channel, the soft operator s i ({z j− }) represents a series of soft structures, containing the soft fields originating from the power-suppressed soft-collinear interaction.In case of the g q channel s i ({z j− }) must contain at least one soft quark, i.e.
where the index i labels the soft structures, and the ellipsis refers to additional structures appearing beyond NLP.Note that, for consistency with the q q channel, we normalise the collinear function such that at tree level it contributes to O(g 0 s ), and we absorb a factor of g s into the soft structures, as can be seen explicitly in Eq. (2.16).Notice also that the collinear gauge invariant building block itself contains a factor of g s , which is compensated by the explicit factor of g −1 s in Eq. (2.15).The matching for the anticollinear leg, with an incoming antiquark, is analogous to Eq. (2.15), and involves the anticollinear function Jī, as for the q q case, and the soft structures given in Eq. (3.45) of [21].After the second matching onto collinear and anticollinear PDF modes has been performed, the derivation of the factorization theorem follows very closely that of the q q channel.Taking into account the factorization of the state |X⟩ = |X c ⟩⊗|X c⟩ ⊗|X s ⟩, the g q matrix element reads where explicit Dirac, Lorentz, and color indices are suppressed.The index k ( k) counts the number of collinear (anticollinear) building block fields in a given current, and we sum over all of the possible currents.The index j ( j ) counts the number of insertions of the subleading power Lagrangians into each collinear (anticollinear) sector, and again we sum over all of the occurrences.
In order to take the square of the matrix element in Eq. (2.17) it proves useful to write it in terms of the momentum-space representation of the various functions, which are obtained by means of Fourier transformation.For the short-distance coefficients C m 1 , m 2 one has (2.18) then we need to introduce the momentum-space representation of the PDF-collinear and anticollinear fields.For instance, for the c-PDF gluon field we use with analogous definition for the c-PDF antiquark field.Furthermore, the Fourier trans-form of the collinear function is given as with an analogous definition for the anticollinear function.Here the set {ω j } denotes the variables with a soft scaling that are conjugate to { z j− } and in the exponents Einstein's summation convention is used.Also, . Making use of these definitions, Eq. (2.17) becomes

.21)
The matrix element can be written in a more compact form by introducing the following amplitude level coefficient functions such that the matrix element reads 2 Let us notice here that, in order to keep equations compact, we use the same symbol to indicate Therefore, when z − appears as argument of functions, as in Gm2 i ( {t k }, u; { z j− } ), or in scalar products z − • k, we refer to the vector z µ − , while when z − appears in exponents or integration measures, such as in e −iωj zj− or dz j− , we refer to the scalar , respectively.
× dg dḡ e i(n + pa)g e −i(n We can now insert this expression into the definition of the hadronic tensor W ρµ in Eq. (2.7), with an equivalent expression for the complex conjugate matrix element, and obtain the factorized expression at general subleading powers for the cross section given in Eq. (2.6).To this end one still needs to identify the PDF-collinear and anticollinear matrix elements respectively with the gluon and anti-quark parton distribution functions.For the former we use the following relation from appendix A of [65] 3 while for the anticollinear matrix element we use Let us remark that, although we do not make indices explicit in this general derivation, it is understood that the indices appearing in Eq. (2.24) are absorbed by the collinear functions.Integrating over the auxiliary variables g, ḡ, g ′ , ḡ′ and the momenta n + p a , n + p b , n + p ′ a , n + p ′ b , after some elaboration we obtain the g q-channel Drell-Yan crosssection as defined in Eq. (2.5): This is the result for the general form of the power-suppressed g q-induced partonic crosssection in the z → 1 limit.The notation with bars (¯) and tildes ( ) is used here in the same way as in the derivation of the q q-induced partonic cross-section of [21].They refer to the anticollinear direction and objects with dependence on the coordinate variables respectively.Also, the contributions from the complex conjugate amplitude are denoted with a prime ( ′ ) symbol.In the last line of equation (2.26) we have introduced the generalised multi-local soft function for the g q-channel.It is given by Let us also recall that, in the same fashion as the q q-channel factorization theorem in [21], the result in Eq. (2.26) is formally valid in d-dimensions.

Factorization at next-to-leading power
Let us now specialise the factorization theorem in Eq. (2.26) to next-to-leading power.
As discussed above, near threshold the production of an off-shell photon from an initial g q state involves at least the emission of a soft quark into the final state.Soft quarks can arise from the factor J s (0) in the expansion of the vector current, but this term is at least of O(λ 3 ), therefore any such contribution is beyond NLP.Hence, the only remaining possibility to achieve a power suppression is provided by time-ordered products with insertions of the Lagrangian terms ξq .Up to NLP only the two terms with i = 1, 2 can appear.The L (2) ξq Lagrangian insertion, however, can also be dropped, as the corresponding amplitude would have to be interfered with a leading power amplitude in order to yield an O(λ 2 ) power suppressed cross-section, and such contribution vanishes.This leaves solely the contribution due to L (1) ξq multiplied with the leading power J A0 current.Therefore, using the Fourier transforms in Eqs.(2.20) and (2.19), at NLP the collinear matching in Eq. (2.15) reduces to where we have written indices explicitly: α and γ are Dirac indices, η is a Lorentz index, a, f and A are a fundamental and adjoint color indices respectively.This is the analogue of Eq. (3.23) in [21] for the q q-channel at NLP.Here we do not sum over the soft structures, as at NLP there is only one, originating from L ξq , given by the term written in Eq. (2.16).With explicit indices it reads In order to simplify the next-to-leading power factorization formula as much as possible, we make use of generic properties of the collinear function G η,A ξq;γα,f a (n + p, n + p a ; ω) in the matching equation (2.28).Firstly, as in case of the q q channel, the collinear function must be proportional to the delta function in the collinear momenta, δ (n + p − n + p a ), since the kinematic set-up does not allow for threshold collinear radiation into the final state.Therefore the incoming c-PDF momentum is the same as the outgoing threshold collinear momentum.Further simplification arises from the fact that L (1) (z) does not contain explicit factors of the position, such as for instance n − z, thus there are no momentumspace derivatives acting on the collinear momentum delta function, as it happens for the collinear function J 1 in case of the q q channel, cf.Eq. (3.40) in [21].Concerning the color structure, we note that G η,A ξq;γα,f a (n + p, n + p a ; ω) carries one adjoint color index A and two fundamental color indices f a, therefore we can extract a T A f a color generator and transfer it into the definition of the soft function.The collinear function also carries a single Lorentz index η and two Dirac indices γ, α.From the matching equation in (2.28) we see that the Lorentz index is contracted with a ⊥ structure.Therefore, in the collinear function a γ η ⊥ must appear, as the only other possible single Lorentz index carrying structures are n η ± , which would vanish upon contraction with ÂPDF A c⊥η .Based on these considerations, we expect the collinear function G η,A ξq,γα,f a to have the following structure at all orders: where we introduced the scalar collinear function G ξq .At tree level the factor / n − arises due to the collinear quark propagator from the point tn + to z.At higher orders the function G η,A ξq,γα,f a (n + p, n + p a ; ω) is given in terms of loops corrections to the collinear quark propagator, therefore it is constrained by gauge-and reparametrization invariance to have the form on the right-hand side of Eq. (2.30).This can be seen explicitly by noticing that additional factors of γ ⊥ cannot appear, because they would be contracted with the collinear momentum, and thus would be power suppressed.The  n − always appears due to the tree level quark propagator, the rules above can be used to systematically reduce any loop correction to the form in Eq. (2.30).We will show this explicitly at one loop in section 3, and this structure was verified at two loops in [52].
In order to simplify further Eq.(2.26), we note now that the time-ordered product at NLP, both in the amplitude and the complex conjugate amplitude involves the leading power current J A0 , therefore the structure Γ m 1 ,m 2 ρ in Eq. ( (2.31) In the following we will absorb the factor of / n −σα /4 into the definition of the soft function, which we give below.
The last simplification which can be applied when considering Eq. (2.26) up to NLP concerns the phase space integration.As discussed for the q q channel in [21], also the kinematic factors in Eq. (2.26) need to be expanded consistently to NLP: this implies that the matrix element squared expanded to NLP needs to be integrated against the LP expansion of the phase space, while the LP matrix element needs to be multiplied by the phase space expanded up to NLP.These two contributions are identified respectively as "dynamical" and "kinematic" terms: ∆ NLP (z) = ∆ dyn NLP (z) + ∆ kin NLP (z).In case of the g q channel the matrix element squared starts at NLP, therefore there is no kinematic contribution: ∆ g q(z)| NLP = ∆ dyn g q (z)| NLP . (2.32) Taking the LP approximation of the phase space implies the following replacement in Eq. (2.26): ξq .Applying all these simplifications to Eq. (2.26) we finally get where the hard function is also given by the LP approximation of the LP short-distance coefficient squared: ), and the soft function S(Ω, ω, ω ′ ) is given by The L ξq insertion in the amplitude is at position z − , whereas in the conjugate amplitude we place the same insertion at position z ′ − .The conjugate variables to these coordinatespace variables are ω and ω ′ respectively.A graphical representation of the factorization theorem in Eq. (2.34) is given in Fig. 2. We note that the factorization formula in Eq. (2.34), in the same way as the results for the q q-channel Drell-Yan cross-section, is formally valid in d-dimension.The objects appearing in the factorization formula, G ξq (x a n + p A ; ω), G * ξq (x a n + p A ; ω ′ ), and S g q(Ω, ω, ω ′ ), should not be treated as renormalized objects, because the convolution linking the collinear and soft functions must be performed in d-dimensions.We stress that the factorization formula in Eq. (2.34) is valid at the level of the partonic cross section.In particular, in the context of SCET, it was argued in [7] that Glauber modes give scaleless contributions at partonic level and therefore Glauber fields are not included in the construction of the effective theory.In the next two sections we compute the collinear and soft functions, respectively to O(α s ) and O(α 2 s ), which is necessary to achieve NNLO accuracy for the invariant mass distribution.

Collinear functions
We now proceed to calculate the collinear function G ξq (x a n + p A ; ω) up to one loop, by means of the matching equation (2.28).Before continuing, let us notice that an equivalent collinear function has been defined in the context of H → gg in [52], which has been calculated up to two loops.We provide here an independent calculation and check that our result for the one-loop collinear function defined in Eq. (2.28) is indeed equivalent to the one considered in [52].In order to set-up the calculation it proves useful to introduce the short-hand notation Tγf (t) ≡ i d 4 z T χ c,γf (tn + ) L (1) for the left-hand side of (2.28).We then consider its Fourier transform and rewrite the matching equation (2.28) in momentum space as follows ×G η,A ξq;γα,f a (n + q, n + p a ; ω) Starting from this definition, we extract the perturbative collinear functions by considering suitable partonic matrix elements of the operator matching equation.In this case, the relevant matrix element involves an incoming c-PDF gluon and an outgoing soft quark ⟨q + (k)|...|g(p)⟩.We compute both sides of the matching equation with the leading power decoupled Lagrangian, considering the soft fields as external.Hence, we only need the soft matrix element ⟨q + (k)|s ξq;α,a (z − )|0⟩ at tree level.Similarly, for the c-PDF matrix element ⟨0|A PDF A c⊥η (un + )|g(p)⟩ the loop corrections are scaleless and ⟨0|A PDF A c⊥η (un + )|g(p)⟩ contributes only at tree level.The calculation of the left-hand side of Eq. (3.3) is done using the momentum-space Feynman rules given in appendix A of [59].Before moving to the description of the actual calculation, we point out that this calculation is far more straightforward compared to the case of the corresponding collinear functions in the q q channel, see [21].This is because, as discussed above, we need to consider the single soft structure given in Eq. (2.29); furthermore, the Lagrangian insertion L (1) ξq which induces the power suppression does not contain any explicit position variables, which means that, in momentum space, no derivatives on the momentum conserving delta functions appear at subleading power soft-collinear interaction vertices.This reduces drastically the number of terms which appear in each diagram and the number of integrals involved in the calculation.

Tree-level collinear function
We first consider the right-hand side of (3.3) and evaluate the tree-level matrix element where B is the adjoint color index and b is the fundamental color index of the external final state.The c-PDF matrix element in (3.4) evaluates to the following The factor Z g,PDF is the on-shell renormalization factor of the c-PDF gluon field.The soft matrix element on the right-hand side of (3.4) becomes In the second step we have used the form of the soft structure as given in Eq. (2.29).As originally discussed below Eq. ( 4.3) of [21] for the case of collinear functions involving PDF and threshold quarks, the matrix elements in Eqs.(3.5) and (3.6) do not have loop corrections.This is true in Eq. (3.5) because the loop corrections are scaleless, and it also holds for Eq.(3.6) because here the soft field is treated as external.Next, we substitute the matrix elements evaluated in Eqs.(3.5) and (3.6) into (3.4) and find This result constitutes the final expression for the matrix element with the chosen partonic external states of the right-hand side of the matching equation in (3.3).Since Having obtained an expression for the right-hand side of the matching equation (3.4), we now compute the left-hand side, which corresponds to the diagram in Fig. 3.The Lagrangian insertion in Eq. (3.1) can be computed by means of the Feynman rule given in Eq. (A.36) of [59] and we obtain The tree-level value of the on-shell wave function renormalization factor of the gluon field is Z g,c | tree = 1 in the EFT.Since at next-to-leading power in the g q-channel only a single soft structure is relevant, no additional manipulations relating to the use of equation-of-motion identity, or on-shell and transversality conditions are necessary here, in contrast to the calculation of the collinear functions for the q q-channel, see section 4.1 of [21].At this point, we simply compare the result for the left-hand side of the matching equation given in Eq. (3.8) to the right-hand side in Eq. (3.7) and read off the tree-level result for the collinear function Using the decomposition introduced in Eq. (2.30) we can extract the scalar collinear function G ξq (n + p; ω) which appears in the factorization formula in (2.34).Namely, we find where the superscript (0) denotes the tree-level result.

One-loop collinear function
We are now ready to consider the O(α s ) correction.For the right-hand side of the matching we use Eq.(3.7).Here the on-shell wave function renormalization factor is unity to all orders in perturbation theory: Z g,PDF = 1.This holds in dimensional regularization, which we use to treat infrared and ultraviolet divergences, since the loop corrections are scaleless.The situation is the same for the Z g,c factor on the left-hand side of the matching equation, see Eq. (3.8).Given that the calculation is relatively straightforward, we skip the details of the computation and present the results for the one-loop matrix element ⟨q + b (k)|T γf (n + q)|g B (p)⟩ on the left-hand side of the matching equation (3.3), which is given by the sum over the one loop diagrams represented in Fig. 4 (a) -(i).Let us only notice that diagram (e), (g) and (i) are zero, the latter because the collinear loop does not have a scale.Furthermore, diagrams (a) and (b) contribute with a color factor C F , (c) and (d) with a color factor (−C A /2), (h) and (f) are proportional to (C F − C A /2).We obtain where the superscript (1) on the left-hand side indicates that the matrix element is considered at order α s .We match this result to the right-hand side of Eq. (3.7), from which we extract the perturbative matching coefficient, the collinear function, at oneloop order ξq;γα,f b (n + q, n + p; ω) = − Using the decomposition introduced in Eq. (2.30), the scalar collinear function G ξq (n + p; ω) at one loop reads This is the O(α s ) correction to the tree-level collinear function presented in Eq. (3.10), valid to all orders in ϵ.As anticipated it agrees with the result given in [52], see in particular Eqs.
It is interesting to compare the above result for the collinear function appearing in the g q-channel to the collinear functions in the q q-channel, which have been given in their expanded form in equations (4.31) and (4.34) of [21].We note that here, in contrast to J (1) 1 and J (1) 6 , the collinear function G ξq exhibits 1/ϵ 2 poles, and finite logarithms ∝ ln 2 (n + p ω/µ 2 ).At cross-section level, as we will show in section 5 below, these correspond to leading logarithmic contributions appearing in the collinear sector.This fact complicates the adaptation of the resummation treatment developed for the q qchannel [20] to the off-diagonal g q-channel.We provide additional details in section 5.It is noteworthy that the leading pole structure in the above equation is proportional to C F − C A .As discussed in several instances in literature, cf.[26,29,32,[66][67][68], this structure characterizes the conversion of a collinear gluon into a collinear quark by means of an emission of a soft-quark; its factorization into soft and collinear parts typically involves divergent convolution integrals.calculate in this section.To this end it is useful to insert a complete set of states between the T and T products, and use the momentum operator to translate the fields in the anti-time-ordered part.Performing the integration over x 0 in Eq. (2.35) gives rise to an energy conserving delta function.Writing spinor and color indices explicitly one has In what follows we evaluate the soft function at first and second order in α s , which is indicated below by a superscript (1) and (2), respectively.

First order
At one loop the final state contains a single soft quark and the matrix element reads 2) The complex-conjugate matrix element can be derived accordingly.Inserting these into Eq.(4.1) and introducing the notation we easily get

S
(1) which we express in terms of the color factor Tr[

Second order: virtual-real contribution
At second order in α s we need to take into account two different types of corrections: the virtual-real contribution, involving a soft gluon loop in the matrix element (or in the complex conjugate one) and the real-real contribution, involving the emission of a soft gluon in addition to the soft anti-quark.We start by considering the virtual-real contribution.The one-loop matrix element is given by where linear propagators are written such that they carry a small positive imaginary factor +i0 + , and the terms on the right-hand side represent respectively Fig. 6 (f), (c), (a), (b), (d), while (e) is immediately zero.Among these, it turns out that only diagram (c) contributes.We focus on this term and insert the matrix element in Eq. (4.1) with the complex conjugate matrix element taken at lowest order; furthermore, we recall that we need to take into account also the complex conjugate one-loop matrix element times the tree level matrix element, such that the full contribution to the virtual-real part of the soft function reads We evaluate the integrations over k and k 1 following two independent approaches: within the first, we integrate directly the terms appearing in Eq. (4.6); within the second we first reduce such terms to a basis of master integrals, in which case Eq. (4.6) reads where the master integrals Ĵi are defined in Eq. (A.3) and calculated in Eqs.(A.4) and (A.5) of appendix A.1.Using the results given there we get

Second order: real-real contribution
We consider now the final state emission of a soft gluon in addition to the soft-antiquark.The matrix element reads .
In practice we have four contributions: three diagrams where the gluon is taken from one of the three Wilson lines in Eq. (4.9), counting also the Wilson line implicit in q + = Y † + q s , and the fourth where the gluon is taken from the insertion of the (leading power) soft Lagrangian L (0) s (y) = (y)g s / A s (y)q s (y).The soft function is obtained according to Eq. (4.1) where the complex conjugate matrix element can be derived from Eq. (4.9).In Fig. 7 we list the diagrams which give a non-zero contribution.We obtain As for the virtual-real contribution, we evaluate the integrals over k 1 and k 2 both by directly integrating the expression in Eq. (4.10), and by reducing such expression to a basis of master integrals, in which case Eq. (4.10) becomes The master integrals Îi are defined and evaluated in appendix A.2.After substituting the expressions for the master integrals in Eqs.(A.14)-(A.17)we finally obtain Eqs. (4.8) and (4.12) provide the complete result for the two-loop soft function.

Asymptotic limits
With the results of the NNLO soft function calculation at hand, we now have the opportunity to study the asymptotic limits of this object.This is useful, because endpoint divergences arise exactly in these limits, thus understanding the emergent asymptotic structure of the soft function is useful for any future study of endpoint refactorization.Indeed, as shown for instance for the off-diagonal "gluon" thrust in [32], the asymptotic limits of the soft (and collinear) functions are the objects one needs to subtract from the full functions, in order to define matrix elements free of endpoint divergences.After this procedure has been completed, it is then possible to implement the standard resummation procedure, based on the renormalization-group evolution of the subtracted functions.
At NLO the soft function is relatively simple, and inspecting Eq. (4.4) we see that an endpoint divergence arises for ω, ω ′ → 0: In particular, notice that, because of the factor δ(ω−ω ′ ), one has effectively a dependence on a single variable ω.This is analogous to the situation in [32] where the relevant asymptotic limit is ω, ω ′ → ∞, however identical simplification occurs, namely, the presence of a Dirac Delta function reduces the functional dependence of the soft function to a single ω variable (see e.g.Eq. ( 27) in [32]).The calculation of the soft function at two loops in sections 4.2, 4.3 gives us the opportunity to explore its asymptotic limits for the first time beyond NLO.Indeed, we find a more involved structure of endpoint singularities, which arise for ω → 0, ω ′ → 0 separately, as well as for ω ′ − ω → 0. To be more specific, let us write the soft function as follows: where for simplicity we have dropped the dependence on Ω and the scale µ.In Eq. (4.14) we classify the various contributions in terms of their domain of integration.The term Ŝ(2) (ω) corresponds to the first contribution of the double-real correction in Eq. (4.12): and the domain of integration for this term is the same as the one of the NLO soft function.We then have S (2A) (ω, ω ′ ) and S (2B) (ω, ω ′ ), which correspond respectively to the first and second term of the virtual-real contribution in Eq. (4.8): Last, we have the double real contribution of Eq. (4.12) not proportional to δ(ω − ω ′ ): The integration domains of the three functions are represented in Fig. 8, with the blue shaded areas representing the endpoint regions where singularities occur.In this regard, let us notice that a further divergence for ω → ∞ would be present for the factors S (2A) (ω, ω ′ ) and S (2B) (ω, ω ′ ) taken separately, but it cancels in their sum.In the boundary regions of Fig. 8, taking the limit ω, ω ′ → 0 with ω ∼ ω ′ we have These limits provide valuable information that will be useful to set up a refactorization procedure such as the one devised in [32].This analysis goes beyond the scope of this paper, and we leave it for future work.

Comparison to fixed order results
We can now use the one loop collinear and the two loop soft functions, which appear in the factorization theorem Eq. (2.34), to evaluate the partonic cross section up to the second order in perturbation theory.We check the result by comparing the bare cross section with an in-house calculation of the NLP NNLO partonic cross section obtained with the method of regions [69].Next we remove the initial-state collinear singularities via PDF renormalization and compare with the finite cross section available in the literature [70].

First order
The partonic cross section in the g q channel starts at first order in perturbation theory.At this order, taking into account Eq. (3.10), and the fact that H(Q 2 ) = 1 + O(α s ), we have ∆ g q (z)| NLP = 2 dω dω ′ S (1) (Ω, ω, ω ′ ) . (5.1) We replace the soft function in Eq. (4.4) and integrate over ω, ω ′ .The exact result reads ∆ (1) Expanding in ϵ, identifying Ω = Q(1 − z) and setting for simplicity µ = Q we get ∆ (1) where, in order to keep equations compact, we introduced the notation Let us notice here that the one loop soft function in Eq. (4.4) is finite for ω ̸ = 0, therefore the single pole in Eq. ( 5.3) arises from the integration over ω, for ω → 0: we see explicitly that the expansion in powers of ϵ and renormalization before taking the convolution in ω, ω ′ would result in an endpoint divergence.

Second order
At second order we need to take into account three contributions.The first involves the hard function at one loop: ∆ g q (z)| NLP,h = 2H (1) (Q 2 ) dω dω ′ S (1) (Ω, ω, ω ′ ) . (5.5) The integration over the soft function is the same as the one occurring at NLO, and the one loop hard function H (1) (Q 2 ), where superscript (1) indicates order α s result, can be found in Eq. (5.6) of [21] 4 .Integrating over ω, ω ′ , the result with exact scale dependence can be expressed as: ∆ Expanding in powers of ϵ, identifying Ω = Q(1 − z) and setting µ = Q we get ∆ (2) This result is in agreement with the in-house method of regions calculation where the virtual gluon is hard and the emitted quark is soft.
The second contribution involves the collinear function at one loop: ξq (x a n + p A ; ω) ξq (x a n + p A ; ω) S (1) (Ω, ω, ω ′ ) , (5.8) which, taking into account the result for the tree level jet function, simplifies to ∆ (2) ξq (x a n + p A ; ω) + G (1) * ξq (x a n + p A ; ω ′ ) S (1) (Ω, ω, ω ′ ) .(5.9) Inserting the one loop jet function from Eq. (3.13), the one loop soft function Eq. (4.4) and integrating we get Expanding in powers of ϵ, identifying Ω = Q(1 − z) and setting µ = Q we get ∆ (2) This result also agrees with the in-house method of regions calculation, where the virtual gluon is collinear and the emitted quark is soft.As discussed after Eq. (3.14), let us highlight that, in contrast to the q q channel, where the collinear function starts contributing only at NLL accuracy, in this case it contributes already at LL level.Comparing Eqs.(3.14) and (5.11) we see more in detail how this happens.Focusing on the highest pole contribution, which is in direct correspondence with the LL, we see that the collinear function itself contains a 1/ϵ 2 pole, which becomes a 1/ϵ 3 leading pole by means of the convolution with the one loop soft function.In this respect, the collinear contribution in the q q channel contains a 1/ϵ pole (see Eqs. (4.30) -(4.34) in [21]), which is then raised to a 1/ϵ 2 subleading pole by the convolution integration.A physical interpretation of this phenomenon has been put forward in [26,29].In general, the idea is that in a splitting 1 → 2 + 3 with soft 3, the leading pole is related to the difference of the Casimir charge of the energetic particles 1 and 2. In the diagonal q q channel the collinear quark 1 emits a soft gluon 3, continuing as the collinear quark 2, thus there is no change of the Casimir charge along the collinear direction.In the g q channel, a collinear gluon 1 emits a soft quark 3, continuing as a collinear quark, which has a different Casimir charge.Indeed, the collinear function in Eq. (3.14) is proportional to the difference of the Casimir charges C F − C A .Last, the third contribution involves the soft function at two loops: ∆ g q (z)| NLP,s = 2 dω dω ′ S (2) (Ω, ω, ω ′ ) . (5.12) As discussed in section 4, the soft function at two loops receives two contributions, where the additional soft gluon is respectively virtual or real.These two contributions have been given respectively in Eqs.(4.8) and (4.12).For a better comparison with the method of regions, it proves useful to evaluate the two contributions separately.Furthermore, the integration of the term involving the virtual gluon, Eq. (4.8), is comparatively more involved.In particular, the term in the last line of Eq. (4.8) requires some care, because of the singularity at ω = 0, which can be integrated by considering the standard identity 1/(ω − i0 + ) = P (1/ω) + iπδ(ω), where P indicates the principal value prescription.The term in the second line of Eq. (4.8) instead can be integrated over ω and ω ′ without particular issues.In the end we obtain a result for Eq.(5.12) valid to all orders in ϵ: The integration of the double real contribution, Eq. (4.12) is straightforward.Again, we obtain a result valid to all order in ϵ: which upon expansion in ϵ, identifying Ω = Q(1 − z) and setting µ = Q gives ∆ Summing together all terms contributing at NNLO we obtain ∆ This gives our result for the bare partonic cross section at NNLO in perturbation theory.The result agrees with an in-house computation of the bare partonic cross section obtained with the method of regions.This provides us with a first confirmation that the (bare) factorization theorem in Eq. (2.34) is correct, and is not missing any contribution up to O(α 2 s ).

Renormalization
As a final check, we need to compare with the known result in literature [70], for which we need the finite partonic cross section.As discussed in appendix B, the latter is obtained by applying PDF renormalization to Eq. (5.17), which relies on the finiteness of the hadronic cross section: where by "ren" we indicate the finite, renormalized functions.In practice, the finite cross section is obtained by applying a counterterm to the bare cross section.We provide these counterterms as a function of the Altarelli-Parisi splitting kernels in Eq. (B.8).At first order in perturbation theory, the explicit calculation gives qg /(ϵ) in the second line of Eq. (B.8).Proceeding in a similar way and evaluating explicitly the counterterm in the third equality of Eq. (B.8) we have: and in the end we get ∆ (2)

Summary
In this work, we derived the factorization formula for the g q-channel of the Drell-Yan production process at general powers in the threshold expansion.Specifying the result to the first non-trivial power, which is at next-to-leading power in (1 − z) compared to the leading q q channel, we arrived at a compact formula comprising of a LP hard function, two NLP collinear functions on each side of the cut and one generalized soft function.The simplicity of this result is in contrast with the more involved structure of the q q case, where four separate soft functions contribute along with their corresponding collinear functions [21].The first main result of this work is the NLP factorization formula given in Eq. (2.34).
In sections 3 and 4 we calculated the collinear and soft functions which appear in the factorization formula.In order to validate this formula to NNLO accuracy, hard and collinear functions are needed up to O(α s ) while the soft function up to O(α 2 s ).The hard function which appears in Eq. (2.34) is the LP hard function which is known in the literature.In section 3, we calculated the NLP collinear function with an external PDFcollinear gluon which we have defined in an analogous way to the collinear functions with an external PDF-collinear quark appearing in [21].Due to its relation to the radiative jet functions appearing in h → gg amplitude, the radiative corrections of this object are known up to O(α 2 s ) [52].We found agreement for the O(α s ) contribution which is given in Eq. (3.13).
In section 4, we calculated the soft function which contains soft quark insertions on both sides of the cut.We carried out the computation by utilising state-of-the-art fixed-order loop integral methods such as reduction to master integrals and application of differential equations.The O(α s ) result for the soft function is given in Eq. (4.4), and two-loop soft function results are presented in Eq. (4.8) and Eq.(4.12) for the one-real one-virtual and double-real contributions, respectively.We have retained all-order ϵ dependence throughout the calculation of the NLO and NNLO contributions to the soft function, which enabled us to study the structure of its asymptotic limits in section 4.4.In particular, we noticed a relatively simple structure at NLO, where singularities arise for ω → 0.More involved structure emerges at NNLO which has implications for resummation beyond LL using refactorization procedures such as the ones developed in [32,44].
Using the calculations in sections 3 and 4, we tested the validity of the factorization formula derived in section 2 through a comparison of our results against in-house expansion-by-regions calculations and known results in the literature.We evaluated the convolutions over the ω, ω ′ variables obtaining results valid in exact d-dimensions at O(α s ) and O(α 2 s ).We then expanded the final results to d → 4, carried out PDFfactorization of the collinear divergences appearing in initial states, which is standard for fixed order calculations, and we compared with [70], finding agreement.
In summary, the bare factorization theorem derived here, along with the higher-order calculations of the functions appearing in our formula and the study of the asymptotic limits of the soft function, will serve as a spring board for future investigations of this channel and an eventual desirable four dimensional renormalization prescription.and the last two belong to topology C and D Î5 (Ω, ω) ≡ ÎC (0, 0, 0, 1, 1, 1, 1), Î6 (Ω, ω) ≡ ÎD (0, 0, 0, 1, 1, 1, 1), (A.13) respectively.Î1 can be easily extracted from the results in [51], we find (A.14) Î2 was not explicitly required for the calculations carried out in [51], hence we evaluate it for the first time here by employing the differential equation method.Î2 appears in a system with Î1 and the differential equation can be solved to all orders in ϵ.The ϵ-dependent integration constant can be fixed by matching to the version of this integral where ω and ω ′ are integrated over first.We find that the constant is zero to all-orders in ϵ and the final integral reads (A.17)

B PDF renormalization
In this appendix we briefly set our notation for PDF renormalization.Let us start from Eq. (2.3), and write it in the form (z) + P (0) gg ⊗ ∆ bare (1) g q (z) 6 Notice that by symmetry we take q = q in the indices of Γ ab and P ab .
+ P (0) qg ⊗ ∆ bare (1) q q (z) − β 0 ∆ bare (1) g q (z) , (B.8) where the last term originates from expressing the bare cross section as an expansion in terms of the renormalized strong coupling constant.Near x → 1 we can expand the Altarelli-Parisi splitting kernels in powers of 1 − x, and up to NLP we need ) where f a/A (x a ) and f b/B (x b ) are the usual parton distribution functions (PDFs) of partons a, b, with momentum fractions x a , x b ; σab (z) represents the partonic cross section for the process a(p a )b(p b ) where Ω = Q(1 − z).As a consequence of this approximation, we can set ⃗ x = 0 in the argument of the position space soft function.Last, we observe that we can carry out a further simplification by setting the anticollinear function to its LP expression, Jī (n − p,x b n − p B ; ω) = δ(n − p − x b n − p B ), given that the suppression to NLP is already provided by the two insertions of L(1)

Figure 2 :
Figure 2: Graphical representation of the factorization theorem in Eq.(2.34), valid for the g q channel in Drell-Yan near threshold at NLP.In this picture red lines represent soft fields, blue lines represent collinear fields, and green lines are anti-collinear fields.The purple circles denoted by "A0" and "A0 * " represent the leading power short-distance coefficient C A0 (Q 2 ) and its complex conjugate, respectively.The blue ovals denoted by G, G * represents the (momentum-space representation of the) collinear functions, introduced in the matching equation Eq.(2.15).The red double lines represent the soft Wilson lines and the red filled circles are insertions of the soft quark building blocks.Last, ω and ω ′ are the convolution variables between the collinear functions and the soft function.

Figure 3 :
Figure 3: Tree-level EFT diagram in the g q-channel.

Figure 4 :
Figure 4: One-loop diagrams contributing to the collinear function.

Figure 5 :
Figure 5: Diagram contributing to the soft function at first order (NLO in α s ).

Figure 6 :
Figure 6: Diagrams contributing to the virtual-real part of the soft function at second order in α s .The part to the left (right) of the cut corresponds to the time-ordered (antitime-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 in this figure stands for the soft quarks and the Wilson lines contained in q + = Y † + q s .

Figure 7 :
Figure 7: Diagrams contributing to the real-real part of the soft function.The part to 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 in this figure stands for the soft quarks and the Wilson lines contained in q + = Y † + q s .

dσ DY dQ 2 = σ 0 a,b 1 τ= σ 0 a,b 1 τ 1 0dy 1 6 ∆
dz dx a dx b δ(τ − x a x b z) f bare a/A (x a , µ) f bare b/B (x b , µ) ∆ bare ab (z, µ), (B.1)where we added a superscript to indicate that all functions are unrenormalized 5 .Due to renormalization invariance of the differential cross section, Eq. (B.1) can be expressed as well in terms of the corresponding renormalized function:dσ DY dQ 2 dz dx a dx b δ(τ − x a x b z) f ren a/A (x a , µ) f ren b/B (x b , µ) ∆ ren ab (z, µ).(B.2)The relation between the bare and renormalized partonic cross section can be obtained by considering the relation between bare and renormalized PDFs:f ren a/A = f bare b/A ⊗ Γ ab , (B.3)where the convolution explicitly readsf ren a/A (x a ) = dy 2 δ(x a − y 1 y 2 ) f bare (y 1 ) b/A Γ ab (y 2 ).(B.4)In turn, the function Γ ab (x) has the following perturbative expansionΓ ab (x) = δ ab δ(1 − x) + Eq. (B.3) into Eq.(B.2) we get the relation∆ bare cd (z) = Γ ac ⊗ Γ bd ⊗ ∆ ren ab (z).(B.7)This equation can then be solved order by order in α s to obtain the renormalized cross section coefficients in terms of the bare ones.Taking into account the UV renormalization and the explicit form of Eqs.(B.5) and (B.6), for the g q channel one has qg (z) + P (0) qq ⊗ ∆ bare (1) g q