The dipole picture and the non-relativistic expansion

We study exclusive quarkonium production in the dipole picture at next-to-leading order (NLO) accuracy, using the non-relativistic expansion for the quarkonium wavefunction. This process offers one of the best ways to obtain information about gluon distributions at small $x$, in ultraperipheral heavy ion collisions and in deep inelastic scattering. The quarkonium light cone wave functions needed in the dipole picture have typically been available only at tree level, either in phenomenological models or in the nonrelativistic limit. In this paper, we discuss the compatibility of the dipole approach and the non-relativistic expansion and compute NLO relativistic corrections to the quarkonium light-cone wave function in light-cone gauge. Using these corrections we recover results for the NLO decay width of quarkonium to $e^{+}e^{-}$ and we check that the non-relativistic expansion is consistent with ERBL evolution and with B-JIMWLK evolution of the target. The results presented here will allow computing the exclusive quarkonium production rate at NLO once the one loop photon wave function with massive quarks, currently under investigation, is known.


I. INTRODUCTION
The partonic structure of hadrons and nuclei in the limit of high collision energies, or equivalently small momentum fractions x, is poorly constrained by existing experimental data. It is believed that at high enough energies the properties of small-x gluons are dominated by gluon saturation, i.e. the dominance of nonlinear interactions in the gluon field. In order to fully understand the behavior of small x gluons, a variety of different experimental measurements is needed. Of particular importance here is exclusive quarkonium production mediated by real or virtual photons. Such measurements are currently made in ultraperipheral heavy ion collisions [1] at the LHC and at RHIC. Exclusive measurements will also be an important part of the program at a future electron-ion-collider [2]. Exclusive quarkonium production is an important process for several reasons. As an exclusive process it depends on the gluon density quadratically and is thus more sensitive to nonlinearities than inclusive cross sections. Exclusive processes can, depending on exactly what final state of the target one measures, be sensitive to separately the average and the fluctuations of the gluon density in the target [3][4][5][6]. On the other hand, the heavy quark masses cut away nonperturbative long distance contributions and make the use of a weak coupling framework safer than for light quark processes [7,8].
The dipole picture of DIS [9][10][11][12][13] (a specialization of the light cone perturbation theory framework of [14] to the DIS process) provides a convenient framework to study deep inelastic scattering at high energy. In particular one expresses both inclusive and exclusive cross sections in terms of the same fundamental quantity, the dipole scattering amplitude, which gives this picture more predictive power than collinear factorization. With light quarks several recent advances have taken calculations of inclusive [15][16][17] and diffractive [18,19] observables to NLO accuracy, and it would be important to do the same for heavy quark cross sections.
At leading order the physical picture of exclusive scattering in the dipole picture is the following [20]: a virtual photon fluctuates into a quark-antiquark pair that interacts with the nucleus elastically with a cross-section σ qq . The resulting dipole can, later on, recombine into a quarkonium state. In this framework, the information about quarkonium is encoded in its light-cone wave function (LCWF), which encodes the overlap of the quarkonium state with various eigenstates of the free light cone QCD Hamiltonian.
In recent phenomenological applications (following e.g. [21]) one commonly uses phenomenological parametrizations for the LCWF, such as the the "boosted gaussian" or "gaussLC". A more model-independent approach would be to exploit the fact that heavy quarkonium is a non-relativistic system. In such a system the typical three-momentum of quarks is of the order of mv,where m the heavy quark mass 1 and v the typical speed of the heavy quarks around the center of mass (which is much smaller than the speed of light). The mass of the heavy quarkonium state M HQ is such that the binding energy is small: M HQ − 2m ∼ mv 2 m. Another important energy scale is Λ QCD which marks the transition between perturbative and non-perturbative physics in QCD.
Using an effective field theory approach non-relativistic QCD (NRQCD) [22] factorization formulas for some processes can be proven. The NRQCD approach has been quite successful, although there is some tension with polarization related observables in charmonium (for a more recent review see [23]). Diffractive quarkonium production has been studied in the non-relativistic limit in both a covariant [24] theory approach and in a light cone formalism [20]. Also velocity expansion corrections have been discussed in [25], although they are a NNLO effect in α s (mv). However, NLO corrections in α s (m) due to radiative corrections to the computation in [20,24] are still missing.
In this paper, we are going to compute the LO relativistic corrections to quarkonium light-cone wave function in the light-cone gauge such that it can be used to obtain NLO predictions in the dipole picture. This will allow to compute the radiative NLO corrections to the leading nonrelativistic result. We are going to check that using our results for the quarkonium light-cone wave function we can obtain the well-known literature result of the NLO correction to quarkonium decay to e + e − [26]. We will then show that the light-cone distribution amplitude (that can be obtained in the limit of small transverse coordinate) fulfils the ERBL evolution equation [27,28] (a NRQCD study of the same quantity can be found in [29][30][31]). We are going to show that using this wave function it is possible to compute NLO corrections to exclusive quarkonium production by analyzing the divergence structure and checking that it is consistent with B-JIMWLK evolution [32][33][34][35][36][37][38][39][40][41][42] of the target. A more complete phenomenological analysis will be possible when the NLO corrections to the photon wave function with massive quarks, currently under investigation, become available.
The paper is organized as follows. In the next section, we review the specific formulas of the dipole approach needed for this computation and fix the notation. In section III we are going to define the procedure to obtain the non-relativistic wave-function and its corrections and we are going to compute the relativistic contribution. In section IV we are going to show that we can use the previous results to obtain the light-cone distribution amplitude and we are going to check that it satisfies the Efremov-Radyushkin-Brodsky-Lepage (ERBL) evolution equation. Section V provides a strong cross-check of our results by applying them to the computation of the radiative corrections to S-wave quarkonium decay into leptons. In section VI we study the divergence structure of exclusive quarkonium production in both the dilute and in the non-dilute limit. Finally, we conclude in section VII.

