Exclusive photoproduction of vector meson at next-to-leading order from Color Glass Condensate

The exclusive photoproduction of vector mesons ($J/\psi$ and $\phi$) are investigated by taking into account the next-to-leading order corrections in the framework of Color Glass Condensate. We confront the next-to-leading order modified dipole amplitude with the HERA data finding good agreement. Our studies show that the $\chi^2/d.o.f$ from leading order, running coupling and collinearly improved next-to-leading order dipole amplitudes are 2.159, 1.097, and 0.932 for the elastic cross section, and 2.056, 1.449, and 1.357 for differential cross section. The outcomes indicate that the higher-order corrections have a significant contribution to the vector meson productions and the description of the experimental data is dramatically improved once the higher order corrections are included. We extend the next-to-leading order exclusive vector meson production model to LHC energies by using the same parameters obtained from HERA. We find that our model can also give a rather good description of the $J/\psi$ and $\phi$ data in proton-proton collision at 7 TeV and 13 TeV at LHCb experiments.


I. INTRODUCTION
Perturbative Quantum Chromodynamics (pQCD) predicts that the gluon density inside a hadron grows rapidly with increasing energy (or decreasing Bjorken-x) and saturates eventually at sufficiently high energies, forming a new state of high density gluonic matter called Color Glass Condensate (CGC) [1]. The rapidity evolution of the CGC matter is known to be described by Balitsky-JIMWLK 1 equation [2][3][4][5] whose mean field version is called Balitsky-Kovchegov (BK) equation [6,7]. One of the hallmark of the BK equation is the geometric scaling. It was found that the experimental data of the total cross section of the electron-proton deep inelastic scattering (DIS) at HEAR in small x (x < 0.01) region shows a geometric scaling behavior [8], which give a strong evidence of the CGC theory. However, a study based on the DGLAP evolution also shows geometric scaling behavior [9]. It is hard to distinguish which one (CGC or DGLAP) is the dominant mechanism to govern the evolution of the partonic system. To get more evidences to support the CGC mechanism, a lot of studies were carried out during the past years. On the one hand, a series of improved QCD evolution equations were proposed, such as running coupling BK (rcBK) equation [10,11], and the full next-to-leading-order (NLO) BK equation [12]. On the other hand, the CGC theory has been used to describe the experimental observables, like the proton structure function and the differential cross section for vector meson production, both from inclusive processes [13][14][15][16][17] and exclusive processes [18][19][20] which may offer more evidences for gluon saturation phenomenon.
In the field of CGC studies, the investigation of exclusive photoproduction of vector meson is especially important since it is very sensitive to the small x gluons, thus it can offer an unique approach to probe the gluon saturation [21,22]. In particular, the quarkonia such as J/ψ and φ are of great interest because they can explore not only the perturbative but also non-perturbative regimes. In recent years, these mesons have being investigated experimentally and theoretically. On the experimental side, the exclusive J/ψ and φ photoproductions have been measured by H1 and ZEUS collaborations at HERA [23][24][25][26]. For higher energies, the LHCb collaboration at LHC has released the exclusive J/ψ production data in proton-proton (pp) collisions at √ s = 7 TeV and √ s = 13 TeV [27,28], which enter even small x region (x ∼ 10 −6 ) and provide high precision experimental data to test the gluon saturation physics.
On the theoretical side, the pioneer study of gluon saturation by using diffractive DIS at HERA based on the Mueller's dipole model [29] can be traced back to two decades ago, in which the Geolec-Biernat and Wusthoff (GBW) model was firstly proposed to search the saturation effect [30]. From then onwards, a lot of efforts were devoted to investigate the phenomenon of gluon saturation via diffraction in DIS. A dipole saturation model, which includes an impact parameter dependent, was developed to describe the differential diffractive J/ψ production data at HERA [18,31]. It was shown that the t−distributions of differential cross sections are sensitive to saturation phenomena. An investigation in Ref. [32] showed that a good description of the diffractive DIS data was obtained by combing the dipole model with the Good and Walker picture. It was found that the diffractive observables can help to discriminate in an unique way between the predictions of different models in the saturation region. Based on the framework of the BK equation at non-zero momentum transfer, the authors in Ref. [33] used the momentum transfer q instead of the impact parameter b in the saturation scale to devise an elegant formulism which is particularly convenient for comparisons between theoretical calculations and experimental data since the data are directly measured as a function of t = −q 2 . In order to investigate whether the diffractive photoproduction of vector meson is a sensitive probe of gluon saturation, a systematic study of the vector meson production was performed with two impact parameter dependent models [34], the IP-Sat [31] and b-CGC [35]. The results further confirm that the t-distribution of differential cross sections of vector meson productions provides an unique method to discriminate among saturation and non-saturation models due to appearance of a pronounced dip in the t-distribution [34]. The aforementioned formulism was extended to study the vector meson productions in proton-proton and nucleus-nucleus collisions at LHC energies [20], which demonstrated that gluon saturation models can give a good description of the experimental data qualitatively. However, all the above mentioned saturation models for the vector meson production are based on the leading-order (LO) dipole amplitude which is inspired by the non-linear BK evolution equation at leading logarithmic accuracy in pQCD, which are insufficient for direct applications to phenomenology. It has been found that the evolution speed of the dipole amplitude resulting from the LO BK equation is too fast to give a precisely description of the HERA data, like proton structure functions [36,37].
Over the past decades, it has been shown that the higher order corrections have a significant contribution to the leading order BK equation. It has been shown that the LO evolution kernel is modified by the running coupling effect, which leads to the rcBK equation [10]. We would like to note that although kernel of the evolution equation is modified, the rcBK equation has the same structure as the LO BK equation. Our studies on the analytic solution to the rcBK equation has been demonstrated that the quadratic rapidity dependence in the exponent of the S-matrix in the LO case is replaced by the linear rapidity dependence once the running coupling correction is included, which indicate that the evolution speed of the dipole scattering amplitude is significantly suppressed by the running coupling effects [38,39]. The full NLO corrections, which includes quark loop (running coupling) and gluon loop contributions, to the evolution equation were calculated by Balitsky and Chirilli in Ref. [12], they found that the kernel and structure of the evolution equation are changed by the full NLO effects. Note that the full NLO BK equation is unstable due to a large double transverse logarithm [40], it can be stabilized by the resummation of the double logarithms leading to a collinearly-improved (ci) BK equation [41]. To see the influence of the full NLO corrections on the dipole scattering amplitude, we analytically solved the full NLO BK equation in the saturation region [42]. Our result shows that the rapidity evolution of the dipole scattering amplitude is still suppressed by the full NLO effects, but the evolution speed is rebound as compared to the running coupling case due to a compensation effect made by gluon loops. Furthermore, our recent studies of the dependence of the dipole scattering amplitude on the running coupling prescriptions have been shown that the argument of the coupling has a great impact on the dipole amplitude [43]. We find that the rapidity evolution speed of the dipole amplitude is significantly slowed down by the smallest dipole size running coupling prescription. Some of the above mentioned NLO theories of the high-energy scattering have been directly applied to phenomenology in the inclusive process. The authors in Refs. [44,45] used the dipole amplitude resulting from the rcBK equation to fit the inclusive small x HERA data. They obtained a rather good description of the data since the running coupling effect significantly slows down the growth of the dipole amplitude with increasing energy. Soon after the collinearly-improved BK equation was established, several groups confronted the resummed equation with recent HERA data [36,46,47], they found that the fit is rather successful and shows good stability up to virtualities as large as Q 2 = 400GeV 2 for the exchanged photon. Although the NLO corrections have been proved to be significant in the inclusive process, they are almost not applied to the vector meson production in the exclusive process in the literature. Based on the aforementioned applications of the NLO effects in inclusive process, we are confident to believe that the higher order corrections are also important in the diffractive vector meson production process. In this work, we find that the χ 2 /d.o.f for elastic cross section (2.159) and differential cross section (2.056) in the leading logarithmic approximation are greatly improved after including one of the NLO corrections, running coupling (1.097 for the elastic cross section and 1.449 for differential cross section). This outcome indicates that the higher order corrections play an important role in the quantitative description of the diffractive vector meson production data. To ensure that the outcome is confident, we study the diffractive vector meson production with the dipole amplitude resulting from the collinearly-improved NLO BK equation. The results show that our theoretical calculations are good in agreement with the experimental measurement with χ 2 /d.o.f 0.932 for the elastic cross section and 1.357 for differential cross section. We extend our model which includes higher order corrections, to study the diffractive vector meson productions at LHC energies. We find that the model can give a rather good description of the J/ψ data from 7 TeV and 13 TeV proton-proton peripheral collisions. The predictions of the diffractive φ production are provided with our model for 7 TeV and 13 TeV proton-proton peripheral collisions at LHC, see Table.V. One can find that the model with LO dipole amplitude gives a larger total cross section, which would not be favored by the data, than the one obtained by NLO dipole amplitude. Again, we see that the high order corrections suppress the evolution of the dipole amplitude.

