Gauge-invariant TMD factorization for Drell-Yan hadronic tensor at small x

The Drell-Yan hadronic tensor for electromagnetic (EM) current is calculated in the Sudakov region s≫Q2≫q⊥2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ s\gg {Q}^2\gg {q}_{\perp}^2 $$\end{document} with 1Q2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \frac{1}{Q^2} $$\end{document} accuracy, first at the tree level and then with the double-log accuracy. It is demonstrated that in the leading order in Nc the higher-twist quark-quark-gluon TMDs reduce to leading-twist TMDs due to QCD equation of motion. The resulting tensor for unpolarized hadrons is EM gauge-invariant and depends on two leading-twist TMDs: f1 responsible for total DY cross section, and Boer-Mulders function h1⊥\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {h}_1^{\perp } $$\end{document}. The order-of-magnitude estimates of angular distributions for DY process seem to agree with LHC results at corresponding kinematics.

The differential cross section of DY process is determined by the product of leptonic tensor and hadronic tensor. The hadronic tensor W µν is defined as (1.1) where p A , p B are hadron momenta, q is the momentum of DY pair, X denotes the sum over full set of "out" states and J µ is either electromagnetic or Z-boson current. In this paper I consider only the case of electromagnetic current, the Z-boson case will be studied in a separate publication. For unpolarized hadrons, the hadronic tensor W µν is parametrized by 4 functions, for example in Collins-Soper frame [18] where X, Z are unit vectors orthogonal to q and to each other (their explicit form is presented in Sect. 8.2). Conventionally, the analysis of hadronic tensor (1.1) in the Sudakov region q 2 ≡ Q 2 q 2 ⊥ is performed by using TMD factorization. For example, functions W T and W ∆∆ can be represented in a standard TMD-factorized way [9,19] is the TMD density of a parton f in hadron A with fraction of momentum x A and transverse momentum k ⊥ , D f /B (x B , q ⊥ − k ⊥ ) is a similar quantity for hadron B, and coefficient functions C i (q, k) are determined by the cross section σ(f f → µ + µ − ) of production of DY pair of invariant mass q 2 in the scattering of two partons. There is, however, a problem with Eq. (1.3) for the functions W L and W ∆ . The reason is that while W T and W ∆∆ are determined by leading-twist quark TMDs, W L and W ∆ start from terms q ⊥ Q and ∼ q 2 ⊥ Q 2 determined by quark-quark-gluon TMDs. The power corrections ∼ q ⊥ Q were found in Ref. [20] more than two decades ago but there was no calculation of power corrections ∼ q 2 ⊥ Q 2 until recently. Also, the leading-twist contribution is not gauge invariant. 1 It is well known from DVCS studies that check of EM gauge invariance sometimes involves cancellation of contributions of different twists (see e.g. [21][22][23][24][25][26][27]) so the fact that we need power corrections to check q µ W µν = 0 should not come as a surprise. Still, the absence of gauge invariance may cause discomfort in practical applications of TMD factorization.
In a recent paper [28] A. Tarasov and the author calculated power corrections ∼ q 2 ⊥ Q 2 to total DY cross section production which are determined by quark-quark-gluon operators. In this paper I present the result of calculation of symmetric part of W µν (q) for unpolarized hadrons at large s Q 2 q 2 ⊥ relevant for DY experiments at LHC. The method of calculation is based on the rapidity factorization approach developed in Refs. [28,29]. The calculations will be performed in the leading order in perturbation theory, first at the tree level and then in the double-logarithmic approximation for coefficient functions C i (q, k). In this paper I consider only the production of leptons by virtual photon and leave the case of Z-boson production for future publication.
To find all functions in Eq. (1.2) we need to have gauge-invariant expression for W µν in terms of TMDs. As noted above, only W T and W ∆∆ come from leading-twist quarkantiquark TMD while two other structures come from higher-twist quark-antiquark-gluon TMDs. Fortunately, in the leading order in N c the latter are related to the former by QCD equations of motion ( [28], see also Ref. [20]). Moreover, in the small-x region x A , x B 1 all structures can be expressed by just two leading-twist TMDs -f 1 (x, k ⊥ ) (responsible for the total cross section) and h ⊥ 1 (x, k ⊥ ) (the Boer-Mulders function [30]). The results for four functions in Eq. (1.2), presented in next Section, are of the type of Eq. (1.3) with TMDs f 1 (x, k ⊥ ) and/or h ⊥ 1 (x, k ⊥ ) and tree-level coefficient functions constructed of q and k ⊥ .
The paper is organized as follows. In section 2 I present the resulting gauge-invariant expression for W µν up to 1 Q 2 terms which is calculated in the rest of the paper. In section 3 the TMD factorization is derived from the rapidity factorization of the double functional integral for a cross section of particle production. In section 4 I explain the method of calculation of power corrections based on approximate solution of classical Yang-Mills equations. Using this method, DY hadronic tensor for small x is calculated in Sections 5, 6, and 7. Section 8 contains results of calculations and order-of-magnitude estimate of angular coefficients of DY cross section. The matching of obtained TMDs and coefficient functions C i in the double-log approximation is discussed in Sect. 9 and the last section 10 is devoted to conclusions and outlook. The necessary technical details can be found in appendices.

Gauge-invariant hadronic tensor
To set up the stage, in this Section I present the final result for tree-level DY hadronic tensor. It is determined by two leading-twist TMDs: the function f f 1 (x, k ⊥ ) responsible for the total DY cross section and Boer-Mulders time-odd function h ⊥ 1 (x, k ⊥ ) (the explicit definition of these functions is presented in the Appendix 11.2). The result reads The gauge-invariant structures W F µν and W H µν are given by where g ⊥ µν and g µν are transverse and longitudinal parts of metric tensor (the explicit definition is given by Eq. (3.2)). It is easy to check that q µ W F µν = 0 and q µ W H µν = 0. As we will see below, in some of the structures the corrections to Eq. Nc times some other higher-twist TMDs discussed in Ref. [28]. It should be also noted that W F part coincides with the result obtained in Refs. [31,32] using parton Reggeization approach to DY process [33].
In the rest of the paper I will derive the above equations and discuss their accuracy. Let me mention upfront that since the approximations made in Eq. (2.1) are x A , x B 1 and q 2 ⊥ Q 2 x A x B s, I hope that the results of this paper can be used for studies of DY process at LHC with Q 2 ∼ 100GeV or less. 2 Last but not least, the derivation of the above equations is lengthly so the readers interested in final formulas for structures W i and the discussion of approximations can go directly to Sect. 8.