II. THE DIPOLE APPROACH
We are going to perform the computations using light-cone coordinates, for a given vector p µ they are defined as which implies As a consequence, the momentum integration measure takes the form and the scalar product of two vectors in these coordinates is It is useful to define the light cone vectors such that n ·n = 1 , Figure 1. Representation of the leading order amplitude for exclusive quarkonium production in the dipole model. A virtual photon splits into a heavy quark-antiquark pair that later interacts with the nucleus (represented by a square box). After this interaction, the pair forms a quarkonium state (represented by a blob).
We also find it sometimes useful to separate from four-momenta in the on-shell part, denoted by a hat such thatp In light cone perturbation theory the off-shellness is included in only the ligh-cone energy (− orn µ -component) of the four vectorp In order to apply the dipole approach to exclusive quarkonium production, we are going to follow the discussion in Ref. [21]. We start with eq. (12) there which corresponds to the leading order contribution, in the limit of no target recoil t = 0: where the sub-index T or L refers respectively to a transverse or longitudinal polarization of the photon and vector meson 2 . Here Ψ γ * and Ψ HQ are respectively the light-cone wave functions of the photon and quarkonium, depending on the transverse coordinate separation r ⊥ and p + momentum fractions z, 1 − z of the quark and antiquark in the meson. The properties of the gluon field of the target are encoded in the "dipole cross section" σ qq which is, more properly speaking, twice the imaginary part of the forward elastic scattering amplitude of the quark-antiquark dipole with the target. Importantly for the predictive power of the framework, the same quantity appears in other cross sections, both for DIS and other high energy scattering processes. A graphical representation can be found in fig.  1. The previous formula actually assumes that only the subspace of the Fock space qq is relevant for quarkonium light-cone wave function. A more general formula would involve a sum over all the possible Fock states of the quarkonium: 2 The interaction with the target is eikonal, thus the polarization of the meson is the same as that of the photon Exclusive quarkonium production is a multi-scale process. In the previous formula, we identify many different energy scales. Firstly, Q is the virtuality of the photon and is encoded in its wave-function. The saturation scale of the target Q s is, in the previous formula, hidden inside the r-dependence of the amplitudes σ qq and σ qqg . Physically it represents the typical size of the dipoles selected by the target. Dipoles much smaller than 1/Q s interact weakly due to "color transparency". For dipoles larger than 1/Q s the growth of the amplitude is limited by unitarity, i.e. gluon saturation, and ceases to compensate for the suppression from the photon and quarkonium wavefunctions. Finally, inside Ψ HQ there are hidden many scales related to quarkonium physics, namely the scales m, mv , mv 2 and Λ QCD , which were discussed before. In this work, we are going to assume that Q s , m mv, Λ QCD and we are going to explore the physical consequences of this regime.
Details on how to obtain the light-cone wave function of quarkonium in this approximation will be given in the next section, however, we can explain the physical implications here. The quarkonium wave function is dominated by quarks that are moving with a non-relativistic velocity. This implies that it extends over distances of order 1 mv and it is not very sensitive to variations at smaller distance scales. Another consequence of the non-relativistic nature of quarkonium is that at leading order Ψ HQ only has support for values of z around 1 2 . Therefore where φ qq is the leading order light cone wave function that only takes into account non-relativistic contributions and λ = z − 1 2 . Extending this idea to higher orders, the quarkonium wave function for short distances In this formula the index n and m represent the particle content; for higher Fock states the single two-particle relative coordinates z, x ⊥ are replaced by the appropriate variables. The notation Ψ HQ refers to the full light-cone wave function of quarkonium, while the φ's are the wavefunctions restricted to the case in which all particles are nonrelativistic (therefore λ 1). Here f (z) is just a test function written to represent the fact that the equality is only valid when integrating for z. At this stage, we are ignoring spin degrees of freedom in the notation, but we will be more explicit about this in the following sections. Note that each power of ∇ M acting on φ gives a suppression of order v (in this equation ∇ is a three-vector, such that ∇ = ∇ ⊥ + i2mλ). Due to the non-relativistic nature of quarkonium the coefficients start with just a heavy quark-antiquark-state which is, to leading order, nonrelativistic: with transitions to higher Fock states entering at higher orders in perturbation theory.
In this paper, we are going to argue that in order to compute exclusive quarkonium production at next-to-leading order (NLO) only C 0 qq←qq and C 0 qqg←qq at order α s are needed and we are going to compute them. The fact that the contributing Fock states at NLO are the qq and qqg ones is similar to other NLO calculations in the dipole picture (e.g. [15][16][17][18][19]). The additional feature in the case of quarkonium is the way how both are related to a common nonrelativistic bound state wave function. Taking into account these two Fock states we can write the cross-section as These terms can be reorganized in such a way that we have an expansion that resembles a typical production process in NRQCD a) b) Figure 2. Comparison between how the degrees of freedom are separated in NRQCD and in the dipole model for the particular case in which the heavy quark pair emits a soft gluon before crossing the target. a) In NRQCD degrees of freedom are separated according to their virtuality. In this picture, red particles have a virtuality larger than m 2 v 2 (we are assuming, like in the rest of the text, that Qs mv) and their influence would be encoded in the hard coefficient. The rest of the particles are represented in blue and will influence the NRQCD matrix elements. b) In the dipole model degrees of freedom are separated according to time. Particles before crossing the target belong to the wave function of the photon and we represent them here in red while particles after the target are encoded in the wave function of quarkonium and are coloured blue. The method that we are implementing in this manuscript combines the two schemes by separating relativistic from non-relativistic degrees of freedom in the wave function of quarkonium.
In this equation, the third line only encodes what in NRQCD would be called soft physics (scales smaller than m), however, the second line receives contributions from both hard and soft physics. For example, there will be a soft physics contribution in the photon wave-function whenever a soft gluon is emitted before crossing the target. Therefore the analogy with NRQCD is not complete. This is illustrated with an example in fig. 2. Note also that the dipole picture is formulated in light-cone gauge and therefore we cannot, in principle, use gauge invariance arguments that are common in NRQCD and that will impose relations between the different terms in the expansion. For example, the fact that derivatives acting on the low energy matrix elements have to be combined with gauge fields to form covariant derivatives connects states with a different number of particles in such a way that we would expect that C 1 qq←qq is constrained by C 0 qq←qqg , however since we are using a setting that is not gauge invariant we are not going to use this kind of relations (for a more complete discussion see [25], and also the discussion about genuine and kinematical twist in [43]).
Let us now discuss the power counting of this computation and which terms do we need to consider if we want to achieve a NLO result. In this problem radiative corrections in the hard part will be suppressed by powers of α s (µ)| µ≥m , each additional ∇ m acting on φ would amount for a suppression of order v ∼ α s (mv). From now on we will ignore, for purposes of order estimation, the scale at which α s is evaluated. With this in mind, we should in principle take into account the first radiative correction in the hard part and the first velocity correction in the soft part. However, the first velocity correction only enters at NNLO. This is because due to rotational invariance the contribution from k = 1 and k = 0 (or vice-versa) in Eq. (20) is zero. At LO the only contributing Fock state is n = m = qq (and correspondingly in the complex conjugate amplitude denoted by primes). At NLO n and m can also have the value qqg. However, we can also check that the case m (or m ) equal to qqg is further suppressed. The reason is that the soft gluon has to be emitted either from the photon wave function (before the target) or by the relativistic part of the quarkonium wave function. In both cases, it would correspond to the emission of a soft gluon by a particle with a virtuality of order m 2 or higher, and this process is always suppressed by powers of v. In summary, to obtain a NLO result we only need to consider the following formula Here we have neglected gradient corrections (higher values of k, k in Eq. (19)), since in the strict weak coupling power counting, valid for mv Λ QCD the quark velocity is v ∼ α s (mv), and thus v 2 α s (m). It is worth keeping in mind that for practical applications in the case of charmonium it might happen that numerically v 2 can be comparable to α s (m), making the gradient corrections important. Nevertheless, it is always true that the NLO radiative corrections are bigger or of the order of the first corrections due to the velocity expansion and, therefore, a requisite of any improvement over the LO result.
The power counting discussed here is, to our knowledge, not clearly expressed in the literature, where typically only the LO limit is considered. Thus Eq. (21) is the first main result of this paper. To make this statement quantitative, we will next proceed to calculate at NLO the coefficient functions C 0 qq←qq and C 0 qqg←qq . We will then show that the cross section formula (21) gives finite results for any σ qq and σ qqg that fulfil B-JIMWLK evolution. However, it is also interesting to discuss the relation with collinear factorization [28,44] that can be used when Q is larger than any other scale in the problem. In that formalism the cross-section can be understood as the convolution of a hard function with the light-cone distribution amplitude (LCDA), which in the light-cone gauge corresponds to the light-cone wave function in the limit r ⊥ → 0. The LCDA must fulfil the ERBL evolution [27,28]. In [29][30][31] the LCDA and ERBL evolution were studied within the formalism of NRQCD. We are going to check that LCDA that we can obtain from our light-cone wave function fulfils ERBL evolution.

