Saturation effects in exclusive ρT,L meson electroproduction

We use recent results for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \gamma_L^{*}\to {\rho_L} $\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \gamma_T^{*}\to {\rho_T} $\end{document} impact factors, computed in the impact parameter representation within the collinear factorization scheme, to get predictions for the polarized cross-sections σT and σL of the diffractive leptoproduction of the ρ meson at high energy. In this approach the helicity amplitude is a convolution of the scattering amplitude of a color dipole with a target, together with the virtual gamma wave function and with the first moments of the ρ meson wave function (in the transverse momentum space), given by the distribution amplitudes up to twist 3 for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \gamma_T^{*}\to {\rho_T} $\end{document} impact factor and up to twist 2 for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \gamma_L^{*}\to {\rho_L} $\end{document} impact factor. Combining these results with recent dipole models fitted to DIS data, which include saturation effects, we show that the predictions are in good agreement with HERA data for photon virtuality (Q2) larger than typically 5 GeV2, without free parameters and with a weak dependence on the choice of the factorization scale, i.e. the shape of the DAs, for both longitudinally and transversely polarized ρ meson. For lower values of Q2, the inclusion of saturation effects is not enough to provide a good description of HERA data. We believe that it is a signal of a need for higher twist contributions in the ρ meson DAs. We also analyze the radial distributions of dipoles between the initial γ* and the final ρ meson states.


JHEP11(2013)062
and ρ L(T ) −meson states 19 6.1 The radial distribution of the γ * L → ρ L transition 24 6.2 The radial distribution of the γ * T → ρ T transition 28 6.3 Comparison with the radial distributions obtained from models of the ρ meson wave function. 32 1 Introduction In the high energy limit, exclusive processes and more particularly the diffractive leptoproduction of vector mesons, provide a nice probe to study the hadronic properties. From the experimental side, data have been extracted in a wide range of center-of-mass energies, from a few GeV at JLab to hundreds of GeV at the HERA collider. The kinematical range which is at the heart of the present paper is the large energy in the center-of-mass of γ * p, denoted W , for which the HERA collaborations H1 and ZEUS measured ρ-meson electroproduction starting from the early period of HERA activity [1,2] till now, with an increasing precision leading to a complete analysis [3,4] of spin density matrix elements, JHEP11(2013)062 polarized and total cross-sections describing the hard exclusive productions of the ρ and the φ vector mesons V in the process These matrix elements and polarized cross-sections can be expressed in terms of helicity amplitudes T λ V λγ (λ γ , λ V : polarizations of the virtual photon and of the vector meson). The ZEUS collaboration [3] has provided data for different photon virtualities Q 2 , i.e. for 2 < Q 2 < 160 GeV 2 , 32 < W < 180 GeV ( |t| < 1 GeV 2 ), while the H1 collaboration [4] has analyzed data in the range 2.5 < Q 2 < 60 GeV 2 , 35 < W < 180 GeV (|t| < 3 GeV 2 ) . The high virtuality of the exchanged photon allows the factorization of the amplitude into a hard sub-process described within the perturbative QCD approach and suitably defined hadronic objects, such as the dipole-nucleon scattering amplitude or the vector meson wave functions and distribution amplitudes (DAs) [5][6][7]. HERA data are thus interesting observables to test the properties of these non-perturbative objects such as the saturation dynamics of the nucleon or the transverse momentum dependence of the vector meson wave functions.
On the theoretical side, three main approaches have been developed. The first two, a k T -factorization approach and a dipole approach, are applicable at high energy, W ≫ Q ≫ Λ QCD . They are both related to a Regge inspired k T -factorization scheme [8][9][10][11][12][13][14], which basically writes the scattering amplitude in terms of two impact factors: one, in our case, for the γ * − ρ transition and the other one for the nucleon to nucleon transition, with, at leading order, a two "Reggeized" gluon exchange in the t-channel. The Balitsky-Fadin-Kuraev-Lipatov (BFKL) evolution, known at leading order (LLx) [15][16][17][18] and next-to-leading (NLLx) order [19][20][21][22], can then be applied to account for a specific large energy QCD resummation. The dipole approach is based on the formulation of similar ideas although not in k T but in transverse coordinate space [23,24]; this scheme is especially suitable to account for nonlinear evolution and gluon saturation effects. The third approach, valid down to W ∼ Q, was initiated in [25] and [26]. It is based on the collinear QCD factorization scheme [27,28]; the amplitude is given as a convolution of quark or gluon generalized parton distributions (GPDs) in the nucleon, the ρ-meson DA, and a perturbatively calculable hard scattering amplitude. GPD evolution equations resum the collinear quark and gluon effects. The DAs are also subject to specific QCD evolution equations [29][30][31].
Though the collinear factorization approach allows us to calculate perturbative corrections to the leading twist longitudinal amplitude (see [32] for NLO), when dealing with transversely polarized vector mesons related to higher twist contributions, one faces endpoint singularity problems. Consequently, this does not allow us to study polarization effects in diffractive ρ-meson electroproduction in a model-independent way within the collinear factorization approach. An improved collinear approximation scheme [33] has been proposed, which allows us to overcome end-point singularity problems, and which has been applied to ρ-electroproduction [34][35][36][37].
In this study, we consider polarization effects for the reaction (1.1) in the high energy region, s = W 2 ≫ Q 2 ≫ Λ 2 QCD , working within the k T -factorization approach, where JHEP11(2013)062 the helicity amplitudes can be expressed 1 as the convolution of the γ * → ρ impact factor Φ γ * (λγ )→ρ(λρ) (k 2 , Q 2 ) with the unintegrated gluon density F(x, k 2 ) density, which at the Born order is simply related to the nucleon impact factor Φ N →N (k 2 , Q 2 ), where k is the transverse momentum of the t−channel exchanged gluons. An important outcome of the k T −factorization is that for the transverse amplitude the end-point singularities are naturally regularized by the transverse momenta of the t−channel gluons [38][39][40]. At large photon virtuality, the γ * − ρ impact factors can be calculated in a model-independent way using QCD twist expansion in the region k 2 ≫ Λ 2 QCD . Such a calculation involves the ρ-meson DAs as nonperturbative inputs. The calculation of the impact factors for Φ γ * L →ρ L , Φ γ * T →ρ L is standard at the twist-2 level [41], the next term of the expansion being of twist 4, while Φ γ * T →ρ T was only recently computed [39,40] (for the forward case t = t min ), up to twist-3, including two-and three-parton correlators, which contribute here on an equal footing.
One remark is in order about the above Fock expansion. We want to emphasize the fact that the twist expansion of the operators is not the perturbative expansion of the hard part of the amplitude. Indeed, the quark antiquark gluon intermediate state contribution does not contain a power suppressed hard part in g compared to the quark anti-quark contribution. Thus, despite the fact that we include an additional gluon in the Fock state, this contribution is still a leading order one, since the additional power of g is absorbed inside the non-perturbative correlators. The dipole which is involved in the scattering process is a pure lowest-order color-anticolor state, undressed by gluons. Thus the question of the magnitude of the correction beyond leading order is beyond the scope of the present paper.
In a previous study [42], we used the results [40,41] for the Φ γ * L →ρ L and Φ γ * T →ρ T impact factors and a phenomenological model [43] for the proton impact factor. It was pointed out that the region k 2 ≫ Λ 2 QCD gives the dominant contribution to the helicity amplitudes while the soft gluon (k 2 < 1 GeV 2 ) contribution cannot be neglected. The soft gluon contribution to the amplitudes involves the interaction of large size color dipole configuration |r| (|r| ≡ 1/|k|) in the fluctuations of the probe and saturation effects could then play an important role.
Following this idea, we have shown in ref. [44] that the helicity amplitudes, expressed in impact parameter space and then computed in the collinear factorization scheme, factorize into the dipole cross-section and the wave functions of the virtual photon combined with the first moments of the ρ meson wave functions parameterized by the DAs, given by the twist expansion up to twist 3 for the production of a ρ T and to twist 2 when producing a ρ L . The results of ref. [44] link the k T −factorization approach, in particular the results of [40], with the calculations performed in refs. [45][46][47][48][49] within the dipole approach. The main difference between our present approach and the one of refs. [45][46][47][48][49] is that instead of using the light-cone wave functions φ(z, r) which in practice need to be modeled, the amplitude of our approach involves the DAs which parameterize the first moments of the wave functions.

