Power corrections to TMD factorization for Z-boson production

A typical factorization formula for production of a particle with a small transverse momentum in hadron-hadron collisions is given by a convolution of two leading-twist TMD parton densities with cross section of production of the final particle by the two partons. For practical applications at a given transverse momentum, though, one should estimate at what momenta the power corrections to TMD factorization formula due to higher-twist operators become essential. In this paper we calculate the first power corrections to TMD factorization formula for Z-boson production and Drell-Yan process in high-energy hadron-hadron collisions. At the leading order in $N_c$ power corrections are expressed in terms of leading-twist TMDs by QCD equations of motion.


Introduction
A typical analysis of differential cross section of particle production in hadron-hadron collisions at small momentum transfer of the produced particle is performed with the help of TMD factorization [1][2][3][4][5][6][7][8][9][10]. However, the question of how small should be the momentum transfer in order for leading power TMD analysis to be successful cannot be resolved at the leading-power level. The sketch of the factorization formula for the differential cross section is [1,11] dσ where η is the rapidity, q is the momentum of the produced particle in the hadron frame (see ref. [1]), D f /A (x, z ⊥ , η) is the TMD density of a parton f in hadron A, and σ(f f → H) is the cross section of production of particle H in the scattering of two partons. The common wisdom is that when we increase transverse momentum q 2 ⊥ of the produced hadron, at first the leading power TMD analysis with (nonperturbative) TMDs applies, then at some point power corrections kick in, and finally at q 2 ⊥ ∼ Q 2 , where Q 2 = q 2 , they are transformed into so-called Y-term making smooth transition to collinear factorization formulas. In this paper we try to answer the question about the first transition, namely at what q 2 ⊥ power corrections become significant.
In our recent paper [12] we calculated power corrections ∼ q 2 ⊥ Q 2 for Higgs boson production by gluon-gluon fusion. The result was a TMD factorization formula with matrix elements of three-gluon operators divided by an extra power of m 2 H . In this paper we calculate power corrections ∼ q 2 ⊥ Q 2 for Z-boson production which are determined by quarkquark-gluon operators. In the leading order Z-boson production was studied in [13][14][15][16][17][18][19][20][21]. The interesting (and unexpected) result of our paper is that at the leading-N c level matrix elements of the relevant quark-quark-gluon operators can be expressed in terms of leading power quark-antiquark TMDs by QCD equations of motion (see ref. [22]). The method of calculation is very similar to that of ref. [12] so we will streamline the discussion of the general approach and pay attention to details specific to quark operators.
The paper is organized as follows. In section 2 we derive the TMD factorization from the double functional integral for the cross section of particle production. In section 3, which is central to our approach, we explain the method of calculation of power corrections based on a solution of classical Yang-Mills equations. In section 4 we find the leading power correction to particle production in the region s Q 2 q 2 ⊥ . In section 5 we perform the order-of-magnitude estimate of power corrections and in section 6 present our result for power corrections to the Drell-Yan cross section. The necessary technical details and discussion of subleading power corrections can be found in appendix.