TMD factorization from rapidity factorization
As was mentioned in the Introduction, to find the TMD formulas of Eq. (1.3) type I use the rapidity factorization approach to developed in Refs. [28,29]. Let me quickly summarize basic ideas of this approach. The sum over full set of "out" states in Eq. (1.1) can be represented by a double functional integral where J µ = flavors e fψf γ µ ψ f is the electromagnetic current. In this double functional integral the amplitude X|J µ (0)|p A , p B is given by the integral over ψ, A fields whereas the complex conjugate amplitude p A , p B |J µ (x)|X is represented by the integral overψ,Ã fields. Also, Ψ p ( A(t i ), ψ(t i )) denotes the proton wave function at the initial time t i and the boundary conditionsÃ(t f ) = A(t f ) andψ(t f ) = ψ(t f ) reflect the sum over all states X, cf. Refs. [34][35][36].
We use Sudakov variables p = αp 1 + βp 2 + p ⊥ , where p 1 and p 2 are light-like vectors close to p A and p B so that p A = p 1 + m 2 s p 2 and p A = p 1 + m 2 s p 2 with m being the proton mass. Also, we use the notations x • ≡ x µ p µ 1 and x * ≡ x µ p µ 2 for the dimensionless light-cone coordinates (x * = s 2 x + and x • = s 2 x − ). Our metric is g µν = (1, −1, −1, −1) which we will frequently rewrite as a sum of longitudinal part and transverse part: Throughout the paper, the sum over the Latin indices i, j, ... runs over two transverse components while the sum over Greek indices µ, ν, ... runs over four components as usual.
Following Ref. [29] we separate quark and gluon fields in the functional integral (3.1) into three sectors (see figure 1): "projectile" fields A µ , ψ A with |β| < σ p , "target" fields B µ , ψ B with |α| < σ t and "central rapidity" fields C µ , ψ C with |α| > σ t and |β| > σ p , see Fig. 1. 3 Our goal is to integrate over central fields and get the amplitude in the factorized "Central" fields < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > "Projectile 00 fields : | | < p < l a t e x i t s h a 1 _ b a s e 6 4 = " v D / 0 s 3 t 0 P 1 M C V Q L z h W f z d U f c I j I = " > A A A C F n i c b V B N S 8 N A E N 3 U 7 / p V 9 e h l s U g 9 l U T 8 x o P o x W M F Y w t N a D f b S b u 6 m 4 T d j V B i / R V e / C t e P K h 4 F W / + G 7 e 1 h 1 p 9 M P B 4 b 4 a Z e U H C m d K 2 / W X l J i a n p m d m 5 / L z C 4 t L y 4 W V 1 S s V p 5 K C S 2 M e y 1 p A F H A W g a u Z 5 l B L J B A R c K g G N 2 d 9 v 3 o L U r E 4 u t T d B H x B 2 h E L G S X a S I 1 C O f O k w M 1 m R c b X Q D X j U C r d 4 5 A B b 6 n e E b 7 z A t D k 7 t h T r C 1 I I 2 k U i n Z 5 1  y E U U P 6 A m 9 o F f r 0 X q 2 3 q z 3 n 9 a c N Z x Z Q 7 9 g f X w D g s C f v w = = < / l a t e x i t > p A < l a t e x i t s h a 1 _ b a s e 6 4 = " B R X J d i T P e u s l g x 0 H c l W E l X l A 4 k M = " > A A A B 6 X i c b V B N T w I x E J 3 i F + I X 6 t F L I z H x R H a N X 9 x Q L x 4 x u k I C G 9 I t X W j o d j d t 1 4 R s + A l e P K j x 6 j / y 5 r + x w B 4 Q f c k k L + / N Z G Z e k A i u j e N 8 o 8 L S 8 s r q W n G 9 t L G 5 t b 1 T 3 t 1 7 1 H G q K P N o L G L V C o h m g k v m G W 4 E a y W K k S g Q r B k M b y Z + 8 4 k p z W P 5 Y E Y J 8 y P S l z z k l B g r 3 S f d q 2 6 5 4 l T P H L d 2 7 u I 5 4 k y B 3 Z x U I E e j W / 7 q 9 G K a R k w a K o j W b d d J j J 8 R Z T g V b F z q p J o l h A 5 J n 7 U t l S R i 2 s + m p 4 7 x k V V 6 O I y V L W n w V J 2 f y E i k 9 S g K b G d E z E A v e h P x P 6 + d m v D S z 7 h M U s M k n S 0 K U 4 F N j C d / 4 x 5 X j B o x s o R Q x e 2 t m A 6 I I t T Y d E o 2 B H f x 5 b / E O 6 n W q s 7 d a a V + n a d R h A M 4 h G N w 4 Q L q c A s N 8 I B C H 5 7 h F d 6 Q Q C / o H X 3 M W g s o n 9 m H X 0 C f P 8 V t j a c = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " B R X J d i T P e u s l g x 0 H c l W E l X l A 4 k M = " > A A A B 6 X i c b V B N T w I x E J 3 i F + I X 6 t F L I z H x R H a N X 9 x Q L x 4 x u k I C G 9 I t X W j o d j d t 1 4 R s + A l e P K j x 6 j / y 5 r + x w B 4 Q f c k k L + / N Z G Z e k A i u j e N 8 o 8 L S 8 s r q W n G 9 t L G 5 t b 1 T 3 t 1 7 1 H G q K P N o L G L V C o h m g k v m G W 4 E a y W K k S g Q r B k M b y Z + 8 4 k p z W P 5 Y E Y J 8 y P S l z z k l B g r 3 S f d q 2 6 5 4 l T P H L d 2 7 u I 5 4 k y B 3 Z x U I E e j W / 7 q 9 G K a R k w a K o j W b d d J j J 8 R Z T g V b F z q p J o l h A 5 J n 7 U t l S R i 2 s + m p 4 7 x k V V 6 O I y V L W n w V J 2 f y E i k 9 S g K b G d E z E A v e h P x P 6 + d m v D S z 7 h M U s M k n S 0 K U 4 F N j C d / 4 x 5 X j B o x s o R Q x e 2 t m A 6 I I t T Y d E o 2 B H f x 5 b / E O 6 n W q s 7 d a a V + n a d R h A M 4 h G N w 4 Q L q c A s N 8 I B C H 5 7 h F d 6 Q Q C / o H X 3 M W g s o n 9 m H X 0 C f P 8 V t j a c = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " B R X J d i T P e u s l g x 0 H c l W E l X l A 4 k M = " > A A A B 6 X i c b V B N T w I x E J 3 i F + I X 6 t F L I z H x R H a N X 9 x Q L x 4 x u k I C G 9 I t X W j o d j d t 1 4 R s + A l e P K j x 6 j / y 5 r + x w B 4 Q f c k k L + / N Z G Z e k A i u j e N 8 o 8 L S 8 s r q W n G 9 t L G 5 t b 1 T 3 t 1 7 1 H G q K P N o L G L V C o h m g k v m G W 4 E a y W K k S g Q r B k M b y Z + 8 4 k p z W P 5 Y E Y J 8 y P S l z z k l B g r 3 S f d q 2 6 5 4 l T P H L d 2 7 u I 5 4 k y B 3 Z x U I E e j W / 7 q 9 G K a R k w a K o j W b d d J j J 8 R Z T g V b F z q p J o l h A 5 J n 7 U t l S R i 2 s + m p 4 7 x k V V 6 O I y V L W n w V J 2 f y E i k 9 S g K b G d E z E A v e h P x P 6 + d m v D S z 7 h M U s M k n S 0 K U 4 F N j C d / 4 x 5 X j B o x s o R Q x e 2 t m A 6 I I t T Y d E o 2 B H f x 5 b / E O 6 n W q s 7 d a a V + n a d R h A M 4 h G N w 4 Q L q c A s N 8 I B C H 5 7 h F d 6 Q Q C / o H X 3 M W g s o n 9 m H X 0 C f P 8 V t j a c = < / l a t e x i t > p B < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > "Target" fields : |↵| < t < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > Figure 1. Rapidity factorization for DY particle production form, i.e. as a product of functional integrals over A fields representing projectile matrix elements (TMDs of the projectile) and functional integrals over B fields representing target matrix elements (TMDs of the target). In the spirit of background-field method, we "freeze" projectile and target fields and get a sum of diagrams in these external fields. Since |β| < σ p in the projectile fields and |α| < σ t in the target fields, at the tree level one can set with power accuracy β = 0 for the projectile fields and α = 0 for the target fields -the corrections will be O m 2 σps and O m 2 σts . Beyond the tree level, the integration over C fields produces logarithms of the cutoffs σ p and σ t which match the corresponding logs in TMDs of the projectile and the target, see the discussion in Sect. 9 From integrals over projectile and target fields in the above equation we see that the functional integral over C fields should be done in the background of A and B fields satisfying Combining this with our approximation that at the tree level β = 0 for A,Ã fields and α = 0 for B,B fields, which corresponds to , we see that for the purpose of calculation of the functional integral over central fields we can set In other words, since A, ψ andÃ,ψ do not depend on x * , if they coincide at x * = ∞ they coincide everywhere. Similarly, since B, ψ B andB,ψ B do not depend on x • , if they coincide at x • = ∞ they should be equal. Summarizing, we see that at the tree level in our approximation It is well known that in the tree approximation the double functional integral (3.5) is given by a set of retarded Green functions in the background fields [37][38][39] (see also appendix A of ref. [29] for the proof). Since the double functional integral (3.5) is given by a set of retarded Green functions in the background fields A and B, the calculation of the tree-level contribution toψγ µ ψ in the r.h.s. of Eq. (3.5) is equivalent to solving YM equation for ψ(x) and A µ (x) with initial condition that the solution has the same asymptotics at t → −∞ as the superposition of incoming projectile and target background fields. The hadronic tensor (1.1) can now be represented as whereÔ µν (q, x;Â,ψ A ;B,ψ B ) should be expanded in a series inÂ,ψ A ,B,ψ B operators and evaluated between the corresponding (projectile or target) states: if where c m,n are coefficients and Φ can be any of A µ , ψ orψ with appropriate Lorentz indices. We get then As we will demonstrate below, the relevant operatorsΦ A andΦ B are quark and gluon fields with Wilson-line type gauge links collinear to either p 2 for A fields or p 1 for B fields.

