Twist decomposition of Drell-Yan structure functions: phenomenological implications

The forward Drell--Yan process in $pp$ scattering at the LHC at $\sqrt{S}=14$ TeV is considered. We analyze the Drell--Yan structure functions assuming the dominance of a Compton-like emission of a virtual photon from a fast quark scattering off the small $x$ gluons. The color dipole framework is applied to perform quantitatively the twist decomposition of all the Drell--Yan structure functions. Two models of the color dipole scattering are applied: the Golec-Biernat--W\"{u}sthoff model and the dipole cross section obtained from the Balitsky--Fadin--Kuraev--Lipatov evolution equation. The two models have essentially different higher twist content and the gluon transverse momentum distribution and lead to different significant effects beyond the collinear leading twist description. It is found that the gluon transverse momentum effects are significant in the Drell--Yan structure functions for all Drell--Yan pair masses $M$, and the higher twist effects become important for $M \lesssim 10$ GeV. It is found that the structure function $W_{TT}$ related to the $A_2$ angular coefficient and the Lam--Tung observable $A_0 -A_2$ are particularly sensitive to the gluon $k_T$ effects and to the higher twist effects. A procedure is suggested how to disentangle the higher twist effects from the gluon transverse momentum effects.


Introduction
The Drell-Yan process is a classical probe of the proton structure and of the strong interactions in hadron collisions [1]. The experiments operating at the Large Hadron Collider have already detected large statistic of the Drell-Yan dileptons and measured the differential Drell-Yan cross sections as functions of several kinematical variables [2][3][4][5][6][7][8][9]. In particular the angular distributions of the dileptons were measured at the Z 0 peak that allow the determination of the Drell-Yan structure functions at the mass close to the Z 0 boson mass [6,9]. The Drell-Yan measurements when extended to the low mass region, M < 10 GeV, may be used to provide unique information about parton densities in the proton at a very low x, at or below x = 10 −5 [10]. In this kinematic region of the low M and very small x the higher twist corrections may affect significantly the Drell-Yan cross section and hence the parton density function determination from the data, see e.g. Ref. [11,12]. Therefore, in order to achieve the highest precision of the parton density function determination the higher twist corrections should be taken into account in the analysis. Unfortunately not much is known about the higher twist contributions to the proton structure. The subleading twist 4 corrections are represented by a set of independent operators whose matrix elements have not been measured yet, see e.g. Refs. [13,14]. Therefore in the estimates of the higher twist corrections to proton scattering cross sections it is still necessary to relay on models. An approach to model of the higher twist effects in high energy scattering at small x was proposed on the basis of the Golec-Biernat-Wüsthoff (GBW) saturation model [15], that provides an efficient unified picture of the high energy scattering down to very low scales where multiple scattering and higher twist effects are expected to contribute. The framework for extraction of the twist components from the GBW model was formulated for the DIS at HERA in [16,17], then further developed and applied to the diffractive DIS at HERA [18] and the forward Drell-Yan cross sections [11,12]. Within this framework an evidence of the higher twist corrections to DDIS structure functions was found in [18] and recently, a related approach revealed an evidence of the higher twist corrections to the proton structure functions at small x [19].
A framework that is capable to provide a QCD guideline for extending theoretical analysis beyond the twist 2 collinear approach, is the k T factorization formalism [20][21][22][23][24]. The treatment of the forward Drell-Yan process within the k T factorization was initially proposed in [25] in the color dipole representation [26], and then further developed and applied to data analysis in numerous papers, see e.g. Refs. [27][28][29][30][31]. Later on also the momentum representation of the forward and general Drell-Yan process were elaborated in detail [32,33]. The dipole formulation of the forward Drell-Yan scattering was used to obtain the twist decomposition of the Drell-Yan cross section integrated over the dilepton angular distribution [11]. In the latter analysis the GBW of the QCD dipole cross section was assumed. In a recent paper [12] we prepared the theoretical framework to extend this type of twist analysis to all the Drell-Yan structure functions. We also discuss there in more detail earlier estimates of the higher twist content of the Drell-Yan cross section performed in Refs. [34][35][36][37][38][39][40].
In the present paper we apply the results of our earlier paper [12] to perform quantitative estimates of the higher twist contributions to the Drell-Yan structure functions based on the GBW saturation model [15]. Moreover we derive the analytic formulae for the twist decomposition of the forward DY structure functions assuming the Balitsky-Fadin-Kurayev-Lipatov (BFKL) pomeron exchange [21,22]. The two approaches assume an essentially different dynamics of multiple hard scattering and have an essentially different twist content. In the GBW model a simple eikonal picture of multiple scattering is applied corresponding to the resummation of independent single exchanges. This leads to the higher twist amplitudes strongly enhanced by inverse powers of x at small x. The BFKL pomeron exchange amplitude emerges as a QCD result obtained from the resummation of the leading logarithms of x. The BFKL pomeron exchange implicitly carries higher twist contributions which, however, are power suppressed by the positive powers of x at small x [41] in striking contrast to the eikonal picture. Also, in the BFKL amplitude the multiple gluon ladder exchanges leading to higher twist contributions are correlated, whereas they are not correlated in the eikonal picture.
The two considered pictures of multiple gluon exchange differ also in the k T shapes of the gluon transverse momentum distribution (TMD). The GBW model leads to a narrow, quasi-collinear gluon k T distribution with the width scale given by the saturation scale, that is O(1 GeV), and the BFKL evolution generates a wide, power-like transverse momentum distribution with the asymptotic (at a very small x) positive anomalous dimension of 1/2, leading to ∼ 1/k T behavior of the gluon TMD F(x, k 2 T ). Since both the initial parton k T and the higher twist effects influence the Drell-Yan structure functions, it is desirable to analyze and disentangle them. The two considered models are particularly useful for this purpose as the GBW model introduces sizable higher twist effects and very small gluon k T , whereas the BFKL exchange generates gluons with large k T but it leads to very small higher twist corrections at small x.
The color dipole description of the Drell-Yan process incorporates small x resummation effects and a multiple scattering resummation (within a model). In particular the Drell-Yan description with BFKL amplitudes may be related to an analysis of small x effects in the DY scattering performed in Ref. [42]. The Drell-Yan process, however, receives large perturbative QCD corrections also in the limit of q T M and from soft gluon radiation near the partonic threshold energy. They are not included in the standard form of the dipole models. In particular the corrections coming from the small q T region of the cross section introduce at all orders n of the perturbative expansion terms enhanced by double logarithms ∼ α n s log 2n−1 (q 2 T /M 2 ), and also subleading logarithmic corrections [43]. Furthermore, it was shown by Collins, Sterman and Soper (CSS) [43] that in the relevant small q T region these corrections may be resummed or parameterized by a universal nonperturbative transverse momentum dependent parton distribution at very small q T . It was proven that the contribution of the resummed corrections of this type cancel after q T integrations [44], but it is essential for the correct description of the Drell-Yan q T -dependent cross section at small q T . Recently the problem of joint resummation of small x effects and the transverse momentum logarithms was addressed [45,46] providing important results for analyses of the q T dependent DY distribution at small x. The scheme for the joint resummations of transverse momentum logarithms and threshold correction is also available (see e.g. [47][48][49]) but it is expected to have less impact on the small x cross section. To summarize, the dipole model and BFKL predictions for the q T dependent DY cross sections need an improvement by the CSS resummation but the effects of this resummation cancel in the q T integrated cross section.
The paper is organized as follows. In section 2 one can find basic definitions of the Drell-Yan structure functions W j and of two models, the BFKL exchange, and the GBW saturation model. In the next section the procedure of the twist decomposition is discussed. The numerical predictions can be found in section 4. These predictions are presented in terms of the dimensionless structure functions A i commonly used for data presentations. Section 5 contains definitions of the structure functionsW i integrated over the lepton pair transverse momenta. It also contains a discussion of their twist decomposition which requires some attention due to apparent singularities. In section 6 numerical results are presented in terms of theÃ i and invariant λ i structure functions. The conclusions are given in section 7.