TMD factorization from functional integral
We consider Z-boson production in the Drell-Yan reaction illustrated in figure 1: h A (p A ) + h B (p A ) → Z(q) + X → l 1 (k 1 ) + l 2 (k 2 ) + X, (2.1) where h A,B denote the colliding hadrons, and l 1,2 the outgoing lepton pair with total momentum q = k 1 + k 2 . The relevant term of the Lagrangian for the fermion fields ψ i describing coupling between fermions and Z-boson is (s W ≡ sin θ W , c W ≡ cos θ W ) where sum goes over different types of fermions, and coupling constants g V i = (t L 3 ) i − 2q i s 2 W and g A i = (t L 3 ) i are defined by week isospin (t L 3 ) i of the fermion i, see ref. [23]. In this paper we take into account only u, d, s, c quarks and e, µ leptons. We consider all fermions to be massless. The differential cross section of Z-boson production with subsequent decay into e + e − (or µ + µ − ) pair is where we defined the "hadronic tensor" W (p A , p B , q) as (2.4) As usual, X denotes the sum over full set of "out" states. It should be mentioned that there is a power correction coming from the leptonic tensor term ∼ q µ q ν . However, if we consider quarks to be massless, the only effect of the q µ q ν term comes from the (square of) axial anomaly which has an extra factor α 2 s , and such two-loop factor is beyond our tree approximation.
The sum over full set of "out" states in Eq. (2.4) can be represented by a double functional integral 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. [24][25][26]. 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 , and the notations x • ≡ x µ p µ 1 and x * ≡ x µ p µ 2 for the dimensionless lightcone coordinates (x * = s 2 x + and x • = s 2 x − ). Our metric is g µν = (1, −1, −1, −1) so that p · q = (α p β q + α q β p ) s 2 − (p, q) ⊥ where (p, q) ⊥ ≡ −p i q i . 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. [12] we separate quark and gluon fields in the functional integral (2.5) into three sectors (see figure 2): "projectile" fields A µ , ψ A with |β| < σ a , "target" fields B µ , ψ B with |α| < σ b and "central rapidity" fields C µ , ψ C with |α| > σ b and |β| > σ a and get Our goal is to integrate over central fields and get the amplitude in the factorized 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 |β| < σ a in the projectile fields and |α| < σ b 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 N σas and O m 2 N σ b s , where m N is the hadron's mass. Beyond the tree level, one should expect that the integration over C fields will produce the logarithms of the cutoffs σ a and σ b which will cancel with the corresponding logs in gluon TMDs of the projectile and the target. The result of integration over C-fields has the schematic form and e S eff represents a sum of disconnected diagrams ("vacuum bubbles") in external fields. As usual, since the rapidities of central C fields and of A, B fields are very different, the result of integration over C fields is expressed in terms of Wilson-line operators made form A and B fields.
After integration over C fields the amplitude (2.5) can be rewritten as 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 (2.7) we can set In other words, since A, ψ andÃ,ψ do not depend on x * , if they coincide at x * = ∞ they should coincide everywhere. Similarly, since B, ψ B andB,ψ B do not depend on x • , if they coincide at x • = ∞ they should be equal. Next, in ref. [12] it was demonstrated that due to eqs. (2.10) the effective action S eff (A, ψ A ,Ã,ψ A ; B, ψ B ,B,ψ B ) vanishes for background fields satisfying conditions (2.9). 2 Summarizing, we see that at the tree level in our approximation It is well known that in the tree approximation the double functional integral (2.11) is given by a set of retarded Green functions in the background fields [27][28][29] (see also appendix A of ref. [12] for the proof). Since the double functional integral (2.11) is given by a set of retarded Green functions (in the background fields A and B), calculation of the tree-level contribution tō ψγ µ ψ in the r.h.s. of Eq. (2.11), is equivalent to solving YM equation for ψ(x) (and A µ (x)) with boundary conditions such that the solution has the same asymptotics at t → −∞ as the superposition of incoming projectile and target background fields. The hadronic tensor (2.8) can now be represented as 3 (2.12) 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ψ) then As we will demonstrate below, the relevant operators are quark and gluon fields with Wilsonline 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 (2.12) we need to calculate the functional integral (2.11) in the background of the fields (2.10). Since we integrate over fields (2.10) afterwards, we may assume that they satisfy Yang-Mills ) and similarly for B fields. It is convenient to choose a gauge where A * = 0 for projectile fields and B • = 0 for target fields. The rotation from a general gauge (Feynman gauge in our case, see footnote on p. 5) to this gauge is performed by the matrix Ω(x * , x • , x ⊥ ) satisfying boundary conditions

