Complete calculation of exclusive heavy vector meson production at next-to-leading order in the dipole picture

Exclusive production of transversely polarized heavy vector mesons in deep inelastic scattering at high energy is calculated at next-to-leading order accuracy in the Color Glass Condensate framework. In addition to the first QCD correction proportional to the strong coupling constant $\alpha_s$, we systematically also include the first relativistic correction proportional to the heavy quark velocity squared $v^2$. When combined with our previously published results for longitudinal vector meson production at next-to-leading order accuracy, these results make phenomenological calculations of heavy vector meson production possible at the order $\mathcal{O}(\alpha_s v^0, \alpha_s^0 v^2)$. When applied to $\mathrm{J}/\psi$ and $\Upsilon$ production at HERA and at the LHC, a good agreement between the next-to-leading order calculations and experimental data is found. Additionally, we demonstrate that vector meson production can provide additional constraints compared to structure function analyses when the nonperturbative initial condition for the Balitsky-Kovchegov evolution equation is extracted.


Introduction
In Quantum Chromodynamics (QCD), emission of small momentum fraction x gluons is preferred, which renders parton densities very large when probed at small momentum fraction x in high-energy collider experiments [1,2]. Consequently, parton densities eventually become so large that the smallness of the QCD coupling α s is compensated by the gluon density, and non-linear dynamics starts to dominate in the hadron wave function [3].
An especially powerful probe of non-linear QCD dynamics is given by exclusive vector meson production in deep inelastic scattering (DIS) experiments. The exclusivity of the process requires that, at leading order, two gluons are exchanged with the target hadron at the amplitude level. Thus the cross section approximately scales as gluon density squared [4] (at next-to-leading order in collinear factorization the relationship is less direct, see Ref. [5]) An additional advantage is that only in exclusive scattering processes it is possible to determine the total transverse momentum transfer ∆, which is a Fourier conjugate to the impact parameter, and as such the spectra are sensitive to the spatial distribution of color charge in the target color field [6].
In the next decade, the Electron-Ion Collider in the US [31,32] and other potential nuclear DIS facilities [33,34] will provide vast amounts of precise data on exclusive vector meson production over a wide kinematical domain. To take advantage of these recent and future developments that provide a unique access to non-linear QCD dynamics at small x, it is important to develop theoretical calculations to the comparable level of accuracy.
To describe QCD in the high energy (and density) regime, we employ the Color Glass Condensate (CGC) effective field theory approach [35,36]. In this formulation the color field of the target is written in terms of Wilson lines that describe an eikonal propagation of partons in the color field, resumming multiple interactions. The purpose of this work is to present the first next-to-leading order (NLO) calculation of transversely polarized exclusive heavy vector meson production cross section at high energy within the CGC framework.
At high energy the scattering process is conveniently described in the dipole picture where the (virtual) photon splits into a quark-antiquark dipole long before the interaction with the target (see discussion in Sec. 2). The dipole then interacts with the target and eventually forms the bound state. In order to develop the CGC calculations to NLO accuracy, all ingredients (the virtual photon and heavy vector meson wave functions, and the dipole-target scattering amplitude) are needed at this order in perturbation theory. In recent years there has been a rapid progress in the field to achieve this. The evolution equations at NLO, describing the center-of-mass energy dependence of the dipole-target scattering amplitude, are derived and solved in Refs. [55][56][57][58][59][60] (and a subset of higherorder corrections are resummed in [61][62][63][64][65]). The initial condition for the perturbative evolution is fitted to HERA structure function data at NLO accuracy in Ref. [66] (see also Refs. [67,68] for an NLO calculation of the proton color charge correlations at moderate x that can potentially be used to initialize the evolution). The NLO light-front wave function for a virtual photon was first derived in the massless quark limit in Refs. [69][70][71][72] and recently the results with finite quark masses have also become available [73][74][75]. The NLO wave functions exist also for heavy [76] and light [77,78] vector mesons. In addition to structure functions and exclusive processes, NLO calculations for dijet production in DIS and hadronic collisions [79][80][81][82], and inclusive particle production in proton-nucleus collisions [83][84][85][86][87][88][89][90][91] are becoming available.
The present paper completes the calculation of exclusive heavy vector meson production at next-to-leading order that we initialized in previous publications. First, the relativistic corrections suppressed by the squared quark velocity ∼ v 2 α 0 s were determined in Ref. [92]. Later, in Ref. [93] we calculated the next-to-leading order corrections ∼ v 0 α s to longitudinally polarized heavy vector meson production. This paper presents the calculation of transversely polarized heavy vector meson in virtual photon-target scattering at the order α s v 0 , and demonstrates how the NLO corrections and the relativistic corrections can both be included consistently. This development enables us to present the first calculation of exclusive J/ψ production at the order O(α s v 0 , α 0 s v 2 ), and comparisons with the HERA and LHC data are presented in this paper.
This manuscript is structured as follows. The exclusive vector meson production process in the dipole picture is first presented in Sec. 2. The next-to-leading order QCD corrections for the transverse heavy vector meson production are calculated in Sec. 3. The implementation of relativistic (velocity) corrections is discussed in Sec. 4 before presenting numerical results in Sec. 5 and conclusions in Sec. 6.
2 High-energy scattering in the dipole picture