II. EXCLUSIVE PHOTOPRODUCTION OF VECTOR MESON AT NLO IN THE DIPOLE FORMALISM
In this section we give a brief review of the formalism for exclusive vector meson photoproduction in dipole model. We firstly introduce the dipole model for calculation of the vector meson productions at non-zero momentum transfer in the CGC framework. We then present the evolution equations of the dipole amplitude which is a key ingredient in the dipole model. The vector meson wavefunctions, which are also a portion of the dipole model, are given in the last part of this section.

A. Exclusive photoproduction of vector mesons in dipole model at non-zero momentum transfer
In terms of the dipole model, the vector meson production in an exclusive diffractive γ * p → V p scattering can be viewed as three separated subprocesses [48], as shown in Fig. 1. The first subprocess is the formation of a dipole (a quark-antiquark pair) derived from a virtual photon fluctuation. The second subprocess is the interaction between the dipole and the proton via exchanging gluons. The last subprocess is the recombination of the outgoing quark-antiquark pair to produce a final vector meson. Therefore, the scattering amplitude of the diffractive process can be factorized into three ingredients: the photon wave function, the dipole-proton scattering amplitude and the vector meson wave function. Putting all ingredients together, one can write the imaginary part of the scatting amplitude for a vector meson production as where z is the longitudinal momentum fraction of the incoming photon carried by quark, x is the Bjorken variable, and Q 2 is the photon virtuality. The variable q denotes the momentum transfer, whose relationship with the squared momentum transfer is t = −q 2 . The remainder two dimensional vector r and b are the transverse size of the quark-antiquark dipole and the impact parameter, respectively. Ψ is the wave function of the incoming photon, which can be accurately calculated by QED and is well known in the literature [49,50]. Ψ * V denotes the final vector meson wave function, unlike the photon wave function, it has various prescriptions as we shall discuss later at the end of this section. (Ψ * V Ψ) T,L represent the transverse and longitudinal overlap function between the photon and vector meson, respectively. We would like to note that Eq.(1) is a scatting amplitude containing only the forward component. To get the nonforward scattering amplitude, one can multiply the forward wave functions by a phase factor exp[±i(1 − z)r · q/2] as it was done in Ref. [51]. Using this approach and assuming that the S-matrix is purely real (or the amplitude is purely imaginary), the scatting amplitude can be written as where T (x, r, b) = 1 − S(x, r, b) describes the scattering amplitude between the dipole and proton, which contains all the basic information about the strong interactions between the dipole and proton. By taking into account the corrections from the real part of the scattering amplitude and the skewness effect, the differential cross section of an exclusive vector meson photoproduction can be written as [18]: where β is the ratio of real part to imaginary part of the scattering amplitude and the factor (1 + β 2 ) is to include the correction from the missing real part of the scatting amplitude due to the amplitude, A γ * p→V p T,L , in Eq.(2) only considering the contribution from the imaginary part. The skewness effect factor R g is derived from the fact that the momentum fraction of the exchanging gluons between the proton and dipole legs can be different. The parameters associated with these two corrections can be expressed by the imaginary part as follows: As we know that the dipole-proton scattering amplitude comes from the solution to the evolution equations, like the IIM model [52] inspired by LO BK equation. In most cases, the impact parameter dependence is disregarded in the BK equation, since it was found that the dipole amplitude develops a power-like b behaviour, called Coulomb tails, which render an unphysical results, i.e. the total cross section violation of the Froissart unitarity bound. To avoid the above mentioned difficulty, a general strategy is to build an impact parameter independent dipole amplitude inspired by the BK equation, then a model is employed to include the impact parameter dependence, such as two typical models IP-Sat [31] and b-CGC [35]. In this work, we use almost the same scheme as just mentioned but with impact parameter independent dipole amplitude resulting from numerical solution to the LO, rc, ci BK evolution equations. We introduce the impact parameter via multiplying the numerical dipole amplitude with a Gaussian b dependence. In view of the advantage of the method 2 , which was proposed in Ref. [33] by Marquet, Peschanski and Soyez (MPS), in study the t-distribution of differential cross sections of photoproduction of vector mesons, we shall follow the MPS strategy in this study. Following Ref. [33], the dipole-proton scattering amplitude can be rewritten in terms of the momentum transfer q instead of the impact parameter b by using the Fourier transform Substituting Eq. (6) into Eq. (2), the scattering amplitude for γ * p → V p exclusive diffractive process becomes For the Fourier-transformed dipole-proton scattering amplitude T (x, r, q), we adopt a generalized formalism where the factor e −Bq 2 comes from the nonperturbative effects, R can be interpreted as the radius of proton and N (r, x) is an impact parameter independent dipole amplitude. We would like to point out that B and R are free parameters in our fit, which shall be determined by fitting to HERA data.