2)
3 As discussed in ref. [12], there is a subtle point in the promotion of background fields to operators.
When we calculate O as the r.h.s. of eq. (2.11) the fields ΦA and ΦB are c-numbers; on the other hand, after functional integration in eq. (2.5) they become operators which must be time-ordered in the right sector and anti-time-ordered in the left sector. Fortunately, as we shall see below, all these operators are separated either by space-like distances or light-cone distances so all of them (anti) commute and thus can be treated as c-numbers. 4 As we mentioned, for the purpose of calculation of integral over C fields the projectile and target fields are "frozen".
where A * (x • , x ⊥ ) and B • (x * , x ⊥ ) are projectile and target fields in an arbitrary gauge and [x • , y • ] A * z denotes a gauge link constructed from A fields ordered along a light-like line: and similarly for [x * , y * ] B• z . The existence of matrix Ω(x * , x • , x ⊥ ) was proved in appendix B of ref. [12] by explicit construction. The relative strength of Lorentz components of projectile and target fields in this gauge was found in ref. [12] / Here m ⊥ is a scale of order of m N or q ⊥ . In general, we consider W (p A , p B , q) in the region where s, Q 2 q 2 ⊥ , m 2 N , while the relation between q 2 ⊥ and m 2 N and between Q 2 and s may be arbitrary. So, for the purpose of counting of powers of s, we will not distinguish between s and Q 2 (although at the final step we will be able to tell the difference since our final expressions for power corrections will have either s or Q 2 in denominators). Similarly, for the purpose of power counting we will not distinguish between m N and q ⊥ so we introduce m ⊥ which may be of order of m N or q ⊥ depending on matrix element.
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 (2.5) in the background fields of the strength given by eqs. (3.4).

Approximate solution of classical equations
As we discussed in section 2, the calculation of the functional integral (2.11) over C-fields in the tree approximation reduces to finding fields C µ and ψ C as solutions of Yang-Mills equations for the action As we discussed above, the solution of eq. (3.6) which we need corresponds to the sum of set of diagrams in background field A + B with retarded Green functions, see figure 3. The retarded Green functions (in the background-Feynman gauge) are defined as  where and similarly for quarks.
Hereafter we use Schwinger's notations for propagators in external fields normalized according to where we use space-saving notation d −n p ≡ d n p (2π) n . Moreover, when it will not lead to a confusion, we will use short-hand notation The solution of eqs. (3.6) in terms of retarded Green functions gives fields C µ and ψ C that vanish at t → −∞. Thus, we are solving the usual classical YM equations 5 5 We take into account only u, d, s, c quarks and consider them massless.
with boundary conditions 6 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. [12], for our case of particle production with q ⊥ Q 1 it is possible to find the approximate solution of (3.11) as a series in this small parameter. We will solve eqs. (3.11) 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 (3.14). The first step is the calculation of the linear term for the trial configuration (3.14). The quark part of the linear term has the form where are operators in external zero-order fields (3.14). Here we denote the order of expansion in the parameter The gluon linear term is The explicit form of gluon linear terms L (0)a µ and L (1)a µ is presented in eq. (3.26) from our paper [12]. For our purposes we need only the leading term L With the linear terms (3.15) and (3.18), a couple of first terms in our perturbative series are It is convenient to fix redundant gauge transformations by requirements Ai(−∞•, x ⊥ ) = 0 for the projectile and Bi(−∞ * , x ⊥ ) = 0 for the target, see the discussion in ref. [30].
for quark fields and (3.20) for gluon fields (in the background-Feynman gauge). Next iterations, like Ψ [3] (x) and µ (x), will give us a set of tree-level Feynman diagrams in the background fields A µ + B µ and ψ A + ψ B .
Let us consider the fields in the first order in perturbative expansion: Here α, β, and p ⊥ are understood as differential operators α = i ∂ ∂x• , β = i ∂ ∂x * and p i = i ∂ ∂x i . Now comes the central point of our approach. Let us expand quark and gluon propagators in powers of background fields, then we get a set of diagrams shown in figure 3. The typical bare gluon propagator in figure 3 is .

