Exclusive heavy vector meson production at next-to-leading order in the dipole picture

We calculate exclusive production of a longitudinally polarized heavy vector meson at next-to-leading order in the dipole picture. The large quark mass allows us to separately include both the first QCD correction proportional to the coupling constant $\alpha_s$, and the first relativistic correction suppressed by the quark velocity $v^2$. Both of these corrections are found to be numerically important in $\mathrm{J}/\psi$ production. The results obtained are directly suitable for phenomenological calculations. We also demonstrate how vector meson production provides complementary information to structure function analyses when one extracts the initial condition for the energy evolution of the proton small-$x$ structure.


I. INTRODUCTION
Deep inelastic scattering (DIS) processes enable precision studies of proton and nuclear structure thanks to the pointlike structure of the electron probe. The vast amount of accurate measurements in electron-proton collisions at HERA has allowed for a detailed determination of the partonic structure of the proton. In particular, a rapid increase of the gluon density towards small momentum fraction x has been observed [1,2]. This rise can not continue indefinitely without violating unitarity, and nonlinear dynamics should eventually start to affect the proton or nuclear structure at small x.
In order to describe the hadron structure at high densities where nonlinear dynamics should be included, an effective theory of quantum chromodynamics known as the Color Glass Condensate (CGC) has been developed, see Refs. [3][4][5] for a review. Despite the success of the CGC framework in describing high-energy scattering measurements, there is no solid evidence that nonlinear saturation effects are visible at current collider energies. In order to access more pronounced nonlinear effects, there are concrete plans to construct an Electron-Ion Collider in the US [6][7][8] with similar longer term plans at CERN [9,10] and in China [11]. As the parton densities are enhanced by approximatively A 1/3 in heavy nuclei, the nonlinear dynamics should be more easily accessible by replacing the proton by a heavy nucleus in these future facilities.
Exclusive vector meson production is an especially powerful tool in probing hadron structure at small-x. In exclusive processes there cannot be a net color charge transfer in the process which, at leading order, requires an exchange of two gluons at the amplitude level. This renders the cross section to be approximatively sensitive to the squared gluon density [12], and nonlinear dynamics is expected to be pronounced in exclusive processes * Electronic address: heikki.mantysaari@jyu.fi † Electronic address: jani.j.penttala@student.jyu.fi with heavy nuclear targets, see e.g. Ref. [13]. An additional benefit of exclusive processes is the fact that only in these events it is possible to measure the total momentum transfer to the target, which is the Fourier conjugate to the impact parameter. This enables studies of the Generalized Parton Distribution Functions (GPDs) [14,15] and the spatial structure of protons and nuclei with event-byevent fluctuations [16,17].
In this Letter we present the first calculation of the exclusive heavy vector meson production at next-to-leading order (NLO) in photon-proton collisions. The advantage of heavy mesons such as the J/ψ is that the mass scale renders the process perturbative even at low virtualities Q 2 , enabling perturbative studies also in ultra peripheral collisions at RHIC and at the LHC where the photons are approximatively real. We include both the first QCD correction ∼ α s , and the first nonrelativistic correction ∼ v 2 . We focus on J/ψ production, in which case both contributions can be expected to be parametrically important as for the charm quark velocity one can estimate v 2 ∼ α s [65]. These developments are crucial to enable precision studies of nonlinear dynamics in exclusive scattering processes in the CGC framework. In this work we consider the case where the photon is longitudinally arXiv:2104.02349v3 [hep-ph] 31 May 2022 polarized, but note that it will be possible to extend the results to the transverse photon case in the future.