Structure functions in Drell-Yan processes
The DY helicity structure functions W j are defined through the formula for the differential cross section [50,51] where M is the lepton pair invariant mass, q T -the transverse momentum of the virtual photon and x F -its Feynman parameter. (θ, φ) are the polar and azimuthal angles of the lepton momentum vector in the dilepton c.m.s. frame. The frame orientation is not unique, and the most common frame choices are the Collins-Soper frame [52] and the Gottfried-Jackson frame [53]. In this paper we apply the description of the forward Drell-Yan process in the color dipole formulation. This formulation [25] assumes the dominance of the Compton-like partonic channel in which the fast collinear quark q, coming from one of the protons scatters off the color field of the other proton by single or multiple virtual gluon g * exchanges. At the lowest order the perturbative partonic channel for the forward Drell-Yan process with an intermediate γ * is qg * → q γ * → q l − l + , where l − and l + denote the produced leptons. The detailed description of the kinematics, and the relevant diagrams may be found in Ref. [12]. Hence, in the k T factorization approach within the color dipole picture defined in [25,27] one can show that in the Gottfried-Jackson frame the DY structure functions may be expressed as, where z is the longitudinal momentum fraction of the initial state quark taken by the virtual photon, ℘(x F /z) is a collinear parton distribution function andσ(s) is a color dipole -proton cross section in the Mellin representation. The parameter Q 0 is the Mellin transform scale. The leptonic impact factorsΦ i calculated in [12] can be found in Appendix A.1. The integration contour in the complex plane is taken as C = (−1/2−i∞, −1/2 +i∞). We consider two models for the description of the color dipole cross section: 1. The GBW model [15] in the form: where Q 0 is the saturation scale, Q 2 0 (x) = (x 0 /x) λ GeV 2 . The Mellin transform of the GBW dipole cross section w.r.t. the ρ 2 is equal toσ GBW (s) = −σ 0 Γ(s). The parameter values of the original GBW model [15] are applied: λ = 0.288, σ 0 = 23.03 mb, and the value ofx 0 = 6.08 · 10 −4 is chosen to be two times larger than the original GBW value x 0 = 3.04 · 10 −4 obtained from the description of the DIS data [15]. This modification was introduced because for the DIS data description the x variable in the dipole cross section was set to the threshold value of the gluon x g , that is to x = Q 2 /W 2 (where W in the proton-γ * collision energy), while in our DY description with q T -dependence the gluon x g is derived from the exact kinematics instead of using its threshold value, see Eq. 2.7. To be more specific, in the k T -factorization picture of the DIS at small x with the exact kinematics, the value of gluon x g in the γ * g → qq process depends on the virtual photon Q 2 and on the mass M X of the produced qq state: x g (Q 2 + M 2 X )/W 2 , (see e.g. [54] for the detailed discussion). Since the typical mass of the produced partonic qq state from the γ * fragmentation is M X ∼ Q, the approximate value of the gluon x g in the DIS is significantly larger than its threshold value Q 2 /W 2 and may be approximately estimated as x g 2 Q 2 /W 2 = 2x. So, if one treats the GBW dipole cross section as a function of true gluon x g instead of its threshold limit x, as we do for the q T -dependent DY cross section, it is necessary to rescale the model parameter x 0 to approximately 2x 0 . Then at given values of the observed parameters Q 2 and W 2 , the dipole cross section expressed through x g 2x and the rescaled parameter 2x 0 is the same as the dipole cross section expressed through x and x 0 . Our default choice for the figures in the paper is the parameterx 0 , but as the described treatment of the gluon kinematics is only approximate, we shall also explicitly display the sensitivity of selected observables to the choice of between the original GBW value x 0 and the rescaled valuex 0 of the dipole cross section.
2. The BFKL dipole cross section based on the solution of the leading-order (LO) BFKL equation [22] with the GBW input at a chosen value of x: where χ(s) is the LO BFKL characteristic function expressed through the digamma function ψ. The cross section parameter σ 0 = 2πR 2 p where R p is an effective radius of the proton which emerges after integration of an imaginary part of the forward dipole-nucleon scattering amplitude over the impact parameter b. The rapidity evolution length Y is given by where the value of gluon x g follows from the kinematics of the forward Drell-Yan process in the qg * → qγ * channel: with S denoting the invariant mass squared of the pair of colliding proton beams. In equation (2.4) we adopted the eikonal form of the initial condition for the BFKL evolution, coming the Golec-Biernat-Wüsthoff (GBW) model [15]: The model parameters were set by the fit to the DIS data [41] and readᾱ s = 0.087, Q 0 = 0.51 GeV, σ 0 = 17.04 mb. In our analysis we shall perform the twist decomposition of the DY structure functions obtained with the two model cross sections. Prior to that however, we test the reliability of the description by comparing the predictions for the forward Drell-Yan cross sections integrated over the lepton angles to the LHCb data. The results are shown in Fig. 1 as functions of the DY pair mass M for both dipole cross section models. In the cross section calculations the kinematical cuts were taken from [55] with an approximate treatment of the lepton transverse momentum (to be precise -when we impose the LHCb cuts on the lepton transverse momenta, q T is neglected in the transverse momentum balance of the leptons. This has a negligible effect on the results for M > 10 GeV). Note that in the region of the Z 0 mass, at M 90 GeV the Drell-Yan cross section is dominated by the Z 0 boson production that is not included in our analysis, hence the data point at M = 90 GeV is not well described. Except of this region, and for M > 10 GeV where the approximate treatment of the cuts may be neglected, the GBW agrees well with the data. The theoretical uncertainty of the GBW predictions due to the choice of the x 0 parameter of the dipole cross section is indicated as the vertical error bars, which are however within the GBW point size for all the points. The upper values correspond tox 0 and the lower ones to the standard GBW x 0 . As seen from the figure, the sensitivity to the choice of x 0 is found to be small. The M shape obtained with the BFKL model is consistent with the data but the overall normalization is slightly overestimated. The overall normalization, however, cancels in the analysis of the relative twist content so the small discrepancy of this parameter does not prohibit using this model in the twist analysis. Note that in the q T integrated cross sections the effects of the CSS resummation cancel [44].