B. The dipole evolution equations
A key ingredient to calculate the differential cross section is the dipole-proton scattering amplitude. It is known that almost all the past studies on the differential cross section of vector meson production in the framework of CGC were hovered on the LO level in the literature [53,54]. Although the LO dipole amplitude can describe the diffractive vector meson production experimental data at HERA under certain uncertainties [33,55], the precision of the model has to be improved to distinguish the dynamic mechanism of the CGC evolution from the DGLAP evolution, since the DGLAP formulism can also give a good description of the data [34]. Indeed a lot of efforts have been made to improve the accuracy of the CGC theory by including other higher order corrections, such as quark loops [10,11], gluon loops [12] and pomeron loops [56]. It was found that the running coupling effects dramatically slow down the evolution of the gluon system, which give a good description of the latest data from HERA on reduced cross sections [45]. Similarly, the direct numerical solution of the full NLO BK equation is also shown that it slows down the evolution [40]. Based on the significance of the NLO corrections, we would like to extend LO vector meson production formalism to the NLO in this work. One shall see in the next section that the descriptions of the experimental data are dramatically improved once the NLO corrections are included.
The LO BK equation describes the evolution of a quark-antiquark (with a quark at x ⊥ and an antiquark at y ⊥ ) dipole with the rapidity Y by the emission of a soft gluon. In large N c limit, it can be written as with the evolution kernel whereᾱ s = α s N c /π. Here z ⊥ denotes the transverse coordinate of emitted gluon in the evolution. In Eq.(9) we have used the notation r = x ⊥ − y ⊥ , r 1 = x ⊥ − z ⊥ and r 2 = z ⊥ − y ⊥ to denote the transverse size of parent and the new daughter dipoles, respectively. The BK equation is obtained at leading logarithmic approximation, it has been found that it is insufficient when confronting with experimental data [44,45,57]. Therefore, a lot of efforts have been made to improve the understanding of the dipole's evolution at NLO accurary.
The first improvement to the LO BK equation was done by including quark loops. After resumming α s N f to all orders, one can get an evolution equation with running coupling corrections [10,11], which is called as rcBK equation. The rcBK equation is given by with a modified evolution kernel The numerical solution of the rcBK equation has been obtained by Albacete et al. [44,45], they found the proton structure function can be well described under this evolution equation. However, the quark loops corrections are not the only source of the higher order corrections, the complete NLO corrections should also include the contributions from gluon loops and the tree gluon diagrams with quadratic and cubic nonlinearities [12]. Considering all these contributions, one can get the full NLO BK evolution equation where the kernels are In Eq.(13) we have used the notation r 1 = x ⊥ − z ⊥ , r 2 = y ⊥ − z ⊥ , and r 3 = z ⊥ − z ⊥ to denote the transverse size of dipoles. From Eq. (14) we can see there is a double logarithmic term ln r 2 1 r 2 ln r 2 2 r 2 in the evolution kernel, which renders the full NLO BK equation unstable [40]. The solution can turn to a negative value for some region due to the double logarithmic term. So, one needs to make a resummation of these double logarithms under the double logarithmic approximation (DLA) as it has done by Iancu et al. in Ref. [41]. When this resummation is applied to the full NLO BK equation, the double logarithmic term is removed from kernel K 1 , and the resummation will modify kernel K 1 by multiplying it with kernel with ρ = ln r 2 1 r 2 ln r 2 2 r 2 . In addition to the double logarithmic term, the single transverse logarithms (STL) will also generate a large logarithmic corrections to the evolution equation as shown in Ref. [36]. The effect of the single transverse logarithm resummation will also modify kernel K 1 by multiplying it with kernel with anomalous dimension A 1 = 11 12 . By resumming the large single and double transverse logarithms as in Ref. [36], the collinearly-improved version of BK evolution equation reads where the collinearly improved kernel in the first integration term becomes Note that the Eqs.(9), (11) and (19) shall be numerically solved, and their solutions shall be used as dipole amplitudes to calculate the elastic and differential cross sections in the next section.

