Exclusive vector meson production and small-x evolution

The process of exclusive elastic vector meson production in deep inelastic scattering is investigated within the dipole model framework supplemented by the small x evolution. The dipole-proton amplitude is obtained from the nonlinear Balitsky-Kovchegov evolution equation with impact parameter dependence. This dipole amplitude is used to compute the differential cross section for exclusive production of J/Ψ, ϕ, and ρ vector mesons. These numerical calculations are compared with the wide range of experimental data from HERA. Good agreement between the experimental data and the calculations is found.


Introduction
Exclusive diffractive vector meson production in deep inelastic scattering, γ * p → V p, has been extensively studied at high energies at HERA electron(positron)-proton collider. Both H1 [1][2][3] and ZEUS [4][5][6][7] collaborations performed detailed measurements of this process in a wide kinematic range. In this process the final state consists of the elastically scattered proton accompanied by a vector meson. It provides important information about the electron-proton interaction at small x as well as the size of the proton in impact parameter. It can be also used to constrain the form of the vector meson wave functions. The main conclusions from these measurements can be summarized as follows (for concise reviews of the experimental results see for example [8,9]). The cross section for the elastic vector meson production increases with the energy for all values of the (minus) photon virtuality Q 2 . This dependence can be phenomenologically parametrized as ∼ W δ where W is the energy of the virtual photon-proton system in DIS. The exponent δ depends strongly on the scales involved in the process, namely on the Q 2 and the mass M V of the vector meson. The exponent δ shows a marked increase with both M V and Q 2 scales. For example, for ρ production the variation of δ with Q 2 is from 0.2 for photoproduction up to ∼ 0.8 for values of Q 2 ∼ 35 GeV 2 . This is consistent with the picture of the soft-hard transition from small to large scales. The value of the δ in the photoproduction region, and at low Q 2 , is consistent with the soft exchange, while in the high Q 2 region perturbative gluon exchange is the dominating process. For low values of Q 2 multiple interactions between the quark-antiquark system and a proton also lead to smaller values of δ, see for example [10]. For J/Ψ production the variation of δ with Q 2 is negligible in the measured range (from photoproduction to electroproduction with Q 2 ∼ 30 GeV 2 ), and its value is JHEP01(2013)001 about ∼ 0.8, indicating that the heavy mass of the J/Ψ vector meson provides the hard scale in this process.
The measurement of the momentum transfer t in this process provides the information about the transverse size of the interaction region, see for example [11]. These measurements indicate a strong dependence of the differential cross section on t, which can be parametrized by an exponential exp(−B D |t|). The experimental data show that the value of B D decreases from about ∼ 10 − 12 GeV −2 at lowest values of (Q 2 + M 2 V ) to about 5 GeV −2 at scale Q 2 + M 2 V ∼ 10 GeV 2 , and stays constant for larger scales. This can be interpreted as follows: the value of parameter B D is closely related to the transverse size of the interaction region which is a combination (or more precisely, a convolution) of the sizes set by the vector meson and the proton. At low values of Q 2 + M 2 V it is the first size that prevails, and as a result the value of B D strongly varies with this scale. At larger values of Q 2 + M 2 V the data indicate that the vector meson size is much smaller than the typical size of the proton, (or more precisely the radius of the gluon density in the proton), and it is this 'gluonic' size of the proton which is dominating B D and thus visible at larger values of Q 2 + M 2 V . This value is about 5 GeV −2 or equivalently ∼ 0.6 fm, and it is smaller than the electromagnetic size of the proton. The data also indicate a slow variation of the slope B D with the increasing energy of the interaction W . The slope can be parametrized by the following Regge inspired formula B D = B 0 + 4α ′ IP ln(W/W 0 ), with the fitted value of α ′ IP ∼ 0.115 ± 0.018(+0.008)(−0.015) [4] and α ′ IP ∼ 0.164 ± 0.028 ± 0.030 [2] for J/Ψ photoproduction (parameter W 0 = 90 GeV). This shows that the size of the gluon distribution in impact parameter space is increasing with rising energy, which is indicative of Gribov-type diffusion of gluons in transverse space.
In this paper we analyze elastic diffractive vector meson production using the dipole model approach for small x physics. The novel aspect of our analysis is that the dipole scattering amplitude is obtained from the numerical solution to the impact parameter dependent Balitsky-Kovchegov (BK) equation [12][13][14]. Numerous studies, which utilized the dipole model, were performed up to date whose goal was a description of the exclusive diffractive vector meson production at small x, see for example [10,[15][16][17][18][19][20][21]. As mentioned above, the t-dependence of the vector meson production cross section is particularly interesting as it is directly related to the transverse size of the interaction region. In the framework of the dipole model this is encoded in the dipole-proton scattering amplitude and its dependence on the impact parameter. The increase of the t slope, B D , with the energy W should be directly coupled with the increase of the transverse interaction area caused by the diffusion of the dipoles in coordinate space. This diffusion in transverse coordinate space is present in the dipole branching Monte-Carlo [21][22][23][24] as well as the Balitsky-Kovchegov equation when the impact parameter dependence is taken into account [25,26]. The main goal of the present work is to use the numerical solution to the BK equation, which includes the full impact parameter dependence, and verify its compatibility with the experimental data on the exclusive vector meson production in deep inelastic lepton-proton scattering.
The outline of the paper is the following, in the next section we recall the formalism of exclusive vector meson production in the dipole model, in subsection 2. the basic facts which pertain to the solution of the BK equation with impact parameter dependence, and discuss modifications due to confinement. In subsection 2.2 additional modifications are outlined, such as the non-perturbative modification of the photon wave function and inclusion of the skewed effect into the initial gluon distribution, as well as inclusion of the real part of the scattering amplitude. This is followed by the description of the model for the vector meson wave function which is used for the calculation. In section 3 we present our results, which include the comparison of the theoretical calculation with a wide range of experimental data from HERA on ρ, φ and J/Ψ production. We compute the cross section integrated over the momentum transfer t and investigate its W and Q 2 dependence. Then we compute the differential cross section with t dependence in bins of W and Q 2 . The experimentally measured scale and energy dependence of the slope B D is also compared with theoretical calculations. Finally, in the last section we state the conclusions and present an outlook.

Exclusive vector meson production in the dipole model
The dipole model [27,28] is a very useful tool in evaluating many processes at small values of x. One of the advantages of this approach is the possibility of including the multiple parton scattering effects. It has been originally formulated for the description of deep inelastic lepton-proton (or nucleus) scattering at small x. In this picture, utilizing the leading logarithmic approximation in x, the incoming electron emits a virtual photon which fluctuates into a quark-antiquark pair (a dipole). This color dipole then subsequently interacts with the parton constituents of the nucleon, as is illustrated in figure 1. The interaction of the dipole pair with the target is given by the scattering amplitude N . This quantity is non-perturbative in principle, however its small x evolution can be found from the BK equation. The qq pair is characterized by a dipole size which is defined as a separation distance of the color charges x 01 = x 0 − x 1 (where x 0 and x 1 are the positions of the q andq in transverse space). 1 The transverse momentum of the quarks in 1 In this paper we shall denote vector quantities in bold, otherwise they should be read as magnitudes of the associated vector. Also, alternatively we will be also using here the notation for the dipole size to be r = x01 and impact parameter b = |x 0 +x 1 | 2 .

JHEP01(2013)001
the dipole is on the order of ∼ 1 x 01 where large dipoles correspond to the infra-red region and need to be regulated. The scattering amplitude N (r, b; Y ) contains all the information about the dynamics of the strong interaction. It depends on the dipole size r, on the impact parameter of the dipole with respect to the target b and the rapidity Y . In the following analysis the full dependence of the scattering amplitude on the impact parameter b will be taken into account.
The structure functions F 2 and F L for the proton can be evaluated using the following standard formulae in the dipole picture in the transverse coordinate representation and The dipole-proton cross section σ dip can be obtained from the scattering amplitude by integrating over the impact parameter b Since the amplitude N is dimensionless and the integration is over the impact parameter, the dipole cross section σ dip has obviously a dimension of the area. We see therefore that although the inclusive quantities are sensitive to the size of the interaction area, the details of the impact parameter profile are not directly accessible through this process. The quantities Ψ(r, Q 2 , Y ) T /L are the photon wave functions. They describe the dissociation of a photon into a qq pair and can be calculated from perturbation theory. The photon wave function has the following form for the case of transverse photon polarization and for longitudinal polarization In the above equationsQ 2 The dipole picture can be also used to compute diffractive processes. Here, we are interested in the process of the exclusive, diffractive production of the vector-meson γ * p → V p ′ . The amplitude for this process is schematically illustrated in figure 2. The virtual photon still fluctuates into qq pair, which then interacts with the proton and a vector meson is formed, which is measured in the final state. The proton scatters elastically, Figure 2. Schematic representation of exclusive vector meson production in the dipole model at small x. A virtual photon fluctuates into a qq pair and interacts with the target(proton). After the interaction the vector meson V is formed which is measured in the final state. The proton scatters elastically with some momentum transfer t.
its 4-momentum in the initial state is p and in the final state is p ′ . The formula for the amplitude for this process reads where h(h) is the helicity of quark (antiquark) and Ψ V h,h is the vector meson wave function. ∆ is the 2-dimensional momentum transfer related to the Mandelstam variable t = −∆ 2 .
The differential cross section for the process is given by The amplitude N (x, r, ∆) can be related to the scattering amplitude N (x, r, b) introduced earlier, the amplitude in the impact parameter representation through the appropriate 2-dimensional Fourier transform (2.8) In this notation the dipole cross section, (compare (2.3)), is which is the expression for the optical theorem for scattering of dipoles. This process, through its dependence on the momentum transfer t, offers a unique possibility of constraining the impact parameter profile of the dipole scattering amplitude.
Formulae (2.6) and (2.8) were original expressions derived under the assumption that the dipole size is much smaller than the proton. In ref. [29], a correction due to the finite size of the dipole was calculated. It was shown that in the non-forward case, ∆ = 0, the amplitude can be written in the similar form as above with the modification of the (2.8) to include the exponential factor exp(−i(1 − z)r · ∆) in the following way (2.10)

JHEP01(2013)001
This modification was included in the calculation [10] and it was shown that it has a nonnegligible effect on cross sections, especially on the values of the B D slope which controls the t-dependence as a function of the scale Q 2 + M 2 V .

Dipole scattering amplitude from impact parameter dependent BK evolution
The dipole-proton scattering amplitude N (r, b; Y ) at high values of rapidity Y (or small x) is found from the solution to the nonlinear integro-differential Balitsky-Kovchegov (BK) evolution equation [12][13][14]30]. The BK evolution equation can be represented in the following form: In the above equation we used the shorthand notation for the arguments of the ampli- ; Y ) which depends on the two transverse positions x i and x j and on the rapidity Y . The branching kernel K(x 01 , x 12 , x 02 ; α s , m) depends on the dipole sizes involved and contains all information about the splitting of the dipoles. In addition, it depends on the running coupling α s . The way the strong coupling runs will be specified later in this work. We have also indicated that the kernel depends on the infra-red cutoff m which we impose in order to regulate large dipoles.
Eq. (2.11) is a differential equation in rapidity and hence suitable initial conditions need to be specified at some initial value of rapidity Y = Y 0 . As in the previous work [31] we are choosing to use the initial condition in the form of the Glauber -Mueller parametrization with (most of) the parameters equivalent to those used in ref. [10] In formula (2.12) the function xg(x, η 2 ) is the integrated gluon density function and T (b) is the density profile of the target in transverse space with the extension set by the parameter B G . The integrated gluon density in (2.12) was also taken from fits performed in [10]. Scale parameter in the gluon density is set to be η 2 = µ 2 0 + C 2 r 2 with parameters µ 0 and C = 2 set to obtain the best description of the data. The values of these parameters are given in table 1. We use (2.12) as the initial condition at Y 0 = ln 1/x 0 , x 0 = 10 −2 and evolve the amplitude with the BK equation to obtain the solution at lower values of x < x 0 . We also note that the initial condition (2.12) depends only on the absolute values of the dipole size and impact parameter. A nontrivial dependence on the angle between vectors r and b is not present in the initial condition, instead being dynamically generated when the initial condition is evolved with the BK equation.
The BK equation was solved numerically by discretizing the scattering amplitude in terms of variables (log 10 r, log 10 b, cos θ), where θ is the angle between the impact parameter JHEP01(2013)001 b and the dipole size r. The amplitude N (r, b, cos θ) was placed on a grid with dimensions 200 r × 200 b × 20 θ . More details about the techniques and the properties of this solution can be found in refs. [25,26].
The running of the coupling in the BK kernel is a next-to-leading effect which has been evaluated in [32,33]. In this calculation we utilize the prescription from ref. [32] which is of the form (2.14) In the above equation, we use and n f is the number of active flavors. The µ parameter effectively freezes the coupling at large dipole sizes at α s,fr = 1 b ln[Λ −2 µ 2 ] . In our simulations we used µ = 0.52 GeV as an infra-red regulator for the strong coupling, Λ QCD = 0.246 GeV, and the rescaled strong coupling is defined asᾱ s = αsNc π . The evolution equation (2.11) has been derived in perturbation theory and it does not include any effects of confinement. The branching kernel is power-like in the dipole sizes and it allows for the splitting of a parent dipole into pair of arbitrarily large daughter dipoles. As a result the impact parameter dependence of the amplitude contains Coulomblike power tails in impact parameter, which need to be regulated. This has to be done by cutting off the large dipole sizes in the branching kernel. This is done by including a mass parameter, m, which accounts for the effect of color confinement. In previous works we have tested several prescriptions and found the simulation with the simple cutoff using the theta functions gave best description of the inclusive data. We have therefore used this scenario for the calculation in this paper. We stress that including the mass parameter into the evolution and in the initial condition cuts of the large dipole sizes in the dipole scattering amplitude and regulates the tails in the impact parameter due to the fact that both are strongly correlated in the evolution. The modified kernel which has been used for the calculation of a solution has therefore the form of the kernel (2.14) with theta functions where m is the mass regulator. The value of the parameter m has been fitted to obtain the best description of the data. Note that the similar procedure has been utilized in refs. [21,23], with the value of the parameter r max = 1 m to be equal around 3 GeV −1 . The formula (2.12) is used as an initial condition for the BK evolution for dipoles with sizes smaller than the cutoff 1/m, for dipoles larger than the cutoff the initial condition is set to zero.