Twist expansion for the helicity structure functions
The forward Drell-Yan helicity structure functions (2.2) can be written as The twist analysis of the forward DY cross section assuming the GBW dipole cross section was performed analytically in the preceding paper [12]. Below we perform an analogous twist decomposition using the BFKL model of the color dipole cross section using the method proposed in [41]. Hence, in order to perform the twist decomposition, we close the contour C with a left semicircle without changing the value of the integral. The integral over the closed contour is proportional to the sum of residues at the enclosed singularities. Hence, we express this integral as a sum of integrals around the singularities, which in this case are at the negative integers. The singularity at s = −n is identified with the twist-2n contribution to the amplitude. Therefore the cross section may be decomposed in the following way: where 2n corresponds to twist-2n term and singularity in s = −n. If the BFKL dipole cross section is assumed then essential singularities appear at s = −n. A procedure to evaluate the corresponding residues was described in [41]. The essential singularities are enclosed by circles with the radius → 0, and the angle θ parameterizes the position on the circle corresponding to the singularity at s = −n. where and is a regular function of in the limit of → 0 and t = log(M 2 /Q 2 0 ). Note that the terms which generate the essential singularities of the BFKL cross section coming from the exponentiated poles ∼ 1/ of the BFKL characteristic function were explicitly isolated. The coefficients h (2n) j can be expanded into an infinite series in , where the arguments q T , z, and Y of the series coefficients a (2n)j m are suppressed After substitution of (3.7) into (3.4) and integration over the angle θ one gets, what combined with (3.1) gives where I m is the modified Bessel function of the first kind. The first few coefficients a   [56].
There are several possible definitions of the Drell-Yan structure functions which are used for data presentation. The dimensionless structure functions A i [52] that can be directly related to coefficients of the lepton angular distribution in the DY pair center of mass frame: are ratios of the structure functions W j and W tot = W T + W L /2. To assess higher twists effects we calculate also exact (sum of all twists) structure functions by evaluation of integral (2.2) numerically.
In Fig. 2 we show a comparison between the exact BFKL results, the exact GBW and the twist 2 GBW components for A 0 , A 1 and A 2 . We do not show the separated BFKL twist 2 results, as they cannot be distinguished from the exact results, similarly to the case of the DIS analysis presented in [41]. The suppression of the higher twist contribution in the BFKL approach is explicitly illustrated in Fig. 5. It clearly follows from Fig. 2 that there is a substantial difference between the exact GBW and BFKL predictions. The difference between the models predictions follows mostly from the fact that the BFKL approach takes into account the transverse momentum parton distribution which is strongly limited in the GBW case (exponentially dumped). On the other hand the GBW model predicts sizable contributions from the higher twist terms in contrast to the BFKL expectations dominated by the leading twist term. Therefore, both models may serve as good benchmarks for the competition between the transverse momentum distributions and higher twist effects. An important observable in the analysis of subtle QCD effects beyond the leading twist collinear approximation, is the Lam-Tung observable A LT = A 0 −A 2 . In the collinear QCD framework A LT = 0 up to the next-to-next-to leading order. Hence A LT is considered to be a good probe of the higher twist effects and parton transverse momentum effects which do not compete here with the leading twist collinear contributions. In the left pannel of Fig. 3 the Lam-Tung observable A LT = A 0 −A 2 is shown as a function of the transverse momentum of the lepton pair at M 2 = 20 GeV 2 . In the GBW model the twist 2 contribution is consistent with A LT = 0, that is for the leading twist GBW the Lam-Tung relation is preserved. This follows explicitly from the analytic expressions for the twist expansion of the forward DY structure functions. It is expected as the color dipole description of the forward DY process is based on the partonic diagrams with the topology of the NLO contribution, and the GBW model leads to an almost collinear gluon distribution. At q T 5 GeV the deviations of A LT from zero appear in the GBW model, that are driven by the twist 4 term between q T 3 GeV and q T 5 GeV. Below q T 3 GeV twist 4 contribution is not sufficient and even the higher twist contributions become relevant. Within the exact GBW model, the Lam-Tung relation is broken at the level of 0.1 -0.2 for q T < 3 GeV due to the higher twist effects. Hence, within the GBW model of the color dipole cross section, the higher twist effects are clearly visible at lower values of the lepton pair transverse momentum.
The pattern following from the BFKL scattering amplitudes is different. Significant violation of the Lam-Tung relation occurs at all of the probed q T range. Recall that the BFKL amplitudes lead to negligible higher twist contributions. Hence, already the leading twist contribution of the BFKL model strongly violates the Lam-Tung relation. This breaking effect is a consequence of the wide transverse momentum distribution of the virtual gluons coming from the BFKL evolution (the more detailed study of the transverse momentum effects in the Lam-Tung relation breaking, taking into account also the g * g * channel, was performed in Ref. [33]). In fact, the gluon transverse momentum effects in the BFKL scattering amplitude lead to stronger breaking of the Lam-Tung relation than the higher twist effects in the GBW model down to q T = 2 GeV. At lower q T , however, the higher twist effects from the GBW and the gluon transverse momentum effects from the BFKL are similar.
In the right panel of Fig. 3 we show the Lam-Tung observable in both models for a very low mass of the DY pair, M 2 = 5 GeV 2 . At this mass we find a similar pattern to the case of M 2 = 20 GeV 2 , with slightly enhanced both the higher twist and the gluon transverse momentum effects. It follows from Fig. 3 that for the Lam-Tung observable at low masses, lowering the mass does not lead to relative amplification of the higher twist corrections w.r.t. the parton transverse momentum effects. Hence, even at lower DY pair masses, even at low q T , and very large energies and in the forward kinematics, the effects of parton transverse momentum and of the higher twists are expected to have a similar contribution to the Lam-Tung relation breaking. At higher masses and at larger q T the higher twist terms become small, but the possible effects of parton k T may still stay sizable, see e.g. Ref. [33]). These results indicate that the observation of the higher twist contributions in A LT requires a good control of the parton transverse momentum effects.