Power counting for background fields
As we discussed in previous section, to get the hadronic tensor in the form (3.6) we need to calculate the functional integral (3.5) in the background of the fields (3.4). Since we integrate over fields (3.4) afterwards, we may assume that they satisfy Yang-Mills equations and similarly for B fields. 5 It is convenient to choose a gauge where A * = 0 for projectile fields and B • = 0 for target fields. (The existence of such gauge was proved in appendix B of Ref. [29] by explicit construction.) The relative strength of Lorentz components of projectile and target fields in this gauge was found in ref. [29] / Here m ⊥ is a scale of order of m or q ⊥ . As discussed in Refs. [28,29], our rapidity factorization (3.8) is applicable in the region where s, Q 2 q 2 ⊥ , m 2 N , while the relation between q 2 ⊥ and m 2 and between Q 2 and s may be arbitrary. Correspondingly, for the purpose of counting of powers of s, we do not distinguish between s and Q 2 so our power counting will be correct at any Bjorken x. The distinction will come at a later time when we specify to small x and disregard 1 s in comparison to 1 Q 2 in final expressions for TMDs and/or coefficient functions. Similarly, for the purpose of power counting we will not distinguish between m and q ⊥ so we introduce m ⊥ which may be of order of m or q ⊥ depending on matrix element. 4 As we mentioned, for the purpose of calculation of integral over C fields the projectile and target fields are "frozen". 5 Since we are dealing with tree approximation and quark equations of motion, it is convenient to include coupling constant g in the definition of gluon fields.

Note also that in our gauge
•i are field strengths for A and B fields respectively. Thus, to find TMD factorization formula with power corrections at the tree level we need to calculate the functional integral (3.1) in the background fields of the strength given by eqs. (4.2).

Approximate solution of classical equations at
As we discussed in section 3, the calculation of the functional integral (3.5) over C-fields in the tree approximation reduces to finding fields C µ and ψ C as solutions of Yang-Mills equations for the action The solution of eq. (4.4) which we need corresponds to the sum of set of diagrams in background field A + B with retarded Green functions, see figure 2. The sum of tree diagrams with retarded Green functions gives fields C µ and ψ C that vanish at t → −∞. Thus, we are solving the usual classical YM equations 6 (4.6) 6 We take into account only u, d, s, c quarks and consider them massless. In principle, one can include with boundary conditions following from C µ , ψ C t→−∞ → 0. These boundary conditions reflect the fact that at t → −∞ we have only incoming hadrons with A and B fields.
As discussed in Ref. [29], for our case of particle production with q ⊥ Q 1 it is possible to find the approximate solution of (4.5) as a series in this small parameter. One solves Eqs. (4.5) iteratively, order by order in perturbation theory, starting from the zero-order approximation in the form of the sum of projectile and target fields and improving it by calculation of Feynman diagrams with retarded propagators in the background fields (4.8).
Let me now explain how the parameter m 2 ⊥ /s comes up in the rapidity-factorization approach (for details, see Ref. [29]). When we expand quark and gluon propagators in powers of background fields, we get a set of diagrams shown in figure 2. The typical bare gluon propagator in figure 2 is . (4.9) In the tree approximation, the transverse momenta in tree diagrams are determined by further integration over projectile ("A") and target ("B") fields in eq. (3.1) which converge on either q ⊥ or m N . On the other hand, the integrals over α converge on either α q or α ∼ 1 and similarly the characteristic β's are either β q or β ∼ 1. Since α q β q s = Q 2 q 2 ⊥ , one can expand gluon and quark propagators in powers of After the expansion (4.10), the dynamics in the transverse space effectively becomes trivial: all background fields stand either at x or at 0. Note that in this statement is solely a consequence of Q 2 q 2 ⊥ and does not rely on small-x approximation.

Power expansion of classical quark fields
Now we expand the classical quark fields in powers of s (the corresponding expansion of classical gluon fields is presented in Ref. [29], but we do not need it here). As demonstrated in Ref. [28], expanding it in powers of p 2 ⊥ /p 2 we obtain where and dots stand for terms subleading in q 2 ⊥ Q 2 and/or α q , β q parameters (hereafter we assume the small-x approximation α q , β q 1 in all calculations). In this formula and similarly for 1 β±i . For brevity, in what follows we denote ψ A . Let us estimate the relative size of corrections Ξ in Eq. (4.12) at small x. As we will see, 1 α and 1 β transform to 1 αq and 1 βq in our TMDs so if α q ∼ β q ∼ Q √ s (recall that we assume that the DY pair is emitted in the central region of rapidity). For example, the correction ∼ [ψ A γ µ Ξ 2 ][ψ B γ ν Ξ 1 ] will be of order of 7 5 Hadronic tensor at s Q 2 q 2

⊥
In general, our method is applicable for calculation of power corrections at any s, Q 2 q 2 ⊥ , m 2 N . However, the expressions are greatly simplified in the physically interesting case s Q 2 q 2 ⊥ which is considered in this paper. As we noted above, we take into account only hadronic tensor due to electromagnetic currents of u, d, s, c quarks and consider these quarks to be massless. It is convenient to define coordinate-space hadronic tensor multiplied by N c 2 s (and denoted by extra "check" mark) as followsW The reader may wonder why there are no corrections ∼ For future use, let us also define the hadronic tensor in mixed representation: in momentum longitudinal space but in transverse coordinate space After integration over central fields in the tree approximation we obtaiň and similarly for J µ B and J µ BA . Here A, B|O(ψ A , A µ , ψ B , B µ )|A, B denotes double functional integral over A and B fields which gives matrix elements between projectile and target states of Eq. (3.8) type.
The leading-twist contribution to W µν (q) comes only from product J µ , while power corrections may come also from other terms like J µ A (x)J ν B (0). We will consider all terms in turn.