JHEP11(2013)062
Our main point here is that one can assume the dominant physical mechanism for production of both longitudinal and transversely polarized mesons to be the scattering of small transverse-size quark-antiquark and quark-antiquark-gluon colorless states on the target. This picture is justified at high center-of-mass energy of the γ * p system, since in this kinematics it is assured that the fluctuation of the γ * into a qq pair and then its interaction with the nucleon are clearly separated in time. We here consider a frame such that the γ * is almost at rest, while the proton is strongly boosted. Our inclusion of both qq and qqg Fock states is therefore not related to the higher order Fock states produced by the high-energy dynamics, which are included in the proton wave function (they are thus part of the dipole-nucleon cross-section). Indeed, the additional gluon which we include is not soft with respect to its quark or antiquark emitters. Therefore we here consider that the fluctuation of the γ * into a qq and qqg arises much before the interaction of these Fock states with the proton. The validity of this dipole factorization is thus only based on the large value of center-of-mass energy of the γ * p system. The potential issues due to the small values of k, which could lead to a diffusion into the non-perturbative domain, is known to be controlled by the saturation effects, which we take as granted, provided by the models which we use which detailed mechanism is not the subject of the present paper.
This picture allows to calculate corresponding helicity amplitudes in a model-independent way, using the natural light-cone QCD language -twist-2 and twist-3 DAs. There, a comment is in order. Indeed, within k T -factorization, the description of the impact factor for produced meson based on the QCD collinear approach requires a modification of the usual twist counting due to the off-shellness of the t-channel partons. For that, we assume that the only large scale is provided by Q 2 . Thus, the virtualities of t−channel gluons involved in the k T -factorization approach are considered to be large, being typically of the order of the photon virtuality Q 2 . Therefore, when here we say "up to twist 3" we only mean twist counting from the point of view of the collinear factorization of the produced ρ-meson, and not of the whole amplitude of γ * p → ρp . This means that the twist counting considered here should be understood only based on non-local operators involved in the definition of meson DAs.
Note that in [38] a related approach, based on the idea that color dipoles of small sizes (r ∼ 1/Q) interact with the target, was used to expand the dipole scattering amplitude and the ρ−meson wave function around small r. Our approach is different since it is based only on the twist expansion of the ρ meson, governed by corrections around the dominant light-cone direction provided by this produced meson. Therefore our framework allows for the inclusion of the whole r−dependence of the dipole-target scattering amplitude while the ρ meson wave function is expanded up to a given twist and parameterized by DAs. Saturation effects, which are beyond the small r approximation, can thus be taken into account in our study. This paper is organized as follows. In section 2, we introduce the impact factor representation of the helicity amplitudes as well as the kinematics of the process. In section 3, we first recall some results for the γ * L → ρ L and γ * T → ρ T impact factors computed in momentum space respectively up to twist 2 [41] and twist 3 [40] accuracies using the collinear approximation to parameterize the soft part associated to the production of the ρ meson JHEP11(2013)062 by distribution amplitudes (DAs), calculated in ref. [50,51]. We recall then the expression of these impact factors in the impact parameter space according to the results of ref. [44] allowing to decouple the dipole-nucleon scattering amplitude from the amplitude of production of dipoles in the initial (γ * (λ γ )) and final (ρ(λ ρ )) states. We terminate section 3 by expressing helicity amplitudes and polarized cross-sections in terms of the dipole crosssection. In section 4, we give a brief review of the main properties of the models for the dipole cross-section [52][53][54] or the proton impact factor [43] that we use in our study. We compare our predictions with the data of HERA [3,4] in section 5 and we obtain a good agreement. In this context we discuss the role of higher twist corrections for small Q 2 values. Finally, we analyze the radial distribution of dipole intermediate states involved between the virtual photon and the ρ meson, and discuss the role of the saturation models on the specific example of the Golec-Biernat and Wüsthoff saturation model [52]. We also compare our radial distribution with the overlap of the γ * and ρ meson wave functions obtained in the approach of the dipole models [46,47,55], where the ρ meson wave function is modeled and the parameters are fitted to HERA data.
2 Helicity amplitudes of the hard ρ meson leptoproduction in the impact factor representation In the impact factor representation at the Born order, the amplitude of the exclusive process as illustrated in figure 1. The γ * (λ γ ) → ρ(λ ρ ) impact factor Φ γ * (λγ )→ρ(λρ) is defined through the discontinuity of the S matrix element for γ * (λ γ ; q)g(k) → g(k − ∆)ρ(λ ρ ; p ρ ) as where κ = (k + q) 2 . In eqs. (2.1) and (2.2) the momenta q and p ρ are parameterized via Sudakov decompositions, in terms of two light-like vectors p 1 and p 2 such that 2 p 1 .p 2 = s, as where Q 2 = −q 2 ≫ Λ 2 QCD is the virtuality of the photon being the large scale which justifies the use of perturbation theory, and m ρ is the mass of the ρ meson. Here −t min denotes the minimal value of −t . The nucleon impact factor Φ N →N in eq. (2.1) cannot be computed within perturbation theory, and is related 2 at Born order to the unintegrated gluon density F(x, k). In the forward limit ∆ ⊥ = 0, the helicity amplitudes read,  Figure 1. Impact factor representation of the γ * N → ρ N scattering amplitude. The impact factors Φ γ * (λγ )→ρ(λρ) (k, Q, µ 2 F ) and the nucleon impact factor vanish at k → 0 or k → ∆, which guarantees the convergence of the integral in eq. (2.4) on the lower limit. 3 The computation of the γ * → ρ impact factor is performed within collinear factorization of QCD. The dominant contribution corresponds to the γ * L → ρ L transition (twist 2), while the other transitions are power suppressed. The γ * L → ρ L and γ * T → ρ L impact factors were computed a long time ago [41], while a consistent treatment of the twist-3 γ * T → ρ T impact factor has been performed only recently in ref. [40]. It is based on the light-cone collinear factorization (LCCF) beyond the leading twist, applied to the amplitudes γ * (λ γ )g(k) → g(k − ∆)ρ(λ ρ ), symbolically illustrated in figure 2. Each of these scattering amplitudes is the sum of the convolutions of a hard part (denoted by H and H µ for two-and three-parton contributions, respectively) that corresponds to the transition of the virtual photon into the constituents of the ρ meson and their interactions with off-shell gluons of the t channel, and a soft part (denoted by Φ and Φ µ ) describing the hadronization of the constituent partons into the ρ meson. The technique is to perform a Taylor expansion of the hard part around the dominant light-cone direction given by the ρ meson momentum. Up to twist 3, the Taylor expansion is truncated to the first order in ℓ ⊥ for the quark antiquark intermediate state contribution. For the quark antiquark gluon intermediate state contribution such expansion is not necessary. Thus, the Taylor expansion of the hard part associated to the quark antiquark contribution up to twist 3 reads