(3.22)
Since we do not consider loops of C-fields in this paper, the transverse momenta in tree diagrams are determined by further integration over projectile ("A") and target ("B") fields in eq. (2.8) 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 The explicit form of operators 1 α+i , 1 β+i , and After the expansion (3.23), the dynamics in the transverse space effectively becomes trivial: all background fields stand either at x or at 0. The formula (3.21) turns into expansion where 1 α and 1 β are understood as 1 α+i and 1 β+i respectively. One may question why we do not cut the integrals in eq. (3.24) to |α| > σ b and |β| > σ a according to the definition of C fields in section 2. 7 The reason is that in the diagrams like figure 3 with retarded propagators (3.24) one can shift the contour of integration over α and/or β to the complex plane away to avoid the region of small α or β. It should be mentioned, however, that such shift may not be possible if there is pinching of poles in the integrals over α or β. For example, if after the expansion (3.23) we encounter 1 (α+i )(α−i ) , the expansion was not justified since actual α's in the integral are ∼ p 2 ⊥ s and hence the field was misidentified: we have a propagator of B-field rather than of C-field. Fortunately, at the tree level all propagators are retarded and the pinching of poles never occurs. In the higher orders in perturbation theory Feynman propagators in the loops cannot be replaced by retarded propagators so after the expansion (3.23) we can get 1 (α+i β)(α+i β ) . In such case the pinching may occur so one needs to formulate a subtraction program to get rid of pinched poles and avoid double counting of the fields.
Note that the background fields are also smaller than typical p 2 ∼ s. Indeed, from eq.
7 Such cutoffs for integrals over C fields are introduced explicitly in the framework of soft-collinear effective theory (SCET), see review [31]. 8 The only exception is the fields B•i or A * i which are of order of sm ⊥ but we saw in ref. [12] that effectively the expansion in powers of these fields is cut at the second term.

Power expansion of classical quark fields
Now we expand the classical quark fields in powers of s (the expansion of classical gluon fields is presented in eqs. (3.35)-(3.38) in ref. [12]). From the previous section it is clear that the leading power correction comes only from the first term displayed in eq. (3.19). Expanding it in powers of p 2 ⊥ /p 2 as explained in the previous section, we obtain In this formula and similarly for 1 β±i . From now on we will denote ψ A It is easy to see that power counting of these quark fields has the form As to quark fields Ψ (1) , we present their explicit form in appendix 8.1 and prove in appendix 8.3 that their contribution is small in the kinematic region s Q 2 .
4 Leading power corrections at s Q 2 q 2 ⊥ As we mentioned in the introduction, our method is relevant to 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 we consider in this paper. 9 As we noted above, we take into account only u, d, s, c quarks and consider them massless. The hadronic tensor takes the form After integration over central fields in the tree approximation we obtain and similarly for J µ B and J µ BA . Hereafter we use notationγ depending on quark's flavor. The quark fields are given by a series in the parameter m 2 ⊥ s , see eqs. (3.27) and (8.2), where Ψ can be any of u, d, s or c quarks. 11 Accordingly, the currents (4.4) can be expressed as a series in this parameter, e.g.
Bs . (4.5) The leading power contribution comes only from product J µ AB (x)J BAµ (0) (or J µ BA (x)J ABµ (0)), while power corrections may come from other terms like J µ A (x)J Bµ (0). We will consider all terms in turn. 10 We denote the weak coupling constant by e/sW and reserve the notation "g" for QCD coupling constant. 11 As we mentioned, we will need only first two terms of the expansion given by eqs. (3.27) and (8.2). region of rapidity. Note that since we want to calculate the leading power corrections, hereafter we 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.

Leading contribution and power corrections from
As to terms ∼Ψ A (0), they can be decomposed using eq. (3.27) as follows: First, let us consider the leading power term coming from the first term in the r.h.s. of this equation.

Leading power contribution
As we mentioned, the leading-power term comes from J and similarly for other matrix elements (summation over color and Lorentz indices is implied).
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 our operators are separated either by space-like intervals or light-like intervals so they commute with each other.
In a general gauge for projectile and target fields these expressions read (see eq. (3.2)) and similarly for A|ψ From parametrization of two-quark operators in section 4.2.1, it is clear that the leading power contribution to W (q) of eq. (4.1) comes from the product of two f 1 s in eq. (4.13) and (4.15). It has the form [32] All other terms in the product of eqs. (4.13) and (4.15) give higher power contributions 12 so they can be neglected at Q 2 s. Similarly, the contribution of two matrix elements in eq. (4.17) is ∼ m 2 ⊥ s in comparison to W lt (q) so it can be neglected as well.