C. The wavefunctions for vector meson
Another ingredient to compute the differential cross section for vector meson production is the overlap function (Ψ * V Ψ) T,L , which depends on the quark momentum fraction z, the dipole transverse size r and the photon virtuality Q 2 . The overlap function has various prescriptions, such as boosted Gaussian, Gauss-LC and DGKP [18]. Although it has been shown in Ref. [33] that for an identified meson not all the overlap functions can give equally good description of the experimental data, namely a meson has its own favorite wavefunction. We focus on studying the higher order effects for vector meson production in this study. So, we shall use an unified formalism of wavefunction for different mesons to have a better insight into the higher order effects. The overlap function between the photon and the vector meson has transverse and longitudinal components and can be written as [18] In Eqs. (21) and (22), φ(r, z) is the scalar function. In our study the boosted Gaussian scalar functions are used since it works well for both light and heavy mesons [58]. In boosted Gaussian formalism, the scalar functions are given by The variable in the Bessel functions in Eqs. (21) and (22) The values of the parameters M V , m f , N T,L , and R T,L in the above equations are given in Table I. It is worth to note that the longitudinal component is ignored in most studies due to its small contribution [54,59]. It is safe in a very small photon virtuality regime, like quasi-real photoproduction. However, the longitudinal component can give a significant contribution in a large photon virtuality region as we have discussed in Ref. [60]. So, the longitudinal component is included in this study since we confront with data at variety photon virtualities.