JHEP11(2013)062
where H Γ α is the hard sub-process projected on the Fierz structure Γ α (see section 3 below for details). In order to link the collinear factorization around the dominant light-cone direction with the dipole picture we express the hard part in terms of its transverse Fourier transform such that the Taylor expansion reads The term "iℓ · r" turns into a transverse derivative into the non-local correlators involved in the soft part after integration by part. The quark antiquark contribution to the impact factor reads where the subscript µ F stands for the factorization/renormalization scale, under which the transverse momenta of the partons are integrated. The Fourier transforms of the non-local operators are parameterized by the set of twist 2 and twist 3 DAs defined in appendix A.
Here we want to pinpoint the fact that the present extended collinear approximation, when expressed in the coordinate space, keeps the whole r ⊥ dependency of the hard part while it expands the non-local correlators around small (z ⊥ < 1/µ F ) relative coordinate between the partonic constituents of the ρ−meson.
The DAs needed to parameterize the ρ meson productions involved in the results of [40] for the impact factors with ∆ = 0 are: the twist 2 DA ϕ 1 (y; µ 2 F ) associated to the production of the longitudinal ρ meson, the two-parton twist 3 DAs ϕ T 1 (y; µ 2 F ), ϕ T A (y; µ 2 F ) and the three-parton twist 3 DAs B(y 1 , y 2 ; µ 2 F ) and D(y 1 , y 2 ; µ 2 F ) that parameterize the production of a transversely polarized ρ meson. 4 Following [42,44], we also use the combinations, We recall in appendix A some basics features of the chiral even twist 2 and twist 3 DAs present in this approach; more details can be found on the DAs in ref. [40]. We recall also in appendix B the dependence of the DAs on the collinear factorization scale µ F which is driven by the renormalization evolution equations given in ref. [50].

JHEP11(2013)062
3 Helicity amplitudes and polarized cross-sections 3.1 Impact factors γ * L,T → ρ L,T In the Sudakov basis, the longitudinal and transverse polarizations of the photon are 5 For t = t min the same parametrization will be used for the ρ-meson polarization with obvious substitutions Q 2 → −m 2 ρ and Q → m ρ . Let us introduce p and n, two light-cone vectors such that p ρ ≈ p at twist 3 and p · n = 1. The polarization of the out-going ρ meson is denoted by e * . For further use, we denote R * ⊥α = ε αλβδ e * λ ⊥ p β n δ . We use the following convention for the transverse euclidean polarization vectors e ± as in eq. (3.1) We have recently shown in [44] that these impact factors, expressed in the impact parameter space, read, and , (3.4) where we have respectively denoted in the case of two-parton exchange as y andȳ = 1 − y the longitudinal momentum fractions of the quark and the antiquark, and in the threeparton exchange y 1 ,ȳ 2 = 1 − y 2 and y g = y 2 − y 1 the longitudinal momentum fractions of the quark, the antiquark and the gluon. The transverse displacement vectors of the color dipole configurations are denoted as r for the interacting dipole (two-and three-parton Fock component) and r ′ for the spectator dipole (present only in the three-parton Fock component). Note that |r| in the case of the two-parton component is the transverse size of the quark anti-quark colorless pair. The functions denoted are the Fourier transforms in the transverse plane of the two-parton component hard parts H and the three-parton component hard parts H µ (illustrated in figure 2), projected on the appropriate Fierz structures Γ µ . 5 In ref. [40] we took ǫ ± = ∓ i √ 2 (0, 1, ±i, 0), which we change here for consistency with the usual experimental conventions [59].