III. THE WAVE FUNCTION OF HEAVY QUARKONIUM
The two body part of the wave function of quarkonium can be determined using the Bethe-Salpeter equation. In the non-relativistic limit, the relation between the light-cone wave function and the wave function that can be obtained in a potential model was studied in [45]. Here we use this result as a starting point. Following [21], we use the convection that the light cone wave function is normalized such that where in this formula we have written explicitly the sum over polarizations. We are going now to analyze in more detail the LO light-cone wave function of a vector meson. Following [45] we get in the spinor matrix space (instead of helicity space) wave function whereφ is a function normalized as it is usually done in the context of potential models and λ is a constant that we are going to fix such that eq. (22) follows from eq. (24). We can transform the last expression to light cone helicity space using the prescription in [28], obtaining The normalization condition then imposes that λ = √ 32m 3 so we finally obtain for the longitudinal component while for the transverse component, we get The previous functions fulfil the Bethe-Salpeter equation represented in fig. 3 which can equivalently be written as a Schrödinger equation [45]. Now we want to go beyond the non-relativistic approximations and consider the influence of relativistic quarks and gluons with an energy of the order of the heavy quark mass. We can differentiate three cases: 1. Components of the wave function that correspond to the probability amplitude of having a relativistic particle inside the quarkonium state.
2. Contributions (from relativistic particles in loops) resulting in a wave function renormalization of a nonrelativistic particle of which the probability amplitude to be in a quarkonium state is considered. An example of such a contribution is shown in fig. 4.