Parametrization of leading matrix elements
Let us first consider matrix elements of operators without γ 5 . The standard parametrization of quark TMDs reads The trivial but important point is that any f (x, k ⊥ ) may have only logarithmic dependence on Bjorken x but not the power dependence ∼ 1 x . Indeed, at small x the cutoff of corresponding longitudinal integrals comes from the rapidity cutoff σa, see the discussion in section 2. Thus, at small x one can safely put x = 0 and the corresponding logarithmic contributions would be proportional to powers of αs ln σa (or, in some cases, αs ln 2 σa, see e.g. ref. [33]). Also, a more technical version of this argument was presented on page 12. for quark distributions in the projectile and for the antiquark distributions. 13 The corresponding matrix elements for the target are obtained by trivial replacements p 1 ↔ p 2 , x • ↔ x * and α ↔ β: Matrix elements of operators with γ 5 are parametrized as follows: (4.17) The corresponding matrix elements for the target are obtained by trivial replacements p 1 ↔ p 2 , x • ↔ x * and α ↔ β similarly to eq. (4.16). 13 In an arbitrary gauge, there are gauge links to −∞ as displayed in eq. (4.11).
Finally, for future use we present the parametrization of time-odd TMDs and similarly for the target with usual replacements p 1 ↔ p 2 , x • ↔ x * and α ↔ β.
Note that the coefficients in front of f 3 , g ⊥ f , h and h ⊥ 3 in eqs.

Power corrections from
The terms in eq. (4.7) proportional to Ξ fields are First, as we demonstrate in appendix 8.3.1, the terms in the second, third, and fourth lines lead to negligible power corrections ∼

Let us start with the term ψ
Next, separating color-singlet contributions we get Using equations (8.9), (8.10), and (8.14) from appendix 8.2 we can rewrite eq. (4.22) as For forward matrix elements we get and similarly for other Lorentz structures in eq. (4.23). The corresponding contribution of the r.h.s of eq. (4.23) to W (α q , β q , x ⊥ ) takes the form 14 Note that for unpolarized hadrons Also, it is easy to see that the last line of eq. (4.23) gives zero contribution. Indeed, let us consider the first term in the r.h.s. of this equation. Since 14 After specifying the projectile and target matrix elements the "A" and "B" labels of the fields become redundant.
and similarly for other Lorentz structures in eq. (4.28). Similarly to eq. (4.25), we get the contribution to W (α q , β q , x ⊥ ) in the form In section 4.3.2, we demonstrated that the matrix elements of quark-antiquark-gluon operators in eqs.
where quark↔antiquark (α q ↔ β q ) term comes from x ↔ 0 contribution in eq. (4.19). As we will demonstrate later, the power corrections which reduce to the leading-power TMDs come with the leading power of 1 Nc in the large-N c approximation -all other power corrections are ∼ 1

Parametrization of matrix elements from section 4.3.1
In this section we will demonstrate that matrix elements of quark-antiquark-gluon operators from section 4.3.1 can be expressed in terms of leading-power matrix elements from section 4.2.1. Let us start with matrix element (4.24) which can be rewritten as (see ref. [22]) Using QCD equations of motion (3.1) we can rewrite the r.h.s. of eq. (4.32) as where we used parametrizations (4.13) and (4.17) for the leading power matrix elements. Now, the second term in eq. (4.33) contains extra α q with respect to the first term, so it should be neglected in our kinematical region s Q 2 q 2 ⊥ and we get For the corresponding antiquark distributions we get The corresponding target matrix elements are obtained by trivial replacements x * ↔ x • , α q ↔ β q and / p 2 ↔ / p 1 .
Next, let us consider Using QCD equation of motion and parametrization (4.18), one can rewrite the r.h.s. of this equation as Again, only the first term contributes in our kinematical region so we finally get (4.40) By complex conjugation we obtain For corresponding antiquark distributions one gets in a similar way (4.42) The target matrix elements are obtained by usual replacements x * ↔ x • , α q ↔ β q and / p 2 ↔ / p 1 .