Twist expansion of the integrated helicity structure functions
It may be useful to study experimentally the forward Drell-Yan structure functions integrated over q T . They are defined in the following way, and their Mellin representation takes the form: where: The rapidity evolution length was defined in (2.6) and for x g we use threshold (q T → 0) value of (2.7), x g = M 2 /(Sx F ) . The expressions for the Mellin transformed impact factor H i (s) were derived in [11,12] and are listed in Appendix A.2. The formulae for the twist decomposition of the q T -integrated forward DY structure functions assuming the GBW dipole cross section were given in [12]. Below we derive the twist expansions for the q Tintegrated DY structure functions with the BFKL exchange using the procedure described in Sec. 3. Hence we insert the explicit Mellin representation of the BFKL cross section into the above formula for the integrated structure functions, (5.4) wheret = ln(4M 2 /Q 2 0 ). Following the steps described in the previous section one gets the expression for the twist-2n component of the structure functions: where the coefficientsh j in the q T -integrated case discussed in section 3. For j = L, T T, LT , after the integration over θ, the leading twist terms are given by,W are the expansion coefficients of the functions where the dependence ofã (2n)j m on the variables z and Y was suppressed. The first two leading coefficientsã  . This follows from the fact that the most leading coefficients in the expansion correspond to the double logarithmic approximation of the BFKL exchange which coincides with the collinear approximation results. As it is expected, the Lam-Tung relation is broken by the non-leading coefficientsã (2)j m , m > 0, that correspond to the BFKL effects beyond the double logarithmic limit. Also interesting is to note that the leading twist 2 coefficient a (2)LT 0 = 0 for the structure functionW LT . As the analytical expression for H LT is not known, the expansion coefficientã A more sophisticated procedure is necessary to obtain the twist 2 component ofW T and the twist components of all the structure functions beyond twist 2, because a logarithmic divergence of the form ln(1−z) occurs in the integrals corresponding to the twist components at z → 1. A treatment of such apparent singularities was developed in [11] and we follow that procedure -see Appendix C for the details.