3.
Contributions that can be encoded as a correction to the kernel K in fig. 3. An example is shown in fig. 5.
We will not discuss contributions of type 3. Their effect can be encoded in a redefinition of the potential (or equivalently, the kernel in the Bethe-Salpeter equation). Hence, in the case of production and decay processes this information is hidden in the value of the non-relativistic wave function at the origin. Contributions of type 2 are important but they just renormalize multiplicatively the non-relativistic wave function by a constant equal to the wave = K Figure 3. Bethe-Salpeter equation that gives the non-relativistic wave function. In this picture all particles are non-relativistic, therefore, the kernel K is just the Green function of the scattering of two non-relativistic heavy quarks.  function renormalization 3 . At leading order they correspond to the diagram in 4 which we compute it in appendix A, giving In this last expression, we have regulated transverse momentum in dimensional regularization such that the dimensions of the transverse space is D −2. The +-component of the momentum has been regulated using an infrared cut-off x 0 in the momentum fraction of the gluon (with respect to the initial quark). This is a similar regularization scheme as used in the NLO LCPT calulations in [15][16][17]. The wave function renormalization that we obtained is similar to the one presented in [46] but not exactly equal. However, our results agree on the double-logarithm contribution. Contributions of type 1 can be computed by considering diagrams in which all relativistic particles annihilate to form non-relativistic ones. They will be proportional to the non-relativistic wave function at the origin. The reason for this simple structure is that in diagrams in which only relativistic particles are involved the only energy scale that appears is m. Since no momenta are parametrically of order mv or mv 2 , there are no diagrams enhanced by inverse powers of the velocity v, which would require a resummation. Furthermore, since all relativistic processes are short distance processes from the point of view of non-relativistic wave function, it is a good approximation (up to additional powers of v) to consider that the non-relativistic particles are created at the same point. In summary, contributions of types 1 and 2 will be encoded in the functions C k n←m that were introduced in the previous sections. The difference is that while contributions of type 2 give rise to corrections that are just constants, contributions of type 1 will have a non-trivial dependence on r ⊥ and z. On the other hand, contributions of type 3 will be encoded in the non-relativistic wave function.
For the computation at hand, we are interested in the components of the light-cone wave function in which we have two relativistic quarks or a relativistic quark, a non-relativistic one and a hard gluon. Let us focus first of the former case. Its LO contribution is represented by the diagram in fig. 6. Details of the computation will be given in appendix B, the result from this contribution gives where p is the momentum exchanged by the gluon and the light cone gauge gluon polarization sum is Note also that the function φ i qq (λ, 0) has to be understood as φ i qq (λ, r ⊥ = 0). For future manipulations, it might be useful to write this result in coordinate space. To do this we make use of the formulas in Appendix C to rewritē and Using the formulas of Appendix D and the expressions for the spinor matrix elements in Appendix E it is possible to find the result in coordinate space where τ = 2m z − 1 2 r ⊥ . The general expression (34) can be further simplified if we distinguish the cases of longitudinal and transverse polarization. For the longitudinal polarization state we obtain The one-loop wave function for the transverse polarization state is Combining the results of eq. (29) and (30), we obtain that This is the first major result of this section. An explicit expression for the structures in (37) in coordinate space and for the meson polarization states can be easily obtained using (35) and (36). The result (37) represents one of the two NLO corrections to the heavy quarkonium wavefunction, needed for the NLO calculation of exclusive quarkonium production. In the non-relativistic limit the polarization of the meson is reflected only in the spin state of the quarks, and the spatial part of the wavefunction becomes rotationally symmetric. We can utilize this simplification and integrate the wavefunction over the polar angle, which results in simpler expressions. The general azimuthally symmetric part of the wavefunction is Decomposing this into meson polarization states we get for the longitudinal polarization while for the transverse case, we obtain It is also possible to write a compact expression in coordinate space for the polar angle average of the coefficients = Figure 7. LO contribution to the component of the wave function with one relativistic heavy quark, one non-relativistic heavy antiquark and a hard gluon.
C 0 qq←qq (z, p ⊥ ; λ 1 , λ 2 ; λ 1 , λ 2 ). For the case of longitudinal polarization, the result is and for the transverse polarization Now let us move to the second contribution, where we have a relativistic quark, a non-relativistic one and a hard gluon. This is represented at leading order in the diagram of fig. 7. We need to introduce some notation in order to write the result. The longitudinal momentum fractions of the relativistic quark, the gluon and the non-relativistic antiquark are respectively z RQ = 1−x 2 + λ, z G = x 2 and zQ = 1 2 − λ. From the nonrelativistic nature of the "spectator" antiquark it follows that λ 1. We shall consequently not write out the integration limits for λ, since it is understood to be dominated by a short interval around the origin. The transverse momenta are respectively P ⊥,RQ = p ⊥ 2 − q ⊥ , P ⊥,G = p ⊥ 2 + q ⊥ and P ⊥,Q = −p ⊥ , where again for the nonrelativistic antiquark p ⊥ m. Using this we get dλ 4π or equivalently dλ 4π In order to write the result in coordinate space we define r ⊥ as the Fourier conjugate of the nonrelativistic momentum p ⊥ and l ⊥ as the one of the relativistic one q ⊥ , leading to It is straightforward to obtain C 0 qqg←qq from the previous results as A cross-check of the results of this section can be obtained by checking that the normalization of the light-cone wave-function is not changed by radiative corrections. To order α s this is ensured if the square of the diagram in fig.  4 cancels with the square of the diagram in fig. 7. An easy way to check this is to start from eq. (43), after doing the square and performing the trace over gamma matrices one immediately obtains an equation compatible with (3.28) in [46], which is equal to −δZ.