JHEP11(2013)062
The computations of the hard parts lead to the following generic expressions, where the function is the scattering amplitude of a quark-antiquark color dipole with the two t−channel gluons, putting apart the color factor T r(t a t b ) = δ ab /2, with a and b color indices and N c the number of colors. The functions ψ are respectively the amplitudes of production of a ρ meson from a quark-antiquark (quark-antiquark gluon) system produced far upstream the target in the fluctuation of the virtual photon. These functions are computed up to twist 3 in the collinear approximation in ref. [44] and they contain information about the relevant color dipole system that interacts with the target. The functions ψ can be expressed in terms of the virtual photon wave functions Ψ where h = ± 1 2 ,h = ± 1 2 denote respectively the helicities of the exchanged quark and anti-quark, and of the combinations of DAs of the ρ meson φ In the two-parton approximation, these functions φ ρ L (h,h) (y; µ 2 F ) and φ ρ T , (λρ) (h,h) (y, r; µ 2 F ) parameterize the moments of the wave functions of the ρ meson, i.e. the first terms of the Taylor expansion of the wave functions at small r. Note that in the approach of ref. [38] (h,h) (y, r; µ 2 F ) are replaced by the Taylor expansion for small r of the modeled wave functions of the vector mesons. Finally, the functions ψ (3.14)

JHEP11(2013)062
The function ψ with three-parton components reads where the function F γ * T describes the fluctuation of the transversely polarized photon into a quark-antiquark-gluon color singlet. The function F γ * T can be expressed in terms of the longitudinally polarized photon wave function 2Nc . Note, that in the large N C limit F γ * T simplifies,

JHEP11(2013)062
where we have identified the dipole cross-section as defined in ref. [56] σ(x, r) = N 2 c − 1 4 Note that we can separate the T 11 as the Wandzura-Wilczek (WW) [60] contribution and the genuine contribution, The formulas (3.21), (3.22) allow us to combine various models of the scattering amplitude of a dipole on a nucleon with the results obtained by twist expansion of the γ * → ρ impact factor. At t = t min the contribution to the cross-sections σ L and σ T are respectively coming from the helicity amplitudes T 00 and T 11 , The t−dependency is expected to be governed by non-perturbative effects of the nucleon, which can be phenomenologically parameterized by an exponential dependence of the differential cross-section This results in the polarized cross-sections The b(Q 2 ) slope has been measured by ZEUS and H1. We will use here quadratic fits of the b(Q 2 ) slope data of ref. [4] to determine the cross-section. In the following section, we will briefly present the dipole models we shall use to compare our predictions with HERA data. Note that the dipole cross-section models which we will use are determined from DIS data using the overlap of the virtual photon wave functions at leading order where only the quark antiquark intermediate Fock state is involved, while the higher Fock states, such as the quark antiquark gluon intermediate state, are neglected. This is consistent with the fact that our approach is leading in the strong coupling g, as mentioned in the introduction.

JHEP11(2013)062 4 Dipole models
In the dipole picture, the DIS cross-section reads [23,24] The low-x saturation dynamics of the nucleon target was first introduced in refs. [52,61] by Golec-Biernat and Wüsthoff (GBW) model to describe the inclusive and diffractive structure functions of DIS, which inspired many phenomenological descriptions of DIS HERA data [47,55,[62][63][64][65][66]. In this model the dipole cross-section readŝ and it involves three independent parameters {σ 0 , x 0 , λ p }. One can see from eq. (4.1) that since the wave functions are peaked at r ∼ 1 Q , the domain in which the saturation effects are significant is given by To make contact with photoproduction, it is customary [52] to make the following modification of the definition of the Bjorken variable x where m f is an effective quark mass which depends on the flavor f and of the model used to fit the data. This modification is adopted in the following.
In the GBW saturation model, light quark masses are taken to be 0.14 GeV in order to get a good fit of the photoproduction region. The inclusion of the charm contribution, with m c = 1.5 GeV has also been performed in ref. [52]. The normalization σ 0 results from the integration over the impact parameter b, assuming that the b−dependence of the dipole amplitude factorizes as is related to the density of gluon within the target, leading tô This model provides a good description of inclusive and diffractive structure functions for the values of the parameters presented in table 1.

JHEP11(2013)062
Fits The small x evolution of the dipole cross-section can be modeled as in the refs. [52], or computed numerically from the running coupling Balitsky-Kovchegov (rcBK) equation [67]. Indeed, the x−dependence is driven at small x by perturbative non-linear equations, the Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) equation [68][69][70][71][72][73] and the Balitsky-Kovchegov (BK) equation [74,75]. The solutions of the two equations are not significantly different and the BK equation is simpler to solve numerically than the JIMWLK equation. The LO-BK equation cannot be used to obtain the low x evolution as it predicts a growth of the saturation scale much faster than the growth expected from the analysis based on phenomenological models. It was shown in [53,54] that taking into account only the running coupling corrections of the evolution kernel allows to get the main higher order contributions, and to solve the discrepancy between the growths of the saturation scale.
A remark is now in order, in view of the use of the dipole for the exclusive non forward process, which is studied in this paper. We could in principle use a dipole model taking into account skewness effects. Indeed, the typical values of the p 2 components of the two t−channel gluons are rather asymmetric, one being much larger than the other one. This skewness dependence can be implemented along approaches of refs. [76,77]. However, the way how to include skewness effects, based on the Shuvaev transform, is still under discussion, see e.g. refs. [78,79], since it requires certain analytical properties of conformal moments of GPDs which are generally not satisfied. Such a treatment seems therefore to be rather considered as an powerful ansatz, to be confronted with data. On the other hand, in the context of saturation, the inclusion of skewness has only been performed using a Glauber-Mueller treatment [47,80]. To the best of our knowledge, the only available approach using the Balitsky Kovchegov equation has only been performed in ref. [81], and uses the Shuvaev ansatz in the initial condition at low Q 2 . We would like to comment on the fact that kinematically, the notion of skewness does not exist in the Balitsky-Kovchegov approach. One can however argue that the skewness effect is mainly important for the last stage (i.e. the coupling with the quark loop), and can be neglected in the fan-diagrams which underlies the BK equation. Anyway, due to the various limitation mentioned above, we will not take into account these skewness effects in the present study.
A numerical solution of the rcBK equation was obtained in refs. [82,83] based on initial conditions inspired by the GBW model [52] and the McLerran-Venugopalan (MV) model [84]. We will denote these numerical solutions for the dipole scattering amplitudes as the Albacete-Armesto-Milhano-Quiroga-Salgado (AAMQS)-model.
where β 0,n f = 11 − 2 3 n f , Λ n f is the QCD scale and C is one of the free parameters of the model. As usual, the scales Λ n f are determined by the matching condition α s, and an experimental value of α s . In eq. (4.9), the scale Λ QCD is identified with Λ 3 .
The parameters are fitted to the experimental data for the inclusive structure function of DIS and give a good description for the data of the longitudinal structure function In the following parts we will focus mostly on the sets (a) and (b) using the GBW initial condition respectively denoted as AAMQSa and AAMQSb, as the effects due to different initial conditions do not involve significant changes in the numerical results of present study. We present in tables 2 and 3 values of the parameters of the fits obtained in ref. [83]. For completeness, we will also display predictions using the Gunion and Soper model (GS-model) [43]. This model was used in our first phenomenological study of ref. [42], and it assumes that the hadron impact factor takes the form where A and M are free parameters that correspond to the soft scales of the proton-proton impact factor. As it was discussed in [42], the ratios of helicity amplitudes are well described JHEP11(2013)062 Fits  for M ≃ 1 GeV and the result is not very sensitive to this parameter. Note that this impact factor indeed vanishes when k → 0 or ∆ − k → 0 in a minimal way. Such a model was the basis of the dipole approach of high energy scattering [85] and used successfully for describing deep inelastic scattering (DIS) at small x [86]. The dipole-proton scattering amplitude in the two-gluon exchange approximation, computed in details in appendix C, resulting from the GS-model for the proton impact factor is:

Comparison with the HERA data
Our main results are the polarized cross-sections σ T and σ L of the process (1.1), that we compare with the data of the H1 collaboration [4]. In the following plots, experimental errors are taken to be the quadratic sum of statistical and systematical errors. The collinear factorization scale µ F only appears in the scattering amplitudes through the DAs and the coupling constants and unless specified we will assume that µ F depends on the virtuality Q 2 as The figures 3(a) and 3(b) show the full twist 3 predictions respectively for σ T and σ L , using the different dipole models. The results of the predictions for σ T (figure 3(a)) and σ L ( figure 3(b)) are in agreement with the data for values of Q 2 larger than approximately 7 GeV 2 for σ T and 7 GeV 2 for σ L depending on the considered dipole model. The results based on the AAMQSa model are giving a slightly better description than with the AAMQSb although both sets give good quality fits of DIS data. The results obtained by GBW model are close to the one obtained by the AAMQSa model; both are in good agreement with the data for Q 2 above 7 GeV 2 . Let us insist on the fact that the agreement for Q 2 ≥ 7 GeV 2 of our predictions with data is non-trivial since all free parameters of the   • the WW contribution, only involving the WW solutions of the DAs.
• the asymptotic (AS) contribution (for µ F → ∞), involving the asymptotic solutions of the DAs. In this limit the genuine contribution vanishes and the WW DAs can be expressed as functions of the asymptotic DA ϕ 1 (y) = 6yȳ.
The results of the predictions for σ T (figure 4(a)) and σ L ( figure 4(b)) are in agreement with the data for values of Q 2 respectively larger than Q 2 min,T ∼ 6.5 GeV 2 and Q 2 min,L ∼ 5 GeV 2 which confirms that the amplitude factorizes into a universal color dipole scattering amplitude and that the truncated twist expansion of the ρ meson soft part is justified.

JHEP11(2013)062
We remind that in our approach the skewness effects are not taken into account. The agreement for large Q 2 values of our predictions with the data suggests that the skewness effects are small for the kinematical range of HERA data we analyzed.
The two scales Q 2 min,T and Q 2 min,L are close to each other. We interpret this fact as an indication that the discrepancy between data and our predictions at low Q 2 are mainly due to the higher twist contributions to the impact factor from the meson structure rather than an effect of the saturation dynamics of the nucleon which should be well described at this scale by the saturation models. Let us mention, again that these saturation models are known to fit very well inclusive DIS as well as diffractive DIS data (for GBW) at these low Q 2 values.
The saturation scale, given by Q S = 1/R 0 (x) is of order 1 GeV in the kinematics of HERA. Since our predictions are only consistent with data in the region Q 2 > Q 2 min,L(T ) > Q 2 S , this limitation do not allow us to access the domain Q 2 S Q 2 where saturation effects can be essential.
The predictions are dominated by the WW-contribution and are not very sensitive to the choice of the collinear factorization scale. This will be further discussed in section 6.
An estimation of the error on the cross-sections caused by the error bars on the b−slope measurements is obtained by fitting the upper and lower bounds of the b−slope as shown in figure 5, and then by computing the predictions based on these fits as shown in figure 6. Note that we have assumed that the longitudinal b L and the transverse b T slopes are equal to the b−slope of the total cross-section. This assumption is supported by H1 data where the measurements of the difference b L −b T for Q 2 = 3.3 GeV 2 and Q 2 = 8.6 GeV 2 are much smaller than the b−slope value. Let us also emphasize that in this approach we compute the polarized differential cross-sections in the limit t = t min ≈ 0 where only the s-channel helicity conserving (SCHC) amplitudes T 00 and T 11 are non-zero. The contributions of other helicity amplitudes are encoded in the phenomenological t−dependence given in eq. (3.28) and it turns out that data for the total differential cross-section are dominated by a t-region of very small values, with a typical spread given by the scale 1 b ≈ 1 6 GeV 2 . We now compare our predictions with the data for the total cross-section σ, given by the sum σ = σ L + σ T according to ZEUS convention in ref. [3] or σ = εσ L + σ T following H1 notation [4], where ε is the photon polarization parameter 6 ε ≃ (1 − y)/(1 − y + y 2 /2) .