Leading-twist contribution and
where quark fields are given by Eq. (4.12). As we mentioned above, in Ref. [28] it is demonstrated that terms neglected in the r.h.s. lead to power corrections ∼ which are much smaller than if DY pair is emitted in the central region of rapidity). Note that since we want to calculate the leading power corrections, we can substitute Q 2 with Q 2 . In the limit s Q 2 q 2 ⊥ this change of variables can only lead to errors of the order of subleading power terms. 8 As to terms ∼Ψ 1 (x)γ µ Ψ 2 (x)Ψ 2 (0)γ ν Ψ 1 (0), they can be decomposed using eq. (4.12) as follows: 8 Except for the leading-twist term where the difference between Q 2 and Q 2 matters.
where the square brackets mean trace over Lorentz and color indices. First, let us consider the leading-twist term coming from the first term in the r.h.s. of this equation.

Leading-twist contribution
As we mentioned, the leading-twist term comes from J µ where all parentheses in the r.h.s. are color singlet. As usual, after integration over background fields A and B we promote A, ψ A and B, ψ B to operatorsÂ,ψ. A subtle point is that our operators are not under T-product ordering so one should be careful while changing the order of operators in formulas like Fierz transformation. Fortunately, all operators in the r.h.s of Eq. (5.7) are separated either by space-like intervals or light-like intervals so they commute with each other. From parametrization of two-quark operators in section 11.2, it is clear that the leadingtwist contribution to W µν (q) comes from 9 . The corresponding leading-twist contribution to to W µν (q) has the form [40] In a general gauge for projectile and target fields these matrix elements read and similarly for other operators.
Let us discuss other terms proportional to different TMDs in parametrizations in Sect. 11.2. To this end, we write down terms from Eq. (2.3) that we are looking for in Sudakov variables: Here zero in the third term means that the contribution of order one is actually absent. As discussed in Sect. 11.2, all TMDs considered here can have only logarithmic dependence on Bjorken x (≡ α q or β q ) but not the power dependence 1 x . It is easy to see that other quark-antiquark TMDs give contributions to W µ (q) which look like terms in Eq. (5.11) but without extra 1 αq and/or 1 βq so they are power suppressed in low-x regime s Q 2 . Let us also specify the terms which we do not calculate. Roughly speaking, they correspond to terms in Eq. (5.11) multiplied by m 2 ⊥ Q 2 or by either α q or β q . Our strategy in the next sections is to compare a certain term inW µν to terms in Eq. (5.11), and, if it is smaller, neglect, if it is of the same size, calculate.