Sixth line in eq. (4.19)
In this section we consider after Fierz transformation (cf. eq. (4.8)). After separation of color singlet contributions which can be rewritten as where we again used formulas (8.9), (8.11), and (8.14) from appendix 8.2.
Next, it is easy to see that 1 α and 1 β in eq. (4.45) give 1 αq and − 1 βq : where Γ is any of the Dirac matrices in eq. (4.45).
The corresponding contribution to W (α q , β q , x ⊥ ) takes the form where we have used the fact that for the unpolarized hadron. Similarly, so the corresponding contribution to W (α q , β q , x ⊥ ) is Using parametrizations (4.52) and (4.53) from appendix 4.3.4 we obtain the contribution of the 6th line in eq. (4.19) in the form where quark↔antiquark (α q ↔ β q ) term comes from x ↔ 0 replacement, cf. eq. (4.31).

Parametrization of matrix elements from section 4.3.3
In this section we present parametrization of matrix elements from section 4.3.3. Similarly to eqs. (4.40)-(4.42) we define and similarly for the target matrix elements. Note that unlike two-quark matrix elements, quark-quark-gluon ones may have imaginary parts which we denote by functions with tildes. By complex conjugation we get and similarly for the target matrix elements. For completeness, let us present the structure of gauge links in an arbitrary gauge, for example:

Power corrections from
Power corrections of the second type come from the terms In appendix 8.3.4, we will demonstrate that terms ∼ Ψ (1) are small in our kinematical region s Q 2 q 2 ⊥ .
As we prove in appendix 8.3.2, the leading power correction comes from last two lines in eq. (4.56). We will consider them in turn.

Last two lines in eq. (4.56)
Using eq. (3.27) and separating color-singlet matrix elements, we rewrite the sixth line in eq. (4.56) as where we used eqs. (8.10) and (8.11). For the forward matrix elements Similarly, for the seventh line in eq. (4.56) using eqs. (3.27) and (8.12) one obtains Using eq. (4.58) one obtains the contribution to W (α q , β q , x ⊥ ) in the form − g 2 e 2 (a 2 + 1) Here we used the fact that the last term in eq. (4.60) ∼ 2a gives no contribution since same as in eq. (4.27).
Next, using parametrizations (4.66) from the next section we obtain the contribution of the 6th and 7th lines in eq. (4.56) in the form where α q ↔ β q contribution comes as usually from the (x ↔ 0) term in eq. (4.59).

Parametrization of TMDs from section 4.4.1
We parametrize TMDs from section 4.4.1 as follows By complex conjugation we get Target matrix elements are obtained by usual substitutions α ↔ β, / For completeness let us present the explicit form of the gauge links in an arbitrary gauge:

Results and estimates
Combining eqs. (4.12), (4.31), (4.51), and (4.64) we get the leading term and first power corrections to W (q) in the kinematic region s Q 2 q 2 ⊥ in the form where the momentum of the produced Z-boson is q = α z p 1 + β z p 2 + q ⊥ .
For completeness, let us present our final result in the transverse coordinate space Note that in the leading order in N c power corrections are expressed in terms of leading power functions f 1 and h ⊥ 1 . To estimate the order of magnitude of power corrections, one can assume that 1 Nc is a good parameter and leave only first term in the r.h.s. of eq. (5.1): Next, eq. (5.3) is a tree-level formula and for an estimate we should specify the rapidity cutoffs for f 1 's and h ⊥ 1 's. As we discussed in section 2, the rapidity cutoff for f 1 (α z , k 2 ⊥ ) is σ a and for f 1 (α z , k 2 ⊥ ) σ b , where σ a and σ b are rapidity bounds for central fields. Since we calculated only tree diagrams made of C-fields we have σ a = β z and σ b = α z in eq. (5.1). 15 Next, power corrections become sizable at q 2 ⊥ m 2 N where we probe the perturbative tails of TMD's f 1 ∼ 1 34]. So, as long as m 2 N k 2 ⊥ Q 2 we can approximate (up to logarithmic corrections). Similarly, for the target we can use the estimate In general, we should integrate over C-fields in the leading log approximation and match the logs to the double-log and/or single-log evolution of TMDs.
as long as m 2 N k 2 ⊥ Q 2 . Substituting this to eq. (5.1) we get the following estimate of the strength of power corrections for Z-boson production Here we used the fact that due to the "positivity constraint" h ⊥ . 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 term is determined by the logarithmical region q 2 ⊥ k 2 ⊥ m 2 N and the second by Q 2 k 2 ⊥ q 2 ⊥ . By this estimate, the power correction reaches the level of few percent at q ⊥ ≥ 20 GeV. Of course, when q 2 ⊥ increases, the correction becomes bigger, but the validity of the approximation q 2 ⊥ Q 2 1 worsens. Moreover, we have ignored all logarithmic (and double-log) evolutions which can significantly change the relative strength of power corrections.

