Twist expansion of Drell-Yan structure functions in color dipole approach

Forward Drell-Yan process at the LHC probes the proton structure at very small Bjorken-$x$ and moderate hard scales. In this kinematical domain higher twist effects may be significant and introduce sizeable corrections to the standard leading twist description. We study the forward Drell-Yan process beyond the leading twist approximation within the color dipole model framework that incorporates multiple scattering effects. We derive the Mellin representation of the forward Drell-Yan impact factors for fully differential cross-sections. These results combined with the color dipole cross-section of the saturation model are used to perform twist expansion of the Drell-Yan structure functions at arbitrary transverse momentum $q_T$ of the Drell-Yan pair and also of the structure functions integrated over $q_T$. We also investigate the Lam-Tung relation, find that it is broken at twist 4 and provide explicit estimates for the breaking term.


Introduction and conclusions
The forward Drell-Yan (DY) processes at the LHC are expected to provide the most sensitive measurements of parton densities in the proton down to very small x ≃ 10 −6 . At the LHCb experiment the DY lepton-antilepton pair may be measured down to invariant mass M of about 2.5 GeV, so the related parton density to the scale µ 2 ≃ 6.25 GeV 2 [1]. This kinematic region has been never probed before. It extends HERA measurements of proton structure at moderate scales towards small parton x by about two orders of magnitude. In fact, the forward Drell-Yan process is a unique tool to explore this region. Hence, it is mandatory to acquire deep theoretical understanding of the process in QCD.
The region of moderate scales and very small x is sensitive to interesting QCD effects. The standard DGLAP description of parton densities evolution in the proton may be significantly affected by small-x resummation effects [2,3] and higher twist contributions related (but not identical) to multiple scattering corrections. Therefore, the forward Drell-Yan process may be used as a sensitive probe of these effects. On the other hand, good theoretical understanding of higher twist corrections to the Drell-Yan cross-section is necessary to extract the standard twist-2 parton densities with higher precision and reduced uncertainties.
Thus, in this paper we address the problem of higher twist effects in the forward Drell-Yan processes.
The forward Drell-Yan process has multiple advantages as the probe of the proton structure and QCD dynamics. The presence of hard electromagnetic probe allows for effective application of perturbative QCD. From the experimental side, the kinematic variables of the final state composed of a lepton-antilepton pair may be measured with good precision, giving access to multiple differential distributions. In particular the Drell-Yan pair angular distributions are determined by four invariant structure functions, T i , i = 1, . . . , 4, describing proton interactions with a virtual photon which mediates the lepton pair production [4,5]. All structure functions T i may be decomposed into twist-series using the Operator Product Expansion (OPE), in which the leading twist-two contributions may be computed using standard parton densities and the (unknown) higher twist terms are suppressed by negative powers of the process hard scale. This power suppression, however, may be compensated in the region of moderate scale and very small-x by rapidly growing higher twist matrix elements, so that the higher twist contributions are expected to become important below µ 2 ∼ 30 GeV 2 [6]. 1 Hence the four DY structure functions carry enriched information on higher twists: higher twist hadronic matrix elements are coupled to four different coefficient functions. This gives opportunities for broader and more detailed tests of the higher twist description and provides tools for more efficient isolation of higher twist corrections. In particular, it follows from the famous Lam-Tung relation [5] that certain combination of the DY structure functions vanishes at twist-2 (up to next-to-next-to leading order corrections), and therefore its deviation from zero is a sensitive probe of higher twist effects.
As yet the higher twist contributions to the proton structure and to Drell-Yan structure functions are poorly known from experiment and the rigorous theoretical description of higher twist terms within QCD is highly involved. A treatment of the higher twist contribution in DY scattering within collinear QCD proposed in Ref. [7], see also [8,9], still requires modeling of higher twist matrix elements. A practical way to circumvent these obstacles is to adopt, as a first approximation, eikonal or Glauber-Mueller picture where multiple scattering in QCD is a product of independent single scatterings. This approach was implemented e.g. in the very successful Golec-Biernat-Wüsthoff (GBW) saturation model [10].
The GBW model was proved to provide an efficient and accurate unified picture of multiple HERA processes at small x: deeply inelastic scattering (DIS), diffractive DIS, elastic vector meson production, deeply virtual Compton scattering [10,11]. In the high energy limit in QCD, relevant for the forward Drell-Yan processes at the LHC, the DY scattering amplitudes may be computed using the high-energy factorization (k T -factorization framework [12]) and the transverse position space. This leads to the 'color dipole picture' [13] of the forward DY process, proposed by [14,15] in which the QED-QCD partonic amplitudes of virtual photon production are combined with the color dipole cross-section, that needs to be fitted and/or modeled (beyond the leading order contribution at the leading twist). The dipole cross-section is constrained by HERA data, and the fit results may be applied to predict DY cross-sections, including higher twist effects [6,[14][15][16][17], see also Ref. [18]. Thus, we shall estimate higher twist effects in the structure functions of the forward DY processes using the color dipole approach and the GBW saturation model. One should stress that at the leading twist the dipole approach is consistent with the standard collinear picture results in the high energy limit, up to NLO, but it provides in addition an estimate of multiple scattering and higher twist effects.
In this paper the cross-section decomposition into its twist components is carried out in the Mellin representation for the color dipole sizes. This decomposition method was initially proposed in Ref. [19], further developed in Ref. [20] and then applied to the total forward DY cross-section [6] and diffractive DIS [21]. In this framework twist contributions are related to complex singularities (poles or branch points) in the Mellin plane. We follow the technique of the total DY cross-section twist decomposition proposed in [6], but extend the results to differential cross-sections in the lepton angles and the DY pair transverse momentum q T . In more detail we compute the Mellin transforms of all forward DY impact factors at given q T and combine the results with the color dipole cross-section to get q Tdependent twist decomposition of all four DY structure functions. Here we restrict ourselves to the simplest GBW eikonal model of the color dipole cross-section but it is straightforward to combine the obtained Mellin representation of the impact factors with other descriptions or parameterizations of the color dipole cross-section.
The findings of this paper may be summarized as follows. The key novel result are the Mellin representations of the forward DY impact factors for all the DY structure functions at arbitrary transverse momentum q T . These impact factors follow directly from perturbative QCD. The impact factors are then combined with the GBW color dipole cross-section to get analytic form of the twist expansion of the DY structure functions. These results are model dependent but exhibit some generic features driven by the perturbative part, like e.g. saturation of the Lam-Tung relation at twist 2, presence or absence of hard scale logarithms. We obtain results both for helicity structure functions and for invariant structure functions. We find that the Lam-Tung relation holds at twist 2 but it is broken at twist 4, and the breaking term is leading in perturbative QCD at this twist. Thus, the Lam-Tung combination of the DY structure functions is a promising observable for experimental measurements of higher twist effects. In the GBW model of the color dipole cross-section the higher twist terms are power-enhanced with decreasing x (besides the generic suppression by negative powers of the hard scale), at twist τ one has T (τ ) Also q T -integrated structure functions are derived. Interestingly enough, this integration leads to emergence of a twist 3 contribution in the LT helicity structure function. Results of this paper provide necessary tools for a forthcoming experiment oriented analysis of the LHC potential to measure higher twist components of the proton structure within the color dipole approach, and prove that the applied twist decomposition method is effective. Figure 1. Dominant diagrams for the forward Drell-Yan process in the t-channel helicity frame (in which the target is at rest), see the text for the notation.

Kinematics and notation
We consider the high energy proton-proton collision with a lepton-anti-lepton pair, l + l − = e + e − or µ + µ − , in the final state, p(P 1 )p(P 2 ) → l + l − X in which the pair is produced in the fragmentation region of one of the protons and the leptons four-momenta l + and l − are measured. At the leading order in QED this process is mediated by a virtual photon γ * (q) , with four-momentum q = l + + l − , and the virtuality q 2 = M 2 > 0 is the lepton pair invariant mass squared. It is convenient to introduce also κ = l + − l − . The proton projectiles four momenta are P 1 and P 2 and they are near light-like, in the center of mass system (c.m.s.) of the pp pair P 1 ≃ ( where invariant collision energy squared S = (P 1 + P 2 ) 2 is much greater than the proton mass squared, m 2 p . We define light-like components of the momenta as p ± = p 0 ± p z , where the z axis is given by the beam direction in the c.m.s. From now on we shall use the light-cone coordinates for four-vectors, v = (v + , v − ; v T ). The Sudakov decomposition of four momenta will be employed, p i = α i P 1 + β i P 2 + p ⊥ , where p ⊥ is the p four-momentum component in the plane perpendicular to P 1 and P 2 . For the DY virtual photon we have q = α q P 1 + β q P 2 + q ⊥ , with q ⊥ = (0, 0; q T ) and q 2 ⊥ = − q 2 T . The forward region is defined by the condition β ≫ α, and the photon (or the DY pair) rapidity y = 1/2 log(β/α). At the LHC the backward region α ≫ β is, in fact, fully symmetric to the forward region, so for simplicity we shall restrict our discussion to β ≫ α.
At the leading order of QED and QCD at parton level the DY hard subprocess is q 1 (p 1 )q 2 (p 2 ) → γ * (q) → l + l − , where q 1 (q 2 ) is a quark (anti-quark) coming from one of the protons (the other proton) that carries four momentum p 1 (p 2 ). At the NLO QCD the partonic subprocesses include also quark (anti-quark)-gluon contributions, q 2 (p 2 )g(k) → q 2 (p ′ 2 )γ * (q) → q 2 l + l − . In the high energy limit p 1 = x 1 P 1 + p 1⊥ and p 2 = x 2 P 2 + p 2⊥ , where x i are the parton x-variables, see Fig. 1.
In the forward region at the LHC one has x 2 ≫ x 1 , x 2 is a sizable fraction of the longitudinal proton momentum P + 2 , say x 2 ∼ 0.1 and x 1 is very small, down to x 1 ∼ 10 −6 . In this region the quark q 1 distribution function (d.f.) is strongly dominated by the seaquarks and the valence quark contribution to q 1 d.f. may be neglected. At small x the evolution of the sea-quark (or anti-quark) distribution function of the target proton is driven by the gluon evolution. In more detail, due to suppression of quark propagation over large rapidity distance, the dominant contribution to sea-quark d.f. comes from diagrams where the sea-quark emerges from the gluon in the last splitting of the QCD evolution. Therefore, in the k T factorization framework the dominant diagrams with the gluon exchange in the t-channel are given in Fig. 1. The topology of the diagram at the right-hand side coincides with the topology of the LO collinear contribution with the sea-quark emerging from the gluon at the last splitting. The diagram at the left-hand side of Fig. 1 contributes to the NLO qg partonic subprocess in the collinear picture.
The diagrams shown in Fig. 1 are the basis of the color dipole approach to the forward DY process. Thus, one considers the virtual photon emission by a fast quark coming from p(P 2 ) in the scattering off the proton p(P 1 ) by a small-x gluon exchange, q( The effective color dipole emerges here through the interference of amplitudes of the virtual photon emission before and after the quark scattering off the target proton. The effective color dipole size corresponds to a displacement of the quark position in the transverse space due to the γ * emission. In the color dipole formulation one relies on k Tfactorization (high-energy factorization) approach, where the gluon transverse momentum does not vanish and the standard gluon distribution function xg(x, µ 2 ) gets replaced by a On the other hand, the fast quark q(p 2 ) carries a rather large x 2 , so its transverse momentum may be neglected, p 2 is saturated by p + 2 .
Having defined the relevant partonic channel and diagrams, we complete fixing the notation: let the gluon-x be the Feynman x of the virtual photon (or DY pair). The fast quark distribution function in proton p(P 2 ), taken in the collinear limit, is denoted by ℘. Finally, it is convenient to analyze the process in the helicity basis, so we denote the helicities of the incoming and outgoing quark (q(p 2 ) and q(p ′ 2 )) by λ 1 and λ 2 respectively, and polarizations of γ * by σ. Since the photon is virtual, it has three polarization states.

Structure functions
The standard description of the differential DY cross-section employs so-called helicity structure functions W L , W T , W T T , W LT [4,5] (see Appendix A). In this approach one factorizes leptonic and hadronic degrees of freedom by contracting both hadronic and leptonic tensors with virtual photon polarization vectors (PPVs). The leptonic tensor reduces to a distribution of lepton angles Ω = (θ, φ) in lepton pair center-of-mass frame while the result of contraction of hadronic tensor with different PPVs are the W -structure functions. The differential DY cross-section is then given by the formula: (3.1) The form of W -structure functions depends on arbitrary choice of axes (which defines PPVs) in lepton pair center-of-mass frame. In this paper we perform calculations in a frame with the Z axis anti-parallel to the target's momentum and the Y axis orthogonal to the reaction plane (in [4] this frame is called the t-channel helicity frame). In order to avoid the helicity frame dependence one introduces invariant structure functions, T i . They are defined as coefficients of hadron tensor decomposition [4]: The invariant DY structure functions are related to the helicity structure functions in the t-channel helicity frame in the following way, see Appendix A for the derivation. The DY helicity structure functions in any helicity frame may be expressed through the invariant structure functions and the explicit formulae for several standard frames may be found e.g. in Ref. [4].
In the following we shall re-derive the DY helicity structure functions in the k Tfactorization approach. Thus, the scattering amplitudes of the fast quark will be computed within the high energy limit of QCD and the corresponding helicity dependent cross-sections will be represented in terms of the impact factors. In order to account for multiple scattering effects we shall introduce the color dipole cross-section. Next, the Mellin representations of the impact factors and of the helicity structure functions will be given.

Forward Drell-Yan impact factors
In the framework applied the forward DY cross-section takes the following form [14][15][16], where the helicity dependent γ * impact factors arẽ the leptonic tensor in the helicity basis reads and ℘(x F /z) is a collinear parton distribution function for the projectile. The amplitudes A (σ) λ 1 ,λ 2 ( q T ) of the virtual photon emission with polarization σ and transverse momentum q T are given in Fig. 1. In the target rest frame P 1 = (m p , m p ; 0) and we choose a standard set of polarization vectors, In the lepton c.m.s. where l + = − l − we define the spatial axes (X, Y, Z) through the γ * polarization vectors, The leptonic helicity tensors, L σσ ′ are then expressed through a set of three scalar products where Ω = (θ, φ) are the standard angles of the spherical coordinate system in the lepton CM frame.
With the chosen set of γ * polarization vectors the DY γ * emission amplitudes take the form where we suppressed the dependence of A (±) λ 1 ,λ 2 ( q T ) on z and k T . The inverse Fourier transforms of the amplitudes to the transverse position space read which is convinient to rewrite as Using representation (3.15) of the amplitudes, with the dipole eikonal factor, 1−e −iz k· r , one may express the forward DY cross-section (3.4) in terms of the color dipole cross-section [14,15],σ where we introducex g being the value of x g at the process threshold. In what follows we suppress thex g dependence ofσ(r). At the leading order the dipole cross-section is in one-to-one correspondence with the unintegrated gluon density, f (x g , k 2 T ), and their inverse relation reads where ∇ 2 is the Laplace operator in two transverse dimensions.
Using the last formula one can rewrite (3.4) as Substituting the inverse Fourier transform of (3.14) into (3.21) and integrating over d 2 k T one obtains the impact factors in the transverse position representation, Finally, let us parameterize the DY structure functions W i in terms of new functions Φ i , Functions Φ σσ ′ and Φ i play a similar rôle to virtual photon wave functions of the color dipole model for DIS, but for the DY also off-diagonal amplitude products in the helicity basis contribute. The functionsΦ σσ ′ are related by (3.21) to standard impact factors in momentum space, so we shall use for them the term impact factors as well, specifying Φ i functions as leptonic impact factors.

Mellin representation of the impact factors
Mellin representation is a useful tool in analysis of QCD amplitudes, in particular of their twist structure. In order to find this representation for the forward DY process we start introducing a Mellin transform of the dipole cross-section, Using (3.28) we rewrite (3.23) as: is the Mellin representation of the leptonic impact factor. The leptonic impact factors Φ i (q T , r, z) are linear combinations of Φ σσ ′ (q T , r, z) (see formulae (3.24)-(3.27)) where Φ σσ ′ are given in (3.22). After integration (3.31) over d 2 r and d 2 r 1 d 2 r 2 we get the following results for Mellin transforms of the leptonic impact factors,Φ The above formulae are one of the main results of this paper. Inserted into (3.30) they allow for efficient analysis of the forward DY cross-section in terms of double integrals over z and s. Let us stress that the obtained Mellin forms of the impact factors do not depend on the chosen form of the color dipole cross-section, they follow directly from perturbative (leading order) QCD amplitudes in the high energy limit.

Mellin representation of the integrated helicity structure functions
For a description of inclusive measurements of the DY processes one applies q T -integrated helicity structure functions,W Clearly, the results depend on the chosen set of virtual photon polarization vectors, in particular one has to account carefully for possible correlation of the polarization vectors with q T .W i can be found by the direct integration of the helicity impact factorsΦ i (q T , s, z) (given by (3.32)-(3.35)) over d 2 q T and then use (3.30). In the t-channel helicity frame the integration over the azimuthal angle is trivial and gives 2π, and the remaining integrals can be performed using a new variable q T /η z . One obtains: where functions χ 1,2 are define by their integral representations (3.41) These expressions provide a compact representation of the q T -inclusive forward DY crosssections differential in the angles of DY leptons in the t-channel helicity frame.

Twist expansion of DY structure functions and Lam-Tung relation
With the Mellin representations of the leptonic DY impact factors derived in the previous section, it is possible to evaluate the DY structure functions and perform their twist analysis. For this one needs, however, to specify the dipole cross-section,σ( ρ). Within OPE approach to hard scattering in QCD,σ( ρ) includes a tower of matrix elements with increasing twist. Currently the higher twist components of the proton structure at small x are not known from experiment. Also theoretical analysis within perturbative QCD is not capable to predict the non-perturbative inputs for evolution of higher twist operators from the first principles. Therefore, current theoretical estimates of higher twist effects must rely on certain assumptions and models. So far two main QCD-inspired approaches to dipole cross-sections beyond leading twist were used in higher-twist analysis: the eikonal Glauber-Mueller scattering cross-section (used in the GBW saturation model) and the Balitsky-Fadin-Kuraev-Lipatov (BFKL) and Balitsky-Kovchegov (BK) dipole crosssections, obtained within small-x resummations in QCD [2,3,22,23]. With both forms of the dipole cross-sections one can describe well the existing data on the proton structure from deeply inelastic scattering (DIS) at small x, see e.g. Ref. [11], but as found in Refs. [20,24], the inclusive DIS data probe the proton structure efficiently only at the leading twist. Hence, current theoretical predictions of higher twist effects are uncertain. Still, it is important to estimate these effects within different models in order to design measurements probing higher twist effect and to make optimal use of the resulting data. In particular, measurements of higher twist contributions in the forward DY structure functions may be used to discriminate between models of higher twist / multiple scattering effects and provide essentially new information about proton structure and QCD evolution beyond the leading twist.
As a first step towards understanding the higher twist structure of the forward DY scattering we consider the simplest case of the eikonal dipole cross-section used in the GBW approach. Using this cross-section we perform an explicit analytic twist decomposition of the forward DY structure functions. This test case exhibits already some interesting features, clearly visible due to simple analytic form of the results. A more detailed numerical analysis of the forward DY process at the LHC, taking into account also the BFKL / BK dipole cross-section will be performed in a forthcoming paper [25].

Twist decomposition of helicity structure functions W i
The twist decomposition of the forward DY structure functions W i may be performed using the Mellin representation (3.30), (3.32)-(3.35) supplemented by the GBW dipole cross section, that takes the form [10], and its Mellin transform readsσ(−s) = −σ 0 Γ(−s).
In order to define properly the twist content of the DY structure functions one should specify the scale that plays the rôle of the OPE hard scale. The forward DY structure functions (3.30), however, depend on two hard scales, M and q T , as visible from the form of the impact factors (3.32)-(3.35), so the choice is ambiguous. In order to avoid this ambiguity we take advantage of the fact that in the OPE series the powers of the hard scale µ are accompanied by the corresponding opposite powers of the hadronic (soft) scale Λ, so that twist-τ terms enter with scale ratio (Λ/µ) τ . In the color dipole formulation the soft scale enters only through the dipole cross-section in accordance with the requirement that the soft scale should be related to the target hadron structure. The GBW dipole cross-section carries a unique soft scale -the saturation scale Q 0 (x g ). In order to match the OPE series structure we identify Λ = Q 0 and define the twist-τ terms as those that scale in Q 0 as Q τ 0 . 2 . The problem of the ambiguous hard scale may be understood from the OPE perspective. In OPE the forward DY structure functions depend on non-perturbative matrix elements of hadronic operators and perturbative coefficient functions. The hadronic operators depend on the hadronic scale Λ and the factorization scale µ, and the coefficient functions depend on M , q T , and µ. Hence, if the factorization scale µ is identified with one of hard scales, M or q T , then in the DY structure functions computed within the color dipole approach dependencies on M and q T coming from the matrix element and from the coefficient functions are entangled and, in general, the scaling in M or q T cannot be used to isolate the terms with definite twist. On the other hand, the procedure based on determination of terms with definite positive powers of Q 0 leads to a unique definition of the twist expansion in the adopted approach.
The twist decomposition of the DY structure functions reduces, therefore, to determination of their Q 0 power series expansion. This may be done with a technique described in Refs. [6,19,20], based on analytic structure in s of the Mellin representation. One closes the contour C of the inverse Mellin transform (3.30) with a semicircle at complex infinity for Res > 0. The integrands are analytic in s except of isolated singularities so one can express the inverse Mellin integrals as sums of contributions coming from the singularities in complex s and these contributions carry definite Q 0 scaling, hence the twist. For the GBW dipole cross-section the only singularities in the complex plane of the Mellin variable s in integrands (3.30) are simple poles coming from Γ(−s), contained in the dipole cross-section, so the evaluation of the s-integrations by taking the residues is straightforward. Integrals over z in (3.30) are convergent for Re s > 0. The obtained leading twist contributions to the forward DY structure functions are the following, (4.5) It is clearly visible from the above equations, that the structure functions have the following leading behavior at small values of the photon transverse momentum, q T ≪ M : W L,T T /σ 0 ∼ O(q 2 T /M 2 ). Therefore, there is a hierarchy of contributions in which W (2) T dominates over W (2) LT and finally come W (4.9) In general, only the even twist contributions are found in the structure functions. Evaluation of twist τ > 4 contributions may be easily done as well.
The obtained results on twist content of the structure functions exhibit certain common features. Clearly, at twist τ one has a suppression by negative powers of the scale and an enhancement due to growth of the saturation scale Q 0 (x g ) ∼x −λ g with the decreasingx g , This small x g -enhancement is specific to eikonal color dipole model and not necessarily holds in other models. One sees as well that the two hard scales of the problem, M and q T , enter in combinations dependent on the term and none of it can be identified as the unique main hard scale of the problem. Moreover, the denominators of the integrands are proportional to powers ofM 2 = M 2 (1 − z) + q 2 ⊥ . Thus, in the region z → 1 and q T → 0 the scaleM → 0 and one finds soft singularities. We will discuss the treatment of this problem in Sec. 4.3.
The terms of twist expansions of the helicity structure functions (4.2)-(4.5) and (4.6)-(4.9), are valid in the t-channel helicity frame only. However, these results may be used to obtain the twist expansion of the invariant structure functions, T i . Twist-2 contributions to the invariant structure functions are given in Appendix B.

Twist expansion of integrated structure functions
It is easy to see that expressions for twist 2 and twist 4 components of the DY structure functions derived in the previous sections cannot be simply integrated over d 2 q T . The reason is a divergence at q T → 0, z → 1, that follows as a consequence of the form of denominators in the integrands, ∝ (q 2 T + M 2 (1 − z)) n . For example, the q T integration in W T leads to a divergent z-integral: (4.10) In order to define properly the twist expansion of the integrated structure functions one has to perform the q T integration of (3.37)-(3.40) prior to evaluating the inverse Mellin integral and carrying out the twist decomposition. Then, the integration over d 2 q T dz introduces additional poles at positive integer values of the Mellin variable s, and double and single poles in s are found in the Mellin representation of the integrated DY structure functions. The problem with z → 1 limit of the integrands in twist expansion of integrated forward DY structure functions was found and solved in Ref. [6]. In short, in this prescription one assumes a power series structure of the structure function integrands in the variable 1 − z (modulo logarithms) at z → 1 and treats analyticly terms of this series that would naïvely lead to divergent integrals of twist components. The corresponding integrals lead to additional poles in s, ∼ 1/(s − τ ) at twist τ . The s integration of these terms may be then performed analyticly. The additional s-singularity increases the order of the twist poles by one, which results with an additional logarithm ln M 2 /Q 2 in the structure functions. The remaining part of the structure function integrands left after separation of the apparently singular terms, lead to convergent integrals that can be directly evaluated by numerical integration. Thus, following the prescription of Ref. [6] we get for the leading twist, Interestingly enough, after the q T integration a non-zero twist 3 contribution emerges for W LT ,W This twist 3 emergence occurs in the perturbative part due to the presence of the √ 1 − z factor in (3.40). It comes from the singular region, z → 1, q T → 0. Since no twist 3 singularities are present in the color dipole cross-section, this contribution comes from a single pole in s and carries no ln(M 2 /Q 0 ). The twist 3 contributions to the other structure functions vanish,W T T = 0. It means, in particular, that the twist 3 contribution affects only the lepton angular distribution and not the total cross-section.
Expressions for twist 4 contributions in the integrated structure functions take the following form, The obtained results for twist components ofW (τ ) T andW (τ ) L may be directly compared to the results of Ref. [6] after an integration of the lepton angles in the DY cross-section. This comparison was performed and full agreement was found. Contributions of two other structure functionsW T T andW LT to the DY cross-section vanish after the integration over the lepton angles so they did not appear in the analysis of Ref. [6] and the twist decomposition of q T -integrated DY structure functionsW T T andW LT is a novel result.
In accordance with Ref. [6] the obtained expression for twist-τ components of the DY integrated structure functions take a general form, with numerical coefficientsÃ (τ ) i given in explicit analytic form andB (τ ) expressed as integrals over z. BothÃ (τ ) andB (τ ) depend on the fast quark density ℘(x, M 2 ), in the projectile proton. The convergent integralsB (τ ) may be obtained by numerical integration. In general termsB (τ ) have no enhancement by ln(M 2 /Q 2 ), so they are subleading, however they may be numerically important at moderate M 2 .

Lam-Tung relation
One of the main goals of this paper is to prepare the ground for experimental measurements of the higher twist contributions in the forward DY scattering. Optimal observables for this purpose should have a suppressed leading twist contribution. One of such observables is given by the well known Lam-Tung relation [5]. Thus, within the parton model the following relation between structure functions was proven [5]: Using (3.4) this can be rewritten as the equation for W i [16]: This relation holds in perturbative QCD at twist 2 up to next-to-next-to-leading order correction, see e.g. [16].
Our results are consistent with the Lam-Tung relation. Indeed, from (4.2) and (4.4) one sees that W T T = 0. Moreover, as in Refs. [8,9,16] we find a breakdown of the Lam-Tung relation at twist 4, and its integrated version, The breakdown term carries ln(M 2 /Q 2 0 ) so the Lam-Tung relation breaking at twist 4 enters at the leading order in QCD. Hence, the Lam-Tung combinations (4.19) and (4.20) of the forward DY structure functions exhibit enhanced sensitivity to higher twist effects and so they are promising observables for experimental finding of higher twist effects at small x.

Discussion and outlook
The estimates of higher twist contributions to the Drell-Yan structure functions given in this paper are based on eikonal GBW color dipole model cross-section. Although this model provides an efficient unified description of the small x DIS data down to photoproduction limit, diffractive DIS and exclusive vector production, its twist content was not tested experimentally yet. In particular, HERA data are consistent with leading twist description, except of the kinematic edge of very small x and moderate scales and also with Balitsky-Fadin-Kuraev-Lipatov (BFKL) and Balitsky-Kovchegov (BK) cross-sections, that have a different twist composition than the GBW model. Therefore, our results for twist expansion of the forward DY structure functions are model-dependent predictions. Such approach is justified by lack of higher twist measurements at small x, and the need to estimate the LHC potential to resolve higher twist effects. On the other hand, clearly, for this latter goal one needs also to provide predictions dedicated for the LHC within other reasonable schemes, e.g. within the BFKL / BK approximation. This paper is a necessary step towards such a broader analysis of higher twists effects in the forward DY scattering at the LHC. It summarizes key theoretical results needed for higher twist extraction from the color dipole picture. In particular, the Mellin representation of the DY impact factors is a novel model independent result, following directly from perturbative QCD. This means in turn that certain features of the higher twist structure, following from properties of the impact factors, are generic, for example the breakdown of the Lam-Tung relation at twist 4. Also presence or absence of ln(M 2 /Q 2 0 ) enhancement factor in twist components of various structure functions is generic. It should be also stressed that the obtained Mellin forms of the DY impact factors may be also used within the BFKL formalism for small x resummations. The next step in our program of theoretical analysis of the LHC higher twist will be a data oriented study in which more models of the dipole cross-section will be considered and suitably refined to reproduce the existing forward DY data. Thus, we leave explicit numerical predictions for the higher twist corrections to forward DY scattering at the LHC to the second part of the analysis to be presented in a forthcoming paper [25].