(5.2)
We show in figures 7(a) and 7(b), the AS, WW and Total contributions to the total crosssection σ as function of Q 2 for fixed averaged W . The predictions are larger than the data for Q 2 smaller than approximately 7 GeV 2 , as expected from the results of σ L,T . The figures 8(a)), (8(b)), (9(a) and 9(b) show the W dependence of the cross-section for several values of Q 2 . We see again a good agreement between the predictions and the data for Q 2 approximately above 6 GeV 2 for H1 data and 8 GeV 2 for ZEUS data, taking into account the uncertainty on the b−slope. Finally, our analysis also provides predictions for the ratios R and for the spin density matrix element r 04 where x 11 = |T 11 |/|T 00 | . H1 and ZEUS measurements of R and r 04 00 = σ L /σ as functions of |t| confirm this weak dependence on |t|. Based on H1 data, we can estimate [42] the correction to the ratio r 04 00 due to the amplitude T 01 , for the t−range of H1, to be below 1%. The results are shown in figure 10 for the ratio R and in figure 11 for the spin density matrix element r 04 00 , using AAMQSa model.
for the amplitude T 00 and W gen 11 (y, r; µ 2 F , Q 2 ) = ψ γ * T →ρ T gen (qq) (y, r; Q, µ 2 F ) +  for the Total, the WW and the genuine contributions to T 11 . The probability amplitudes W λρλγ permit in turn to define the radial distributions P λρλγ of the interacting dipole, as where N λρλγ are normalization factors N λρλγ = ∞ 0 dr r dy W λρλγ (y 1 , y, r; µ 2 F , Q 2 ) .  Expressed in terms of these functions, the scattering amplitudes read The probability amplitude for a dipole of size r to scatter on the nucleon is then proportional to P λρλγ (r, Q 2 , µ 2 F )σ(x, r), which justifies the inclusion of the factor r in (6.5).   Below we will also use the rescaled radial distributions P λρλγ (λ, µ 2 F ),

JHEP11(2013)062
which depend on r and Q only through the variable λ = r Q and we choose to put µ 2 F = µ 2 F (Q 2 ), see eq. (5.1). Note that the rescaled asymptotic distributions, P   This change of variable leads to the formulas The average value of a function f (y) depending of the longitudinal fraction of momentum y carried by one of the partons, will be estimated by f (y) λρλγ = 1 N λρλγ dr dy f (y) r W λρλγ (y, r; µ 2 F , Q 2 ) . (6.12) 6.1 The radial distribution of the γ * L → ρ L transition The distributions P 00 (r, Q 2 , µ 2 F ) and P (AS) 00 (r, Q 2 ) ≡ P 00 (r, Q 2 , ∞) are close to each other, as it is shown in figures 12(a) and 12(b), which indicates that the distribution P 00 (r, Q 2 , µ 2 F ) is not sensitive to µ 2 F . We can then restrict the study of P 00 (r, Q 2 , µ 2 F ) by considering only P (AS) 00 (r, Q 2 ), which is also simpler for the analytic treatment. At first glance we see that the distributions are peaked around r ∼ 1.3 Q and consequently the peak moves to the right and the distribution becomes wider as Q 2 decreases. Note that the dependency ofσ(x, r) with respect to Q 2 , which can be seen in figure 12, only occurs through the dependency of R 0 (x), according to eq. (4.4).
In figure 13 we show the Total and AS rescaled radial distributions P 00 (λ, µ 2 F (Q 2 )) and P (AS) 00 ). This last one reads P (AS) 00 The average value of λ estimated with P  Figure 11. Predictions for r 04 00 vs W and Q 2 compared respectively with H1 [4] data (figures (a)) and ZEUS [3] data (figures (b)), using the AAMQSa-model. About half of the dipoles are contained in the region 1 < λ < λ (AS) 00 , the peak of the distribution being at λ peak ∼ 1.3. The typical transverse scale µ = yȳ Q 2 entering the wave functions overlap can be estimated using eq. (6.12), The choice of the factorization scale µ F (Q 2 ) given by eq. (5.1) is then a good approximation of the transverse dynamical scale µ (AS) 00 involved in the process. The dipole scattering amplitude plays the role of a filter that selects dipoles of λ > λ Sat. (Q 2 , W ) = 2R 0 (x) Q. Note that the critical saturation line eq. (4.5) is given by   (r, Q 2 ) vs the size r of the interacting dipole, for Q 2 = 1 GeV 2 (12(a)) and Q 2 = 10 GeV 2 (12(b)), and the dipole cross-sectionσ(x, r) normalized by the factor 4σ 0 for W = 50 GeV and W = 150 GeV. where in accordance with eq. (4.6). In the kinematics of HERA, the energy in the center of mass W varies roughly from 50 GeV to 150 GeV, leading to the following bounds for the two values Q 2 = 1 GeV 2 and Q 2 = 10 GeV 2 , as one can see in figure 13 (plotted in logarithmic scale). The large difference between N λ>λ Sat. (1,90) and N λ>λ Sat. (10,90) indicates that the integrand of T 00 shown in figure 14, is very sensitive, when Q 2 varies between 1 GeV 2 and 10 GeV 2 , to the overlapping of the dipole cross-section bandwidth (λ > λ Sat. (Q 2 , W )) and the radial dipole distribution P (AS) 00 (λ); we are then probing with a high accuracy the quality of the shape of the dipole cross-section.
At high Q 2 the tails of the distributions plays a dominant role. The tail of the distribution corresponds to the region where the integrand of the radial distribution can be approximated by an exponential fall, where typically (6.23) In the case Q 2 = 1 GeV 2 , the bandwidth of the dipole cross-section mostly overlaps with the peak of the distributions as λ Sat. (1, 90) ∼ λ peak , while it only overlaps with the tail of the distribution when Q 2 = 10 GeV 2 , λ Sat. (10, 90) ∼ λ tail . Figure 14 shows the normalized integrand of T 00 . This summarize our discussion on the respective role of the radial distribution and of the dipole cross-section. This integrand is peaked near the saturation radius r ∼ 2R 0 (x). Comparing the cases Q 2 = 10 GeV 2 and Q 2 = 1 GeV 2 , we see that this peak is moving to the right as Q 2 decreases, going through the bandwidth of the dipole cross-section.
Contrary to the γ * L → ρ L transition, we see that the dependence on µ 2 F is quite strong if one compares P 11 (λ, µ 2 F ) to P (AS) 11 (λ). We cannot then restrict ourselves to only study the asymptotic case.
It is interesting to estimate the average λ obtained with the different distributions P 11 , P The average value of λ estimated with the WW distribution is, The effect of the term involving a 2 (µ 2 F ) in the r.h.s. of eq. (6.26) is under 4%, which indicates that P W W 11 (λ, µ 2 F ) ∼ P AS 11 (λ, µ 2 F ) is a good approximation. The computation of the average values of λ for the all the different contributions are given in table 4. 7 The tilde is to differentiate the contributionsP (λ) and the dipole cross-sectionsσ(x, λ) (black dot-dot-dashed lines) at W = 90 GeV 2 , for Q 2 = 1 GeV 2 (thick lines) and Q 2 = 10 GeV 2 (thin lines).
The results in table 4 show that the γ * T → ρ T transition is more sensitive to saturation effects than the γ * L → ρ L transition as λ 11 is about twice larger than λ 00 . Indeed it means that more dipoles are produced in the bandwidth of the dipole cross-section by the radial distribution P 11 (λ, µ 2 F ) than P 00 (λ, µ 2 F ). As a consequence the polarized crosssection σ T should be a more sensitive observable to probe features of the saturation regime than σ L .
The transverse dipole scale associated to the WW contribution, using eq. (6.12), is which is not so far from the values of the function µ 2 that is used here. The tail of the distributionP W W 11 (λ, µ 2 F ) can be defined as, is approximated by √ yȳ AS 11 ≈ 0.37. The genuine contribution P gen 11 (λ, µ 2 F ), which vanishes in the limit µ 2 F → ∞, depends strongly on the factorization scale µ 2 F (Q 2 ), as one can see in figure 15. The distribution P gen 11 (λ, µ 2 F ) can be split into two contributions, whereP gen (qq) 11 (µ 2 F ) ∼ 0.5 and is not very sensitive to µ 2 F . The choice of µ 2 F (Q 2 ) ∼ Q 2 is then a good choice for this contribution.
The contributionP gen (qqg) 11 (λ, µ 2 F ) have several transverse scales which are µ 1 , µ 2 , µq g , µ qg and µ qq defined by eqs. (3.18), each corresponding to a dipole configuration involving two of the three partons available in the process. In order to estimate these transverse scales, we first evaluate the average fraction of momentum carried by the quark where we assume that the transverse scale values are roughly approximated by using the average fractions of longitudinal momentum in eqs. (3.18). These values are evaluated at µ 2 F (1 GeV 2 ). Other values of µ 2 F have also been used, leading to approximately the same results. We note that the function µ F (Q 2 ) ≈ Q 2 is close to the µ 1 and µ qq values, and of the same order of magnitude than µ qg but in principle, one should adapt the choice of the factorization scale to the relevant transverse scales at stake for each part of the process.
The WW and the genuine contributions to the radial distribution P 11 (λ, µ 2 F ) are of the same order of magnitude for Q 2 ∼ 10 GeV 2 , and the genuine contribution becomes even more important for Q 2 = 1 GeV 2 . At Q 2 = 10 GeV 2 , as However the fact that the genuine contribution is important even at large Q 2 , indicates that the three-parton exchange between the γ * T and ρ T states is important in such relations as the normalization eq. (6.36) of the ρ meson wave function [7,87], and the electronic decay width eq. (6.37) [7,88], where the exchange of only two-parton is assumed. Indeed, the r.h.s. of eq. (6.37), if one expands at large Q 2 the ρ meson wave function around r = 0 , is the WW approximated result, which therefore misses the genuine contributions arising from three-parton correlators, which can have a significant effect even for large Q 2 values, see figure 15. These relations are usually used to constrain the parameters of the wave function models of the ρ meson, assuming that the meson is solely constituted of a quark and an antiquark, which then consist in neglecting the higher Fock state contributions like the genuine contribution.
cross-section bandwidth (λ Sat. (1, 90) < λ gen 11 ) and gives an important contribution to the integrand of T 11 . When increasing Q 2 , as shown in figure 16 where we display the product P 11 (r, Q 2 , µ 2 F )σ(x, r) for Q 2 = 1 GeV 2 and Q 2 = 10 GeV 2 , we can note that the difference between the AS and the Total results is a consequence of the genuine contribution growing when Q 2 decreases. It is the convolution with the dipole cross-section which washes-out the effect of these genuine twist-3 contributions.
6.3 Comparison with the radial distributions obtained from models of the ρ meson wave function.
It is instructive to compare shapes of the radial distributions P 00 and P 11 used in our analysis with those used in two other approaches which involve the overlap of the virtual photon wave functions and of the ρ meson wave functions: • the "Boosted Gaussian" (BG) model [46], • the "Gaus-LC" model [55].
We use here the convention and the parameter values of ref. [47], which for completeness are shown in table 6. The scalar parts of the wave functions are given by, The overlaps with the virtual photon wave function are, with δ = 0 for the Gaus-LC model and δ = 1 for the BG model. The radial distributions thus read, h,h (y, r) , (6.43) where the factors N L,T normalize the distributions P L,T (r). Comparing the distributions at Q 2 = 1 GeV 2 (figures 17(a) and 17(c)) and Q 2 = 10 GeV 2 (figures 17(b) and 17(d)) we see that at large Q 2 , in the bandwidth of the dipole cross-section, our distributions are converging with the distributions obtained from the Gaus-LC and BG models.
In the γ * T → ρ T case in figures 17(c) and 17(d), we see that the BG and Gaus-LC models are closer to the distribution P 11 than to the asymptotic distribution. When compared to other distributions, the asymptotic distribution P AS 11 is shifted to the right, thus selecting larger dipole sizes.
In the γ * L → ρ L case, in figure 17(a) we see that the distributions from the Gaus-LC and BG models are not close to our predictions, indicating that the higher twist corrections are presumably more important at small Q 2 than in the γ * T → ρ T transition. These qualitative remarks remains the same when using the AAMQS dipole models, and we expect that they are model independent.