II. EXCLUSIVE SCATTERING AT HIGH ENERGY
At high energies, it is convenient to describe exclusive vector meson production in the dipole picture in the frame where a highly energetic photon scatters off the target proton or nucleus and forms a vector meson. The exclusive nature of the process requires that there is no net color exchanged in the process. In this frame, the partonic Fock states of the virtual photon, e.g. |qq dipole and |qqg at NLO, have a long lifetime compared to the timescale of the interaction. These partons propagate eikonally through the color field of the target at fixed transverse coordinates x i and pick up Wilson lines V (x i ) in the fundamental (quarks) or adjoint (gluons) representation. The scattering amplitude can be conveniently written as a convolution of the photon and vector meson wave functions Ψ γ * and Ψ V , and Wilson line operators in the mixed longitudinal momentum fraction z i , transverse coordinate x i space. At high energies the imaginary part dominates, and the scattering amplitude at NLO in the case where the transverse momentum transfer to the target vanishes can be written as Here x 0 , x 1 and x 2 are the quark, antiquark and gluon transverse coordinates respectively, and z 0 , z 1 and z 2 are the fractions of the photon plus momentum carried by these partons. The calculations are done in the target rest frame where the photon plus momentum is chosen to be large. The wave functions depend implicitly on quark and gluon helicities, and a summation over the helicity states is also implicit. We consider coherent scattering processes, where the target proton remains in the same quantum state (i.e. no target breakup). The coherent vector meson production cross section reads The Wilson line operator describing the qq dipole-target scattering known as the dipole amplitude N 01 is where x 0 and x 1 are the quark and antiquark transverse coordinates. Here the average refers to the average of the target color charge configurations. Using the Fierz identity, the operator for the qqg-target scattering can be written as (see e.g. Ref. [50]) The transverse momentum transfer |∆| ≈ √ −t is the Fourier conjugate to the impact parameter. Consequently, an accurate calculation of the cross section differentially in |t| would require a realistic description of the impact parameter dependence in the small-x evolution equation. Such equations exist and have been solved in the literature [29,[66][67][68], but come with the price that one has to model confinement scale effects that suppress long range Coulomb tails. In this work, we want to perform a rigorous NLO calculation and limit our analysis to the t = 0 case in which the diffractive scattering amplitude is only sensitive to the dipole-target scattering amplitude integrated over the impact parameter, which we take to satisfy an impact parameter independent small-x evolution equation. The necessary ingredients in the NLO calculation are the photon and vector meson wave functions and the dipole scattering amplitude, all at NLO level. The wave functions have been recently calculated in the literature, first in the massless quark limit in Refs. [50,51,55], and recently the heavy quark contributions have become available [54,56]. These wave functions are reviewed in Sec. III. The dipole amplitude N 01 = N 01 (Y ) whose energy (or rapidity Y ) dependence is given in terms of the BK equation is also available at NLO accuracy. The BK equation requires a non-perturbative input that can be taken to describe the dipole-target scattering amplitude at initial x, typically around x ∼ 0.01. This perturbative evolution then predicts the dipole amplitude at smaller x (higher energies). In this work we use a BK evolved dipole scattering amplitude with the initial condition fitted to the HERA structure function data [1,2] at NLO accuracy in Ref. [69] (including only light quarks), using the public codes from Ref. [70].

III. VECTOR MESON PRODUCTION AT NEXT TO LEADING ORDER
In order to calculate exclusive longitudinal quarkonium production at next-to-leading order accuracy, light front wave functions for the longitudinally polarized virtual photon Ψ γ * and vector meson Ψ V are needed at this order in α s (see also Ref. [71] for a discussion of negligibly small polarization changing contributions).
The virtual photon light front wave function at NLO accuracy in the case of massive quarks has recently been calculated in Ref. [54]. In our NLO calculation, we need the wave function describing the photon fluctuation to the |qq Fock state at NLO, Ψ qq γ * , and the tree level result for the formation of the |qqg state, Ψ qqg γ * . Detailed expressions can be found from Ref. [54], and we also ex-plicitly show these results in Appendix A 1, equations (A2) and (A3).
For the heavy quarkonium, a systematic method to include both the higher order QCD corrections and the relativistic corrections has been derived in Ref. [56]. In this approach, the quarkonium wave function is written in terms of the coefficients C k n←m defined as are the meson and heavy quark masses, and φ m is the leading order light front wave function which is generally nonrelativistic in the case of heavy quarkonium. The transverse separation between the quark and the antiquark is r, and z is the fraction of the meson plus momentum carried by the quark. Here k = (k 1 , k 2 , k 3 ) is to be understood as a multi-index, and the sum goes over all positive values of the single indices k i , and As we neglect contributions ∼ α s v 2 , we need the quarkonium wave function Ψ n V at NLO in the nonrelativistic limit, which corresponds to the |k| = k 1 + k 2 + k 3 = 0 case. The |k| > 0 terms are suppressed by quark velocity, and are used to calculate relativistic corrections in Sec. IV. The variables m and n refer to the partonic Fock states, and at NLO in the nonrelativistic limit the non-zero coefficients are C (0,0,0) qq←qq and C (0,0,0) qqg←qq . These have been calculated in Ref. [56] and are explicitly shown in Appendix A 2 in Eqs. (A20) and (A21).
With m = qq, the leading order wave function for the longitudinally polarized quarkonium is φ qq h 0 h 1 with the helicity structure given by δ h 0 ,−h 1 . Here h 0 and h 1 are the quark and antiquark helicities. The coefficients C k n←m depend implicitly on the helicities of all the partons in state n, denoted by h i , as well as on h 0 and h 1 , and on parton colors as shown explicitly in Appendix A 2. Summation over helicities and colors is implicit in Eq. (5). In particular, the helicity structure of the Ψ qq V wave function becomes δ h0,−h1 at lowest order in v, where h 0 and h 1 are the quark and antiquark helicities. The helicity flip contribution (where the quark helicities are the same) would only appear in the longitudinal quarkonium wave function Ψ qq V at the order v 4 [72] which is not included in our calculation. Consequently we can also drop the helicity flip component Ψ h.f from the longitudinal photon wave function at NLO.
Using the NLO wave functions available in the literature, we calculate longitudinally polarized exclusive heavy vector meson production amplitude −iA at nextto-leading order in the nonrelativistic limit. In this calculation there are both ultraviolet (UV) and infrared divergences that cancel at the level of the total NLO amplitude. First, both the real qqg and virtual qq contributions contain ultraviolet divergences that cancel in the sum. In practice, we use dimensional regularization and subtract this UV divergence from the real contribution, and add it to the virtual contribution which renders both terms finite. As the UV subtraction term can contain an arbitrary finite piece, division of the NLO contributions between the "real" and "virtual" parts is not unique. In this work, we follow the UV subtraction scheme developed in Refs. [50,54] which differs slightly from the one used in Ref. [56] and results in simpler expressions.
In addition to UV divergences, the NLO amplitude is singular in the infrared where the gluon plus momentum is very small. In order to cancel this divergence, we include two contributions as discussed in Ref. [56]. First, the leading order wave function φ qq h 0 h 1 is divergent at NLO, with the soft gluon divergence regulated by an infrared regulator α. The leptonic decay width Γ(V → e − e + ), however, is finite and connects the wave function and the infrared regulator: [56] Here e is the elementary charge and e f the fractional charge of the quark. In practice the wave function φ qq can be written in terms of the decay width and the infrared regulator, and the 1/α divergence will cancel when combined with the virtual NLO contribution. Additionally, the real gluon emission contribution is also singular in the soft gluon limit. This contribution can be absorbed in the target BK evolution. The final result for the scattering amplitude at nextto-leading order reads is the impact parameter. Detailed expressions for the NLO contributions K NLO qq and K qqg are shown in Appendix A 3, Eqs. (A24) and (A26). This corresponds to the "unsubstracted scheme" discussed e.g. in Refs. [53,69]. Following the same terminology, we refer to the second term in Eq. (7) as the virtual "dipole" contribution (denoted by NLO dip later), and the third term as the real contribution (NLO qqg ). As discussed above, the division of the NLO corrections between these two terms is not unique. The dipole amplitudes are evaluated at evolution rapidities Y 0 , Y dip and Y qqg that are discussed in detail below.

Evolution equations and rapidities
The integral of K qqg over z 2 in Eq. (7) is singular in the limit z min → 0. The singular part is related to the rapidity evolution of the dipole amplitude as can be seen by writing out the singularity explicitly: We can recognize from this the leading-order Balitsky-Kovchegov equation in integral form. It corresponds to the evolution over ln 1 2zmin units of projectile rapidity Y , defined as Y = ln k + /P + . We recall that we work in the frame where the incoming photon has a large plus momentum q + and the gluon plus momentum reads k + = z 2 q + . The target plus momentum P + is obtained as P + = Q 2 0 /(2P − ), where the transverse momentum scale of the target is taken to be Q 2 0 = 1 GeV 2 following [69]. The photon-proton center of mass energy squared is W 2 = 2q + P − .
When the singular part in Eq. (8) is combined with the term K LO qq (Y 0 ), one obtains the leading order contribution but with the dipole amplitude evolved from rapidity Y 0 to rapidity using the LO BK equation at fixed coupling. This evolution is part of the actual leading order contribution, as the BK evolution resums α s ln 1/x contributions, that at high energy are of the order 1, to all orders. In this work we use the dipole amplitudes obtained as a result of the NLO fit to HERA structure function data [69] where Y 0 = 0, and we also use the same running coupling prescription as in Ref. [69] in numerical analysis. We can write the leading order scattering amplitude as We note that what actually is the leading order contribution is not unique. One could as well define it as the sum of the lowest order contribution evaluated at Y 0 and the singular part of NLO qqg . These two definitions match at fixed coupling if no higher order corrections were included in the BK evolution [53]. In the NLO DIS fit of Ref. [69] applied in this work, modified versions of the BK equation that include the most important higher order corrections, in addition to the running coupling effects, were used and consequently Eq. (10) can not be obtained from Eq. (7). The definition of what is considered as the leading order contribution matters, because when calculating the cross sections an interference between the leading order and the genuine next-to-leading order contributions is needed, in which case the NLO correction is obtained as −iA L − (−i)A L LO . Let us then determine the lower limit of the z 2 integral z min which controls the amount of small-x evolution. The applicability of the eikonal approximation requires that the invariant mass of the qqg system M qqg satisfies M 2 qqg W 2 . There is some freedom in determining how strong ordering is required, and resulting differences at the cross section level will again formally be of higher order in α s . In this work we use the same convention as in Ref. [69] and require M 2 qqg < W 2 , which gives In the NLO qqg term the rapidity at which the dipole amplitudes are evaluated depends on z 2 , and using again Y = z 2 q + /P + we get We note that the total evolution range probed in the NLO dip contribution is exactly Y dip discussed above. For consistency, we choose to evaluate the dipole amplitude in the NLO dip term at the same rapidity Y dip as the leading order contribution. In Ref. [69] the virtual corrections to the structure functions were evaluated at Y = ln 1/x bj , which in our case would correspond to the rapidity Y incl dip = ln 1/ is the fraction of the target longitudinal momentum transferred to the meson. Although it is not exactly consistent to use a different scheme to set the evolution rapidity in the structure function fit and in the application of these fit results in exclusive vector meson production, here we choose to apply the more natural choice for the evolution rapidity. The difference between these two choices is formally of higher order in α s . In Ref. [69] the fits are performed using initial conditions parametrized at both at Y 0,BK = 0 and at Y 0,BK = 4.61, in which case there is no evolution in the region 0 < Y < Y 0,BK . The evolution equations that approximate the full NLO BK [43] in the fits are the so-called KCBK, ResumBK and TBK equations (following the terminology of Ref. [69]) derived in Refs. [73][74][75][76]. The "kinematically constrained BK equation" (KCBK) [75] is obtained by explicitly enforcing the required time ordering between the subsequent gluon emissions in the evolution. This procedure effectively resums corrections that are enhanced by two large transverse logarithms ∼ α s ln in the evolution, and the same double logarithms are also resummed in Ref. [73]. When additional contributions enhanced by single transverse logarithms ∼ α s ln 1/(x 2 ij Q 2 s ) (where Q s is the saturation scale of the target) are also resummed following Ref. [74] one obtains the evolution equation referred to as the Re-sumBK equation. The third evolution equation (TBK) refers to the BK equation where the evolution rapidity η ("target rapidity") is related to the fraction of the total longitudinal momentum of the target. When using the fit result that is written in terms of the target rapidity η in the impact factors written in terms of the (projectile) rapidity Y , we apply the same shift as in Ref. [69]: η = Y +ln min(1, x 2 01 Q 2 0 ) . For more details of the different evolution equations, we refer the reader to Ref. [69].

IV. RELATIVISTIC CORRECTIONS
As we have parametrically α s ∼ v 2 , it is also interesting to consider the first relativistic corrections of order v 2 at leading order in α s . Using Eq. (5), we note that each term in the expansion corresponds to a correction of order v |k| . The coefficient functions C k qq←qq are straightforward to calculate at leading order in α s as then the wave function gets no loop corrections (and C k qqg←qq = 0), and we can write where δ (z − 1/2) , and (14) φ qq Here r = (r 1 , r 2 ) is the transverse separation of the two quarks and α 0 , α 1 refer to the quark colors. Calculating the production amplitude at order v 2 corresponds to keeping terms with k 1 +k 2 +k 3 ≤ 2. The nonperturbative constants φ qq h 0 h 1 (k 1 , k 2 , k 3 ) can be related to the derivatives of the leading-order rest frame wave function φ RF at the origin as shown in Ref. [72]. This allows us to write (see discussion in Appendix A 2 for more details) With this the order v 2 correction to the production amplitude is: The value for ∇ 2 φ RF (0) for J/ψ can be determined using the nonrelativistic QCD (NRQCD) matrix elements from [65]: The long-distance matrix elements O 1 V and q 2 V are related to the J/ψ wave function and its derivative at the origin and explicitly defined in Ref. [65]. They are determined using m q = 1.4 GeV for the charm mass. In this work, on the other hand, we use the nonrelativistic limit relation m q = M V /2 also when calculating the relativistic correction. This difference is of higher order in v, see also discussion in Ref. [72].

V. NUMERICAL RESULTS
In this section we present numerical results for the exclusive J/ψ production at t = 0 in the kinematics covered by HERA and future EIC and LHeC/FCChe measurements. Unless otherwise stated, we use the KCBK evolved dipole amplitude with the initial condition parametrized at Y 0,BK = 4.61 from Ref. [69]. The qualitative features do not depend on the actual dipole amplitude fit used.
The scattering amplitudes for exclusive J/ψ production at leading and next-to-leading order are shown in Figs. 1 and 2 as a function of center-of-mass energy W and photon virtuality Q 2 . The total NLO amplitude is shown in Eq. (7), and should be compared to the leadingorder result including the small-x BK evolution defined in Eq. (10) and denoted by LO(Y dip ) in the figures. Note that all results in Figs. 1 and 2 are obtained by using the same dipole amplitude N 01 from Ref. [69]. The NLO corrections are found to be sizeable, of the order of ∼ 75%, and depend weakly on W and Q 2 . Only at highest Q 2 values ∼ 100 GeV 2 (where the high scale renders α s smaller) the NLO corrections become slightly less important.
In Figs. 1 and 2 the different contributions to the NLO amplitude are also shown separately. First, the LO(Y 0 ) curve refers to the leading order result with no BK evolution. The virtual NLO correction NLO dip is found to be small and positive (by positive we mean that it has the same sign as the leading order result) at all W and Q 2 . The real contribution NLO qqg includes a leading order part in terms of the BK evolution. In the figures we  show separately the contribution from the BK evolution shown in Eq. (8), and the genuine next-to-leading order correction to it due to the exact gluon emission kinematics included in the full NLO calculation. This NLO correction NLO qqg (no BK) = NLO qqg −NLO qqg (BK) significantly suppresses the effect of the small-x BK evolution as expected. This systematics in the real and virtual corrections is similar to what is observed in case of structure function calculations at NLO in Ref. [53]. However, we emphasize that the division of the NLO corrections between the NLO dip and NLO qqg terms is not unique, see the discussion in Sec. III. In Fig. 3 we show the energy dependence of the J/ψ electroproduction cross section at NLO using different dipole amplitude fits which describe the HERA struc-ture function data approximatively equally well [69]. For comparison, the LO result using the fit with LO BK evolution and leading order impact factors from Ref. [23] is also shown. In the case of the LO BK evolved result the evolution rapidity is chosen as Y = ln 0.01/x P , consistently with the leading-order fit procedure of Ref. [23].
Despite the fact that all dipole amplitudes result in almost identical descriptions of the HERA structure function data, we find that the resulting J/ψ production cross sections can differ by almost a factor of 2. This demonstrates that vector meson production provides complementary information for the extraction of the initial condition for the BK evolved dipole scattering amplitude, as it is sensitive to the dipole-target interaction at different distance scales (see also Ref. [29,30]) compared to total cross section measurements.
It should be noted that the NLO results for the J/ψ production cross section are closer to the LO BK result than one might expect judging from Figs. 1 and 2. This can be understood by noting that the LO BK result in Fig. 3 is calculated using a LO dipole amplitude resulting from a leading-order fit, whereas in Figs. 1 and 2 the same NLO fit result was used for the dipole amplitude in all cases. In particular, the fit parameters obtained from the LO fit include an effective description of some of the higher-order effects. However, we still find that the NLO cross section generically evolves more slowly in W at high energies compared to LO results.
In Fig. 4 the effect of the relativistic corrections is shown. As previously discussed in Ref. [72], the relativistic corrections are significant and decrease the cross section up to ∼ 50% at low photon virtualities, and become insignificant (but non-zero [77]) at large Q 2 . At low virtualities the relativistic correction is more important than the next-to-leading order contribution. However, when comparing the relativistic ∼ v 2 and NLO ∼ α s corrections one has to keep in mind that the leading order BK evolution effectively includes higher order corrections encoded in the fit parameters as discussed above. The relativistic correction is less important when it is added on top of the next-to-leading order result, ∼ 40% at low Q 2 , as we do not include corrections of the order α s v 2 .
The next-to-leading order correction becomes large at high virtualities as can be seen from the lower panel of Fig. 4. We note that both LO and NLO fits provide a good description of the Q 2 dependence of the HERA structure function data at small x. The stronger virtuality dependence at leading order can be again understood to result from the fact that J/ψ production is sensitive to dipole scattering amplitude at smaller length scales compared to structure functions. The small dipole size region is also only weakly constrained by the structure function data when the initial condition for the BK evolution is fitted.
Technically, the dependence on the virtuality is related to the anomalous dimension γ which describes the behaviour of the dipole amplitude at small dipole sizes: N 01 ∼ (x 2 01 Q 2 s ) γ . At leading order the BK evolution results in γ ∼ 0.7 at large rapidities, but as Y ∼ ln 1/x P in the LO fit, at high Q 2 one is actually sensitive to the dipole amplitude close to the initial condition where γ ∼ 1.2 [22,23]. On the other hand, in our NLO setup there is a long evolution at high Q 2 , see Eq. (12). However, the anomalous dimension at asymptotically small dipoles does not actually change when higher order corrections are resummed in the NLO fit. As the NLO fits also result in γ ∼ 1.2 [69] at the initial rapidity, in principle we would expect to see comparable Q 2 evolution speeds in the exclusive J/ψ production. In practice one is not probing the dipole amplitude at asymptotically small dipoles but at x 2 01 ∼ 1/Q 2 , and in the NLO fits γ decreases in the evolution at intermediate dipole sizes [69]. As a result, one finds that the NLO exclusive vector meson cross section decreases more slowly as a function of Q 2 than the leading order case at high virtualities.

VI. CONCLUSIONS
We have calculated, for the first time, exclusive heavy vector meson production at next-to-leading order in the Color Glass Condensate framework. In the calculation we apply the recently derived wave functions for the virtual photon and vector meson including massive quarks. The main result of this work, the scattering amplitude for longitudinal vector meson production at NLO, is Eq. (7). We emphasize that this result is free from any ultraviolet or infrared divergences and suitable for phenomenological applications.
We have numerically evaluated the derived scattering amplitude, using dipole-proton scattering amplitudes recently obtained as a result of an NLO fit to HERA structure function data. We have presented the first numerical calculation of the exclusive J/ψ production cross section at NLO in the CGC framework. As the future Electron-Ion Collider and other nuclear DIS facilities will provide vast amounts of precise vector meson production data in the future, these developments that promote the CGC calculations to the precision level are extremely important.
We have shown that the next-to-leading order corrections to the J/ψ production cross section are significant, although these corrections can partially be captured in leading order calculations by the non-perturbative fit parameters. We also demonstrate that the vector meson production data provides complementary informa-tion compared to structure function measurements. A global analysis including both the reduced cross section and exclusive vector meson production data would be preferable in the future when extracting the initial condition for the Balitsky-Kovchegov evolution of the dipole scattering amplitude. Comparing the NLO ∼ α s correction to the relativistic ∼ v 2 correction, we have observed that especially at low virtualities both corrections are numerically important, with the relativistic correction generically larger.
In the future, we will include the contribution from the transversely polarized virtual photons. This development will enable comparisons with the vector meson production data from HERA [34][35][36] and from the UPC physics program at the LHC [38,39], as well as calculation of precise predictions for the EIC. Extending the calculation from protons to heavy nuclei will also enable precision studies of saturation phenomena in current [40][41][42] and future nuclear DIS experiments.