Power corrections for Drell-Yan process
In this section we consider γ * contribution to the cross section of the Drell-Yan process which is determined by the hadronic tensor where J em µ = e uψu γ µ ψ u + e dψd γ µ ψ d + e sψs γ µ ψ s + e cψc γ µ ψ c is the electromagnetic current for active flavors in our kinematical region.
From the results of the present paper it is easy to extract power corrections to W µ µ . 16 We replace constants a u in eq. (5.1) by e 2 f and remove factors "1" from expressions like a 2 ± 1. One can formally set a u → ∞ inγ µ ≡ γ µ (a u − γ 5 ), divide the result (5.1) by a 2 u , and multiply by e 2 u . After that, we repeat the procedure for other flavors and get Let us present also the large-N c estimate similar to eq. (5.6) Obviously, the relative strength of leading-twist terms and power corrections is the same as for Z-boson production so from our naïve estimate (5.7) one should expect power corrections of order of few percent starting from q ⊥ ∼ 1 4 Q.

Conclusions and outlook
In this paper we have calculated the higher-twist power correction to Z-boson production (and Drell-Yan process) in the kinematical region s Q 2 q 2 ⊥ . Our back-of-the-envelope estimation of importance of power corrections tells that they reach a few percent of the leading-twist result at q ⊥ ∼ 1 4 Q which surprisingly agrees with the same estimate made in ref. [21] by comparing leading-order fits to experimental data.
Of course, we made our estimate without taking into account the TMD evolution, notably the most essential double-log (Sudakov) evolution. One should evolve projectile TMD from σ a = β q toσ a = and match to the result of leading-log calculation of integral over central fields in the rapidity interval betweenσ a andσ b .
To accurately match these evolutions, we hope to use logic borrowed from the operator product expansion. We write down a general formula (2.14) where the coefficient functions c m,n (q, x) are determined by integrals over C-fields and do not depend on the form of projectile or target. To find these coefficients in the first-loop order, we integrate over C-fields in eq. (2.11) with action S C = S QCD (C + A + B, ψ C + ψ A + ψ B ) − S QCD (A, ψ A ) − S QCD (B, ψ B ) but without any rapidity restrictions on C-fields, and subtract matrix elements of the operatorsΦ A (z m )Φ B (z n ) in the background fields A, ψ A and B, ψ B multiplied by tree-level coefficients. Both the integrals over C-fields in eq. (2.11) and matrix elements ofΦ A (z m )Φ B (z n ) will have rapidity divergencies which will be canceled in their sum so what remains are the logarithms (or double logs) of the ratio of kinematical variables (Q 2 in our case) to the rapidity cutoffs σ a of the operatorsΦ A (z m ) and σ b ofΦ B (z n ). Using the above logic we hope to avoid the problem of double-counting of fields which arises when integrals over longitudinal momenta of C-fields got pinched at small momenta (see the discussion in the end of secttion 3.2). The work is in progress.
It should be mentioned that, as discussed in ref. [12], 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 [1]. 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 [36]) equations. Also, since matrix elements of TMD operators with our rapidity cutoffs are UV-finite [37,38], the only UV divergencies in our approach are usual UV divergencies absorbed in the effective QCD coupling.
It is worth noting that recently the treatment of power corrections was performed in the framework of SCET theory (see e.g. refs. [39][40][41]). However, since our rapidity factorization is different from factorization used by SCET, the detailed comparison of power corrections to Z-boson (or Higgs) production would be possible when SCET result for TMD corrections in the form of 1 m 2 Z times matrix elements of quark-antiquark-gluon operators will be available. Let us note that we obtained power corrections for Drell-Yan hadronic tensor convoluted over Lorentz indices. It would be interesting (and we plan) to calculate the higher-twist correction to full DY hadronic tensor. Also, it is well known that for semi-inclusive deep inelastic scattering (SIDIS) and for DY process the leading-order TMDs have different directions of Wilson lines: one to +∞ and another to −∞ [42,43]. We think that the same directions of Wilson lines will stay on in the case of power corrections and we plan to study this question in forthcoming publications.
The authors are grateful to S. Dawson, A. Prokudin, T. Rogers, R. Venugopalan, and A. Vladimirov for valuable discussions. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under contracts DE-AC02-98CH10886 and DE-AC05-06OR23177 and by U.S. DOE grant DE-FG02-97ER41028.