III. NUMERICAL RESULTS
In this section, we use the dipole amplitudes, which come from the numerical solutions to the LO, rc and ci BK evolution equations, to calculate the vector meson productions. Firstly, we give a brief description about numerical method to solve differential equations and experimental data sets used in our fit. Then, we show our theoretical calculations of J/ψ and φ productions and compare them with the experimental data from HERA. Finally, we extend the formalism to LHC energies and make predictions for the rapidity distributions of J/ψ and φ productions in pp collisions at 7 TeV and 13 TeV.

A. Numerical setup and data selection
As well known that the LO, rc and ci BK evolution equations are integro-differential equations. To get the numerical solutions, we can solve them on a lattice. In this study, we discretize the variable r into 256 points (r min = 2.06 × 10 −9 and r max = 54.6). Throughout this numerical study, the unit of dipole size r is GeV −1 . For the rapidity, the number of points are set to 100 with the step size Y = 0.2. This setup can insure that the grid is small enough for our purpose. To perform the numerical simulations, we employ the GNU scientific library (GSL). The main GSL subroutines we have used are Runge-Kutta for solving ordinary differential equations, the adaptive integration for numerical integrals, and the cubic spline interpolation for interpolating data points.
To solve these integro-differential equations, one needs initial conditions. There are several kinds of the initial conditions in the literature, such as GBW [30] and MV [61] models. The Refs. [44] and [45] are used both of them as initial conditions for the rcBK equation in the fit of the reduced cross sections, they showed that the MV initial condition is much more favorable by the experimental data than the GBW model. So, we adopt the MV model as initial condition in this work [61], with γ = 1.13, Q 2 s0 = 0.15 GeV 2 , and Λ QCD = 0.241 GeV. In our analysis, we use Eq.(3) to fit the differential cross section and the elastic cross section for J/ψ and φ productions. The experimental data are taken from ZEUS Collaboration (J/ψ [24], φ [23]) and H1 Collaboration (J/ψ [26], φ [25]). It should be noted that our studies are within the framework of the CGC which is valid in the range x ≤ x 0 with x 0 = 10 −2 . Therefore, the data points with x larger than x 0 are automatically excluded in the data set. In addition, we exclude the data with large error bars at large photon virtuality. After selection, the total number of data are 177 points, which are used in the fit. For the details, the elastic cross section data of J/ψ and φ productions are 58 and 61 points, and the differential cross section data of J/ψ and φ productions are 24 and 34 points, respectively.   To demonstrate the significance of the high order corrections in the description of the HERA data, one needs to compute the vector meson productions with the LO and NLO dipole amplitudes and compare these calculations. In this study, there are only two free parameters B and R, see Eq.(8), other parameters like, x 0 , C 2 , are directly taken from [44]. Tables II and III show Tables II and III, one can see that the NLO descriptions of the vector meson productions are better than the LO case, which indicate that the NLO corrections play an important role in diffractive process. Especially, the χ 2 /d.o.f resulting from the fit to the elastic cross section, σ, there is a large improvement as compared to the LO description once the NLO corrections are included. By Global analysis, one can see that the values of χ 2 /d.o.f calculated from the rc and ci dipole amplitudes are more close to unit than those calculated from the LO amplitude. Figure. 2 shows the elastic cross sections σ for J/ψ and φ productions as a function of photon virtuality Q 2 . The blue, red and black lines represent the results calculating by using the LO, rc and ci dipole amplitudes, respectively (similarly hereinafter). For each of mesons, we consider the experimental data both from H1 and ZEUS Collaborations. For J/ψ, one can see that the higher order dipole amplitudes are good in agreement with the experimental measurement. It can be seen that for φ production it seems that all the dipole amplitudes give a similarly good description of the data in moderate Q 2 , however only the rc and ci amplitudes can give a precise description of the data in low Q 2 . From Fig. 2 it is almost clear that the NLO amplitudes are more favored by the experimental data.
The elastic cross section σ for J/ψ and φ productions as a function of photon-hadron center of mass energies W γp at different photon virtuality Q 2 are shown in Fig. 3. It is shown from the left panels of Fig. 3   that the theoretical calculations from the NLO amplitudes are more consistent with the J/ψ data. From the right panels of Fig. 3, one can see that the LO calculations have rather poor description of the experimental data, while the NLO computations give a relatively good description of the data although the quality is not as good as the J/ψ case, since the experimental data for φ meson have large uncertainties. From Fig. 3, we can see that the NLO calculations have a better agreement with experimental data than the LO BK equation both for J/ψ and φ.
The differential cross section dσ/dt for J/ψ and φ as a function of the squared momentum transfer t at different photon virtuality Q 2 are shown in Fig. 4. From Fig. 4, it seems that the LO and NLO calculations give a similar quality description of the experimental data. It because of the data set is small and with large error bars. However, one can clearly see from the last column in Table III that the χ 2 /d.o.f computed from the NLO dipole amplitudes are much smaller then the ones from LO cases, which indicate that the NLO corrections take an effective role in the diffractive vector meson productions.
Furthermore, it should be noted that there is not significant difference between the description of the experimental data from the rcBK and ci BK calculations. To better interpret the underlying reasons, we have plotted the dipole amplitude, N (r, x), as a function of dipole size, r, for three different rapidities in Fig. 5. From Fig. 5, we find that the difference between the LO and the NLO (rcBK and ci BK) dipole amplitudes is obvious. But, the difference between the amplitudes from rcBK and ci BK equations is tiny up to large rapidities, i.e. Y = 5. While the largest available rapidity at HERA is about Y = 5, therefore it is almost impossible to discriminate the NLO running coupling effect from collinear resummations with current HERA data. This is the reason why we cannot see a remarkable difference between χ 2 /d.o.f resulting from running coupling and resummation improved dipole amplitudes.