Terms coming from
We separate terms in Eq. (5.6) according to number of gluon fields (contained in Ξ's ).
where leading-twist terms without gluons (quark-antiquark TMDs) were considered in previous Section, anď The corresponding contributions to W µν (q) will be denoted W µν , and W (2)c µν , respectively. We will consider these contributions in turn.

Term with Ξ 1
Let us start with the last term in Eq. (6.2). The Fierz transformation (11.1) yields . Note that all colors are in the fundamental representation so e.g.
Promoting A and B fields to operators and sorting out the color-singlet contributions we get 10 It is convenient to treat terms ∼ g µν separately so we defineW Hereafter we omit "hat" notation from from operators: O A,B ≡ Ô A,B for brevity. 10 We will keep different notations Ai and Bi for the projectile and target fields because of the relations (11.10) and (11.14) Let us now estimate this contribution toW µν . First, recall that B i is of order of m ⊥ (more accurately, it will be ∼ q i after the Fourier transformation, see e.g. Eq. (11.42) or Eq. (11.48)). Next, as demonstrated in Sect. 11.3 (see Eqs. (11.30), (11.31)), 1 α in the target matrix element turns to ± 1 αq after Fourier transformation. Due to this fact we will replace 1 α by 1 αq in our estimates, even in the coordinate space. Similarly, for the estimate of the target matrix elements we will replaces operator 1 β by 1 βq wherever appropriate. Now we will demonstrate that three terms in the r.h.s. of Eq. (6.8) are small in comparison to terms listed in Eq. (5.11). The projectile matrix element in the first term in the r.h.s. of Eq. (6.8) brings factor s (see Eq. (11.29)) but the target matrix element can produce only factor x i so the first term is ∼ gµν m 2 ⊥ αqs which is smaller than gµν q 2 ⊥ Q 2 that we have in Eq. (5.11) (and will calculate in the next Section). As to the second term in the r.h.s. of Eq. (6.7), it can be rewritten as The projectile matrix element in the first term in the r.h.s. of this equation brings factor s but, as we discussed above, the target matrix element cannot produce factor s so this As to the second term, converting three γ-matrices in the projectile matrix element to a combination of γ's and γγ 5 's and looking at the parametrization of Sect. 11.2, we see that 2 In addition, as discuss in Sect. 11.2, the target matrix element ψ B i (0)γ µ ψ(x) B knows about p 1 only via the direction of Wilson lines so it can be proportional only to p 1µ p 1 ·p 2 that does not change at rescaling of p 1 . Thus, ψ B i (0) p 2 ψ(x) B is ∼ O(1) and therefore the second term in Eq. (6.9) is even smaller than the first one. Finally, let us discuss the third term in the r.h.s. of Eq. (6.8). If both α and β are transverse similarly to the first term in Eq. (6.9). If both indices are longitudinal, we get The projectile matrix element brings a factor s, but the target one is ∼ O(1) due to the reason discussed above, so this contribution is negligible. Finally, let us consider the case when index α is longitudinal and β is transverse Again, the target matrix element is ∼ O(1) while the projectile one can bring one factor of s as can be seen from parametrization (11.29) by reducing the number of γ-matrices to two. Thus, the contribution (6.12) is negligible and so is the total contribution (6.8).
We geť Let us start with the case when both of the indices µ and ν are transverse. It is easy to see that the power counting for the first term in the r.h.s. of Eq. (6.13) is the same as for Eq. (6.9) so it is small. Also, the estimate of the second term in Eq. (6.13) is similar either to the estimate of Eq. (6.10) or (6.12) so it can be neglected.
Next, let us consider both µ and ν longitudinal. It is easy to see that multiplication of the r.h.s. of Eq. (6.13) by p µ 2 p ν 2 gives zero so there is no term proportional to p µ 1 p ν 1 . The term proportional to p µ 2 p ν 2 has the form It is easy to see that both projectile and target matrix elements are proportional to the first power of s so the resulting estimate is in comparison to the corrseponding term in Eq. (5.11). If one index is p 1 and the other p 2 we get It is easy to see that in all terms the projectile matrix element is ∼ s but the target one is Finally, let us consider the case when one of the indices in Eq. (6.13) is longitudinal and one transverse. For example, let µ be longitudinal and ν transverse, the opposite case will differ by replacement µ ↔ ν. Using the decomposition of g µν in longitudinal and transverse part (3.2) we get The term proportional to p 2µ in the r.h.s. can be expressed using Eq. (11.13) as follows Let us at evaluate two the most important contributions. The first is As we shall see below, due to QCD equations of motion B in the r.h.s. of this equation can be replaced by transverse momentum of the target TMD k ⊥ . Also, 1 α will be replaced by 1 αq so from the parametrizations (11.24) and (11.27) we see that which is of order of fourth term in Eq. (5.11). The second relevant term is where we used formula (11.4) and the fact that for unpolarized protons from parity conservation. 11 Again, 1 α will turn to 1 αq and B can be replaced by k ⊥ for the target, so (6.11) is of order of 11 A rigorous argument goes like that: the matrix element (6.21) can be rewritten as . As demonstrated in Sect. 11.3, A in this formula can be replaced by k ⊥ so the contribution is proportional to matrix element k i ψ (0)iσ•iγ5ψ(x) = k i ij ψ (0)σ•jψ(x) which vanishes as seen from the parametrization (11.29).
Let us demonstrate that the remaining terms in the r.h.s. of Eq. (6.17) are negligible. First, term coming from replacement ψ(0) ⊗ ψ(x) ↔ γ 5 ψ(0) ⊗ γ 5 ψ(x) in Eq. (6.18) vanishes since ψ (x) p 2 γ 5 ψ(0) A = 0 for unpolarized hadrons, see Eq. (11.28). Next, term p 2µ is small because neither projectile no target matrix elements can bring factor s. Last, using Eq. (11.4) we get It is easy to see that neither the projectile nor the target matrix element in the r.h.s. of this equation gives s so these terms can be neglected in comparison to Eqs. (6.19) and (6.22). Thus, the two non-negligible terms in Eq. (6.17) givě Using formulas (11.30), (11.31), (11.33), (11.36), (11.41), and (11.43) for quark-antiquarkgluon operators and parametrizations from Sect. 11.2 we get the contribution to W µν in the form Next we consider the remaining ∼ p 1µ term in Eq. (6.13) which can be rewritten as where again we used formula (11.4). Note that while the matrix elements between projectile states give contributions ∼ s αq k ⊥ , the target matrix elements cannot give s. Indeed, these target matrix elements know aboutp 1 only through direction of Wilson lines so they should not change under rescaling p 1 → λp 1 , see the discussion in Sect. 11.2. Thus, the r.h.s. of Eq.
Thus, the contribution of the first term in Eq. (6.2) to W µν is Let us consider now the second term in Eq. (6.2). The calculation repeats that of the first term so we will indicate here main steps and pay attention to non-negligible terms only. If one of the indices (say, µ) is longitudinal and the other transverse, we get The most important terms are those proportional to p 2µ . Using Fierz transformation and separating color singlets, they can be rewritten as (cf. Eq. (6.17)) After some algebra with γ-matrices this can be transformed to plus terms small in comparison to p 2µ q ⊥ ν αqs . Using Eq. (11.31) we can transform ψ 1 Using QCD equation of motion and other formulas from sections 11.2 and 11.3 one gets so the contribution of Eq. (6.28) is effectively doubled. Again, the term with x ↔ 0 exchange leads to Eq. (6.34) with f 1 ↔f 1 and h ⊥ 1 ↔h ⊥ 1 replacement. Thus, the sum of first and second terms in Eq. (6.2) leads to twice Eq. (6.28) In this Section we calculate the third term in Eq. (6.2).
Again, main contribution correspond to one index (e.g. µ) being longitudinal and the other transverse so we need . Let us consider first the term proportional to p 1µ . Performing Fierz transformation (11.1) and sorting out the color-singlet contributions we get (cf. Eq. (11.31))   while contribution from the last line is even smaller, of order of It remains to prove that the last term in Eq. (6.37) proportional to p 2µ is small. One can rewrite that term similarly to Eq. (6.26) with replacement p 1 ↔ p 2 and (projectile matrix elements) ↔ (target ones). After that, the proof repeats arguments after Eq. (6.26) and one obtains the estimate Similarly, by repeating arguments from Section 6.1.1 with replacement p 1 ↔ p 2 and projectile matrix elements ↔ target ones, one can demonstrate that terms in Eq. (6.36) with µ, ν both longitudinal or both transverse are small in comparison to terms listed in Eq. (5.11). Thus, Using QCD equation of motion and formulas from Appendix, we obtain the corresponding contribution to W µν in the form Same as in previous Section, the term with x ↔ 0 exchange leads to Eq. (6.43) with f 1 ↔f 1 and h ⊥ 1 ↔h ⊥ 1 replacement so we get Repeating arguments from previous Section, it is possible to show that the contribution of the fourth term doubles that of the third term so we get the full contribution of the terms with one quarkantiquark-gluon operators in the form This result agrees with the corresponding 1/Q terms in Ref. [20].
6.2 Term with two quark-quark-gluon operators coming from Ξ 1 and Ξ 2 Let us start with the first term in the r.h.s. of Eq. (6.3). Performing Fierz transformation (11.1) we obtain It is convenient to defineV 3µν to be traceless. In next Sections, we will consider these terms in turn.

Term propotional to g µν
Using Eq. (4.12) and extracting color-singlet contributions one obtainš Let us start with the first term. Using Eq. (11.21) and the fact that ψ (x) A k σ * j − A j (x)σ * k ψ(0) A = 0 (see the footnote 11), we obtain where we used the fact that projectile and target matrix elements in the two last terms in the l.h.s. cannot produce factor of s. Next, consider second term in Eq. (6.51). Using Eqs. (11.17) and (11.21), one can rewrite is as 1 Similarly, from Eq. (11.21) we get the third term in the form Since both projectile and target matrix elements cannot give factor s this contribution is O q 2 ⊥ s in comparison to that of the two first terms. Thus, we geť Next, using QCD equations of motion (11.40), (11.43) and formulas from Appendix 11.2, we obtain the contribution to W µν in the form where replacements f f 1 ↔f f 1 and h ⊥ 1f ↔h ⊥ 1f come from x ↔ 0 term.

Term with TMD's f 1
Separating color-singlet contributions one can rewrite Eq. (6.49) aš We need to consider three cases: both µ and ν are transverse, both of them are longitudinal, and µ is longitudinal and ν transverse (plus vice versa).
In the first case we can use formula (11.18) and geť which gives the contribution to W µν in the form where we again used formulas from Appendices 11.2 and 11.3. Next, if both µ and ν are longitudinal, we geť Using formula (11.20) we rewrite r.h.s. of Eq. (6.57) as followš Since matrix elements in the r.h.s. cannot give factor s, the contribution of this term to W µν is ∼ q 2 ⊥ s times that of Eq. (6.59). Finally, let us consider the case when one index is longitudinal and the other transverse. Using Eq. (11.23) we geť It is clear that ψ A(x) p 2 γ i 1 α ψ(0) A and ψ B(0) p 1 γ i 1 β ψ(x) B bring one factor s sǒ which is 1 αqs or 1 βqs correction in comparison to Eq. (6.46). Thus, the contribution to W µν is given by Eq. (6.59)

1
Let us consider now (the trace will be subtracted after the calculation). Separating color-singlet contributions, we geť First case is when µ and ν are transversě With the help of Eq. (11.4) the first term in the r.h.s. turns to After some algebra, it can be rewritten as where again we used property (6.21). Using QCD equations of motion (11.41), (11.43) and parametrization (11.47) one can write the corresponding contribution to W µν as where we introduced the notation The second term in Eq. (6.67) can be rewritten aš where we used Eq. (11.4). It is clear that neither projectile no target matrix element in the r.h.s. can bring factor s soV in comparison to Eq. (6.70). Next, consider the case when both µ and ν are longitudinal. The non-vanishing terms areV The first two terms in the r.h.s. can be rewritten aš The corresponding contribution to W µν yields where again we used QCD equations of motion (11.41), (11.43) and parametrization (11.47). Next, it is easy to see that the third term in Eq. (6.74) is small in comparison to Eq. (6.76): because neither projectile no target matrix element can bring factor s. Finally, take one of the indices (say, µ) longitudinal and the other transverse. From Eq. (6.66) we geť Using formulas (11.4) this can be rewritten as followš As we discussed above, projectile matrix elements in the r.h.s. like ψ B j (0)σ •i 1 β ψ(x) B can bring factor s but the target matrix elements cannot so the corresponding contribution to W µν is of orderV in comparison to Eq. (6.46).
Next, the sum of Eqs. (6.70) and (6.76) is so subtracting trace we obtain As we will see in Sect. 8, cancellation of terms ∼ g µν proportional to h A in the r.h.s of this equation is actually a consequence of (EM) gauge invariance.
Let us now assemble the contribution of terms (6.47) to W µν . Summing Eqs. (6.56), (6.64), and (6.82) we get µν (q) of Eq. (6.3) we need to add the contribution of the term . Similarly to the case of one quark-quark-gluon operator considered in Sect. 6.1, it can be demonstrated that this contribution doubles the result (6.83) so we get where Q 2 ≡ α q β q s 6.3 Term with two quark-quark-gluon operators coming fromΞ 2 and Ξ 2 Let us start with the first term in Eq. (6.4).
After Fierz transformation (11.1) we obtaiň We will consider them in turn.

Term proportional to f 1f1
Let us start with g µν term in Eq. (6.88).
It is obvious that the target matrix element can bring factor s. On the contrary, as we discussed above, the projectile matrix element cannot produce s since Indeed, since projectile matrix elements know about p 2 only through the direction of Wilson lines, the l.h.s. can be proportional only to factor p 2α p 1 ·p 2 that does not change under rescaling of p 2 . Also, due to Eq. (11.31) ψ 1 β (0) ⊗ 1 β ψ(x) B can be replaced by − 1 Consequently, the r.h.s. of Eq. (6.90) is ∼ g µν in comparison to Eq. (6.84). We getV If the index ν is transverse, the contribution of this equation to W µν is of order of βqs in comparison to Eq. (6.46). For the longitudinal indices µ and ν we geť Similarly to Eq. (6.92), the contribution of the second term to W µν is The corresponding contribution to W µν is obtained from QCD equation of motion (11.44) and formula (11.31) from Appendix 11.3: Let us start with g µν term in Eq. (6.89).
The target matrix element is proportional to s while the projectile one cannot bring s due to Eq. (6.91), so the contribution of the r.h.s of Eq. (6.99) to W µν is of order h.s. of Eq. (6.84) (6.100) similarly to Eq. (6.96). We geť Let us at first consider the second term in this formula: Similarly to Eq. (6.91), projectile matrix elements cannot give factor s so the corresponding contribution to W µν is of order of in comparison to Eqs. (6.46) and (6.84), respectively. We are left witȟ First, note that the two last terms are small, of order of Eq. (6.103), for the same reason as Eq. (6.101) above. As to the first term in r.h.s. of Eq. (6.104), using Eq. (11.7) it can be rewritten aš so the corresponding contribution to W µν takes the form where we used Eqs. (11.31) and (11.45). The full result for W (2b) µν is given by the sum of Eqs. (6.98) and (6.106)

Second term in Eq. (6.4)
Let us start now consider the second term in Eq. (6.4).
After Fierz transformation (11.1) we obtaiň Sorting out color-singlet terms, we get similarly to sum of Eqs. (6.88) and (6.89) Starting from this point, all calculations repeat those of Sections 6.3.1 and 6.3.2 with replacements of p 1 ↔ p 2 , α q ↔ β q and exchange of projectile matrix elements and the target ones. The result is Eq. (6.111) with these replacements so we finally get

Third term with two quark-quark-gluon operators
Let us consider the first term in the r.h.s. of Eq. (6.5). After Fierz transformation it turns to After separation of color singlet contributions Since projectile and target matrix elements can bring s each and 1 α and 1 β convert to 1 αq and 1 βq , the typical contribution of (6.115) to W µν is In Ref. [28] we calculated the sum of these structures corresponding to convolution of µ and ν. In principle, one can repeat that calculation and find contribution to these structures separately. However, since the corresponding matrix elements of quark-quarkgluon operators are virtually unknown, in this paper we we will disregard such 1 N 2 c terms. Thus, the contribution of Eq. (5.6) to W µν (q) is given in the leading order in N c by the sum of equations (6.46), (6.84), and (6.111).

Power corrections from
Power corrections of the second type come from the terms where Ψ 1 and Ψ 2 are given by Eq. (4.12). 12 We get First, let us demonstrate that contributions to W µν from the second to fifth lines in eq. (7.2) vanish. Obviously, matrix element of the operator in the second line vanishes. Formally, and, non-formally, one hadron cannot produce DY pair on his own.
It is easy to see that contributions toW µν from the third and the fourth lines in Eq. (7.2) vanish due to the absence of color-singlet structure. Indeed, let us consider for example the term which obviously does not have color-singlet contribution. Similarly, other three terms in the third and fourth lines in Eq. (7.2) vanish.
Next, let us demonstrate that the contribution of the fifth line in Eq. (7.2) vanishes for the same reason as in Eq. (7.3). Let is consider for example the first term in the fifth line In the appendix 8.3.2 to [28] it is demonstrated that higher-order terms in the expansion Eq. (4.11) (denoted by dots) are small in our kinematical region s where we separated color-singlet contribution in the last line. The corresponding term in W µν is Similarly, the contribution of the second term in the fifth line of Eq. (7.2) will be proportional to δ(β q ) and hence vanish.
Let us now discuss the non-vanishing contributions coming from last two lines in Eq. (7.2). For example, the first term in the sixth line is Separating color-singlet contributions with the help of the formula we get the corresponding term inW µν in the form which is similar to Eq. (6.57) with exception of extra color factor Nc Nc . Consequently, as discussed in Sect. 6.2.2, non-negligible contributions come from transverse µ and ν only. We calculate them in next Section.

Last two lines in eq. (7.2)
In this section we calculate the traceless part of sixth and seventh lines Eq. (7.2). Since we consider only transverse µ and ν, to simplify notations we will call them m and n in this Section. Using eq. (4.12) and separating color-singlet matrix elements with the help of Eq. (7.9), we rewrite the traceless part of sixth and seventh lines in Eq. (7.2) as To save space, hereafter we do not display subtraction of trace with respect to m, n indices but it is always assumed. Using formulas (11.22) we can write down the contribution toW µν from sixth and seventh lines in Eq. (7.2) in the form Let us now consider corresponding matrix elements. It is easy to see that where we used parametrization (11.49). For the target matrix elements, we obtain 14) The corresponding contribution to (traceless) W (α q , β q , x ⊥ ) takes the form where we have recovered the subtraction of trace.
The trace part can be obtained in a similar way. Using Eq. (11.21) one gets The corresponding contribution to trace part of W (α q , β q , x ⊥ ) takes the form which agrees with Eq. (6.2) from Ref. [28] after replacements j 2 = j tw3 2 − ij tw3 2 and j 2 = j tw3 1 + ij tw3 1 . It should be noted that the difference between j 1 and j 2 in traceless vs trace part is due to difference in formulas (11.22) and (11.21).
Thus, the result is the sum of Eqs. (7.15) and (7.17) As we mentioned in the Introduction, in this paper we we will take into account only leading and sub-leading terms in N c and leave the 1 N 2 c corrections discussed above for future publications.
Finally, as proved in Appendix 11.5, we can neglect contributions proportional to the product of quark and gluon TMDs.