Conclusions
We performed phenomenological analysis of experimental data from HERA on ρ meson electroproduction within the approach based on the recently derived [44] impact parameter representation of γ * → ρ impact factor up to twist-3 accuracy. The important feature of this representation consists in the inclusion of contributions coming from both two-and three-partonic Fock states, maintaining a close connection with the dipole model picture. Consequently, it was possible to include in our framework based on collinear factorization, the saturation effects. The inclusion of saturation effects in our approach was done in a standard way by using the existing models for the dipole cross-section, which saturate for large dipole sizes, in the factorized form of the scattering amplitude. Our predictions show that we can get simultaneously good predictions for the polarized cross-sections σ T and σ L . The ability of the model to reproduce the data is the confirmation of the following points:     and AS (red, dashed) radial distributions for the γ * L → ρ L transition (top) and for the γ * T → ρ T transition (bottom), vs r for Q 2 = 1 GeV 2 (left) and Q 2 = 10 GeV 2 (right), as well as the dipole cross-sectionσ(x, r) rescaled by the factor 5σ 0 for W = 90 GeV (black, dot-dashed).
• the factorization of the dipole cross-section in the helicity amplitudes of the electroproduction of the ρ meson works and, as an universal quantity, is the same for T 00 and T 11 , giving the good energy dependence and normalizations of the polarized cross-sections, • the collinear factorization procedure of the ρ meson is justified and works successfully beyond the leading twist.
As expected the model has some limits due to the truncation of the twist expansion. Thanks to HERA data, we have identified the virtualities Q 2 min,L(T ) ∼ 5 GeV 2 were the higher twist corrections become important, which is a motivation to compute impact factors beyond the twist 3 accuracy in order to probe the genuine saturation regime which starts at Q 2 S ∼ 1 GeV 2 .
Other helicity amplitudes could be computed keeping the same approach, they would be useful in the t = t min regime. The kinematics of the impact factor can be also extended to take into account the t−dependence of the impact factors, which would be a test for the dipole models which include the impact parameter dependence, providing a probe of the proton shape [7], in particular through local geometrical scaling [89,90].
Data also exist for φ leptoproduction. In this case quark-mass effects should be taken into account, in particular, because this allows the transversely polarized φ to couple JHEP11(2013)062 through its chiral-odd twist-2 DA. Indeed, as it was pointed out in [42], the fact that the ratio T 11 /T 00 is not the same (after trivial mass rescaling) for ρ and φ mesons points to the importance of this effect. This is also beyond the scope of our present study, but may open an interesting way for accessing chiral-odd DAs.
The next-to-leading order effects -both on the evolution and on the impact factorshould be studied, since it is now known that both may have a important phenomenological effect [91][92][93][94][95][96].