Exclusive vector meson production
The high-energy limit allows us to describe exclusive scattering in a factorized form where different parts of the process can be written independently. We work in a frame where the photon plus momentum q + is very large and it has no transverse momentum. The splitting of the virtual photon and the vector meson formation are described by the (boost invariant) light-front wave functions of the photon (Ψ γ * ) and meson (Ψ V ). At leading order the only contribution comes from the photon splitting into a quark-antiquark dipole. Additional Fock states have to be introduced at higher orders in α s , and at next-to-leading order one has to include a contribution from the photon splitting into a qqg state. The corresponding NLO scattering amplitude for vector meson production at t ≈ −∆ 2 = 0 can be written as Here x i are the transverse coordinates of the quark (i = 0), the antiquark (i = 1) and the gluon (i = 2), and z i are the corresponding fractions of the photon plus momentum. The different helicity and color components of the wave functions are summed over implicitly. The coherent vector meson production cross section then reads [94] dσ The action of the Wilson lines V (x i ) on the quark-antiquark dipole is given by the dipole amplitude N 01 : The dipole-target scattering amplitude N 01 depends on the transverse separation x 01 = x 0 − x 1 , impact parameter b = (x 0 + x 1 )/2 and projectile evolution rapidity Y . The evolution rapidity depends on the photon-nucleon system center-of-mass energy W as discussed in Sec. 2.2. The notation · · · corresponds to the average of the target color charge configurations, which is done at the amplitude level in Eq. (2.2) when calculating coherent (i.e. no target dissociation) vector meson production [94] (see also e.g. Refs. [41,45,95] related to incoherent diffraction and discussion about the averaging procedure). Similarly, interaction of the qqg system with the target is given in terms of the dipole amplitude N 012 which in the mean field limit can be written as [96] 1 As the impact parameter is conjugate to the momentum transfer in the process, the impact parameter dependence of the dipole amplitudes can be connected to the t-dependence of the production amplitude. However, the impact parameter dependence of the dipole amplitude requires additional modeling and an effective description of confinement effects (see e.g. [39,[97][98][99]), and for simplicity we choose to study only the case t = 0 given by Eq. (2.1) where only the impact parameter integrated dipole amplitude constrained by structure function measurements [66] contributes.
In general, the production amplitude depends on both the polarization of the photon λ γ and the vector meson λ V . The polarization mixing λ γ = λ V is heavily suppressed and consequently it is sufficient to consider only the case λ γ = λ V [100]. Vector meson production can then be divided into longitudinal and transverse production, of which longitudinal channel has already been calculated at NLO by us in Ref. [93]. In this paper we complete the NLO production calculation by computing the transverse production case, allowing us to consider total vector meson production.

High-energy evolution
The center-of-mass energy or, equivalently, Bjorken-x dependence of the Wilson lines can be obtained by solving the perturbative JIMWLK [101][102][103][104][105][106][107] evolution equation. In the large-N c limit one can derive from it the BK equation describing the energy (evolution rapidity Y ) dependence of the dipole amplitude N 01 : The kernel K BK describes the probability to emit a gluon at the transverse position x 2 from the quark-antiquark dipole at the coordinates x 0 and x 1 . Including the running-coupling corrections following [53], the kernel reads The BK equation at next-to-leading order, and a numerical solution to it, are available [55,[58][59][60] (as well as the NLO JIMWLK equation [56,57]). In principle it would be consistent to use the full NLO evolution equation when calculating vector meson production at this order in α s . However, the NLO BK equation is numerically demanding due to an extra transverse integral, which is also the reason why there is no initial condition to it fitted to experimental data. In this work we follow Ref. [66] and use the leading-order BK evolution equation combined with different implementations of a resummation of the most important higher-order contributions. These resummations are known to approximate the full NLO BK equation well as shown in Ref. [59,108]. The initial conditions for these evolutions are determined in Ref. [66] by performing a fit to HERA structure function data [1,2]. In our numerical analysis we use the fit results from publicly available codes [109]. The running strong coupling constant in coordinate space is evaluated using the same parametrization as in the corresponding dipole amplitude fits in Ref. [66]. The explicit expression for the running coupling is with Λ QCD = 0.241 GeV, c = 0.2, µ 0 /Λ QCD = 2.5, β 0 = (11N c − 2N F )/3 and N F = 3, and C 2 is a fit parameter determined when the initial condition for the BK evolution is fitted to the HERA data. The three different schemes to include resummation of higher-order corrections into the BK equation equation used in this work are, following the terminology of Ref. [66], KCBK [110], ResumBK [63,64] and TBK [61]. The evolution rapidity in the KCBK and ResumBK equations is related the fraction of the projectile (photon) plus momentum carrried by the gluon: where k + = z 2 q + is the gluon plus momentum and P is the target momentum. The evolution rapidity in the TBK equation is related to the target longitudinal momentum fraction as we will discuss shortly. The KCBK ("kinematically constrained BK equation") is derived in Ref. [110] by requiring the necessary time ordering between the subsequent gluon emissions. This procedure effectively resums corrections that are enhanced by double transverse logarithm ∼ α s ln x 02 x 01 ln x 12 x 01 . The same logarithms are included in the ResumBK ("resummed BK") equation, with the difference that in Ref. [63] a form of the evolution equation which is local in rapidity Y is derived. Additionally, the ResumBK evolution equation further includes a resummation of single transverse logarithms [64] ∼ α s ln 1 to all orders. For explicit expressions for these evolution equations, see Ref. [66]. The third evolution equation used in this work is the TBK equation ("BK equation in target rapidity"), where the evolution rapidity η is expressed in terms of the fraction of the target longitudinal (minus) momentum transferred in the process x P (see detailed discussion in Ref. [61]): where M V is the meson mass. Consequently the TBK evolution can be thought of as evolution in ln 1/x P , whereas the KCBK and ResumBK evolutions written in terms of the projectile rapidity Y are evolutions in ln W 2 [61,66]. When using a solution to the TBK evolution, written in terms of the target rapidity η, in the NLO impact factors calculated in this work that are written in terms of projectile rapidity Y we use the same shift as in Ref. [66]: where the target transverse momentum scale is set to Q 2 0 = 1 GeV 2 . Initial conditions for all these three evolution equations are obtained in Ref. [66] by parametrizing the initial condition and fitting the free parameters to the HERA reduced cross section data. In this work we use the fit results obtained using the "Balitsky + smallest dipole" running coupling scheme. We note that in Ref. [66] only the light quark contribution is included in the NLO structure function calculations. On the other hand, in this work we consider heavy vector meson production, and as such it is not fully consistent to use the fit results from Ref. [66]. However, the main purpose of this work is to derive the cross section at NLO accuracy, and detailed phenomenological comparisons to experimental data should be done later when the initial condition for the BK evolution is determined including the effect of quark masses.
3 Vector meson production at next-to-leading order Next-to-leading order corrections to exclusive vector meson production consist of corrections from perturbative gluons. These can be included by calculating the virtual photon and meson wave functions at proper order in α s such that we have all the corrections at the order α s at the amplitude level. This means that we have to include the O(α s ) loop corrections to the light-cone wave functions Ψ qq γ and Ψ qq V , and also take into account the contribution from the qqg state with the wave functions Ψ qqg γ and Ψ qqg V . The NLO wave function for the transverse photon with massive quarks has been calculated in Refs. [74,75], and the NLO heavy vector meson wave function in the nonrelativistic limit is evaluated in Ref. [76]. These results are applied in this work.
For completeness, we present here the next-to-leading order wave functions that enter our calculations. Our notation follows mostly Refs. [74,75] with the exception that the integration measure is chosen to be i where i goes over the partons of the Fock state corresponding to the wave function. This introduces additional normalization factors z i compared to the photon wave functions presented in [74,75]. Also, we choose to use the conventional dimensional regularization (CDR) scheme for our calculations, which corresponds to the case D s = D in Refs. [74,75].
The wave functions contain divergences that need to be regularized. Ultraviolet (UV) divergences are regularized using dimensional regularization in D − 2 dimensions for the transverse coordinates. Infrared (IR) divergences originating from gluons with zero plus momenta are removed by introducing a cut-off α for the gluon plus momenta, k + 2 > αq + where α > 0 and q + is the plus momentum of the photon. The divergences will cancel in the calculation, and at the end we will take the limit D → 4 and α → 0.

Virtual photon wave function at next-to-leading order
With these conventions, the LO transverse photon wave function for the qq state (at zero photon transverse momentum, q = 0) is [75] Ψ γ * →qq where u(0) and v(1) are spinors corresponding to the quark and antiquark, m q is the heavy where Q 2 is the photon virtuality, and α 0 , α 1 are the color indices of the quark and antiquark. The fraction of the photon plus momentum carried by the quark is z 0 . Repeated indices in the Latin alphabet are summed over in D − 2 transverse dimensions. The functions K ν are modified Bessel functions of the second kind. The NLO correction to the qq wave function is [74,75] where the wave function is written in terms of different form factors. It will turn out that we do not need to know the explicit expressions for all of these form factors, which follows from the fact that at leading order in α s and v the vector meson spin structure is very simple, picking up only parts with an odd number of transverse gamma matrices. In addition to this, the traceless part P i P j P 2 − 1 2 δ ij also vanishes as after taking the gamma matrix traces one gets j λγ i * λ V P i P j P 2 − 1 2 δ ij = 0 which is valid for λ γ = λ V in 2 transverse dimensions. In the polarization mixing case, λ γ = λ V , there would be a non-zero contribution from this term. Thus, only the last term in Eq. (3.2) contributes to vector meson production in the nonrelativistic limit, and the required combination of the form factors reads Here µ is the mass scale coming from dimensional regularization. The functions Ω T V and L are defined as where Li 2 is the dilogarithm function and The function I T VMS (r, z) can be written in the form where the substitution (z ↔ 1 − z) corresponds to the whole expression.
In addition to the qq wave function of the photon, we also need the light-front wave function for the qqg state (again for the photon with zero transverse momentum). This can written as [75] where σ is the gluon helicity and (3.10) The special functions use the following labeling for the subindices (and analogously for the special functions with transverse indices), where (3.14) The special functions I are defined as: Contribution of the leading-order wave function φ qq to the total vector meson wave function.
The coefficient function C qq←qq in terms of Feynman diagrams. The self-energy corrections are not shown. (3.20)

Meson wave function at next-to-leading order
For the heavy vector meson wave function we can use the nonrelativistic expansion developed in Ref. [76]. In this expansion, corrections in α s are included in factors multiplying the leading-order wave function, and corrections suppressed by the heavy quark velocity v appear as derivatives of the leading-order wave function. For a general Fock state n we write this in the following form: where φ m is the leading-order wave function (LOWF) for the Fock state m and the sum goes over all Fock states m. The case m = qq, n = qq is shown schematically in Fig. 1, where the primed indices 0 , 1 correspond to the nonrelativistic quark and antiquark in the LOWF and the non-primed indices 0, 1 correspond to the quark and antiquark in the wave function Ψ qq . It should be noted that the coefficient functions C k n←m depend on colors and helicities of the particles in Fock states n and m, and the sum over these is left implicit. The derivative ∇ is defined in the mixed space as ∇ = (∂ r 1 , ∂ r 2 , (z − 1/2)2m q i), where r i are the components of the transverse separation r, and k = (k 1 , k 2 , k 3 ) is to be understood as a multi-index: The strength of this expansion is that it is now straightforward to include corrections at a given order in α s and v. Higher-order corrections in α s can be calculated perturbatively from Feynman diagrams, and they are defined as parts of the coefficient functions C k n←m . Relativistic corrections in v can be read from the derivatives which amount to a suppression of v |k| .
At next-to-leading order in the nonrelativistic limit, we need the following wave functions [76]: qq←qq is calculated to the order O(g 2 ) and C (0,0,0) qqg←qq to the order O(g) in the strong coupling constant. Only the LOWF for the dominating Fock state qq contributes at this order in the expansion, as soft gluons in the LOWF Fock state would bring additional suppression in velocity (see also discussion in Ref. [76]). Relativistic corrections at leading order in α s are discussed in more detail in Sec. 4.
As the functions C (0,0,0) n←m are fully perturbative (and calculated at NLO accuracy in Ref. [76]), the nonperturbative physics is contained in the constants dz 4π φ qq where h 0 and h 1 are the helicities of the quark and antiquark in the LOWF. The nonrelativistic limit requires that the helicity structure of LOWF is simply δ (h 0 +h 1 )/2,λ V where λ V is the polarization of the vector meson. To write the expressions in a more compact form we extract the helicity and color structure from the LOWF in the coefficients C k n←qq : and write the spin-independent part of the LOWF as dz 4π φ qq (r = 0, z ). The spinors u(0 ), v(1 ) correspond to the nonrelativistic quark and antiquark in the LOWF with z 0 = z 1 = 1 2 . The perturbative coefficients required for the transverse NLO calculation in the nonrelativistic limit are then and The wave function renormalization factor δZ is calculated in the pole-mass scheme and is given by The contribution from an antiquark emitting the gluon is left implicit in Eq. (3.26). Its contribution to the final expression is equal to the quark emitting the gluon, meaning that we can include its contribution by multiplying the result from the quark contribution by two.

Calculation of the next-to-leading order production
The next-to-leading order amplitude (2.1) has contributions from two terms: the quarkantiquark dipole contribution 2 x 0 x 1 The next-to-leading order corrections to the dipole term come from virtual gluon loops. Using the above expressions for the photon and meson wave functions it is possible to calculate these two contributions.
The evaluation of the dipole term yields For simplicity, we have chosen the LOWF to be real. The real emission contribution is where we have set D → 4 wherever possible and b = z 0 x 0 + z 1 x 1 + z 2 x 2 . Both the dipole term (3.29) and real gluon emission contributions (3.31) have divergences in the D → 4 (UV) and α → 0 (IR) limits. However, the UV divergences cancel in their sum. Therefore it is useful to subtract the UV divergent part of the real correction and add it to the dipole term. Note that now the division of the NLO contributions between the two terms is not unique but depends on the chosen UV subtraction scheme. We choose to do this subtraction following the scheme presented in Ref. [96] and used in Refs. [74,93]. In our case, this means that we write: Although the subtraction procedure is not unique, this particular choice has the correct behavior at x 20 → 0 and results in relatively simple expressions. With the subtraction of Eq. (3.32) we can perform the x 20 and z 2 integrals before adding it to the dipole part: (3.34) Adding this to the dipole term we get Note that this expression is now UV finite, allowing us to take the limit D → 4. The UV subtraction also cancelled the dependence on the scale µ introduced by dimensional regularization.
For the real correction, the UV subtracted form is This expression is also UV finite and does not depend on the renormalization scale µ.
The dipole term (3.35) still has an apparent IR divergence coming from the 1 α term. This is related to the fact that the LOWF is also divergent and has to be renormalized. This can be seen explicitly in the NLO equation for the leptonic width which in the nonrelativistic limit is given by [76] As the leptonic width has to be finite, the LOWF has to have an IR divergent part that cancels the 1 α divergence appearing in this equation. One way to account for this IR divergence of the LOWF is to invert Eq. (3.38) and solve the integrated LOWF directly from it, which gives at the order O(α s ) This expression can then used in Eqs. (3.35) and (3.36) to cancel the IR divergence in the dipole term and to connect the nonperturbative integral over LOWF to the leptonic width for which one can use the experimental value in numerical calculations. The dipole part of the amplitude is then 41) and the NLO part, which contains the corrections from virtual gluon loops, is defined as We define this way of renormalizing the LOWF to be the decay width scheme. Note that this renormalization scheme adds the term αsC F 2π ×4 from the equation of the leptonic width to the virtual correction K NLO qq,Γ . We can also renormalize the LOWF in a different way where such an additional term does not appear. This is done by connecting the LOWF φ qq in our regularization scheme to the dimensionally regularized one φ qq DR following Ref. [76]: This also cancels the IR divergence in the virtual correction. Now the dipole part of the amplitude can be written as where the LO part is still given by Eq. (3.41) but the NLO correction is slightly different: We will refer to this renormalization scheme of the LOWF as the wave function scheme.
Here one still has to determine the value of the dimensionally regularized LOWF for numerical calculations. This can likewise be done with the leptonic width, from which we can solve the dimensionally regularized LOWF in the nonrelativistic limit as Note that when the relativistic corrections are taken into account in Sec. 4 the relation (3.46) will be modified. The difference between the decay width and wave function schemes boils down to the location of the NLO correction αsC F 2π × 4 from the equation for the leptonic width, i.e., whether the correction appears in the equation for the virtual correction or as an overall coefficient when solving the LOWF. The difference between the schemes is parametrically of the order O(α 2 s ) which is of higher order than considered here. This difference is also numerically small in realistic kinematics as we will demonstrate in Appendix A.
When choosing the scheme one also has to take into account the running of the coupling constant α s . In the decay width scheme we use the running coupling in the coordinate space α s (x ij ), Eq. (2.7), whereas in the wave function scheme the coupling constant is evaluated at the momentum scale of the decay process. Following Ref. [111], we take this to be the vector meson mass so that the coupling constant in Eq. (3.46) is chosen to be α s (M V ).

Rapidity divergence and the leading-order result
The real correction (3.36) is still IR divergent as can be verified by taking the α → 0 limit for the lower bound of the z 2 integral. This divergence is actually related to the rapidity evolution of the dipole amplitude. The divergent part of the real correction (3.37) is from which we recognize, using Eq. (2.4), the leading-order BK equation (2.5) if the running of the coupling is neglected. Thus at fixed coupling it is possible to combine this divergent part with the leading-order part (3.41) of the production amplitude by taking into account the rapidity evolution of the dipole amplitude. The amount of BK evolution is controlled by the lower limit of the z 2 integral (recall that Y = ln(z 2 q + /P + ) as discussed in Sec. 2.2). In practice one should not take here the α → 0 limit. Instead, the lower limit has to be set to a finite value. The reason is the following: as the invariant mass M qqg of the qqg system goes like M 2 qqg ∼ 1/z 2 at small z 2 , in the limit z 2 → 0 we would have M qqg → ∞. However, in the calculation we employ the eikonal approximation which assumes that M 2 qqg W 2 where W is the center-of-mass energy of the photon-proton system. Therefore we should set the lower limit, denoted by z min from now on, such that the eikonal limit is satisfied. We follow Ref. [66] and choose z min = P + /q + .
The target plus momentum is given by P + = Q 2 0 /(2P − ) where Q 2 0 is again the transverse momentum scale of the target which we have taken to be Q 2 0 = 1 GeV 2 following [66]. The center-of-mass energy of the photon-proton system is then given by N . Using these expressions we get the following condition for z 2 : (3.48) Using this integration limit we find that the z 2 integral of the K sing qqg part corresponds to the evolution of the dipole amplitude from the initial rapidity Y 0 = 0 to the rapidity using the leading order BK equation.
Evaluating the dipole amplitude at this rapidity corresponds to a resummation of large logarithms α s ln W 2 . As parametrically α s ln W 2 ∼ 1, the actual leading-order part of the production amplitude corresponds to the term K LO qq (Y dip ) where we use the rapidity Y dip to evaluate the dipole amplitude: Note that the IR divergent part (3.47) combined with the lowest order part in (3.41) results in Eq. (3.50) (which corresponds to the subtracted scheme of Ref. [70]) exactly only at fixed coupling. As discussed in Sec. 2.2, in this work we use dipole amplitudes that are evolved using running coupling BK equations that include a resummation of most important higher order corrections to all orders and as such approximate the full next-to-leading order BK equation accurately. Consequently the definition of the leading order scattering amplitude is not unique, but we use the definition (3.50) as it naturally includes also the parametrically large resummation effects included in the used BK equations. At next-to-leading order, we choose to evaluate the leading-order part at the initial rapidity Y 0 and let the qqg part take care of the rapidity evolution. This corresponds to the unsubtracted scheme used also in Ref. [66]. For the virtual correction K NLO qq we choose to evaluate the dipole amplitudes at the evolved rapidity Y dip which corresponds to the total evolution range as discussed above (but note that the dependence on the evolution rapidity is formally of higher order in α s ). For the term K qqg in Eq. (3.37) we use the definition Y = ln k + 2 /P + and evaluate the dipole scattering amplitude at the rapidity: Taking the rapidity dependence of the dipole amplitudes into account we can write the final scattering amplitude for the next-to-leading order production amplitude of transversely polarized vector meson in the form Here K LO qq , K NLO qq,Ψ and K qqg are given by Eqs. (3.41), (3.45) and (3.37), and the rapidity values in parentheses correspond to the rapidities at which the dipole amplitudes are evaluated. An analogous result can also be written in the decay width scheme. The initial rapidity is chosen to be Y 0 = 0 following Ref. [66]. The strong coupling constant is evaluated at the distance scale set by the parent dipole, |x 01 | 2 , in the first two terms and by the smallest dipole min{|x 01 | 2 , |x 20 | 2 , |x 21 | 2 } in the last term, consistently with the NLO fit of Ref. [66].
We note that in Ref. [66] the virtual contribution K NLO qq is evaluated at the rapidity Y = ln 1/x bj which in our case would correspond to Y incl = ln 1/x P = Y dip . It is not entirely consistent to use a different evolution rapidity in our NLO calculation of vector meson production and in the fit procedure used to determine the dipole-proton amplitude. However, as in Ref. [93] we choose to use the more natural choice for the evolution rapidity and note that the difference between the usage of the two rapidities is formally of higher order in α s .

Relativistic corrections at leading order
We can use the nonrelativistic expansion (3.21) to include the first relativistic corrections of order v 2 . At leading order in α s , the nonrelativistic expansion reduces to a distributional identity from which it is easy to read off the coefficients C k qq←qq . In general, including terms of order v 2 corresponds to including the terms with |k| = k 1 + k 2 + k 3 ≤ 2. From now on, we will focus on the case where the meson polarization is λ V = +1. The final result, Eq. (4.10), will be the same for both transverse polarizations λ V = ±1. We can then write the relativistic correction at the order v 2 as where we introduced a simplifying notation and r i = (x 01 ) i . Strictly speaking, it is the complex conjugate of Eq. (4.2) that corresponds to vector meson production, but as this quantity is real we have chosen to write it in this form to get rid of the complex conjugates on the meson wave function. Note that we do not have here terms like φ qq +− (1, 0, 1) or φ qq −− (2, 0, 0). The reason for this is that non-dominant spin components also bring additional velocity suppression, giving a total velocity suppression of v |k|+| 1 2 (h 0 +h 1 )−λ V | . This is because in momentum space the meson wave function must have angle-dependence given by (p x ± ip y ) | 1 2 (h 0 +h 1 )−λ V | which then has to be coupled with similar terms when combined with the photon wave function. This explains why the non-dominant spin terms have to go like p meaning that there is an additional velocity suppression of v | 1 2 (h 0 +h 1 )−λ V | . This can be seen explicitly in Ref. [92] where φ λ V =+1 We can simplify Eq. (4.2) by using the spin-parity J P C = 1 −− of the vector meson. Although parity is only a dynamical symmetry in the light-front, one can still use it to derive symmetry relations of the light-front wave function using similar properties such as the socalled mirror parity (see the discussion in Ref. [112] and the references therein). The spin of the meson allows us to write the dependence on the azimuthal angle ϕ ⊥ as φ is the magnetic quantum number and the part φ does not depend on the angle ϕ ⊥ . The C-and P -parities give the requirements . Using these, the relativistic corrections can be simplified as In Ref. [93] where longitudinal production is calculated at NLO, the decay width scheme is used throughout the paper to renormalize the LOWF. This is possible for the longitudinal production even with the relativistic corrections, as it turns out that the leptonic width for a longitudinally polarized vector meson depends only on the fully nonrelativistic part of the LOWF. For transverse production it is no longer possible to use the decay width scheme with the relativistic corrections as then the leptonic width (at leading order in α s ) has the form which also has contributions from relativistic components of the wave function. The leptonic width (4.5) can be calculated using the light-cone perturbation theory (see e.g. Ref. [113]), and Eq. (4.5) is the form one gets without making any assumptions about the meson wave function. In the nonrelativistic limit φ qq ++ ∼ δ(z − 1 2 ), M V = 2m q , Eq. (4.5) reduces to Eq. (3.38) at LO.
For consistency, one should use the same scheme when combining NLO results of transverse and longitudinal production. We choose to use the wave function scheme throughout this paper as then it is possible to quantify the significance of the relativistic corrections without additional complications coming from the scheme dependence. The previously derived longitudinal cross section is presented in the wave function scheme in Appendix B.
Using Eqs. (3.52) and (4.4) requires that we know the nonperturbative constants related to the LOWF. Note that when relativistic corrections are included we cannot use Eq. (3.43) to directly express dimensionally regularized wave function in terms of the leptonic decay width, as in the case of transverse polarization the relativistic corrections contribute to the decay width as can be seen from Eq. (4.5). We instead calculate these constants using the heavy vector meson wave function from Ref. [92] that includes relativistic corrections of order v 2 . This is a convenient choice as this wave function connects the nonperturbative constants in Eq. (3.52) and (4.4) to the universal NRQCD matrix elements at order v 2 . Using this wave function allows us to write where φ RF ( r) is the value of the rest-frame wave function in the position space (see Ref. [92] for the corresponding light-front wave function). One advantage of using this particular choice of the wave function is that it results in the same leptonic width for both the longitudinal and transverse polarizations at the order O(v 2 ), which follows from the spherical symmetry of the wave function in the rest frame. The rest-frame wave function can be connected to the NRQCD matrix elements by [111] φ Numerical values for the NRQCD matrix elements can be obtained from the decay width data. For J/ψ production, we use the values for O 1 V and q 2 V (with their correlated uncertainties) from Ref. [111]. Similarly, for Υ production the matrix elements from Ref. [114] are used. We note that the matrix elements for J/ψ in Ref. [111] are determined using a charm mass m c = 1.4 GeV. On the other hand, in our calculation we use the nonrelativistic value m c = M V /2, where M V is the J/ψ mass, effectively neglecting the quark momentum contribution to the meson invariant mass (similarly the b quark mass is taken to be half of the Υ mass). Consequently, the different mass values result in a difference which is of higher order in quark velocity v, see also the discussion in Ref. [92].
The relativistic corrections to the production amplitude can then be written as The main result of this work, the exclusive heavy vector meson production cross section at the order O(α s v 0 , α 0 s v 2 ) is then Eq. (3.52), to which Eq (4.10) is added.

Numerical results
We show numerical results for the transverse vector meson production amplitude and for the total (longitudinal and transverse) coherent vector meson production cross section in γ * +p scattering at next-to-leading order. We consider separately the fully nonrelativistic limit O(α s v 0 ), and the case where first relativistic corrections are included, O(α s v 0 , α 0 s v 2 ). The longitudinal vector meson production is calculated using the results of Ref. [93] included for completeness in Appendix B.
The numerical calculations are done in the wave function scheme for the renormalization of the LOWF. We note that it would be possible to use the decay width scheme in the nonrelativistic limit, but this would introduce additional scheme dependence when comparing to the case with the relativistic corrections included. The difference between the two wave function renormalization schemes is studied numerically in detail in Appendix A.

Transverse vector meson production amplitude
First we study in detail different contributions to the forward (t = 0) transverse J/ψ production amplitude, Eq. (3.52), in the nonrelativistic limit. This equation is finite and can be directly evaluated numerically. We limit ourselves to the forward production case as we do not want to specify any specific form of impact parameter dependence for the dipoleproton scattering amplitude, and at t = 0 only the dipole-proton amplitude integrated over the impact parameter appears. The dipole-proton amplitude is obtained from Ref. [66] where the assumption is that the impact parameter b integral only results in a constant   factor σ 0 /2 interpreted as the proton transverse area, also determined from the fit to HERA structure function data. Different contributions to the scattering amplitude as a function of center-of-mass energy W and photon virtuality Q 2 are shown in Fig. 2. In these calculations we have chosen to use as the dipole amplitude the fit performed with the KCBK evolution and using an initial evolution rapidity Y 0,BK = 4.61 in Ref. [66]. Note that in Ref. [66] two different initial rapidities for the BK evolution are used (Y 0,BK = 4.61 and Y 0,BK = 0), and the dipole amplitude is frozen in the region Y 0 < Y < Y 0,BK . Results using both of these initial evolution rapidities are shown later in this Section.
The different contributions to the scattering amplitude shown in Fig. 2 are labeled as follows. First, NLO corresponds to the full NLO level scattering amplitude of Eq. (3.52). The leading-order result, obtained using a BK-evolved dipole amplitude evaluated at the rapidity Y = Y dip is labeled as LO(Y dip ) and shown in Eq. (3.50), and LO(Y 0 ) corresponds to the ∼ α 0 s part (first line of Eq. (3.52)) where the dipole amplitude is evaluated at the initial rapidity. The virtual NLO contribution is denoted by NLO dip , but we emphasize that an UV divergence has been cancelled between the real and virtual contributions and as such the division of NLO corrections into the real and virtual parts is not unique. The real gluon emission correction obtained using the UV subtraction scheme used in this work, Eq. (3.36), is shown as NLO qqg (BK) and NLO qqg (no BK). The NLO qqg (BK) term refers to the singular part of the gluon emission contribution (see Eq. (3.47)) which can be included in the BK evolution, and NLO qqg (no BK) corresponds to the remaining pure NLO correction. In this notation the total NLO amplitude can be expressed as The real gluon emission contribution (NLO qqg terms) has a large contribution, which is expected as the BK evolution resumming terms ∼ α s ln 1/x ∼ 1 to all orders should be considered to be part of the leading order result. However, we also find a significant negative NLO correction NLO qqg (no BK) to the BK evolution from the actual NLO calculation where the exact gluon emission kinematics is included. In our UV subtraction scheme the virtual NLO contribution NLO dip is very small. The total NLO correction is significant, about ∼ 50% of the (BK-evolved) LO contribution. These conclusions are valid at all W and Q 2 .
In the longitudinal production case presented in Ref. [93], the NLO corrections were found to be even more significant (∼ 75% of the LO result); however, these results cannot be directly compared as the longitudinal calculation used the decay width scheme for the renormalization of the LOWF as opposed to the wave function scheme. In the decay width scheme the NLO corrections are larger because of the larger scheme-dependent constant in K qq NLO , which is true for both longitudinal and transverse production. It should be noted that this difference in the NLO corrections is compensated by the overall LOWFrelated constant in the amplitude which is smaller in the decay width scheme, bringing the numerical values of the full NLO result in the two schemes closer to each other and thus reducing the scheme dependence.
The term NLO qqg (BK) corresponds to the (leading order) BK evolution, and as such one could also take the leading order scattering amplitude to be LO(Y 0 ) + NLO qqg (BK). As can be seen in Fig. 2, this differs from LO(Y dip ) by roughly a factor of 2. As discussed in Sec. 3.4 this difference would vanish at fixed coupling if the dipole amplitude satisfied the leading order BK equation, but in this work where we use resummed BK evolution equations to approximate the full NLO BK equation it is more natural to use LO(Y dip ) as a leading order amplitude. This choice includes most of the parametrically large resummation corrections to the leading order amplitude and renders the NLO corrections moderate. On the other hand, if one used LO(Y 0 ) + NLO qqg (BK) as a leading order amplitude, the NLO corrections would be dominant and even render the cross section negative at high Q 2 .

Differential cross section at t = 0
Next we calculate the coherent transverse J/ψ production cross section at t = 0. When squaring the scattering amplitude only the genuine NLO corrections are kept, and NNLO contributions proportional to α 2 s are dropped out. In practice, Eq. (3.50) is used as the leading-order amplitude, and when squaring the amplitude its interference with genuine NLO contributions is needed. This NLO correction is obtained from Eq. (3.52) by subtracting the leading-order amplitude (3.50). Similarly, when including relativistic corrections we do not keep the square of the relativistic correction that would be proportional to v 4 (but note that such contribution was included in the numerical results reported in Ref. [93]). Instead, when relativistic corrections are included we add the interference term between the leading-order amplitude, Eq. (3.50), and the v 2 suppressed part of the amplitude, Eq. (4.10).
The differential J/ψ production cross section at t = 0 is shown in Fig. 3a as a function of the center-of-mass energy W , using different fits for the dipole-proton scattering amplitude from Ref. [66]. Results obtained using fits where the initial evolution rapidity is Y 0,BK = 4.61 (or η 0,BK = 4.61 in the case of TBK evolution) are shown as solid lines, and dashed   lines correspond to calculations where fits with the initial rapidity Y 0,BK = 0 (η 0,BK = 0) are used. The LO result is "LO LOBK" which uses the leading-order dipole-proton amplitude from Ref. [42] (we use the fit referred to as "MV e " in [42]) at the rapidity Y = ln 1/x P as this is the rapidity scale used in the LO fit. We see that the NLO corrections reduce the cross section slightly. This is in contrast to what is seen in Fig. 2, and can be explained by the fact that in Fig. 2 the same NLOfitted dipole amplitude was used for both the LO and NLO results. When nonperturbative parameters describing the initial condition for the BK evolution are determined in leadingorder fits such as in Refs. [42,115], they effectively absorb a part of the higher-order contributions.
These results are in line with what has been obtained for longitudinal vector meson production at NLO [93]. The NLO corrections generally change the center-of-mass energy dependence (faster evolution at low W compared to LO, slower or similar to LO at high W ). We also find some deviations between the results with different NLO dipole amplitudes, similarly to the longitudinal J/ψ production case [93]. As all of these dipole amplitudes were fitted to the same HERA structure function data, this deviation shows that vector meson production gives us complementary information to structure function analyses. This will be discussed more in Sec. 5.3.
In Fig. 3b, we show the effect of relativistic corrections at LO and NLO. At small Q 2 , the relativistic corrections reduce the cross section by ∼ 60% at LO and ∼ 40% at NLO and are numerically more important than the next-to-leading order QCD corrections, that in turn have a much larger effect at large Q 2 . The relativistic corrections have a smaller relative effect at NLO than at LO because of the sizable NLO corrections (recall that we do not include O(α s v 2 ) corrections). It should be noted that the relativistic corrections do not vanish at high Q 2 , which is in contrast with longitudinal production where the relativistic corrections are negligible for high photon virtualities [93,116].

Total vector meson production cross section
The transverse vector meson production can be combined with longitudinal production to calculate the total vector meson production cross section σ tot = σ L + σ T . This is a phenomenologically more interesting quantity as most exclusive vector meson production data is measured in terms of the total production cross section.
As discussed above, we only calculate vector meson production at t = 0 (see Eq. (2.1)) as we do not want to specify any particular model for the proton impact parameter profile. In order to obtain results that can be compared with experimental t-integrated cross section measurements, we use the experimentally measured t slopes that allow us to write the γ * + p → V + p cross section as For the J/ψ production the slope b can be written as b = b 0 + 4α ln(W/(90 GeV)), where the experimentally measured values for the J/ψ production are b 0 = 4.15 GeV −2 and α = 0.116 GeV −2 [8]. This allows us to calculate the total t-integrated vector meson production.
The results are shown in Fig. 4 separately for the nonrelativistic case and with the v 2 relativistic corrections, and they are compared to the experimental data measured by H1 [9,14], ZEUS [7,8], ALICE [19,20], and LHCb [22,23] collaborations. However, we emphasize again that the phenomenological analysis here is not fully consistent as we are using the dipole amplitudes extracted from a fit to structure function data where only the light quark contribution is included [66], and a fully consistent setup would require a heavy quark contribution to be included in the structure function calculations also. Consequently, strong conclusions cannot be drawn from these data comparisons.
Keeping this uncertainty in mind, we find that both the W and Q 2 dependence of the experimental data is described reasonably well, especially when the relativistic corrections are included. For the virtuality dependence, the relativistic corrections are important at low Q 2 and the next-to-leading order corrections also modify the dependence on Q 2 slightly, however both LO and NLO results are compatible with the HERA data. Generally, we again find that both the relativistic and NLO corrections can be numerically important and need to be included when considering the J/ψ production.
At W 100 GeV the calculations using dipole amplitudes from BK fits where the evolution starts at the smallest possible evolution rapidity Y 0,BK = 0 (or η 0,BK = 0 in the     : Total exclusive J/ψ production at next-to-leading order, with HERA data from [7][8][9]14], ALICE data from [19,20], and LHCb data from [22,23]. case of TBK evolution) result in W slope which is not compatible with the data. The nextto-leading order corrections also become extremely large, even rendering the cross section negative. As we will demonstrate in Appendix A, the results obtained with Y 0,BK = 0 (η 0,BK = 0) in the low-W region are also sensitive to the wave function renormalization scheme, but this is not the case for the fits with Y 0,BK = 4.61 (η 0,BK = 4.61). We consider this behavior in the low-energy region to be an artifact of the unphysical initial condition obtained in the BK evolution fits in Ref. [66] when the evolution is started at Y 0,BK = 0, in which case there is a long evolution before one enters in the region probed by small-x structure function data. In that case the fit results in unphysical parameters, and especially the anomalous dimension γ is very large at the initial condition 1 . As heavy vector meson production is sensitive to smaller size dipoles than the structure function fitted in [66], similar unrealistically large NLO corrections were not observed in the NLO fit of Ref. [66]. 1 The dipole amplitude behaves as N01 ∼ (x 2 01 Q 2 s ) γ in the dilute region, and large γ corresponds to e.g. negative unintegrated gluon distribution [117].   Consequently, we consider the results obtained with Y 0,BK = 4.61 to be our main numerical results and emphasize that heavy vector meson production data provides additional constraints for the determination of the nonperturbative initial condition for the small-x evolution.
The leading-order result is constant at W < 31 GeV which corresponds to x P > 0.01. This is because in the leading-order calculation the dipole amplitude is evaluated at Y = ln 1/x P and no BK evolution is included in the region Y < ln 1 0.01 in the leading order fit [118] used in this work.
Next we study the longitudinal-to-transverse J/ψ production cross section ratio where e.g. the normalization uncertainty cancels. This is plotted in Fig. 5 where we show the results for the nonrelativistic case (Fig. 5a) and with the relativistic corrections (Fig. 5b), compared to the HERA data [7,9]. Excluding the TBK η 0,BK = 0.00 dipole, the NLO corrections have only a modest effect on this ratio, slightly increasing it in general. In the nonrelativistic limit the differences between the different dipole amplitude fits having initial evolution rapidity Y 0,BK = 4.61 (or η 0,BK = 4.61) are negligible, whereas with Y 0,BK = 0 (η 0,BK = 0) there is some variation. Including the relativistic corrections increases this variation in both cases. The agreement with the experimental data is similar for both the LO and NLO results, and in both cases the ratio seems to be somewhat overestimated. This ratio is sensitive to the form of the wave function (see e.g. Ref [92]), and in particular including the relativistic corrections improves the agreement with the HERA data slightly.
Finally, we consider exclusive Υ photoproduction. As Υ is much heavier than the J/ψ studied above, relativistic corrections become very small and it can be expected to be more sensitive to the next-to-leading order QCD corrections. The Υ photoproduction cross section as a function of the center-of-mass energy W is shown in Fig. 6 and compared with HERA [119][120][121], CMS [122] and LHCb [123] data. The t-integration of the analytic result is done using Eq. (5.2) with the experimentally measured W -independent slope parameter b = 4.3 GeV −2 [124]. The relativistic corrections, calculated using the NRQCD matrix elements from Ref. [114], are indeed small at ∼ 3% level at all W . The next-to-leading  Figure 6: Exclusive Υ photoproduction cross section as a function of center-of-mass energy W compared with HERA [119][120][121] and LHC [122,123] data. Results are obtained using the KCBK evolution equation with initial evolution rapidity Y 0,BK = 4.61.
order contributions are larger, and result in a slower W dependence compared to leadingorder results. Results obtained at leading and next-to-leading order are both compatible with the available data. Again the leading-order result is constant at W < 95 GeV where x P > 0.01.