Results
Assembling Eqs. (5.10), (6.46), (6.84), (6.111), and (7.18) we get the result for W µν (q) that consists of two parts: W µν (q) = W 1 µν (q) + W 2 µν (q) (8.1) The first, gauge-invariant, part is given by where F f and H f are given by Eq. (6.29) and where q µ ≡ α q p 1 + β q p 2 andq µ ≡ α q p 1 − β q p 2 . These are the same expressions as in Eq. ⊥ Q 4 corrections due to difference between Q 2 and Q 2 . It is easy to see that q µ W F µν = 0 and q µ W H µν = 0. Note that q µ W F µν and q µ W H µν are exactly zero without any q 2 ⊥ Q 2 corrections. This is similar to usual "forward" DIS, but different from off-forward DVCS where the cancellations of right-hand sides of Ward identities involve infinite towers of twists [41][42][43] The second part is where H A is given by Eq. (6.71) and These terms are not gauge invariant: q µ W 2 µν (q) = 0. The reason is that gauge invariance is restored after adding terms like m 2 ⊥ Q 2 × Eq.(5.11) which we do not calculate in this paper. Indeed, for example, They are of the same order so one should expect that gauge invariance is restored after calculation of the terms ∼ p µ 2 q ν ⊥ q 2 ⊥ α 2 q βqs 2 which are beyond the scope of this paper. For the same reason we see that all structures in Eq. (5.11) except g ⊥ µν q 2 ⊥ αqβqs and q ⊥ µ q ⊥ ν αqβqs are determined by leading-twist TMDs f 1 and h ⊥ 1 . Sometimes it is convenient to represent hadronic tensor in transverse coordinate space. Introducing (and similarly for target TMDs) we get (the question about rapidity cutoffs for TMDs will be addressed in Sect. 9). Similarly, we can write down W H contribution in coordinate space. For future use, however, it is convenient to define Fourier transform in a slightly different way. Introduce then W 1H µν can be represented as everywhere except h ↔h terms where it is opposite, cf. Eq. (8.9).