C. Predictions for LHC
It believes that the experimental data from LHC offer a peculiar way to test the hadronic structure since the higher energy collision will touch even small-x region. As it was shown in Fig. 5, one can see the difference between the dipole amplitudes from rc and ci BK equations at the larger rapidities, Y > 5 (smaller-x region). So, the predictions for the LHC energies are meaningful as higher precision and rapidity data will be released by LHCb collaboration.
In the high energy proton-proton collisions, there are events involving interacting at large impact parameters where the electromagnetic interaction is dominant. In these photon-induced processes, the two protons keep intact after interaction. For the total cross section, it can be written in terms of a convolution of the equivalent photon flux and the photon-proton production cross section. Therefore, the rapidity distribution for the exclusive vector meson production is given by where y is the rapidity of the produced vector meson, σ γp→V p is the total photon-proton cross section, and ω is the photon energy (ω L = M V 2 exp(−y) and ω R = M V 2 exp(y)). Note that there are two terms in the right hand of the rapidity distribution equation. It is because the photon can be emitted either from left proton or from the right proton.
In Eq. (26), dN dω is the equivalent photon spectrum of the relativistic proton. In the Weiszäcker-Williamsis approximation it can be written as where ξ = 1 + [ (0.71 GeV 2 )/Q 2 min ] with Q 2 min ≈ (ω/γ L ) 2 at high energy limit, √ s is the proton-proton center of mass energy, γ L is the lorentz factor.
Moreover, the total photon-proton cross section σ γp→V p can be integrated from the differential cross section in Eq.(3). The integral over t can be rewritten as follows Using the above formalism and the parameters obtained from fitting the HERA data, we can predict the rapidity distributions for diffractive J/ψ and φ productions in proton-proton collisions at LHC energies. Figures. 6 and 7 show our predictions for the rapidity distributions of exclusive J/ψ and φ in proton-proton collisions at 7 TeV and 13 TeV, respectively. For the LHC kinematics region, there are possible rapidities whose corresponding Bjorken-x can be larger then x 0 for one of the proton but still smaller then x 0 for the other proton. In order to get a smooth curve, we make a linear extrapolation for the dipole amplitudes when x > 0.01. We consider three kinds of dipole amplitudes (LO, rc, ci amplitudes) to calculate the rapidity distributions for exclusive vector meson productions and compare with the released data from LHCb [27,28]. The numerical results in Figs. 6 and 7 show that the NLO dipole amplitudes give a better agrement with experimental data points as expected. For completeness, we present our predictions of the total cross section with different kinds of dipole amplitudes in Tables IV and V. From the tables, one can see that the production rates of the vector mesons (J/ψ and φ) are suppressed by the NLO effect, which satisfy with the theoretical execrations. In conclusion, we have investigated the exclusive vector meson photoproduction for J/ψ and φ at HERA in the framework of color glass condensate. By comparing the results from the rcBK and ci BK equations with those from the LO BK equation, we find that the results from NLO equations are more consistent with experimental data than LO BK equation. We also present our predictions for the rapidity distributions in pp collisions by using parameters obtained from fitting the HERA data. These results indicate that the NLO effects are significant in the calculation of vector meson production at LHC energies. Furthermore, the  higher order corrections considered in this work are part of NLO corrections to the BK evolution equation. As we have studies in Ref. [42] that the rare fluctuations also have a large corrections to the evolution equation once the gluon loop contributions are included into the rcBK equation. Therefore, the exclusive vector meson production with a rare fluctuation corrections are worth to explore in the next work.