Results for the integrated helicity structure functions
The results in this section are presented using the dimensionless coefficientsÃ i : A 0 ( 4 ) -A 2 ( 4 ) B F K L whereW tot =W T +W L /2. Additionally, we introduce coefficients λ 1 , λ 2 and λ 3 which are invariant with respect to rotations in the X − Z plane in the lepton center of mass frame [57]. They read In Fig. 4 the Lam-Tung observableÃ LT =Ã 0 −Ã 2 for the q T -integrated distributions is presented for the exact values of the GBW and BFKL modelsÃ (E) LT and for the twist 4 components of the modelsÃ (4) LT . Recall that in GBW model the twist 4 component provides the leading contribution to the Lam-Tung observable. It follows from the figure that in the GBW model the higher twist contributions inÃ LT become visible below M 2 = 100 GeV 2 . However, the BFKL model predicts sizable violation of the Lam-Tung relation in the integrated DY structure functions in the whole plotted range of the masses (at the leading twist). The origin of the violation may be traced back to the strong parton transverse momentum effects in the BFKL approach. These effects are stronger than the GBW higher twist effects down to M 2 = 4 GeV 2 . Only below that threshold the higher twists prevail. This supports our previous conclusion that disentanglement of the higher twist effects from the parton transverse momentum effects requires a careful analysis and acquiring a good understanding of the parton k T distribution. In BFKL twist 4 is negligible soÃ (4) LT is very close to zero.
In Fig. 5 we illustrate the higher twist contributions to a selected structure functioñ W L using ratios: R (R (4) W L ) is the relative negative contribution of the twist n components with n > 2 (n > 4). The figure shows that the BFKL result is dominated by the leading twist component, and the higher twist components enter at the level of 10 −3 of the dominant twist 2 term for the whole range of M 2 > 1 GeV 2 . In the GBW model the twist 4 correction toW L becomes relevant below M 2 100 GeV 2 , and below M 2 30 GeV 2 all the twist components should be taken into account. Fig. 6 shows comparison between the predictions of the GBW and BFKL models for W i structure functions (the exact values). ForW L ,W T , andW LT the largest differences  between the predictions appear at higher values of M 2 . In this kinematical range the higher twist contributions are suppressed and the shape of the curves is determined by the transverse momentum effects. However, theW T T structure function exhibits the opposite behavior. The largest difference between the BFKL and GBW models appears at the lower values of M 2 . One concludes thatW T T is the particularly sensitive structure function to the higher twist effects. The same pattern is found in the dimensionless integrated structure functionsÃ i . On the left hand side of Fig. 7 we show a comparison between predictions ofÃ i obtained with the GBW and BFKL models. Comparing this plot with the right panel of Fig. 7 one concludes that the higher twist contribution in the GBW model starts being visible already for M 2 300 GeV 2 . Note that the coefficientÃ 2 , that is directly related toW T T , is particularly sensitive to the effects beyond the leading twist collinear approximation. At the level of twist 2 the difference between the BFKL model and the GBW model results may be treated as an approximate measure of the gluon k T effects that are sizable in the BFKL approach and almost negligible in the GBW model. Comparing the higher twist content A 2 in the GBW results (see Fig. 7, the right panel) with the spread between the GBW and BFKL predictions forÃ 2 , one concludes that for M 2 30 GeV 2 the higher twist effects that follow from the GBW model are larger than the gluon k T effects following from the BFKL evolution. ThereforeÃ 2 at low M 2 should be particularly useful observable for an experimental discrimination between the models and for constraining the higher twist contributions.
As discussed in Sec. 2, we introduced a modified value of the original x 0 parameter x 0 = 2x 0 in the model of the GBW dipole cross section. In order to display the sensitivity of the GBW predictions to this parameter variation we compare in Fig. 8 the twist 2 and the exact estimates for the q T -integrated angular coefficients A i obtained with the original GBW value x 0 and the modified valuex 0 . The resulting theoretical uncertainty bands are found to be rather narrow.
In Fig. 9 we show also the results in terms of invariant coefficients λ i , that are valid in the frames with the Y -axis transverse to the beam -DY pair plane [57]. Particularly strong higher twist effects in the GBW approach are found in λ 2 for M 2 100 GeV 2 (see the right panel). Note also since λ φ 1 and 1 − λ 1 = λ 3 /(1 − λ φ ) we have 1 − λ 1 ≈ λ 3 which can be seen on the plot. At the leading twist the GBW model yields λ 1 = 1 and λ 3 = 0 that implies that in this approximation the Lam-Tung relationW L − 2W T T = 0 is satisfied, as expected.
To sum up, the forward Drell-Yan angular coefficientsÃ i in particularÃ 2 for M 2 100 GeV 2 are sensitive to the higher twist effects and the measurements at the LHC can be used to constrain the higher twist contributions.