Four Lorentz structures of hadronic tensor
The four Lorentz structures of hadronic tensor in Collins-Soper frame are given by Eq.
First, let us check the structure corresponding to the total cross section of DY pair production. From Eq. (8.1) we get which agrees with Eq. (6.2) from Ref. [28]. This equation gives the sum of structures W µ µ = − (2W T + W L ).

W L
The easiest structure to get is W L . Multiplying Eq. (8.1) by Z µ Z ν and comparing to Eq.
(1.1) we get Thus, one may say that W L is known at LHC energies at q 2 ⊥ Q 2 as far as f 1 and h ⊥ 1 are known.
Again, we see that W ∆ is expressed via f 1 and h ⊥ 1 with great accuracy.

W ∆∆
Finally, the easiest way to pick out W ∆∆ is to multiply W µν by q ⊥ µ q ⊥ ν /q 2 ⊥ . One obtains from Eq. (1.1) and from Eq. (8.5) Thus, we get This is the only function which has a O 1 Q 2 , leading-N c contribution proportional to twistthree TMD H A not related to leading-twist TMDs by equations of motion. The functions W T , W L , and W ∆ do not have such contributions (although they have such contributions at the 1 Nc level).

Order-of-magnitude estimates
Following the analysis in Ref. [28], let us estimate the relative strength of Lorentz structures W i at q 2 ⊥ m 2 . First, we assume that 1 Nc is a good parameter and leave only terms leading in N c . Second, at q 2 ⊥ m 2 we probe the perturbative tails of TMD's f 1 ∼ 1 44]. So, as long as Q 2 q 2 ⊥ m 2 we can approximate (up to logarithmic corrections). Similarly, for the target we can use the estimate as long as k 2 ⊥ Q 2 . Thus, we get an estimate Note that due to the "positivity constraint" [45] h we can safely assume that the functions f (x) and h(x) defined in Eqs. (8.20) and (8.21) are of the same order of magnitude. Moreover, both theoretical [46] and phenomenological [47,48] analysis indicate that h ⊥ 1 is several times smaller than f 1 so in numerical estimates we will disregard the contribution of h ⊥ 1 .

Power corrections for total DY cross section
Substituting the above approximations to Eq. (8.13) we get the following estimate of the strength of power corrections for total DY cross section [28] where we used estimates (8.22) and the fact that (k, q − k) ⊥ ∼ q 2 ⊥ m 2 . Thus, the relative weight of the leading term and power correction is determined by the factor 1 − 2 (k,q−k) ⊥ Q 2 . Due to Eqs. (8.20) and (8.21), the integrals over k ⊥ are logarithmic and should be cut from below by m 2 N and from above by Q 2 so we get an estimate where we assumed that the first integral is determined by the logarithmical region q 2 ⊥ k 2 ⊥ m 2 N and the second by Q 2 k 2 ⊥ q 2 ⊥ . Taking these integrals to Eq. (8.24) one obtains By this estimate, the power correction reaches the level of few percent at q ⊥ ∼ Q 4 .

Power corrections for W T
Let us now consider estimates described in Sect. 8.3.1 for W T given by Eq. (8.16). At large N c , we can omit the third line so Again, due to q 2 ⊥ m 2 the second term in braces can be neglected and we get Thus, for W T the power correction reaches 10% level at q ⊥ ∼ Q 2 .

Estimate of W L
Again, using estimates from Sect. 8.3.1 one obtains which gives approximately in agreement with Eqs. (8.26) and (8.29). The estimate of the ratio of W L /W T is

Magnitude of W ∆
It is easy to see that W ∆ vanishes if one uses the estimates (8.20) and (8.21). Indeed, with these formulas F (q, k ⊥ ) and H(q, k ⊥ ) are symmetric under replacement k ⊥ ↔ (q − k) ⊥ whereas (q, q − 2k) ⊥ in the integrand in Eq. (8.15) is antisymmetric. Moreover, this vanishing of W ∆ will occur for any factorizable model of TMDs (8.15) vanishes. On the other hand, W ∆ is only ∼ Q ⊥ Q W T so without better knowledge of TMDs it is impossible to tell whether W ∆ is smaller or bigger than, say, W L . Also, if the parameter α q Q Q ⊥ is not negligible, to compare W ∆ and W L one needs to take into account O(α q ) corrections to W ∆ defined by TMDs other than f 1 and h ⊥ 1 .

Estimate of W ∆∆
Let us consider the relative weight of the terms in the r.h.s. of Eq. (8.19). As we mentioned, we assume that 1 Nc is a valid small parameter so we can omit the last J 1 term. Also, it is natural to assume that H f A (q, k ⊥ ) is of the same order of magnitude as H f (q, k ⊥ ) and, since the term with H A is a power correction, it is not unreasonable to neglect this term in the first approximation. Using now estimates (8.22) and the integrals one gets an estimate

Lam-Tung relation
It is easy to see that if one neglects H in Eq. (8.14) the ratio of W L and 2W ∆∆ is approximately It seems like the Lam-Tung relation works better if we move closer to the domain of collinear factorization Q 2 ∼ Q 2 ⊥ m 2 .

Estimates of asymmetries
The differential cross section of DY process is parametrized as 1 + λ cos 2 θ + µ sin 2θ cos φ + ν 2 sin 2 θ cos 2φ (8.36) where Ω is the solid angle of the lepton in terms of its polar and azimuthal angles in the center-of-mass system of the lepton pair. The angular coefficients λ, µ, and ν can be expressed in terms of the hadronic tensor: For an estimate, let us take s=8 TeV and Q=90 GeV so that x A ∼ x B ∼ 0.1 in central region of rapidity. Although we did not include the contribution of Z-boson, we can compare our order-of-magnitude estimates with experimental data at this kinematics [49,50]. Let us take Q ⊥ = 20 GeV so the power corrections ∼ Q 2 ⊥ Q 2 are small but sizable, of order of few per cent. At this kinematics, we obtain As to µ coefficient, as we mentioned, we cannot estimate it since with factorization hypothesis for TMDs it vanishes. Reversing the argument, if µ will be measured to be much smaller than ν, it will be an argument in favor of factorization hypothesis for TMDs f 1 and h ⊥ 1 . Actually, there are experiments at much lower q ⊥ and Q ∼ few GeV which indicate that µ is very small [52].
Last but not least, let us estimate Lam-Tung relation. With our approximation in the above kinematics we get so it seems to be violated at this kinematics. Again, these order-of-magnitude estimates do not include the contribution to DY cross section mediated by the Z-boson.

