Exclusive production of light vector mesons at next-to-leading order in the dipole picture

Exclusive production of light vector mesons in deep inelastic scattering is calculated at next-to-leading order in the dipole picture in the limit of high photon virtuality. The resulting expression is free of any divergences and suitable for numerical evaluations. The higher-order corrections are found to be numerically important, but they can be mostly captured by the nonperturbative fit parameters describing the initial condition for the small-$x$ evolution of the dipole scattering amplitude. The vector meson production cross section is shown to depend only weakly on the meson distribution amplitude and the factorization scale. We also present phenomenological comparisons of our result to the existing exclusive $\phi$ and $\rho$ production data from HERA and find an excellent agreement at high virtualities.


I. INTRODUCTION
Deep inelastic scattering (DIS) is a powerful tool to study the partonic structure of protons and nuclei at high energies. This process has been studied in detail in electron-proton collisions at HERA, where the vast amount of measured data has revealed a rapid increase in the density of gluons with small momentum fraction x [1,2]. The observed increase cannot continue indefinitely without violating unitarity, and as such the saturation effects are expected to dominate the small-x part of the hadron wave function. To describe QCD in this region of phase space where parton densities become of the same order as the inverse of the strong coupling, an effective field theory approach to QCD, called the Color Glass Condensate (CGC), has been developed [3][4][5].
In CGC, the high density of small-x gluons gives rise to nonlinear dynamics that slows down the growth of the gluon density. Despite the success of the CGC-based calculations in describing various high-energy collider experiments [6], there has not been definitive experimental evidence of saturation. To get precise DIS data from the saturation region new experimental facilities have been proposed, such as the upcoming Electron-Ion Collider (EIC) in the US [7][8][9] and a similar collider at CERN [10]. These facilities would allow for DIS measurements with heavy nuclei where the saturation effects are amplified approximately by A 1/3 . To meet the precision of these future experimental studies where nonlinear QCD dynamics is probed, it is necessary to promote the theory calculations in the CGC framework to higher-order accuracy.
One powerful process to probe gluon saturation is exclusive vector meson production as it requires an exchange of at least two gluons with the target. This renders the cross section roughly proportional to the gluon density squared [11] at leading order (but the situation is more complicated at next-to-leading order in a collinear factorization based approach, see Ref. [12]). Another advantage of it is that only in exclusive processes it is possible to measure the momentum transfer squared t in the process. The momentum transfer dependence can be related to the impact parameter dependence via a Fourier transform, providing access to the spatial distribution of nuclear matter in nuclei at high energy [13,14] and to the generalized parton distribution functions [15].
A convenient approach for describing exclusive vector meson production in DIS is the dipole picture where the process can be written in terms of the virtual photon and meson light-front wave functions along with the dipoletarget scattering amplitude [16,17]. The dipole amplitude satisfies perturbative small-x evolution equations, such as the Balitsky-Kovchegov (BK) equation [18,19] which resums large logarithmic contributions ∼ α s ln 1/x. The photon wave function can be calculated perturbatively [20,21], but the meson wave function is instead nonperturbative and therefore requires additional modeling. For heavy vector mesons one can take advantage of the small relative velocity of the quark-antiquark pair in the meson and model it as a fully nonrelativistic bound state [11], with velocity corrections that can be linked to the nonrelativistic QCD (NRQCD) matrix elements [22]. Another possibility is to take the limit of high photon virtuality, Q 2 M 2 V (where M V is the meson mass), where one can make a twist expansion for the process [23,24]. This corresponds to writing the meson wave function in terms of a nonperturbative distribution amplitude, on top of which higher-order corrections can be calculated perturbatively. This is especially suitable for light vector mesons and is a basic assumption in this paper.

A. High energy factorization
The scattering amplitude for exclusive vector meson production at high energy and in the zero squared momentum transfer t = 0 limit can be written in a factorized form and the coherent vector meson V electroproduction cross section can now be obtained as dσ γ * +p→V +p dt Here x 0,1,2 are the quark, antiquark and gluon transverse coordinates and z i denote the fractions of the photon's plus momentum carried by these partons. This factorization is justified at high energy as the lifetimes of the virtual photon qq and qqg Fock states are much longer than the timescales related to the interactions with the target color field. We use the eikonal approximation and describe the interactions with the target in terms of Wilson line correlators. The Wilson line V F,A (x) describes a color rotation of a quark (fundamental representation F ) or a gluon (adjoint representation A) when it propagates through the target, and the relevant correlators read where S 01 = 1 − N 01 . Here we took the mean field limit where the average over the target color charge configurations denoted by · · · factorizes. These correlators satisfy the BK evolution equation discussed in Sec. II C and depend implicitly on evolution rapidity which we will specify later. We only consider forward production in this work even though this framework can also be extended to calculate the momentum transfer dependent cross section. The momentum transfer is the Fourier conjugate to the impact parameter, and thus being able to calculate the cross section at finite momentum transfer is an advantage of exclusive processes as this can be used to do spatial imaging of the hadron structure [13]. On the other hand, this means that calculating vector meson production at t = 0 would require us to implement a model to describe the nonperturbative spatial structure of the proton. As the purpose of this work is to focus on a rigorous NLO calculation of vector meson production cross section we choose not to employ any such modeling and limit our studies to t = 0 where only the dipole amplitude integrated over the impact parameter is required. The same quantity is also probed in structure function measurements that are used to constrain the initial condition for the BK evolution of the dipole scattering amplitude N 01 [37].