Conclusions
We have presented the first calculation for transversely polarized exclusive heavy vector meson production at next-to-leading order accuracy in the Color Glass Condensate framework. The main result of this work is Eq. (3.52), which is the scattering amplitude for the transverse vector meson production at NLO in the nonrelativistic limit. We have also presented how relativistic corrections, which are generally as important numerically as the NLO QCD corrections (in the case of J/ψ production), can be consistently included in the NLO calculation. The corresponding part of the scattering amplitude is given in Eq. (4.10). Combined with the NLO calculation for the longitudinal production presented in Ref. [93], the results of this paper allow for phenomenological studies of heavy vector meson production at next-to-leading order accuracy.
The NLO corrections are numerically significant for both the transverse and longitudinal production amplitude. This is largely compensated by the smaller dipole amplitude in the NLO calculation, making the NLO results mostly in line with the LO production and rendering the NLO corrections generally moderate. When the first relativistic corrections are added the agreement of the coherent J/ψ production cross section with the HERA and LHC data is improved, especially at small photon virtualities where the relativistic corrections are larger than the NLO corrections.
If the NLO cross sections are calculated using dipole amplitude fits from Ref. [66] where there is a long evolution before one enters the region constrained by the small-x structure function data (Y 0,BK = 0 or η 0,BK = 0 fits), the NLO corrections become very large at small center-of-mass energies and even result in negative cross sections. However, we also note that the nonperturbative parameters describing the initial condition in these fits are not physically well motivated. Large NLO corrections observed in this case illustrate how heavy particle production is sensitive to different length scales than structure function calculations and can provide additional constraints when the nonperturbative initial condition for the Balitsky-Kovchegov equation is determined.
Now that the results for both longitudinal and transverse [73][74][75] photon wave functions with massive quarks are available, it will be possible to extend the dipole amplitude fits of Ref. [66] to the massive quark case. This will allow for a consistent phenomenological study of NLO vector meson production at t = 0, which was not possible in this paper. Furthermore, as the impact parameter dependence of the gluon structure is directly related to the t dependence of vector meson production, it would be especially interesting to study t-dependent vector meson production amplitudes. This requires additional modeling for the impact parameter dependence of the dipole amplitude, which is the reason it was not considered in this work.
With these possible future developments in mind, the results presented in this paper can be used for extensive comparisons with heavy vector meson production data from HERA [7,8,14] and from the UPC physics program at the LHC [20,22,22,23], along with making predictions for the future EIC. The results can also be extended from proton targets to heavy nuclei by changing the dipole-target scattering amplitude, which enables studies of non-linear QCD dynamics in heavy nuclei at small-x. This is especially interesting given the existing and future data from ultra-peripheral Pb+Pb collisions at the LHC [24][25][26][27][28][29] and possibilities at the future nuclear DIS experiments.