Next-to-leading quark fields
In this Section we present the explicit expressions for the next-to-leading quark fields Ψ (1) . It is convenient to separate these fields in "left" and "right" components: The next-to-leading term in the expansion of the fields (3.19) has the form: where P i = i∂ i + gA i + gB i , see eq. (3.16). The expressions forΨ should be read from right to left, e.g.
(and 1 α ≡ 1 α+i 1 β ≡ 1 β+i as usual). It is easy to see that the power counting of quark fields has the form (cf eq. (3.29)): The gluon fields A (0) • and A (0) * were calculated in ref. [12]: A ab j B jb (8.5) and their power counting reads

Formulas with Dirac matrices
In the gauge A • = 0 the field A i can be represented as Sorting out the color-singlet contributions 18 we get and therefore where 1 α ≡ 1 α+i , see eq. (3.28). For the forward matrix elements we get (8.19) where Γ is any of γ-matrices with transverse indices. Next, consider dx ⊥ dx * e −iβqx * +i(k,x) ⊥ B|ψ(0)gB i (0)ψ(x * , x ⊥ )|B (8.20) where f tw3 (β q , k 2 ⊥ ) is some function of order one (by power counting (3.4) this matrix element (8.20) is ∼ 1). Also, this function may have only logarithmical singularities in β q as β q → 0 but not the power behavior 1 βq . 19 The corresponding contribution to W (q) of eq. (4.1) is proportional to so it can be neglected in comparison to the contributions ∼ q 2 ⊥ Q 2 W lt (recall that we assume that Z-boson is produced in a central range of rapidity so . In a similar way one can show that the remaining three terms in the second and third lines of eq. (4. 19) give small contributions to W (q).
Next, it is easy to see that the matrix element of the fourth line of eq. (4.19) vanishes. Indeed, let us consider the first term in the fourth line and perform Fierz transformation (4.8): The matrix element B|Â ai (x)Â aj (0)|B for unpolarized hadrons can be proportional either to 2x i x j + x 2 ⊥ g ij or to g ij . Since the former structure does not contribute due to (2x i x j + x 2 ⊥ g ij )γ i/ p 2 γ k ⊗ γ k/ p 2 γ j = 0 (8. 28) we get For the forward target matrix element one obtains dx * e −iβqx * B|Â a i (x)Â ai (0)|B (8.30) where we used parametrization (4.6) from ref. [12]. Since the gluon TMD D g (x B , x ⊥ ) behaves only logarithmically as x B → 0 [38], the contribution of eq. (8.29) to W (q) is of order of First, let us notice that terms likē give zero contribution sinceψ A (x)γ µ ψ A (x) does not depends on x * so dx * e −iβqx * = 2πδ(β q ). Let us consider now the last two lines in the power expansion (4.6) of J µ AB (x)J BAµ (0): B (x)Ψ