IV. LIGHT-CONE DISTRIBUTION AMPLITUDES
In the cases in which the virtuality of the photon is much larger than other dimensionful scales in the problem, in particular the mass of the heavy quark, one can use collinear factorization [28,44]. The information about the bound state in this formalism is encoded in the light-cone distribution amplitude (DA). This coincides with the light-cone wave function in the limit in which the transverse radius r ⊥ goes to zero. This object was already studied in [29][30][31] using the non-relativistic expansion (see also [47]), here we provide an independent computation in light-cone gauge.
Naively, one would expect that one can obtain the light-cone distribution amplitude trivially making the limit r ⊥ → 0 in the results in the previous section. This is not so because taking the limit introduces an ultraviolet divergence that needs to be regularized. This can be done in dimensional regularization using the following prescription Applying this method to eq. (30) one obtains for the longitudinal polarization and for the transverse one We now define the light-cone distribution amplitude as Using the previous results, the light-cone distribution amplitude D(z) can be written as For the longitudinal polarization, this is while for the transverse one the DA is The first lines in Eqs. (52) and (53) are the leading order bare amplitude (apart from the δZ). The NLO corrections involve a divergence proportional to ∼ 1/(D − 4) from the transverse integrals, and the associated logarithm of the regularization scale µ. In the usual fashion these divergences must be absorbed into the bare distribution, leaving the renormalized distribution dependent on the scale µ. This dependence is related to the ERBL evolution equation [28,44]: To check that we recover this equation in the standard form we can take the expressions for the kernels in Eq. (54) from the literature and compute the r.h.s. The comparison of this result with the coefficient of the logarithm of µ in Eqs. (52) and (53) then provides a check of our calculation of the wave function.
The ERBL kernel for the longitudinal polarization is given by while for transverse one it is The [] + notation means that Note that the [] + prescription can be implemented with a cut-off, as we did when computing δZ. For example, We can then check that the light-cone distribution amplitudes that we have computed fulfil the ERBL equation by computing explicitly ∂D i (z) ∂ log µ 2 . In the case of longitudinal polarization while for the transverse one Thus we have verified that we recover a standard ERBL evolution of the light cone wave function. This is a check of our NLO wavefunction for the quark-antiquark Fock state in the short distance limit that is probed in a high virtuality (large Q 2 ) process. We emphasize that the quark momentum scale (or coordinate separation) in the whole wave function (37) is controlled by the quark mass, and this full result can be used more generally than just in the large Q 2 limit.