B. Twist expansion
The meson light-front wave function is highly nonperturbative. For heavy vector mesons one can model the wave function based on the nonrelativistic nature of heavy quarks [22] but this simplification cannot be made for light mesons. On the other hand, the high-virtuality limit Q 2 M 2 V , which is justified for light mesons, can be used to simplify the mathematical description of the process. In this limit transverse momentum scales on the meson side become corrections suppressed by powers of 1/Q 2 , leading to the twist expansion of the meson wave function [23,24]. The leading-twist term then does not depend on the transverse momentum scales of the meson, meaning that only the dependence on the longitudinal momenta remains.
The twist expansion can be explained formally using the virtual photon wave function. The photon wave function is exponentially suppressed in Q 2 r 2 where r is the dipole size, which renders the relevant dipole sizes to be r 2 ∼ 1/Q 2 . We can then do a Taylor expansion for the meson wave function Ψ V (r, z) = Ψ V (0, z) + 1 6 r 2 ∇ 2 r Ψ V (0, z) + . . . = Ψ V (0, z)+O 1 Q 2 (see also Ref. [22]). Thus, only the dependence on the momentum fraction z remains at leading order in the twist expansion. The momentum space equivalent of this is a delta function in terms of the quark transverse momentum k: . This first term in the wave function corresponds to the twist-2 distribution amplitude φ(z) of the meson. In an NLO calculation the distribution amplitude has to be renormalized as we will demonstrate explicitly below, and the scale dependence of the renormalized distribution amplitude is described in terms of the Efremov-Radyushkin-Brodsky-Lepage (ERBL) evolution equation [20,56] which is discussed in more detail in Sec. IV C.
The twist expansion also guarantees that we need the nonperturbative part of the meson wave function only for the qq state. The nonperturbative part for other Fock states, such as qqg, is higher order in twist and can therefore be neglected at high Q 2 [57,58] . This means that the meson wave function for the qqg state can be calculated perturbatively by considering a gluon emission from the qq state, i.e. at high virtualities the Fock state qqg is created through the process V → qq → qqg (see Figs. 1f and 1g).
Another consequence of the high virtuality is that we need to consider only the longitudinal polarization for both the photon and the meson. A polarization flip is highly suppressed in coherent vector meson production such that the meson and photon effectively have the same polarization [59] (see also Ref. [60]). In fact, in the limit of zero momentum exchange t = 0 the polarization flip contribution vanishes exactly in our calculation. In the case of transverse production the leading-twist distribution amplitude is twist 3 [23], meaning that transverse production is suppressed relative to longitudinal by σ T /σ L ∼ M 2 V /Q 2 for high virtualities. Thus, total light vector meson production is given by the longitudinal cross section σ(γ * C. High-energy evolution The dipole amplitude, given by the correlator N 01 = 1 − S 01 , satisfies the perturbative BK equation describing its energy dependence. At leading order the BK equation reads [18,19]: This equation is written in terms of a rapidity variable Y which is discussed in more detail shortly. The kernel K BK describes the probability density for a dipole with transverse coordinates x 0 and x 1 to emit a gluon at the transverse coordinate x 2 . Including the running coupling corrections following Ref. [61], the kernel can be written as where we use the notation x ij = x i − x j . The BK equation effectively resums the contributions α s ln 1/x ∼ 1 from small-x gluons which is necessary for the stability of the perturbative calculations at high energy. When higher-order corrections enhanced by large double transverse logarithms are resummed [62], the NLO BK equation [25] becomes stable and can in principle be used in phenomenological applications [27]. A usual and numerically convenient approach, however, is to include resummations of the most important higher-order corrections to the leading-order BK equation. The leading-order BK equation with such resummations can be used to accurately approximate the full NLO BK equation [27,63]. Several resummation schemes exist, and in this work the following equations are used (we adopt the terminology used in Ref. [37]): KCBK [55], ResumBK [62,64] and TBK [65]. The nonperturbative initial conditions for these evolution equations have been determined in Ref. [37] by performing a fit to the HERA structure function data [1]. Of these, the evolution rapidity in the KCBK and ResumBK equations is the projectile rapidity Y = ln k + P + , where k + and P + are the gluon and target plus momenta. We work in the frame where the photon plus momentum q + is large and the photon has no transverse momentum. The target plus momentum is obtained as P + = Q 2 0 /(2P − ), where the target transverse momentum scale is taken to be Q 2 0 = 1 GeV 2 (note that the photon-nucleon center-of-mass energy reads W 2 = 2q + P − ).
Both KCBK ("kinematically constrained BK") and ResumBK ("resummed BK") involve a resummation of double transverse logarithms ∼ α s ln |x02| |x01| ln |x12| |x01| , with ResumBK also resumming single transverse logarithms α s ln 1 at all orders. In the KCBK equation the double logarithms are resummed by explicitly requiring a time ordering between the subsequent gluon emissions, which results in a non-local equation. The third evolution equation, TBK ("BK equation in target rapidity"), uses the target rapidity η as an evolution variable. This rapidity variable is related to the fraction of the target longitudinal momentum fraction transferred in the scattering process in the frame where the target has a large longitudinal momentum (see Ref. [65] for a detailed discussion). This evolution rapidity corresponds to Consequently, the TBK evolution can be thought as evolution in ln 1/x P whereas the KCBK and ResumBK equations correspond to evolution in ln W 2 . In order to use dipole amplitudes as a function of the target rapidity η in the impact factors written in terms of the projectile rapidity Y , we use the same shift as in Ref. [37]: The BK equation contains a transverse-coordinate dependent coupling constant. We model the running of the coupling in the coordinate space following Ref. [37]: This running coupling approaches a constant value in the infrared region 1/|x ij | Λ QCD , with the constants µ 0 and c controlling its behavior there. The values of these constants are chosen as in Ref. [37]. The constant C 2 is a fit parameter that describes the relation between momentum and coordinate spaces, k 2 = 4C 2 /r 2 , with the expected value C 2 = e −2γ E from Fourier analysis [66,67]. The same coordinate space coupling constant is used when calculating the scattering amplitude, Eq. (1), where the coupling constant is included in the next-to-leading order photon and meson wave functions. As the running coupling prescription (6) can be seen to effectively choose the smallest of the three distance scales x 2 01 , x 2 12 , x 2 02 , when calculating the qqg contribution in Eq. (1) we choose to evaluate the coupling at the scale set by the smallest of the daughter dipoles, as in Ref. [37]. When evaluating the qq term the scale choice is x 01 .

III. LIGHT-FRONT WAVE FUNCTIONS AT NEXT-TO-LEADING ORDER
The NLO corrections to exclusive vector meson production can be calculated in terms of the NLO wave functions for the photon and meson. In this section, we first list the relevant photon light-front wave functions at NLO accuracy calculated in Refs. [29][30][31]. Then, we proceed to calculate the light vector meson wave function at NLO in terms of the twist-2 distribution amplitude, and present the results Fourier transformed to mixed transverse coordinate, longitudinal momentum fraction space.

A. On the regularization scheme
The calculation will be done in two different regularization schemes. The first one is the conventional dimensional regularization (CDR) where the momenta and polarization vectors of all particles are continued to D dimensions. The second one is the four-dimensional helicity (FDH) scheme where the polarization vectors are kept in four dimensions [68,69]. In our case, this amounts to real gluons having two polarization states.
To do the calculations simultaneously in both schemes we follow the notation of Ref. [32]. The dimension arising from the gluon polarization vectors is denoted as D s to distinguish it from the dimension D in the dimensional regularization. The CDR scheme corresponds to the case D s = D, and for the FDH scheme we have D s = 4. Sums over gluon helicities can be calculated as λ i λ j * λ = δ ij (Ds) where the subscript denotes that this Kronecker delta has D s − 2 transverse dimensions. In the sums over spin and Lorenz indices we take D s ≥ D so that the following relations for the Kronecker deltas hold: We will also make use of the following spinor identity [29]: where h = ±1 is the quark helicity and Here the D s -dimensional Levi-Civita tensor has to be understood through the Fierz identity This identity is valid in D s dimensions if there are no more than two Levi-Civita tensors [29], which holds for the calculations considered in this paper.

B. Photon wave function
The photon light-front wave functions in the massless quark case have been calculated in Refs. [29][30][31] and are shown here for completeness. In our notation, an additional factor of 1 2q + i 1 √ zi appears in the wave functions. This additional factor follows from a different choice of the integration measure which we choose to be i . With this choice the leading-order wave function for the virtual photon in the mixed transverse coordinate and plus momentum fraction space is We always work in the frame where the photon transverse momentum is zero. Here e f is the fractional charge of the quark with flavor f , Q 2 is the virtuality of the photon, z i = k + i /q + is the (anti)quark's fraction of the photon plus momentum, and α i and h i are the color and helicity indices. We also use the short-hand notation Quantities corresponding to the quark are denoted with i = 0 and antiquark with i = 1. We note that the last factor, which is equal to 1 at D = 4, is absent in Ref. [29] where the transverse momenta of the observed particles are kept in two dimensions. Here "observed" particles are those which appear as the final state in the wave function, not including soft or collinear particles. In this paper, we choose to evaluate the transverse momenta of the observed particles in D − 2 dimensions, as this is necessary for regularizing the NLO meson wave function. However, this term does not have any contribution to the final cross section where all 1 D−4 divergences have been canceled. In principle, this factor multiplied by K γ * L contributes a finite logarithm term ∼ ln |x 01 |. It however cancels when we perform the UV subtraction in Sec. IV B.
The next-to-leading order correction to the photon wave function can be written as where Here α is the infrared cut-off for the gluon plus momentum fraction and µ is the mass scale for dimensional regularization, and the last term depends on the regularization scheme. The virtual photon wave function for the Fock state qqg can be written as where and Quantities with the subscript 2 correspond to the emitted gluon. The function I i differs from the similar one in Ref. [29] by an additional power 4−D 2 for the variable u, which has the same origin as the last factor in Eqs. (14) and (15).

C. Light vector meson wave function
In this section, we calculate the NLO corrections to the light vector meson light-front wave function. The calculation is done in the limit where the transverse coordinate dependence of the meson leading-order wave function can be neglected. As discussed in Sec. II, this follows from the large photon photon virtuality. This means that we can neglect all mass scales in the meson, allowing us to set the meson mass M V to zero along with the transverse momenta k i of the quark and the antiquark. We work in a frame where both the photon and the vector meson transverse momenta are zero, as we consider forward production. Consequently, at leading order the meson wave function is given by a delta function in the transverse plane: Here φ 0 (z) is the (bare) distribution amplitude of the meson that describes how the meson plus momentum is shared by the two quarks. This wave function is normalized in such a way that it gives the correct decay constant f V given that the distribution amplitude is normalized as The decay constant f V is related to the leptonic width by The wave function (20) describes the probability of the meson to split into a quark-antiquark pair with the flavor f . Here c f is a normalization factor needed for mesons that consist of a superposition of different flavored quark-antiquark states. For example, the ρ meson can be written at leading order as |ρ = 1 The normalization factors are also related to the effective charge fraction of the meson which is defined by e V = f c f e f . We emphasize that in the high-Q 2 limit the dependence on the vector meson type is included in the nonperturbative distribution amplitude φ 0 (z 0 ) (in addition to the normalization factors f V , e V and c f ).
At next-to-leading order, we get perturbative corrections to the meson wave function from Feynman diagrams shown in Fig. 1. Of these, the diagrams 1a and 1b, corresponding to the self-energy corrections of the quark and antiquark, evaluate to zero. This is a consequence of the dimensional regularization used in the calculation, as these diagrams give transverse integrals with no mass scales (in the high-Q 2 limit considered here where we neglect the quark and meson masses and the quark transverse momenta) such as d D−2 k0 To calculate the rest of the diagrams we use the Feynman rules of the light-cone perturbation theory from Ref. [29]. For the diagram 1c this gives: where the identity (11) has been used to simplify the result. The square root factor in the first line comes from our choice for the integration measure, and the quark and antiquark transverse momenta after the gluon exchange are k 0 and k 1 = −k 0 . We use a notation u(0) = u h0 (k 0 ), v(1) = v h1 (k 1 ) for the quark and antiquark spinors, and the primed quantities correspond to the intermediate quark and antiquark whose spins and helicities are summed over (see Fig. 1). The contribution of the diagram 1d is similar to the diagram 1c. An explicit calculation gives the result Eq. (23) with the substitutions z 0 → 1 − z 0 and z → 1 − z : Here we used the symmetry condition φ(z ) = φ(1 − z ) which follows from the parity of the vector meson.
The final contribution to the NLO qq wave function comes from the diagram 1e describing an exchange of an instantaneous gluon between the quark and the antiquark. The contribution from this diagram can be evaluated to give Summing these contributions together, we get the NLO correction to the meson qq wave function: It should be noted that this NLO correction does not affect the normalization (21) of the distribution amplitude. The reason for this is that the decay constant is given by , and this integral vanishes for Eq. (26) in dimensional regularization. We also need the wave function for the qqg state. This is simply given by the sum of the diagrams 1f and 1g, which evaluates to: Note that the momentum conservation implies z 0 + z 1 + z 2 = 1 and k 0 + k 1 + k 2 = 0. These wave functions are presented in the momentum space. For the meson production calculation we need the mixed space wave functions which can be calculated from the momentum space wave functions by a Fourier transform in the transverse plane. The leading-order wave function in the mixed space is given by The NLO correction to the qq wave function is given by and the wave function for the qqg state is where

A. Production amplitude
Having determined the NLO corrections to the meson wave function, we now have all the ingredients to calculate the exclusive light meson production amplitude. We substitute the photon wave functions for the qq (sum of Eqs. (14) and (15)) and qqg (Eq. (17)) states, along with the meson wave functions for the qq (sum of Eqs. (28) and (29)) and qqg (Eq. (30)) states, into Eq. (1) to obtain the production amplitude and keep terms up to O(α s ). The production amplitude can then be divided into the dipole (qq) and real emission (qqg) parts. The dipole part contains the leading-order result and the NLO correction These choices for the impact parameter b follow Ref. [37], but we note that in the t = 0 case the weighting of the coordinates by the momentum fractions z i in the definition of b does not affect the results.
The coherent vector meson V electroproduction cross section (2) can now be evaluated using the scattering amplitude When squaring the amplitude, we keep terms up to O(α s ). However, the iA qqg amplitude also contains a large contribution enhanced by a large logarithm ln 1/z 2 ∼ 1/α s , and as such this contribution has to be considered as being part of the leading order amplitude. This is in practice done by taking into account the BK evolution as we will discuss in more detail in Secs. IV D and V.

B. UV subtraction
The dipole (−iA qq NLO ) and real emission (−iA qqg ) parts of the amplitude are separately UV divergent. However, most of the divergences cancel in their sum. Therefore it is useful to subtract the UV divergent part of the real emission and combine it with the dipole part. The subtracted term is chosen to be where This choice for the UV subtraction term is analogous to the one in Ref. [29] and also what is used when considering heavy vector meson production in Ref. [41]. Unlike in Ref. [29], we choose to include the additional factor Q 2π|x01| D− 4 2 to the UV subtraction to cancel the same factor in the dipole part. The integrals over x 20 and z 2 can be done analytically, which simplifies the UV subtraction term to We then add this to the dipole part, which gives us Here K(z, z ) is the kernel of the ERBL equation [20,56] which describes the scale dependence of the distribution amplitude as we will discuss in Sec. IV C: This form for the ERBL kernel is equivalent to the usual one written in terms of the plus distributions in the limit α → 0. After the UV subtraction, the real emission part becomes finite and reads: where

C. ERBL evolution and the renormalized distribution amplitude
The dipole part −iA qq sub , Eq. (39), still contains a divergence of the form 1 D−4 which is canceled when the distribution amplitude is renormalized. We define the renormalized distribution amplitude φ(z, µ F ) as where µ F is the factorization scale. This choice for the finite terms in the subtraction corresponds to the MS scheme. We note that the distribution amplitude depends on the regularization scheme (FDH or CDR), as in practice it has to be determined from some experimental process for which an NLO calculation also depends on the same scheme choice. In principle the scheme-dependent term ∼ (D s − 4)/(D − 4) in Eq. (39) could be also included in the definition of the renormalized distribution amplitude (42). However, in this work we choose to keep the scheme dependence explicitly visible in the dipole term, Eq. (39). This allows us to straightforwardly quantify the scheme dependence which is shown in Appendix A to be negligible. The renormalized distribution amplitude satisfies the ERBL evolution equation [20,56] ∂φ(z, µ F ) ∂ ln µ 2 where the kernel K(z, z ) is given in Eq. (40). We note that this renormalization does not change Eq. (21) for the normalization of the distribution amplitude as the z-integral over the ERBL kernel vanishes: 1 0 dz K(z, z ) = 0. Next we use Eq. (42) to write the bare distribution amplitude φ 0 (z) in −iA qq sub , Eq. (39), in terms of the renormalized distribution amplitude. We also choose to use the scale dependent renormalized distribution amplitude instead of the bare distribution in the NLO part, as their difference is now formally higher order in α s . This results in the finite expression Similarly we can replace φ 0 (z) by φ(z, µ F ) in the real emission part (41). Let us briefly consider the evolution of the renormalized distribution amplitude. It is useful to write the distribution amplitude in terms of the eigenfunctions f n (z) of the ERBL kernel The eigenfunctions can be written in terms of the Gegenbauer polynomials C n (2z − 1), and the corresponding eigenvalues are given by [20] Writing the distribution amplitude as a sum of the eigenfunctions, the ERBL equation then tells us that the coefficients in the sum depend on the factorization scale: Taking into account the running of the coupling constant as α s (µ 2 , we can solve the evolution of the coefficients a n explicitly [20]: a n (µ F ) = a n ln µ 2 Here Λ QCD = 0.241 GeV and β = (11N c − 2N f )/3 with N f = 3. Values for the coefficients a n at an initial scale are a nonperturbative input for the calculation. These coefficients also depend on the considered vector meson, and should be determined from experimental data. It should be noted that the eigenvalue λ n is zero for the term n = 0 and negative for the n > 0 terms. This means that the first term is actually constant in µ F , and the higher-order terms become suppressed as µ F increases. In the asymptotic limit µ F → ∞ only the first term contributes, and the distribution amplitude then simplifies to φ(z, µ F = ∞) = 6z(1 − z). Here we have also used the fact that the coefficient a 0 of the first term is actually determined by the normalization condition Eq. (21), as the orthogonality of the Gegenbauer polynomials guarantees that only the first term contributes to the normalization, giving us a 0 = 1. It should also be noted that parity conservation demands that the distribution amplitude is invariant under the substitution z ↔ 1 − z, meaning that all terms with n = odd are zero in the sum. We point out that Eq. (48) is divergent for µ F = Λ QCD . In practice, we avoid this singularity by introducing an infrared (IR) cutoff µ F 0 for the ERBL evolution and freeze the distribution amplitude below this scale: φ(z, µ F ) = φ(z, µ F 0 ) for µ F < µ F 0 . We choose the value of the IR cut-off to be µ F 0 = 1 GeV. The dependence on the IR cutoff is quantified in Appendix A.

D. Soft gluon divergence
The real emission part still has an IR divergence from the lower limit α of the z 2 integral. This is related to the emission of soft gluons from the dipole, and to the rapidity evolution of the dipole amplitude. This can be seen by noting that the singular part of the real emission can be written as where the identity [29] has been used. Note that as we do not have an explicit dependence on the infrared cutoff α in the integrands anymore, from now on the lower limit of the z 2 integral is denoted by z min whose value will be discussed shortly. We can recognize the integrand in Eq. (49) as the kernel of the (fixed coupling leading order) BK equation (5). This can then be combined with the leading-order term (α 0 s part of Eq. (44)), and the sum of these two contributions corresponds to using in the leading-order term a dipole amplitude evolved from the initial rapidity Y 0 to the rapidity At finite center-of-mass energy the lower limit z min of the z 2 integral should not be taken to zero. In particular, the invariant mass of the qqg system should be much less than W 2 in order to justify the usage of the eikonal approximation, which imposes the lower limit z min . We follow Refs. [37,41] and choose Here the minimum comes from the kinematic constraint z 0 + z 2 ≤ 1 which guarantees that the dipole does not evolve backwards in rapidity. As we are interested in the high (but finite) energy limit, the minimum is only needed in a small subset of the integration region and in practice the evolved rapidity is For the α s -suppressed terms the dependence on the evolution rapidity is formally of higher order in the coupling constant. Following again Refs. [37,41] we choose to use the same evolution rapidity Y dip when evaluating the nextto-leading order terms in the dipole part, Eq. (44). The z 2 -dependent evolution rapidity used with real gluon emission term is obtained from the definition Y = ln k + P + and can be written as [41] Y qqg = ln z 2 + ln We can now write the scattering amplitude for light meson electroproduction in its full form. It reads where the LO part is and the NLO corrections are for the dipole part and for the real emission. The lower limit for the z 2 integral is given by Eq. (52). This expression is finite and suitable for numerical evaluation. The rapidity scales at which the different dipole amplitudes are evaluated, Y 0 , Y dip and Y qqg , are shown explicitly. In numerical calculations we follow Ref. [37] and take Y 0 = 0.
The NLO correction to the dipole part has a dependence on the regularization scheme given by a term proportional to This regularization scheme dependence is in principle canceled by the regularization scheme dependence of the distribution amplitude at the given order in α s . The distribution amplitude is a nonperturbative quantity that has to be determined from some process where the same regularization scheme dependence should also appear. In this paper, we choose to use the CDR regularization scheme when we show numerical results in Sec. V. However, it will turn out that the regularization scheme dependence is very small even if the same distribution amplitude is used in both schemes, which is a consequence of the fact that for the first term in the Gegenbauer expansion (47) of the distribution amplitude this regularization scheme dependent term vanishes. The regularization scheme dependence of the cross section will be discussed quantitatively in Appendix A.
The dependence on the factorization scale µ F is of higher order in α s , as can be verified by taking into account the ERBL equation. However, as we are keeping terms only to the order α s , the results do have a dependence on the factorization scale. The value of µ F can be chosen in different ways. Eq. (44) suggests the choice µ 2 F = 4e −2γ E /x 2 01 , as with this choice the logarithm multiplying the ERBL kernel K(z, z ) vanishes (we will refer to this term as the "ERBL term"). Note that the factor 4e −2γ E is the same one that appears in the Fourier analysis of the coordinate space running coupling [66,67]. This choice for the factorization scale will be referred to as the r-scheme to emphasize its dependence on the dipole size. In the qqg term we choose to use the smallest dipole size min{|x 01 |, |x 20 |, |x 21 |} for the factorization scale, in accordance with the running of the coupling constant α s . Another possible choice for the factorization scale is to use µ F = Q, which is supported by the fact that the relevant length scales for meson production are Q ∼ 1/|x 01 |, meaning that the logarithm ln Q 2 x 2 01 of the ERBL term in Eq. (44) should also be small in this scheme. This will be referred to as the Q-scheme, and its main advantage is that the hard scale does not depend on the integration variable. In this paper, we will use the r-scheme in our calculations as then the ERBL term vanishes completely. In the Q-scheme, there can in practice be a large contribution from the ERBL term as the dipole amplitude amplifies the contribution of larger dipoles that can have a numerically significant contribution even at moderately large virtualities [70,71]. This scheme dependence is studied in more detail in Appendix A, where it is shown that the factorization scale dependence at the cross section level is a few percent. The result (55) can be compared to the previously calculated NLO light vector meson production from Ref. [53]. In that paper the production amplitude is presented in the momentum space as opposed to the mixed space used in this paper. Comparing these results is non-trivial, as one has to perform complicated Fourier transforms from momentum space to coordinate space to match the results. We have only been able to make comparisons for the dipole part of the amplitude, finding that the result of Ref. [53] matches our result in the CDR scheme apart from differences in the UV subtraction procedure. The real gluon emission part is much more complicated, and so far we have not been able to make actual comparisons of the results.

V. NUMERICAL RESULTS
In this section, we present numerical results for coherent light vector meson electroproduction at next-to-leading order, calculated using Eq. (55). As our default setup, we use the CDR scheme for regularization and r-scheme for the factorization scale µ F . For the distribution amplitude, we choose to keep only the first two terms in the Gegenbauer expansion (47). The reason for this is that the exact values for the higher-order terms are not well known, but estimated to be small [72]. For the ρ meson, the coefficient of the second term has been extracted in many different ways, with relatively large uncertainties [73]. We choose to use the value a 2 (µ F = 1 GeV) = 0.1 which is in agreement with most of the values tabulated in Ref. [73]. We also choose to use this same value for the φ meson, as current analyses suggest that they are of the same order of magnitude [57,72]. As we then use the same distribution amplitude for both mesons, the only difference between ρ and φ production is the decay constant f V which appears as an overall coefficient in Eq. (55). These decay constants can be calculated from the experimental values for the leptonic widths [74] using Eq. (22).
The numerical results are calculated using the dipole amplitude fits from Refs. [37,75] for the different schemes of the BK evolution equation discussed in Sec. II C. We use the fits where the "Balitsky+smallest dipole" running coupling scheme is used, and use both fits with initial evolution rapidities Y 0,BK , η 0,BK = 0 and Y 0,BK , η 0,BK = 4.61 (in which case the dipole amplitude is frozen in the region Y 0 = 0 < Y < Y 0,BK or η 0 = 0 < η < η 0,BK ). In these fits the impact parameter dependence is assumed to factorize and one can replace d 2 b → σ 0 /2, and the proton transverse area σ 0 /2 is a fit parameter which is also determined in Ref. [37].
In Fig. 2 we show different contributions to the exclusive ρ production amplitude at NLO as a function of the center-of-mass energy W (Fig. 2a) and photon virtuality Q 2 (Fig. 2b). The same dipole amplitude, corresponding to the KCBK equation with the initial rapidity Y 0,BK = 4.61 for the BK evolution [37], is used in these figures. Here the leading-order result is denoted by LO(Y dip ) which is calculated from the leading-order part of Eq. (55) with the dipole amplitude evaluated at rapidity Y dip . Using the evolved rapidity Y dip means that the LO(Y dip ) contains the resummation of large logarithms ∼ α s ln 1/x included in the BK evolution. The result LO(Y 0 ) is the leading-order term of Eq. (55) at the initial rapidity Y 0 , and NLO dip is the NLO correction to the dipole term corresponding to Eq. (57). The contribution from the qqg term, Eq. (58), has been divided into two parts: the result NLO qqg (BK) contains only the part corresponding to the BK equation, Eq. (49), and NLO qqg (no BK) is Eq. (58) from which the BK contribution has been subtracted. The total NLO result is then the sum From these plots we see that the contributions from the LO result at the initial rapidity and from the NLO dipole term are small. Both qqg contributions are large, but they mainly cancel each other. These findings are similar to what has been observed in the case of heavy vector meson production [41]. The total NLO correction, the difference between NLO and LO(Y dip ), is large and positive. However, we point out that in these plots the same NLO fitted dipole amplitude was used to calculate all of the results. Consequently, these results only quantify the largeness of the NLO correction terms in Eq. (55) and not the actual difference between the NLO and LO results where the corresponding dipole amplitude fits should be used.
We also note that at fixed coupling the identification LO(Y dip ) = LO(Y 0 ) + NLO qqg (BK) would be exact if the dipole amplitude satisfied the leading-order fixed coupling BK equation. In that case there would be no ambiguity in defining the leading-order amplitude. In our setup this is not the case, and consequently the leading-order amplitude is not uniquely defined. In this work we choose it to be LO(Y dip ) following Ref. [41], as this is the most natural choice   when using a dipole amplitude that satisfies a resummed BK evolution equation. Identifying LO(Y 0 ) + NLO qqg (BK) as a leading order term instead would have maximally a ∼ 20% effect on the calculated cross sections discussed below.
Next we show numerical comparisons to the existing coherent vector meson production data for ρ and φ mesons at (moderately) large Q 2 . The H1 data is from Ref. [59], and the ZEUS data is from Ref. [76] for φ and Ref. [77] for ρ. The results are shown with various different dipole amplitude fits that all give a good description of the HERA structure function data. As discussed above, the NLO results use fits from Ref. [37]. For the leading order, the dipole amplitude used is the "MV e " fit from Ref. [35]. In the leading-order calculation the evolution rapidity is chosen as , consistently with the fit. The differential cross section is proportional to the square of the production amplitude as given by Eq. (2). When calculating the cross section at NLO, we drop the higher-order terms proportional to α 2 s so that we only keep the genuine NLO correction at the cross section level. In Fig. 3 we show the differential cross section for the longitudinal φ and ρ production at t = 0. Here the experimental data is for the total production which is the sum of the longitudinal and transverse channels. However, the longitudinal production dominates at Q 2 M 2 V and therefore it is expected that for high virtualities these data points accurately correspond to the longitudinal case.
In general, we see that both the LO and NLO results describe the H1 data well. The difference between the LO and NLO results is smaller than one would expect based on Fig. 2, as in the leading-order fit the nonperturbative parameters describing the initial condition of the dipole amplitude effectively capture part of the higher-order effects. This difference becomes small at high virtualities, where our approach is expected to be most reliable. The NLO results also give a surprisingly accurate description of the data for smaller values of the photon virtuality where the framework cannot be trusted (we have assumed that Q 2 M 2 V ). We consider this agreement to be accidental for two reasons. First, we have only calculated the longitudinal cross section, leaving out the transverse contribution which is significant at small Q 2 . Secondly, we neglect any dependence on the dipole size in the meson wave function, keeping only the r = 0 case which corresponds to the distribution amplitude. In general, the wave function is expected to be a decreasing function of |r|, meaning that this approximation overestimates the results. These two corrections, which may have significant numerical contributions at small virtualities, affect the result in opposite ways and therefore their total contribution at least partially cancels.
Next we will consider the t-integrated cross sections for which more data exists. To avoid additional modeling for the impact parameter dependence of the dipole amplitude, we evaluate the t-integral by using the following experimental parametrization for the t-dependence of the cross section: Here b is the slope parameter that in general depends on Q 2 and W . It has been measured for both φ and ρ [59,76,77] at different values of the virtuality at W = 75 GeV. The slope parameter can be thought of as the effective transverse area of the meson-target system, and we model its dependence on virtuality and center-of-mass energy by assuming The measured slope parameter b for ρ and φ production as a function of photon virtuality [59,76,77] and a fit to this data.
The W dependence determined from HERA data [59] gives α = 0.12 ± 0.04. The model for the virtuality dependence is chosen for its simplicity and that it approaches a constant value at high Q 2 . Also, the dependence on the virtuality and the center-of-mass energy does not seem to be correlated [59]. We fit the parameters b 0 and b 1 to H1 and ZEUS data at W 0 = 75 GeV, with the fit shown in Fig. 4, and note that the errors on these fitted parameters are significant, which results in ∼ 10% uncertainty in the calculated total cross sections. The virtuality dependence of the coherent φ and ρ production cross sections is shown in Figs. 5 and 6. In Fig. 5, the results are calculated using different dipole amplitudes fitted to the HERA structure function data in Ref. [37], using fits with both choices for the initial evolution rapidities Y 0,BK (η 0,BK in the case of TBK evolution). The H1 collaboration has measured, in addition to the total production cross section, the longitudinally polarized ρ production, which exactly corresponds to the presented theory calculations. In general we find an excellent agreement with the H1 and ZEUS data [59,76,77], except that the φ production cross section is overestimated at low virtualities where our approximations are not justified.
In Fig. 6, we show results obtained using the dipole amplitudes fitted to the structure function pseudodata generated in Ref. [37] such that it includes only the approximative light quark contribution. For comparison, the results calculated     with the dipole amplitudes fitted to the full HERA structure function data are also shown. In the fit process of Ref. [37] only the light quark contribution is calculated, and as such the fit to light quark pseudodata is in principle better motivated than the fit to the full HERA structure function data. On the other hand, the light-quark-only data contains a larger nonperturbative contribution, and the determined parametrizations describing the initial condition of the dipole amplitude are not physically as well motivated. Here we use these light quark fits with the KCBK evolution equation, but different schemes for the BK evolution result in very similar cross sections at all Q 2 . We see that the results calculated with dipole amplitudes fitted to the light quark pseudodata also show a relatively good agreement with the virtuality dependence of the H1 and ZEUS light meson production data. However, the cross sections at large virtualities are somewhat underestimated. This difference in Q 2 dependence between the two fit setups is expected, as the light-quark-only pseudodata is close to the full structure function data at low Q 2 where similar results for other observables are also expected. On the other hand, at high Q 2 the charm contribution on structure functions is significant, and consequently the light quark fit should result in smaller cross sections in this kinematical region, which is exactly what we observe in Fig. 6.
In Figs. 7 and 8, we show the dependence of the integrated cross section on the photon-proton center-of-mass energy W . Again, Fig. 7 shows results obtained with the dipole amplitudes fitted to the HERA data, and Fig. 8 shows results calculated with dipole amplitudes fitted to the light quark pseudodata for comparison. The center-of-mass energy      dependence of the results agrees with the data, although the results with the light quark fit seem to underestimate the data by a constant factor as already seen in Fig. 6. The differences in the results with different schemes for the BK evolution start growing at larger W , as at high W one starts to be sensitive to the region not constrained by the structure function data. This suggests that light meson production data can provide additional constraints when the nonperturbative initial condition for the BK evolution is extracted from experimental data. The dependence on the center-of-mass energy is similar for the HERA and light quark fitted dipole amplitudes. In Fig. 9a we show the dependence of the cross section on the distribution amplitude. We have normalized the results by the cross section calculated using the asymptotic form for the distribution amplitude which corresponds to the case where the higher-order terms in the Gegenbauer expansion (47) vanish, i.e. a n = 0 for n > 0. The case a 2 (1 GeV) = 0.1 corresponds to our default setup, and the cases a 2 (1 GeV) = ±0.2 are estimates for the upper and lower bounds for the coefficient, chosen based on Ref. [73]. The final setup shown has the coefficients a 2 (1 GeV) = −0.054 and a 4 (1 GeV) = −0.022 chosen such that the distribution amplitude matches the Boosted Gaussian wave function parametrization for the ρ meson from Ref. [16] at r = 0 where the higher Gegenbauer terms are neglected. The reason for this choice is that the distribution amplitude should roughly correspond to the wave function at r = 0, and the Boosted Gaussian is a phenomenological wave function that describes well vector meson production at leading order [16]. We see that the dependence on the distribution amplitude is moderate, and maximally ∼ 30% in the considered kinematical domain. The different distribution amplitudes are illustrated in Fig. 9b, both at the initial scale µ 2 F = 1 GeV 2 and after the ERBL evolution up to µ 2 F = 50 GeV 2 using the extreme values for a 2 . While the form of the distribution amplitude depends considerably on the value of a 2 , the effect on the cross section in Fig. 9a is small. In Fig. 9b we also see that as the factorization scale µ F increases the distribution amplitude approaches the asymptotic form, but there is still a significant deviation from the asymptotic shape at µ 2 F = 50 GeV 2 .

VI. CONCLUSIONS
In this paper, we have presented a next-to-leading order calculation of exclusive light vector meson production in the dipole picture using the light-cone perturbation theory. The main result of this work is the scattering amplitude for longitudinally polarized light vector meson production at high Q 2 given in Eq. (55). This amplitude is finite and directly suitable for numerical evaluations. In particular the 1 D−4 divergences have been canceled between the real and virtual diagrams after the ERBL evolution of the renormalized distribution amplitude is also taken into account. The apparent soft gluon divergence is shown to be factorizable into the small-x Balitsky-Kovchegov evolution of the dipole amplitude. While a similar calculation has already been done in the momentum space [53] (using also a different scheme to subtract the rapidity divergence), the results presented in this paper are in the mixed transverse coordinatelongitudinal momentum fraction space where the numerical calculations using existing results for the dipole-target scattering amplitude are straightforward. Comparing the results in the different spaces is cumbersome due to the requirement of calculating complicated Fourier transforms, and thus an explicit comparison of the two results has been made only partially.
We have also calculated numerically exclusive light meson production at NLO and compared the results to the existing HERA data for ρ and φ mesons. The NLO corrections are numerically important, but their effect can be partially captured when the initial condition for the small-x evolution of the dipole amplitude is fitted to the structure function data. Consequently, the differences between LO and NLO results are moderate in the high virtuality region Q 2 M 2 V where our framework is valid. The different schemes used to capture higher-order effects in the small-x evolution result in similar cross sections for vector meson production. Some deviations can be seen in the centerof-mass energy W dependence, which means that the exclusive vector meson production data can further constrain the nonperturbative initial condition for the small-x evolution. Both the Q 2 and W dependencies of the production cross section are in excellent agreement with the HERA data. If a dipole amplitude with an initial condition fitted to the structure function pseudodata that only includes a light quark contribution is used, the experimental cross sections are underestimated at high Q 2 . We also note that there is some overall normalization uncertainty due to e.g. modeling the t dependence of the vector meson production cross section. We additionally left out the commonly used phenomenological corrections (see e.g. [16]) whose role should be further clarified.
Our result for the analytic expression of the production amplitude is presented in two different schemes for regularization in the transverse plane: the CDR and FDH schemes. The regularization scheme dependence is shown to be very small. The results also depend on the choice for the factorization scale µ F , for which we present two different choices, taking this scale to be either a function of the dipole size or of the photon virtuality. The dependence on the factorization scale is also relatively small. Both of these scheme dependencies have numerically small effects because the distribution amplitudes for the ρ and φ mesons are close to the asymptotic form, and the dependence on regularization scheme and factorization scale vanishes in the Q 2 → ∞ limit. The dependence on the exact form of the distribution amplitude, on the other hand, is somewhat larger with effects of up to ∼ 30% in the HERA kinematics at the cross section level for realistic values of the higher-order terms in the Gegenbauer expansion.
The results in this paper are calculated at zero momentum transfer t = 0. Calculating the t-dependence of exclusive vector meson production is also interesting as it allows access to the spatial distribution of the target color field including its event-by-event fluctuations [14,78,79]. This requires additional nonperturbative modeling of the dipole amplitude which we wanted to avoid in this paper. We also note that the dipole amplitudes used in numerical calculations in this paper were fitted to HERA data using only massless quarks, while there is a significant contribution from the massive c quark to the structure functions in HERA kinematics. As the NLO photon wave functions with massive quarks are becoming available [32,33], it will be possible to make a new NLO fit for the dipole amplitude to the HERA data including heavy quark contributions. This is needed for accurate phenomenological comparisons with the HERA data. The results of this work can then be used for predicting exclusive light vector meson production in the future EIC which will also produce data for DIS off heavy nuclei, allowing for precision studies of saturation phenomena. In this Appendix we quantify the dependence of the light vector meson production cross section on the choices for regularization and factorization schemes. The cross sections have been calculated as described in Sec. V, and our default setup is the same. The dipole amplitude used in this Appendix is the KCBK evolved one from Ref. [37] with the initial rapidity Y 0,BK = 4.61.
First, we show the dependence of the cross section on the regularization scheme in the transverse plane. The cross section has been calculated in both the CDR and FDH schemes and their ratio is shown in Fig. 10. This ratio depends on the distribution amplitude, and it is equal to unity in the asymptotic limit µ 2 F → ∞ where only the first term in the Gegenbauer expansion (47) contributes. For this reason we show the ratio with two different distribution amplitudes: one with a 2 (1 GeV) = 0.1 (our standard setup), and one with a 2 (1 GeV) = 0, a 4 (1 GeV) = 0.1. Higher-order terms are set to zero. We see that the dependence on the regularization scheme is very small in both cases, of the order 0.2% at most. Small scheme dependence is expected, as the first dominant term in the Gegenbauer expansion vanishes when one calculates the scheme dependent term in Eq. (44). It should be noted that the ratio does not seem to approach the asymptotic limit in the considered kinematics. This is a consequence of the ERBL evolution with running coupling being extremely slow, and the scheme dependence vanishes if we go to even higher values of virtuality.
The cross section depends on the factorization scale µ F at which the distribution amplitudes are evaluated as discussed in Sec. IV C. In Fig. 11, the cross sections have been calculated evaluating the distribution amplitude in both the r-and Q-schemes using our default choice for the distribution amplitude with a 2 (1 GeV) = 0.1. The factorization scale has also been scaled by factors of 0.5 and 2. These results have been normalized by our default setup (r scheme with µ F = 2e −γ E /r). In the Q-scheme, varying the factorization scale by a factor 2 has a very small effect. For the r-scheme the cross section varies somewhat more, but the variation is still only ∼ 2%. The difference between the r-and Q-schemes is also only a few percent at most, meaning that the dependence on the factorization scale is small. This follows from the fact that the distribution amplitude receives only a small correction from the second scale-dependent Gegenbauer term, and the dominant term is factorization scale independent.
Finally, we show the dependence on the infrared cutoff µ F 0 in Fig. 12 using our default setup (r scheme and a 2 (1 GeV) = 0.1). There is some dependence on the IR cut-off, almost 5% at most with our choice for a 2 . The reason for the cut-off dependence is that the dipole amplitude amplifies the contribution of large dipoles, meaning that dipoles of size 1/r ∼ 1 GeV may have a numerically significant contribution even when Q 2 1 GeV 2 . The dependence on the IR cut-off vanishes exactly in the limit Q 2 → ∞.