Corrections to the dipole scattering amplitude and photon wave function
There are several additional phenomenological corrections which we have included in this calculation. As we will see explicitly they have a non-negligible impact on the calculations. The first of these is to take into account the effect of the skewed gluon distribution. This is necessary as the two gluons exchanged in the production of the vector mesons need not have the same longitudinal momentum fractions x and x ′ . The skewed effect vanishes at small x (in the leading logarithmic limit in ln 1/x), however it can be substantial correction when the energy is not very large. As suggested in ref. [34] this correction can be taken into account by a multiplicative factor on the standard gluon distribution as follows Strictly speaking, in the dipole formalism there is no integrated gluon distribution, but rather the dipole scattering amplitude. However, since the initial condition is taken as the Glauber-Mueller model in the form (2.12), with explicit dependence on the gluon density xg(x, η 2 ), we can implement the correction for skewedness inside the initial condition. Then the evolution to lower values of x is computed according to eq. (2.11).
In addition we have included the correction for the real part of the scattering amplitude. This effect can be taken into account by multiplying the amplitude by ( (2.20) In the above equation β is the ratio of the real to imaginary part of the scattering amplitude. A similar procedure was used previously in the other descriptions of the vector meson production [17,19,34]. Finally, a correction to the photon wave function is necessary in order to modify and enhance the contribution at low Q 2 . The photon light-cone wave function given in eqs.
( 2.4) and (2.5) has been derived in the perturbation theory under the assumption of the presence of hard scale Q 2 , which in turn leads to the small dipole sizes r. This statement is most accurate for the longitudinal contribution, but in the case of the transverse contribution the endpoint singularities in z cause the distribution in dipole sizes to be broader, even for very large scales Q 2 . At low values of photon virtuality Q 2 the corresponding dipoles which contribute to the cross section are very large and non-perturbative. Consequently the cross section with which they interact with the proton is also large. Therefore, at these large dipole sizes the photon has a hadronic component, which is non-perturbative. To account for this effect we use a modification of the photon wave function in a way suggested JHEP01(2013)001 Table 1. Values of the free parameters used in the calculations. Description 'fit' means that the parameter was fitted to the inclusive data, description 'fixed from KMW' means that the parameter was set equal to the value in ref. [10].
in ref. [19] and also used in ref. [21] The constants B, ω and R are the parameters which have to be adjusted to fit the data. This factor provides an enhancement for dipoles which have a hadronic size. The numerical values of the parameters B, ω, R which are used in the calculation are given in table 1.

Vector meson wave function
Several different models for the vector meson wave function exist in the literature, for example see [15][16][17][35][36][37][38][39]. Typically, there are many uncertainties in obtaining the wave functions for the vector mesons, however they are constrained by model independent features. First of all, ψ h,h V has to satisfy the following normalization condition Here, we have neglected possible contributions of the gluon or sea-quark states. In addition, the value of the wave function at the origin is related to the leptonic decay width Γ(V → e + e − ) of the vector meson, given by the following formula Here f V is the coupling of the meson to the electromagnetic current andê V the isospin factor, which is the effective charge of the quarks in units of the elementary charge e: for the ρ meson it is the charge of the combination (uū − dd)/ √ 2, i.e.ê V = 1/ √ 2. Finally, one requires that the mean radius be consistent with the electromagnetic radius of the vector meson.
In this paper we utilize the model proposed in [15,17]. In this approach the information from spectroscopic models is used to constrain the long distance physics. One assumes that the meson is composed of a constituent quark and antiquark which move in a harmonic  oscillator potential. This results in a wave function which has a gaussian dependence on the spatial separation between the quarks. Additionally, this model is supplemented by the short-distance physics driven by QCD exchange of hard gluons between the valence quarks of the vector meson. Finally, a relativization technique has to be applied to the wave function.
The wave function of the vector meson is parametrized as in [15,17]. In this model, the radial part of the wave function φ(r, z) satisfies the following normalization condition for the transversely polarized part and for the longitudinally polarized part The parameter m V is the mass of the vector meson. The function φ(r, z) (L/T ) itself is parametrized by: The masses of the light quarks are taken to be 0.14 GeV and the charm mass is taken to be 1.4 GeV. The two parameters Ψ 0L/T (1S) and R 2 are chosen by taking into account several constraints: the normalization conditions (2.25), (2.24) for the wave function has to be satisfied, and the value of the leptonic decay width must agree with the experimental measurement. Additionally, the mean radius of the meson has to be of the order of a hadronic scale. The normalization Ψ 0L/T (1S) is different for the longitudinal and transverse wave functions of each vector meson. On the other hand the value of R 2 is the same for each polarization but different for each vector meson. This means that the leptonic decay width (2.23) will be matched to the experimental value for the longitudinally polarized vector mesons while it will be slightly different for the transversely polarized ones. This procedure is identical to the one presented in [10] for this model of the wave function. The values of the parameters used for the vector meson wave functions can be found in table 2.

JHEP01(2013)001 3 Comparison with the experimental data
In this section we compare the results of our calculation based on the numerical solution of the BK equation with the impact parameter dependence to the data from H1 and ZEUS on ρ, φ, and J/Ψ production.
As discussed above, the initial condition for the evolution in x is given by eq. (2.12). This formula contains several free parameters, such as the value of B G , the parameters A g , λ g in the initial gluon distribution xg(x, η 2 ) = A g x −λg (1 − x) 5.6 , and the parameters µ 0 and C which enter into the definition of the scale η 2 = C 2 r 2 +µ 2 0 . The initial impact parameter profile T (b) (2.13) is assumed to be gaussian. This profile is substantially changed during the course of the evolution with x. In addition, the kernel contains two scales: the scale µ in the running of the coupling as well as the mass parameter m. All these parameters enter into the calculation of the dipole scattering amplitude which is found by the iterative solution of the eq. (2.11). The computational time that is required to obtain the solution for each set of parameter values is rather long (on the order of 24 hrs on 32 cores). Therefore we have chosen to vary only a subset of these parameters. To be precise, the parameters A, λ g , C, µ 0 were fixed and their values were taken from ref. [10]. The only parameter which was varied in the initial condition eq. (2.12) was B G . The value of B G was adjusted to obtain the best description of the t dependence of the differential cross section dσ/dt. In particular this parameter controls the magnitude of the experimentally determined value B D , which is the slope of the t dependence, i.e. dσ/dt ∼ exp(−B D |t|). In principle, the mass parameter in the kernel m is also a free parameter, but in this analysis we have set it to be equal to m = 1 √ 2B G , thus reducing the number of fit parameters. Lastly, the scale µ was also left as a fitting parameter to the inclusive data.
There are also parameters that are present in the formulae for the cross section for the production of the various vector mesons. These include the quark masses as well as the parameters in the vector meson wave function and non-perturbative correction to the photon wave function. The variation of these parameters does not cost so much computational time therefore a larger number of them were varied. The complete set of parameters used in this paper (with the description whether they are fixed or varied) can be found in table 1.
Using the framework described above we have performed a fit of the calculation to the inclusive data on the structure functions in DIS. The resulting values of the free parameters after the fit are shown in table 1. Note that this was necessary as the details of the procedure presented above are slightly different than that originally used in [31]. Good agreement between the inclusive data and the calculation was found using this procedure.
Before we proceed to the discussion and the comparison of the calculation with the experimental data on vector meson production we show the results of the calculation for the dipole scattering amplitude obtained from the solution to the BK equation. The dipole scattering amplitude is shown in figure 3   BK equation, whereas the dashed lines represent the model b-Sat used in [10] which serves as an initial value for the BK simulation. The calculations are shown for two values of rapidity Y = 4 and Y = 8. Also, initial condition at Y = 0 is shown by the dotted-dashed line which is also given by the b-Sat model. In the case of the dipole size dependence the characteristic cut-off at size 1/m is visible in the case of the solution using the BK evolution equation. The calculation is strongly modified in the vicinity of the cutoff with respect to the b-Sat model, but it converges to this model for small values of the dipole size, i.e. corresponding to the perturbative regime. We recall that the cutoff imposed on the large dipole sizes is necessary and corresponds to including confinement into the calculation. The sharp spike in the amplitude stems from the fact that the cutoff is imposed in the form of the theta functions and from the cutoff in the initial condition. It is possible to impose a different cutoff which will lead to the smoother behavior of the solution in that region. The resulting cross section for the physical observables has a smooth behavior as a function of virtuality Q 2 .
The left plot in figure 3 shows the profile in impact parameter both for BK solution and the b-Sat model. The unique feature of the BK solution is the diffusion in impact parameter with increasing rapidity which is evident in this plot, especially when compared with the b-Sat model which is characterized by the fixed extension in impact parameter. In    σ (nb) Figure 5. Cross section σ(Q 2 , W 2 ) for the vector meson production plotted as a function of (Q 2 + M 2 V ) for ρ elastic production. The experimental data are from [1,3,7].
diffractive electroproduction of vector mesons, where the cross section has been integrated over the momentum transfer t. The experimental data from H1 and ZEUS for ρ [1,3,7], φ [3,6] and J/Ψ [2,5] were used. Figures 5 and 6 shows the cross section for production of ρ, φ and J/Ψ vector mesons as a function of variable (Q 2 +M 2 V ), where M 2 V is the mass squared of the corresponding vector meson. This variable is commonly used instead of Q 2 itself as it provides the scale for the vector meson. One cannot, however, expect that cross sections for different vector mesons will behave identically when this variable is used. We therefore show the cross sections for different species of vector mesons separately. The momentum transfer t has been integrated in the ranges provided by the experimental data. For the data on ρ production from [1] the range in t is |t| < 0.5 GeV 2 and |t| < 1.0 GeV 2 for data from [7]; |t| < 3 GeV 2 for JHEP01(2013)001 σ (nb) Figure 6. Cross section σ(Q 2 , W 2 ) for the vector meson production plotted as a function of (Q 2 + M 2 V ) for φ (left plot) and J/ψ (right plot) elastic production. The experimental data are from [2,3,5,6]. data from [3]. For the data on φ production from [6] the range is |t| < 0.6 GeV 2 and |t| < 3 GeV 2 for data from [3]. Finally for J/Ψ the range of |t| < 1.2 GeV 2 in [2] and |t| < 1.0 GeV 2 for data from [5]. The theoretical curves shown in figures 5, 6 are evaluated for constant energy W where W = 90 GeV for the J/Ψ cross section and W = 75 GeV or W = 90 GeV for the φ and ρ cross sections depending on the data set used.
We have tested the sensitivity of the cross section due to the inclusion of the various corrections discussed in section 2.2. This is shown in two plots in figure 7. In the left JHEP01(2013)001 plot the effect of including the correction due to the skewed effect in the initial condition is shown with data for J/Ψ production. We see that the effect is the substantial change in the normalization of the cross section, without distortion of the Q 2 dependence, and calculation without skewed effect correction completely fails to match the data. We note that even with the substantial variation of the free parameters it would be difficult to match the normalization of the cross section without this additional effect.
In the plot on the right hand side we demonstrate the effect of including the nonperturbative modification of the photon wave function in the form of eq. (2.21). This correction is especially important for low values of Q 2 and light vector mesons. This is to be expected as this correction vanishes for large Q 2 in order to reproduce the perturbative expression for the photon wave function. This will be better illustrated in figure 8, 9 which demonstrate the overlap distribution between the photon and meson wave function. On the other hand, this correction has a negligible effect at large values of Q 2 > 20 GeV 2 and the cross section for J/Ψ production, therefore the calculation is free from this non-perturbative contamination at higher, perturbative scales. We note that this is different than in approach presented in [10] where such non-perturbative correction was not necessary. The need for the additional non-perturbative enhancement in the photon wave function is most likely related to the fact that the dipole scattering amplitude obtained from the BK evolution equation has a confinement cutoff 1/m on large dipole sizes. On the other hand in the approach [10] such cutoff is absent and the cross section gets sizeable contributions from very large dipole sizes. We note that such additional non-perturbative correction in the approach using BK evolution with confinement was needed to obtain good description of the F 2 data at low values of Q 2 , see [31].
To illustrate the effect of the non-perturbative correction included in the photon wave function we examine the overlap between the vector meson wave function and the photon wave function as well as the integrand of the amplitude as a function of the dipole size, see figures 8, 9. The overlap function, see for example [18], is defined as There are two calculations shown in these figures, one which includes the non-perturbative correction to the photon wave function and the second one without it. As expected, in the case of ρ production at small values of Q 2 , there is a moderate enhancement in the vicinity of the cutoff which compensates for cutting the large dipole sizes. For J/Ψ this enhancement is almost negligible, meaning that for the heavier vector mesons and/or larger scales the non-perturbative correction is irrelevant.
The W energy dependence of the cross sections is shown in figures 10, 11, 12 for ρ, φ, J/Ψ respectively. The different curves are plotted in different bins of Q 2 . Overall trend in both cases, i.e. for Q 2 and W dependence, is such that the calculations describe the dependence in Q 2 and W dependence very well for the case of φ and J/Ψ. The data for ρ production are not well described, in particular the normalization in this case is systematically little bit low, especially for lowest values of Q 2 . This region is however the one which is not under perturbative control and some unknown non-perturbative corrections, Overlap function 500 x P with wavefunction correction Overlap function 500 x P without wavefunction correction Integrand 100 x (σ X P) with wavefunction correction Integrand 100 x (σ X P) without wavefunction correction The ratio of the longitudinal to the transverse part of the cross section R = σ L σ T was analyzed as well. This ratio has a significant sensitivity to the exact form of the wave functions used. In figure 13 the calculation is compared with the experimental data as a function of Q 2 . The data sets shown in figures are for very wide bins on W , and therefore we have shown the curves which correspond to the middle value of the bins.
By inspecting the formulae for the transverse and longitudinal cross sections one would think that the ratio should be approximately independent of the energy W . This is not entirely true as the longitudinal and transverse components of the cross sections are sensitive to somewhat different distributions of the dipole size configurations in the photon wave functions. It is well known that the transverse photon wave function takes more contributions from the large dipole sizes, and therefore one should expect that the energy dependence of σ T should be flatter than that of σ L . Experimental data show however that the ratio R does not depend on the W indicating that the large dipole components may be suppressed in the transversly polarized exclusive vector meson production. The overall description is very good with the exception of the φ data where the theoretical curves indicate the growth with the energy which is slightly too fast as compared with the data. Overlap function 500 x P with wavefunction correction Overlap function 500 x P without wavefunction correction Integrand 100 x (σ X P) with wavefunction correction Integrand 100 x (σ X P) without wavefunction correction This trend is similar to what was observed in [10] which indicates a generic feature of the calculations based on the dipole model.

Differential cross section and t-distribution in vector meson production
The t-distribution of the differential cross section for elastic vector meson production is key in unraveling the impact parameter profile of the target at small-x. The t-dependence can be related to the impact parameter dependence via two-dimensional Fourier transform of the amplitude as indicated in eqs. (2.6), (2.8). The differential cross section is usually parameterized as dσ dt ∝ e −B D |t| in bins of Q 2 and W . The dimensionful slope parameter B D thus contains the information on the spatial distribution of the interaction region in the scattering process. Three plots in figure 14 show the dependence of the slope parameter on the variable Q 2 + M 2 V for ρ, φ, and J/Ψ. The theoretical curves follow the trend of the experimental data. We observe that for the ρ production the dependence of B D on Q 2 is well described but the normalization is underestimated, which is most probably related to the lower normalization for the resulting integrated cross section. The decrease of the slope for low values of Q 2 + M 2 V is related to the initial dependence on the size of the vector meson.  Figure 11. W dependence of the vector meson cross section for elastic production of φ. The experimental data are from [3,6].
large values of Q 2 indicates that in this regime the B D indeed characterizes the size of the proton through the interaction with the small probe which is high Q 2 dipole. This characteristic size of the gluon density inside the proton r 2 ∼ 0.6 fm is markedly smaller than the electromagnetic radius which is of the order ∼ 0.8 fm. This indicates that the gluon distribution differs from the spatial extension of the quarks in the proton. In figure 15 we show the same quantity B D but as a function of W for two different values of Q 2 for both J/Ψ and ρ. While the error bars on the experimental data for B D are relatively large, we see that the theoretical curves describe very well the increasing trend of the data, which is especially visible in the bin for lower value of Q 2 (in fact the bin with higher Q 2 is consistent with the flat dependence as well). In the calculation presented, the energy dependence of the slope is naturally obtained from the evolution of the dipole scattering amplitude with the energy. The diffusion of the dipoles in the impact parameter space is what provides the change of B D with energy. Since it is encoded in the BK evolution it naturally leads to the broadening of the impact parameter profile with the energy. The normalization and slope (in energy) of the B D , however, do depend on the two free parameters in the calculation which are not calculated from first principles. The intercept depends on the value of B G which is used in the initial condition provided by Glauber-Mueller formula, eqs. (2.12), (2.13). The slope in energy of B D is, on the other hand, controlled by the value of the mass m in the dipole evolution kernel which cuts off the large dipole sizes. The B D parameter for J/Ψ is very well described both in normalization and in the W dependence. On the other hand the the normalization of ρ is once again low, nevertheless the dependence in W is well described. We note that the value of B G was fitted from the normalization of the B D slope of J/Ψ and that the mass m parameter which is related to the W dependence is correlated in the calculation with B G . The comparison of the calculation with the data shows that perhaps the initial size in the dipole scattering amplitude may not be universal between ρ and J/Ψ production or that the wave function of ρ needs to include additional non-perturbative components which will increase the interaction size for this meson.  Figure 15. Dependence of the slope parameter B D versus W for J/ψ and ρ production. Data are from H1 experiment [2] and [7].
Finally we note that, the differential cross section of J/Ψ production compares very favorably with the H1 data [2] both in t and in W dependance, which is shown in figures 16 and 17.

Conclusions
In this paper we have computed the integrated and differential cross section for exclusive diffractive vector meson production in deep inelastic scattering using the dipole model framework for small x. The dipole -target scattering amplitude was obtained from the numerical solution to the impact parameter dependent Balitsky-Kovchegov equation. We have found that the overall description of the data is very good, meaning that the calculation based on BK equation with impact parameter dependence is able to reproduce all of the features seen in the data. Several comments are in order: 1. For the good description of the experimental data within the framework presented it is important to include additional corrections. Non-perturbative modification of the photon wave function at low values of Q 2 was necessary, which enhanced the cross JHEP01(2013)001  Figure 17. Differential cross section of J/Ψ production for a fixed W in bins of Q 2 as a function of momentum transfer |t|. Calculations were done with W = 100GeV and W = 90GeV. The experimental data are from H1 experiment [2].
sections of the dipoles with sizes of the order of the hadronic scales. This is likely related to the fact that the dipole scattering amplitude in our approach has a nonperturbative cutoff on large dipole sizes which is absent in most other approaches, like for example in [10]. Thus the enhancement of the photon wave function compensates for the absence of the very large dipole sizes in the present approach. We stress that the correction is only necessary for low values of Q 2 and for lighter vector mesons. The description of the J/Ψ is not sensitive to this correction.
2. The skewedness effect was included in the gluon density distribution. This distribution is present in the initial conditions for the small x evolution. This had a substantial impact on the normalization of the resulting cross section and helped to bring the calculations to agreement with the experimental data.
3. The BK equation was modified to include physical confinement effects by cutting off large dipole sizes. The parameter r max = 1 m , which sets the maximal size of the JHEP01(2013)001 interaction, together with the initial proton size control the slope of the differential cross section with respect to t as well as its variation with the energy. The presented calculation shows very good agreement with the experimental data on B D , including its W dependence in the case of J/Ψ. The slope of B D is reproduced for ρ but the normalization remains low. The W dependence is generated dynamically in the dipole evolution. The speed of this increase is controlled by the parameter r max = 1 m which is not calculable from perturbation theory and needs to be adjusted. 4. The calculation presented includes the running coupling in the evolution, but misses other important NLL effects which are known to be non-negligible, see for example [40]. In particular these effects further slow down the growth of the dipole scattering amplitude with the decreasing x. These effects are thus likely to bring the calculation to a better agreement with the data, especially as far as the W dependence is concerned. The analysis which includes these effects is thus left for further investigation.