A Dependence on the wave function renormalization scheme
As discussed in Sec. 4, the renormalization of the LOWF can be done in different ways. In this paper we consider two different renormalization schemes, called the decay width scheme (Eq. (3.39)) and the wave function scheme (Eq. (3.43)). The reason for the different schemes is that the decay width scheme is convenient in the nonrelativistic case as then one can use the same running coupling constant (2.7) as in the rest of the calculation when renormalizing the meson wave function. On the other hand the wave function scheme is necessary when considering the relativistic corrections as in that case the decay width scheme is not possible for transverse production. This makes the NLO cross section dependent on the choice of the wave function renormalization scheme. The choice of the scheme appears parametrically at α 2 s , and is thus of higher order than we consider here, but it can still have an effect on the numerical results. This is what we will study in this Appendix.
In the nonrelativistic case we can choose to use either of these two schemes. The decay width scheme is used as described by Eq. (3.39). For the wave function scheme, we can calculate the dimensionally regularized LOWF in Eq. (3.43) from the nonrelativistic limit of the leptonic width using Eq. (3.46). In that case one has to choose the scale at which to calculate the coupling constant. A natural choice is the mass of the vector meson which in the case of J/ψ evaluates to α s (M J/ψ ) ≈ 0.25 [111] (for Υ this is α s (M Υ ) ≈ 0.18 [114]). The difference between the decay width and wave function schemes is then where the NLO contribution to the decay width appears when calculating the NLO production amplitude. In the decay width scheme it is calculated as part of the virtual correction K NLO qq ; in the wave function scheme it appears when we calculate the value of the dimensionally regularized LOWF.
To quantify the effects of the scheme choice, we have evaluated the NLO differential cross section for exclusive J/ψ production as a function of Q 2 and W in the nonrelativistic case using the two different schemes. The ratio of these cross sections is shown in Fig. 7 for both longitudinal and transverse production. We see that in the calculations where a BK evolution starting at rapidity Y 0,BK = 4.61 (or η 0,BK = 4.61 in the case of TBK evolution) is used there is only a small dependence on the scheme, of the order 10%. On the other hand, if the BK evolution starts at initial rapidity 0 the differences between the two schemes become very large at small center-of-mass energies W . However, we note that as discussed in Sec. 5 the initial conditions for evolutions starting at rapidity Y 0,BK = 0 have unphysical features that do not strongly affect the structure function calculations and the fit process of Ref. [66], but have a large effect here as heavy vector meson production is sensitive to smaller dipole sizes.