A Distribution amplitudes in the LCCF parametrization
The seven chiral-even 8 ρ-meson DAs up to twist 3 are defined through 9 the following matrix elements of nonlocal light-cone operators [40], for two-parton

JHEP11(2013)062
and three-parton correlators where we used the standard notation . DAs are linked by linear differential relations derived from equations of motion and n−independency condition [39,40]. The solutions for ϕ P (y) ≡ {ϕ 3 , ϕ A , ϕ T 1 , ϕ T A } are the sum of the solutions in the so-called WW approximation and of genuine solutions, The WW approximation consists in neglecting the contribution from three-parton operators, thus taking B(y 1 , y 2 ; µ 2 F ) = D(y 1 , y 2 ; µ 2 F ) = 0. Then, ϕ W W P (y) become functions of ϕ 1 (y) only, and their explicit expressions [109] are given by Genuine solutions only depend on {B(y 1 , y 2 ; µ 2 F ), D(y 1 , y 2 ; µ 2 F )} or equivalently on the combinations {S(y 1 , y 2 ; µ 2 F ), M (y 1 , y 2 ; µ 2 F )} defined by eq. (2.8), namely where A(u; µ 2 F ) has the compact form and it obeys the conditions The correspondence between our set of DAs and the one defined in ref. [50] is achieved through the following dictionary derived in ref. [40]. It reads, for the two-parton vector DAs, ϕ 1 (y) = φ (y), ϕ 3 (y) = g The dependence on the renormalization scale µ F of the coupling constants a 2 , ω A {1,0} , ζ A 3 , and ζ V 3 are given in ref. [50]. In appendix B we present both the evolution equations and the values of these constants at µ 2 F = 1 GeV 2 used in our analysis, as well as the dependence on µ F of the DAs.  Table 7. Coupling constants and Gegenbauer coefficients entering the ρ meson DAs, at the scale µ 0 = 1 GeV updated in ref. [51]. Note that in ref. [50] the normalization are such that f V,A 3ρ [50] =

B Evolutions of DAs and coupling constants with the renormalization scale
The parameters entering the DAs at µ 2 0 = 1 GeV 2 are updated 10 in ref. [51] and their evolution equations are given in ref. [50], we recall in table 7 their values for the ρ meson.

D Results using the GBW and AAMQSb models
We present some of the predictions obtained by using the GBW or the AAMQSb model for the dipole cross-section. As expected the results are not so far from the ones obtained with the AAMQSa model. In figure 18 and 19 are respectively shown for the AAMQSb and the GBW model, the polarized cross-sections σ T and σ T . The spin density matrix element r 04 00 predictions using these dipole models are shown in figure 20, for completeness we show also the prediction obtained with the GS model.