V. RADIATIVE CORRECTIONS TO S-WAVE QUARKONIUM DECAY
The decay width of quarkonia to leptons offers an important constraint on the nonrelativistic wavefunction. In this section we will see how the decay width is related to the light cone wavefunction. We will see that our regularization involving a cutoff in the p + momentum modifies the relation between the nonrelativistic wavefunction and the decay width compared to a usual rotationally invariant nonrelativistic potential model calculation.
The decay of quarkonium in a S-wave state into leptons is expressed in terms of the meson decay constant, which is defined as the matrix element of the electromagnetic current operator with the vector meson state [48]. The current operator is local. As a consequence, at leading order, the decay constant is related to the value of the wave function at the origin. This is true both for the light cone wave function and the nonrelativistic wavefunction. Beyond leading order, the relation between the decay constant and the nonrelativistic wavefunction receives a NLO correction [26]. In light cone perturbation theory, on the other hand, the equivalent diagrams become corrections to the light cone wavefunction itself, not to the relation with the decay constant. Equating the decay constants in terms of the light cone and nonrelativistic wavefunctions yields the equation where the sum in n is over all the Fock state components that contribute to the wave-function at the origin. In the case of longitudinal polarization, the only component at NLO is n = qq. For the transverse polarization, there is also a contribution where the overlap between the electromagnetic current operator and the vector meson state receives a correction from the n = qqg Fock state, mediated by an instantanous interaction. The different one-loop corrections are shown diagrammatically in fig. (8).
In the following we shall concentrate, for simplicity, on the longitudinal polarization state. In this case we have, by construction Therefore, we can do the computation using the light-cone distribution amplitude that we computed in the previous section. It is useful to note that Using this on eq. (52) we get We see that, apart from a finite piece that coincides with the result in [26], we also get a divergent term. Its origin is the Coulomb singularity, which would not be present at this order in dimensional regularization. Degrees of freedom with transverse momentum much smaller than m and z − 1 2 close to zero should be included only in φ. However, when we computed eq. (34) we did not include any mechanism to subtract them. Looking at eq. (30) we can see that, the only scale that affects the transverse momentum when computing Ψ qq HQ (z, 0 ⊥ ) is 2m z − 1 2 . Then, because we are using dimensional regularization in the p ⊥ integration, if we choose a cut-off x 0 much bigger than the non-relativistic velocity around the centre of mass v we can make sure of only taking into account relativistic degrees of freedom. However, dλφ(λ, 0) will also depend on x 0 . The condition that must be fulfilled is that the decay constant is independent of the cutoff x 0 : In order to compute how dλφ(λ, 0) depends on the cut-off we can make use of Bethe-Salpeter equation represented in Fig. 3, following the lines of the discussion in Appendix B. Note that from the non-relativistic point of view, z ∼ x 0 or 1 − z ∼ x 0 represent the ultraviolet since the quark relative longitudinal momentum is large. Therefore, we can substitute the kernel by a one gluon exchange between the quark and the antiquark and the non-relativistic part of the wave-function renormalization. Let us begin with the computation of the one gluon exchange part of the r.h.s. of the Bethe-Salpeter equation of Fig. 3, that we name Ψ qq sub , From this expression we can obtain, and therefore The non-relativistic part of the wave-function renormalization is computed in Appendix A. Here we use that Putting the two one-loop corrections pieces together we get We can then see that eq. (65) is fulfilled and the dependence on the cutoff x 0 cancels. Our regularization scheme with a cutoff x 0 is not the same as conventional rotationally invariant dimensional regularization. Cutoff-dependent, not directly measurable, quantities do not need to be the same in both schemes. The nonrelativistic wavefunction at the origin is such a quantity; and indeed the relation between our wavefunction and the dimensionally reduced one of e.g. [26] is where the subscript DR means that this would be the result that we would have obtained if we had regulated all divergences using dimensional regularization. Then the light cone wave function at the origin, i.e. the vector meson decay constant, is related to the regularized nonrelativistic wavefunction in the natural way In conclusion, we recover the result in [26] after taking into account the different renormalization prescription. To our knowledge, this is the first computation of this relation (which can be directly related with the decay into leptons) done in light-cone perturbation theory.
Finally, let us discuss a more phenomenological explanation. We can use the experimental value of the decay into electrons to obtain the wave function at the origin as a function of x 0 The formalism is consistent if we use the previous formula to compute the relation of another observable with the decay width into electrons and we get a result that does not depend on x 0 . In the following section we are going to see that this is indeed the case.