B Longitudinal vector meson production at next-to-leading order
For completeness, we list here the expressions for longitudinal vector meson production at NLO from Ref. [93]. The NLO production amplitude can be divided into similar parts as in the case of transverse production. In the decay width scheme, the production amplitude     : Scheme dependence of the differential J/ψ production cross section at NLO illustrated as a ratio of the cross sections calculated in the decay width (Γ) and wave function (Ψ) schemes in the nonrelativistic limit.
can be written as where K LO qq is defined in Eq. (3.41), The special functions L(γ, z), I (j) , I (k) , I i (j) and I i (k) are defined in Sec. 3.1. The terms K L and Ω L V can be written as − K 0 (ζ) 1 + 2 z − 1 2 1 + 2γ E + 2 ln(m q |x 01 |) + 2 ln with and I V (c)+(d) (z, x 01 ) = m 2 The coefficient C L m in the above expression reads When considering relativistic corrections, one can no longer use the decay width scheme for transverse production. It is then more consistent to use the wave function scheme which works also with the relativistic corrections. In the wave function scheme, the longitudinal production amplitude can be written as where the LOWF φ qq,L DR is now the dimensionally regularized one from Eq. (3.43), and the virtual correction becomes instead of K NLO qq,Γ (Y dip ). The relativistic corrections to longitudinal production can be written as where the expression was simplified using the identity φ qq,L +− (r, z ) = φ qq,L −+ (r, z ) that follows from the spin-parity of the vector meson as discussed in Sec. 4.
The nonperturbative constants related to the LOWF can be written in terms of the rest-frame wave function φ RF ( r) [92], giving us This allows us to write the v 2 relativistic correction in the compact form The value of the rest-frame wave function and its derivatives can be related to NRQCD matrix elements as described in Sec. 4.