Matching for W 1F terms
We need to multiply Eq. (8.9) in coordinate space by M (α q , β q , b ⊥ ; ς p , ς t ). First, recall that does not actually depend on the "rapidity divides" ς p and ς t . However, the differentiation ∂ ∂b i affects evolution equations (9.2). In this case we modify the derivative with respect to b i as follows∂ , It is easy to check gauge invariance: (α q p µ 1 + β q p µ 2 + i∂ µ ⊥ )W 1F µν (α q , β q , b ⊥ ) = 0.

Matching for W 1H terms
First, with our definitions (8.10) the Eq. (11.29) reads and similarly for the target matrix elements. With such definition, the evolution equation and similarly forh i . The reason is that one-loop rapidity evolution forψ(x • , x ⊥ )Γψ(0) in the Sudakov region is the same for all matrices Γ betweenψ(x • , x ⊥ ) andψ ( 0) due to the fact that the "handbag" diagram in Fig. 3c is small and in two other diagrams (as well as self-energy corrections) the matrix Γ betweenψ(x • , x ⊥ ) andψ(0) just multiplies the result of calculation.
x 0 x 0 x 0 Figure 3. Typical diagrams for the rapidity evolution of quark TMD in the Sudakov regime.
Next, to write down the product of m 2 W H µν and the coefficient function we need modified derivatives of h i 's of Eq. (9.6) type: (and similarly forh's) so that∂ i h j will satisfy same evolution equations (9.9) as h j . We get then except h ↔h terms where it is opposite, cf. Eq. (8.9). Let us comment on the choice of "rapidity divides" ς p and ς t in the product M (α q , β q , b ⊥ ; ς p , ς t )f 1 (α q , b ⊥ ; ς p )f 1 (β q , b ⊥ ; ς t ) (and in similar M hh product). As we mentioned in the beginning of this Section, in order to avoid double counting one should write down factorization of the amplitude in projectile, target and central fields at ς p , ς t ∼ 1.
After that, as discussed in Ref. [53], one can use the double-log Sudakov evolution (9.2) until ς p ≥ς p ,ς p ≡ 1 At this point, the result of Sudakov evolution is so in the final result (8.9) one should take at the end of Sudakov evolution (9.9). It should be emphasized that since factor M is universal for (9.14) and (9.15), our estimates of asymmetries in Sect. 8.3.8 are not affected by summation of Sudakov double logs.

Rapidity-only cutoff for TMD
As discussed in Refs. [28,29], from the rapidity factorization (3.8) we get TMDs with rapidity-only cutoff |α| < σ t or |β| < σ p (or with modifications (9.12)). Such cutoff, relevant for small-x physics, is different from the combination of UV and rapidity cutoffs for TMDs used by moderate-x community, see the analysis in two [55][56][57] and three [58] loops. For the tree-level formulas of Sect. 8, this difference in cutoffs does not matter, but if one uses the formulas from Sect. 9 and integrates models for TMDs with Sudakov factor M of Eq. (9.4), one has to relate TMDs with rapidity-only cutoffs to the TMD models with conventional cutoffs. This requires calculations at the NLO level which are in progress.

Conclusions and outlook
Main result of this paper is Eq. (8.1) which gives the DY hadronic tensor for electromagnetic current at small x with gauge invariance at the 1 Q 2 level. The part (8.2), determined by leading-twist TMDs f 1 and h ⊥ 1 , is manifestly gauge invariant. The only non-gauge invariant term at the 1 Q 2 level is Eq. (8.6) with transverse µ and ν which is ∼ q ⊥ µ q ⊥ ν Q 2 times twist-3 TMDs. Also, in the leading-N c approximation the only structure affected by those terms is W ∆∆ , all other structures are calculated up to O q 4 ⊥ Q 4 terms. It is interesting to note that 1 Q 2 terms necessary for gauge invariance are calculated more than than two decades after the calculation of 1 Q corrections in Ref. [20]. It should be mentioned that, as discussed above, our rapidity factorization is different from the standard factorization scheme for particle production in hadron-hadron scattering, namely splitting the diagrams in collinear to projectile part, collinear to target part, hard factor, and soft factor [9]. Here we factorize only in rapidity and the Q 2 evolution arises from k 2 ⊥ dependence of the rapidity evolution kernels, same as in the BK (and NLO BK [59]) equations. Also, since matrix elements of TMD operators with our rapidity cutoffs are UV-finite [60,61], the only UV divergencies in our approach are usual UV divergencies absorbed in the QCD running coupling. For the tree-level result (8.1) this does not matter, but if one intends to use the result like (8.9) with Sudakov logarithms for conventional TMDs with double UV and rapidity cutoffs, one needs to relate our TMDs with rapidityonly cutoff to conventional TMDs. Needless to say, the gauge-invariant tree-level result (8.2) should be correct for TMDs with any cutoffs.
An obvious outlook is to extend these results to the "real" DY process involving Z-boson contributions which are relevant for our kinematics. The study is in progress.
Note that the coefficients in front of f 3 , g ⊥ f , h and h ⊥ 3 in eqs. (11.24), (11.26), (11.28), and (11.29) contain an extra 1 s since p µ 2 enters only through the direction of gauge link so 17 In an arbitrary gauge, there are gauge links to −∞ as displayed in eq. (5.9).
the result should not depend on rescaling p 2 → λp 2 . For this reason, these functions do not contribute to W (q) in our approximation. Last but not least, an important point in our analysis is that any f (x, k ⊥ ) may have only logarithmic dependence on Bjorken x but not the power dependence ∼ 1 x . Indeed, the low-x behavior of TMDs is determined by pomeron exchange with the nucleon. The interaction of TMD with BFKL pomeron is specified by so-called impact factor and it is easy to check that the impact factors for all leading-twist TMDs are similar and do not give extra 1 x factors. The only 1 x may had come from some unfortunate definition of TMD which includes factor s artificially, but from power counting (5.11) we see that all definitions of leading-twist TMDs do not have such factors.

Matrix elements of quark-quark-gluon operators
In this section we will demonstrate that matrix elements of quark-antiquark-gluon operators from section 6 can be expressed in terms of leading-power matrix elements from section 11.2.
First, let us note that operators 1 α and 1 β in Eqs. (4.13) are replaced by ± 1 αq and ± 1 βq in forward matrix elements. Indeed, and Γ can be any γ-matrix. Similarly, The corresponding formulas for target matrix elements are obtained by substitution α ↔ β (and x • ↔ x * ). Next, we will use QCD equation of motion to reduce quark-quark-gluon TMDs to leading-twist TMDs (see Ref. [20]). Let us start with matrix element For corresponding antiquark distributions we get Also, as we saw in Sect. 6.2.3, at the leading order in N c there is one quark-antiquarkgluon operator that does not reduce to twist-2 distributions. It can be parametrized as follows (cf. Eq. (11.43)) (11.47) and similarly for the target matrix elements.

Parametrization of TMDs from section 7.1
We parametrize TMDs from section 7.1 as follows To estimate the magnitude of this contribution, first note that dx * e −iβqx * B|A a i (x)A a j (0)|B (11.53) where we used parametrization (3.26) from Ref. [29]. Since the gluon TMDs D g (x B , x ⊥ ) and H g (β q , x ⊥ ) behave only logarithmically as x B → 0 [61], the contribution of eq. (11.52) to W (q) is of order of m 2 ⊥ βqs m 2 ⊥ Q 2 . (As discussed in Ref. [29], the projectile TMD in the r.h.s. of eq. (11.52) does not give 1 αq after Fourier transformation). Also, this contribution is ∼ 1 Nc with respect to our leading terms. Similarly, all other terms in eq. (11.51) are either m 2 ⊥ βqs or m 2 ⊥ αqs times 1 Nc so they can be neglected. 19