A. Leading order
In this section, we study the cancellation of divergences in photoproduction of quarkonium with longitudinal polarization at NLO. In particular, we will show that the dependence on the dimensional regularization scale µ cancels, and that the dependence on the longitudinal momentum cutoff x 0 can be factorized into the B-JIMWLK-evolution of the target. To simplify the calculation we assume that Q m. This assumption does not modify ultraviolet or collinear divergences and has the advantage that we can use the same photon wave-function that was used when studying massless quarks. Thus our results can to some extent be seen as a reformulation of the result obtained for light vector mesons in a somewhat different language in Ref. [19].
The LO result is given by eq. (15). The LO qq component of the photon wave-function is [21]: where e f is the charge of the corresponding quark. Regarding the cross-section σ qq , its specific form is model dependent and in general phenomenologically realistic models do not allow for an analytical treatment. Instead, we will use here the small r ⊥ (dilute) limit In this paper, we are going to discuss in general the structure of the divergences, in order to check the consistency of the method, and we are going to apply the dilute limit to get analytic results when possible. Let us remind the reader of the LO result, which was computed in [20]. This can be obtained by writing the LO photon wave-function in eq. (15) In the dilute limit, this becomes Using the relation between Q s and the gluon density from [21] and the relation between quarkonium's wave-function at the origin and its decay into electrons we get which agrees with the classic result from [24] (and is four times the result found in [20]). Now we move on to discuss the different corrections that this process gets at NLO, focusing on the longitudinal polarization state at Q m where we can use results from the literature for massless quarks. At NLO, the production process gets contributions from both the qq and the qqg Fock states. In addition, we have to consider separately the corrections to the vector meson and the photon light cone wave functions, although in the large Q 2 limit they behave in a similar way.
In our case, we are going to need the expression for z = 1 The effect of this correction on quarkonium production is the following In the dilute limit, this is C. NLO corrections to the qq component of quarkonium wave-function.
In this subsection, we will compute one-loop corrections to the qq component of quarkonium wave-function in the hard scattering limit. We will first study the divergence structure in the general case, and then perform the computation in the dilute limit for the dipole-target scattering amplitude. Using Eq. (19) the contribution that we are studying in this subsection is First, we look at the dependence on µ. For this, we note that σ qq goes like r 2 ⊥ in the small r ⊥ region while both C 0 qq←qq and (Ψ γ * ) L go like log r ⊥ . On the other hand, for large r ⊥ both wave-functions are exponentially suppressed. Therefore we conclude that the only dependence on µ comes from the wave-function renormalization of the heavy quarks. Therefore (86) . This means that the dependence on µ of the wave-functions of the photon and heavy quarkonium is the same. This should not be surprising, since in the hard scattering limit both come from a perturbative one-gluon exchange between the quarks. Now we focus on the dependence on x 0 . For the longitudinal polarization state where the wavefunction is azimuthally symmetric it is convenent to work again with the angularly integrated wave function We first note that Using this and computing dδZ dx0 we get The first term in the rhs will be cancelled by the dependence of the non-relativistic quarkonium wave-function on x 0 , c.f. Eq. (65). The other terms have the same dependence on x 0 as the contribution of the NLO corrections to the photon wave-function. Now we move to the dilute limit, in which we wish to obtain an analytic expression in the limit Q m. Note that this implies that if Qr ⊥ ∼ 1 then mr ⊥ 1 while if mr ⊥ ∼ 1 then the photon wave-function is exponentially suppressed. Therefore, we can use the quarkonium wave-function in the mr ⊥ 1 limit. As a starting point, we take Eq. (41) It is convenient to separate this into a r ⊥ independent and r ⊥ dependent piece as lim r ⊥ →0 dθ r C 0 qq←qq (z, r ⊥ ; λ 1 , λ 2 ; λ 1 , λ 2 ) long = 8π 2 δ(z − where The final result for the qq component of the vector meson wave function that we obtain is Note that the divergences of this part, both in µ and x 0 , are the same as in the photon wave function (84).
D. Contribution of the qqg Fock state Figure 9. Transverse momentum, transverse coordinate and longitudinal momentum fraction assignments for the qqg Fock state, and the definitions of the coordinate separations r ⊥ , l ⊥ , u ⊥ in terms of these.
We now move to the qqg Fock state, which interacts with the target with the scattering amplitude σ qqg . Note that the NLO correction to the cross section is given by the interference of this NLO correction to the amplitude with the leading order amplitude that only involves a qq state (and of course the complex conjugate). This contribution to the cross section is given by where r ⊥ and l ⊥ are defined as in Eq. (45), with r ⊥ roughly interpreted as the large distance separating the nonrelativistic antiquark from the system of a relativistic quark and gluon, which has a small size l ⊥ . These coordinate assignments are demonstrated in Fig. 9. Since we are working in the hard scattering limit Q M , the wave-function of the photon can be taken from the massless quark calculation in [16], which using our notation and normalization and fixing the antiquark to be non-relativistic is: where and As in the inclusive DIS case [16,17] there are cancellations of divergences between the qq and qqg states. To make these manifest, it is useful to separate the contribution into two different terms, by adding and subtracting a term that involves a qq scattering amplitude Note that the argument of σ qq is the separation between the quark and the antiquark, which is u ⊥ = r ⊥ + l ⊥ 2 for the qqg Fock state wave-function, see Fig. 9.
First, we discuss the possible µ dependence. Note that this is related with the limit l ⊥ → 0. In this limit σ qqg − σ qq also vanishes. Therefore, we only need to focus on the second line in eq. (97). For small values of mxl ⊥ , using the formula in Appendix D, Regarding Ψ qqg γ * , we can approximate and Note that this approximation keeps the correct small l ⊥ behaviour without introducing additional ultraviolet divergences. Therefore, in this limit, we can write Performing the integration over l ⊥ , we obtain Note that since we performed approximations valid in the limit l ⊥ → 0 we can only claim at the moment to capture the ultraviolet behaviour, or in other words, the µ dependence. This is captured by the following equation (103) This is required so that the µ-dependence of the contribution from the qqg Fock state crossing the target shock will cancel against the µ-dependence of quark wavefunction renormalization, which is a contribution where only the qq dipole crosses the shock wave. Now we focus on the dependence on x 0 . In the x → x 0 limit, we can write where Looking at the original expressions for I j (a) and I j (b), we might wonder if by simplifying the pole in k ⊥ we are modifying the x 0 dependence. However, note that for small x 0 the integral is only sensitive to the exact value of the pole in the infrared. The infrared limit of I j (a) cancels with that of I j (b) for x close to x 0 , and therefore this sensitivity is cancelled out. In the previous formula this cancellation is given by the factor 1 − e ik ⊥ u ⊥ . Regarding C 0 qqg←qq , in the x → x 0 limit we can approximate it by Therefore, we get coefficient for the qqg component in the meson wavefunction as Using these formulae we can compute the dependence on the cutoff x 0 separately for the two parts of the decomposition (97) The first term gives Since the l ⊥ integration introduces no divergences (recall that (σ qqg − σ qq ) → 0 when l ⊥ → 0) we can take the exact limit D = 4. Then we get where in the last equality we have used the B-JIMWLK equation [35][36][37][38][39][40][41][42]. The other contribution, depending only on the qq amplitude, that was added and subtracted in the decomposition (97) becomes Finally, we put all the pieces together. Note that until now we have only explicitly shown the computation in the case in which the intermediate quark is relativistic and the antiquark is non-relativistic. However, the opposite case also gives a contribution which is, due to symmetry considerations, exactly equal. Taking this into account we get d dµ and

E. Summary
In this subsection, we check that the sum of all the contributions to the cross section at NLO accuracy has the expected dependence on the cutoff scales µ and x 0 . Let us first remind the reader that the leading order cross section is finite and does not depend on the dimensional regularization scale Combining this with the scale dependence from the photon wave function (83), the vector meson wave function (86) and the one-loop correction to the wave function (111) we see that the dependence on µ cancels d dµ The overall dependence on x 0 must also cancel, when the B-JIMWLK evolution of the target is taken into account. The leading order cross section depends on the cutoff x 0 both through the x 0 -dependence of the dipole cross section, and of the nonrelativistic wavefunction (this resulted from the fact that the decay width must be independent of x 0 , see (65)): Both of these x 0 dependences of the leading order cross section are proportional to α s , and are thus at the same level as the NLO contributions. Relevant NLO contributions are the ones from the photon wave function (i.e. quark wavefunction renormalization), Eq. (83), from the "vertex correction" to the qq component of the meson light cone wave function in (88), and from the qqg Fock state, Eq. (112). All of these together are needed to make the cross section at NLO independent of the cutoff, x 0 , naturally up to terms of higher order in α s : This demonstrates that the leading high energy behavior (i.e. dependence on the cutoff x 0 ) can be factorized into B-JIMWLK evolution of the target. Note that our discussion here has been framed in the language of a simple small-x factorization procedure where the amplitudes σ qq and σ qqg depend on the cutoff x 0 . It is known for other processes (see e.g. [49]) that this can result in an oversubtraction that can give rise to unphysical behavior for the NLO cross section. We expect this problem to require, just as for inclusive hadron production or inclusive DIS [16,[50][51][52], taking the amplitude to in fact depend on the gluon longitudinal momentum fraction (λ G in Eq. (46)) in the appropriate way. This modification is formally of higher order in α s , and does not affect our conclusion that the leading high energy logarithm can be factorized into B-JIMWLK evolution.

VII. CONCLUSIONS
In this manuscript, we have studied the light-cone wave function of heavy quarkonium in the non-relativistic limit with a focus on future applications within the dipole model. In section II we discussed how to combine the dipole model with the non-relativistic power counting. Within this scope, the main results of this work are the NLO correction to the qq component, given in eqs. (41) and (42), and the leading contribution to the qqg component, given in eq. (44) in momentum space and in eq. (45) in coordinate space.
We have performed several cross-checks of these results: • We have checked that the light-cone distribution amplitudes that are obtained from the light-cone wave functions fulfil ERBL equation.
• We have confirmed that, for the longitudinal case, we recover the results of [26] for radiative corrections to S-wave quarkonium decay.
• Finally, we have also checked that using the light-cone wave functions that we found in the longitudinal case to compute exclusive quarkonium production in the limit Q m we get that all divergences vanish at NLO assuming B-JIMWLK evolution of the target.
The importance of our results lies in a future application to the computation of quarkonium exclusive production at next-to-leading order. This would be the heavy quark analogue of the exclusive light meson production calculation of Ref. [19], which is performed in a slightly different formalism but uses the same physical picture of high energy scattering. The heavy quark mass introduces additional complications on one hand, but simplifies the description of the meson on the other, due to the nonrelativistic nature of the bound state. For this the photon light-cone wave function with finite mass effects is needed, this result will appear in the near future, following the massless case [16]. Overall, we believe our work here is a necessary part of a broad effort to increase to next-to-leading order the accuracy of computations in the dipole picture for a variety of processes. Such improvements are urgently needed to fully confront QCD theory in the saturation regime with experimental data from an active program at the LHC and from a future Electron-Ion Collider. Finally, Using the previous results and expanding around D ∼ 4 we can obtain the terms that are needed to compute δZ A (m 2 ) = ig 2 C F 8π 2 m 2 1 D − 4 + 1 2 log m 2 4πµ 2 + γ E 2 (2 log x 0 + 1) + log x 0 + 1 2 + (log x 0 ) 2 , (A19) Finally, we obtain the wave-function renormalization (2 log x 0 + 1) + 2 log x 0 + 2(log x 0 ) 2 and