Conclusions
A study of the forward Drell-Yan cross sections in pp collisions was performed in the k T factorization framework in the color dipole realizations. We assumed the LHC energy √ S = 14 TeV and the forward kinematics -the Feynman x of the DY pair, x F = 0.05. There are two important ingredients which influence the Drell-Yan process beyond the collinear approximation: the higher twist contributions and the transverse momenta of partons. At lower values of the boson invariant mass M 2 and at a very small x, both the effects are significant and of similar magnitude and may compete with each other. The disentanglement of those distinct contributions requires a careful analysis of the Drell-Yan structure functions over a broader range of kinematical parameters. The higher twist corrections are strongly suppressed at large process scales given by the Drell-Yan pair mass M and/or the transverse momentum p T . Thus in that kinematical region of M 10 GeV and/or p T > 10 GeV the higher twist effects may be safely neglected and the parton k T effects may be isolated. The results may be then used to provide necessary input to fit the parton transverse momentum distributions. A good and encouraging example is provided by a recent ATLAS measurement of the Lam-Tung relation breaking at the Z 0 boson peak [9] which may serve as an excellent test for the gluon transverse momentum distributions [33]. With good understanding of the parton k T achieved, the higher twist effects may be probed in detail. Hence the observation of the higher twist effects should be possible at lower values of M 2 , desirably at M 2 < 100 GeV 2 , and very small x < 10 −5 of one of the partons, and at a low p T . For optimal quality of the higher twist determination, the measurements of all the Drell-Yan structure functions should be performed, in a wide range of kinematical parameters. The W T T structure function and the Lam-Tung observable A LT = A 0 − A 2 exhibit particularly high sensitivity to effects beyond the leading twist collinear approximation.

B Coefficients of twist expansion
Here we present a (2n)i m that were obtained by the series expansion of (3.7). Coefficients for m ≥ 1 are not presented as they are lengthy.
(B.1) Notice that the twist 4 terms are suppressed by the e −2Yᾱs factor. In fact, it turns out that contribution from the higher twist terms in the BFKL model is negligible as is seen e.g. from Fig. 5.

C Twist expansion forW T in BFKL model
The twist expansion of the structure functionW T in the BFKL model requires a careful treatment of the z integration at the singular z → 1 limit. It is performed following the procedure proposed in Ref. [11]. Hence we rewrite (5.5) for W (2) T in the following form: After the integration of the first term, ℘(x F ), one has: