Relativistic corrections to the vector meson light front wave function

We compute a light front wave function for heavy vector mesons based on long distance matrix elements constrained by decay width analyses in the Non Relativistic QCD framework. Our approach provides a systematic expansion of the wave function in quark velocity. The first relativistic correction included in our calculation is found to be significant, and crucial for a good description of the HERA exclusive $\mathrm{J}/\psi$ production data. When looking at cross section ratios between nuclear and proton targets, the wave function dependence does not cancel out exactly. In particular the fully non-relativistic limit is found not to be a reliable approximation even in this ratio. The important role of the Melosh rotation to express the rest frame wave function on the light front is illustrated.


I. INTRODUCTION
At large densities or small Bjorken-x, non-linear QCD dynamics is expected to manifest itself in nuclear structure. To describe the QCD matter in this non-linear regime, an effective field theory known as the Color Glass Condensate (CGC) has been developed, see e.g. [1,2]. Diffractive scattering processes at high energies are especially powerful probes of this region of phase space. The advantage in diffractive, with respect to inclusive, scattering is that since no color charge transfer is allowed, even at leading order in perturbative QCD at least two gluons have to be exchanged with the target. Consequently, the cross sections approximatively probe the square of the gluon density [3], and can be expected to be highly sensitive to non-linear dynamics.
An especially interesting diffractive process is exclusive vector meson production in collisions of real or virtual photons with the target, where only one meson with the same quantum numbers as the photon is produced. In these processes only vacuum quantum numbers are exchanged between the target and the diffractive system. Thus the target can remain intact, and the transverse momentum transfer can be used to probe the spatial structure of the target. This momentum transfer is by definition the Fourier conjugate to the impact parameter. As such, it becomes possible to study the target structure differentially in the transverse plane. A particularly important channel is the production of J/ψ mesons. The charm quark is heavy enough to enable a weak coupling description of its elementary interactions. Nevertheless the quark mass is not large enough to make the process insensitive to saturation effects. Also experimentally the J/ψ is relatively easily identifiable and produced with large enough cross sections to be seen.
One major source of model uncertainties in the theoretical description of vector meson production follows from the nonperturbative vector meson wave function. For the J/ψ , a natural first approximation is to treat it as a fully nonrelativistic bound charm-anticharm state, which is the limit taken in the seminal work in Ref. [3]. The calculation of Ref. [36] recovers the same nonrelativistic result in the dipole picture (see also Ref. [37]). Already early on, it has been argued that this nonrelativistic approximation obtains important corrections from the motion of the charm quark pair in the bound state [38,39]. More recently, much of the phenomenological literature on J/ψ photoproduction has used phenomenological light cone wave functions to describe the meson bound state. This has the advantage that the light cone wave function is invariant under boosts in the longitudinal direction, and is thus naturally more suited to high energy collision phenomena. A disadvantage of some recent phenomenological parametrizations has been that they do not fully use the information on the nonperturbative bound state physics, most importantly decay widths, of quarkonium states that are usually analyzed in terms of nonrelativistic wavefunctions.
In recent literature, the applied different phenomenological wave functions result in e.g. J/ψ production cross sections that differ up to ∼ 30% from each other [27,40,41]. This is a large model uncertainty, compared to the precise data that is already available from HERA and the LHC, and especially given that the Electron Ion Col-lider (EIC) [42,43] is in the horizon (and similar plans exist at CERN [44] and in China [45]). The EIC will perform vast amounts of precise DIS measurements over a wide kinematical range, which calls for robust theoretical predictions.
To reduce the model uncertainty related to the vector meson wave function, we propose in this work a new method to constrain the wave function for heavy mesons based on input from the Non Relativistic QCD (NRQCD) matrix elements. These matrix elements capture non-perturbative long-distance physics, and can be obtained by computing the vector meson decay widths in different channels as a systematic expansion in both the coupling constant α s and the quark velocity v. As we will demonstrate, these matrix elements can be used to determine the value and the derivative of the vector meson wave function at the origin. As such, this approach provides more constraints than the phenomenological parametrizations widely used in the literature. In particular, starting from manifestly rotationally invariant rest frame wave functions, one by construction obtains consistent parametrizations of longitudinally and tranversally polarized vector mesons simultaneously, which is not obvious in many light cone approaches.
This manuscript is organized as follows. First, in Sec. II we review how vector meson production is computed in the dipole picture within the Color Glass Condensate framework, and how the cross section depends on the vector meson light front wave function. In Sec. III we first present how to obtain the rest frame wave function in terms of the NRQCD matrix elements, and then show how this is transformed to the light cone by applying the Melosh rotation [46,47]. We compare the obtained NRQCD based wave function to other widely used wave functions that are reviewed in Sec. IV. The numerical analysis including vector meson-photon overlaps and J/ψ production cross sections is presented in Sec. V.

II. VECTOR MESON PRODUCTION IN THE DIPOLE PICTURE
A. Exclusive scattering At high energies exclusive vector meson production in virtual photon-proton (or nucleus) scattering can be described in a factorized form. The necessary ingredients are the virtual photon wave function Ψ λ γ describing the γ * → qq splitting, the dipole-target scattering amplitude N and the vector meson wave function Ψ V describing the transition qq → V . The scattering amplitude reads [48] (note that the correct phase factor coupling the dipole size r to the transverse momentum transfer ∆ is determined in Ref. [49]) Here Q 2 is the photon virtuality, r the transverse size of the dipole, b the impact parameter and z the fraction of the photon light cone plus momentum carried by the quark. The photon polarization is λ, with λ = ±1 referring to the transverse polarization and λ = 0 to the longitudinal one.
In this work will study coherent vector meson V production. The coherent cross section refers to the scattering process where the target proton (or nucleus) remains intact. In this case, the cross section as a function of squared momentum transfer t ≈ −∆ 2 can be written as The dipole amplitude N depends on the longitudinal momentum fraction x P the target loses in the scattering process, which reads ( Here M V is the mass of the vector meson V and m N is the proton mass. The scattering amplitude A T,L is obtained from Eq. (1) by summing over the quark helicities, and in the case of transverse (T) polarization, averaging over the photon polarization states λ = ±1. In Eq. (2) two phenomenological corrections are included following Ref. [48]. First, β = tan πδ 2 is the ratio between the real and imaginary parts of the scattering amplitude. It can be obtained from an analyticity argument as The so called skewedness correction is included in terms of the factor R g , which reads, This correction can be derived by considering the vector meson production in the two-gluon exchange limit, assuming that the two gluons carry very different fractions of the target longitudinal momentum [50]. In this case, the cross section can be related to the collinearly factorized parton distribution functions scaled by the factor R g . In the dipole picture applied here, where the two quarks are color rotated in the target color field and undergo multiple scattering, this limit is not reached. In this work we include both of the real part and skewedness corrections widely used in the previous literature, but emphasize that these numerically large corrections should be used with caution when predicting absolute normalizations for the cross sections. In addition to coherent scattering, one can study incoherent diffraction where the target breaks up, but there is still no exchange of color charge between the produced vector meson and the target remnants. These processes are recently studied extensively in the literature as they probe, in addition to saturation effects [41], also the event-by-event fluctuations of the scattering amplitude resulting from the target structure fluctuations, see e.g. Refs. [51][52][53][54] or Ref. [55] for a review. As the focus in this work is on the vector meson wave function which enters in calculations of both incoherent and coherent cross sections similarly, from now on we only consider coherent scattering here.

B. Virtual photon wave function
The virtual photon splitting to a cc dipole is a simple QED process, and the photon wave function Ψ γ can be computed directly by applying the light cone perturbation theory (see e.g. [56,57]). Using the diagrammatic rules of light front perturbation theory and the conventions used in Refs. [58,59], the wave function can be written as Here e f is the fractional charge of the quark (in this work we consider only charm quarks with e f = 2/3), k, k and q are the quark, antiquark and photon momenta respectively, e = √ 4πα em and h andh refer to the quark and antiquark light front helicities [60]. The factor √ N c is included to obtain a squared wave function proportional to the number of colors N c . The spinors which are the eigenstates of light front helicity read, in the Lepage-Brodsky convention u h (k) = 1 v h (k) = 1 where the four-component helicity spinors readχ h=+1 = 1 √ 2 (1, 0, 1, 0) T andχ h=−1 = 1 √ 2 (0, 1, 0, −1) T , and α T = (γ 0 γ 1 , γ 0 γ 2 ). We use the light cone variables defined as In the light cone gauge, in which ε + = 0, the photon polarization vectors read where and Q 2 = −q 2 . The wave function can be evaluated by substituting the polarization vectors and explicit expressions for the spinors in Eq. (6) and setting the photon transverse momentum q to zero. It is convenient here to define a wave function in terms of the momentum fraction z and pull out a factor 4π. This should be done so that probability is conserved, so that we can write Ψ λ γ,hh (z, k) = 4πq + Ψ λ γ,hh (k + , k). In momentum space, the wave functions read The wave function in the mixed transverse coordinate, longitudinal momentum fraction space entering in the vector meson production cross section (1) is then obtained by performing a Fourier transform The mixed space wave function for the longitudinal polarization is Similarly, for the transverse photon with λ = ±1 the wave function reads We note that the these wave functions agree with those derived in Ref. [59] using the same convention, except for the overall sign in case of transverse polarizations which does not affect any of our results. On the other hand, when compared to the widely used wave functions reported in Ref. [48], the relative sign between the mass term and z terms in the λ = +1 case is different.
We emphasize that the quark light cone helicity structure above does not exactly correspond to the spin structure in the rest frame of the meson (there is no rest frame for the spacelike photon). In particular, when transformed to the meson rest frame, there are both S and D wave contributions in both longitudinally and transversely polarized photons. The transformation between the light front wave function expressed in terms of the quark light front helicities, and the rest frame wave function in terms of the quark spins is discussed in Sec. III B. We will discuss the decomposition of light cone wave functions, including the virtual photon one, into the S and D wave components in more detail in Appendix A.

C. Dipole-target scattering
The dipole-target scattering amplitude N in Eq. (1) is a correlator of Wilson lines, corresponding to the eikonal propagation of the quarks in the target color field. In principle, it satisfies perturbative evolution equations describing the dependence on momentum fraction x P , the so called JIMWLK equation [61][62][63][64][65][66][67], or the BK equation [68,69] that is obtained in the large-N c limit. These perturbative evolution equations, combined with a nonperturbative input obtained by fitting some experimental data, can in principle be used to evaluate the dipole amplitude at any (small) x P . This has been a successful approach when considering structure functions in DIS or inclusive particle production in hadronic collisions, see e.g. Refs. [31-35, 70, 71].
In diffractive scattering considered here one explicitly measures the transverse momentum transfer ∆, which is the Fourier conjugate to the impact parameter. Consequently, the dependence on the transverse geometry needs to be included accurately in the calculation. However, perturbative evolution equations generate long distance Coulomb tails that should be regulated by some non-perturbative physics in order to avoid unphysical growth of the cross section [72]. There have been attempts to include effective confinement scale contributions in the BK and JIMWLK evolutions and use the obtained dipole amplitudes in phenomenological calculations of e.g. vector meson production [73][74][75][76] (see also [77]). As the main focus of this work is in vector meson wave functions, we apply a simpler approach and use the so called IPsat parametrization to describe the dipole-proton scattering amplitude.
The IPsat parametrization [78] consist of an eikonalized DGLAP-evolved [79][80][81][82] gluon distribution, combined with an impact parameter b dependent transverse density profile. The advantage of this parametrization is that it matches perturbative QCD result in the dilute (small dipole size |r|) limit, and respects unitary in the saturation regime. The dipole amplitude in the IPSat parametrization reads where the proton transverse density profile is assumed to be Gaussian: with B = 4 GeV −2 . The initial condition for the DGLAP evolution is obtained by fitting the HERA structure function data [83][84][85][86], and the fit results in an excellent description of the total reduced cross section and the charm contribution [87]. The scale choice is µ 2 = C/r 2 + µ 2 0 , with the parameters C and µ 0 , among with the DGLAP initial condition, are determined in the fit performed in Ref. [87] (see also [88]).
Following Ref. [78] (see also [87]), the dipole-proton scattering amplitude can be generalized to coherent scattering in the dipole-nucleus case as (21) This estimate is valid in case of large nuclei, assuming that the dipole size |r| is not very large, which is the case in heavy vector meson production. Here T A (b) is the Woods Saxon distribution integrated over the longitudinal coordinate, with the normalization d 2 bT A (b) = 1. The nuclear radius used here is In order to calculate vector meson production, it is still necessary to determine the vector meson wave function. It can not be computed perturbatively, and consequently there are many phenomenological parametrizations used in the literature. The main goal of this paper is to obtain the meson wave function in a systematic expansion in quark velocities given by the NRQCD approach. We will also discuss, for comparison, some other wave function parametrizations in Sec. IV.

III. LIGHT CONE WAVE FUNCTION FROM NRQCD
NRQCD is an effective field theory describing QCD in the limit where quark masses are large, or v = p/m is small, where p is e.g. quark momentum and m is the quark mass. In this approach, it becomes possible to factorize cross sections into universal long distance matrix elements and perturbatively calculated process dependent hard factors.

A. Vector meson wave function in the rest frame
The J/ψ decay width in the NRQCD approach is written as an expansion in the quark velocity v [89]. At lowest order in v, the decay width is only sensitive to the long distance matrix element O 1 J/ψ , which itself is determined by the value of the (renormalized) wave function at the origin. At next order, one finds a contribution proportional to the long distance matrix element q 2 J/ψ which is suppressed by a relative v 2 . This matrix element is sensitive to the derivative of the wave function at the origin.
In this work we follow Ref. [90], where these matrix elements are determined. There, a subset of higher order (in v) contributions to the decay width including higher powers of ∇ 2 are resummed to all orders following Ref. [91]. As a result, the J/ψ decay width in the leptonic channel can be written as Here, e q = 2/3 is the fractional charge of the charm quark and M V is the J/ψ mass. At this order in v, the J/ψ is a pure S wave state, and its wave function can be factorized into a spin part and a scalar part. We will discuss the spin and angular momentum structure in more detail later. The extraction of the matrix elements that we use [90] has been done in a calculation that includes both velocity and α s corrections, such as in (22). Here, on the other hand, we will be using the light cone wave functions in a leading order calculation of cross sections, including only velocity corrections to the wavefunction. In a strict NRQCD power counting sense in α s , the α s corrections could be considered more important. Although steps have been taken to take them into account in the dipole picture exclusive cross section calculations [92] (see also recent work in a different formalism [93]), fully including them in the cross section is not yet possible at this point since the full photon to heavy quark pair wave function is not known to one loop accuracy. Thus we will leave a computation that includes also the perturbative α s calculations to future work, and continue with our focus on the velocity corrections to the wave function here.
Since our cross section calculation does not include pure α s corrections, taking the wave function to be given by just the operator O 1 J/ψ in (22) would lead to an inconsistent treatment of the α s corrections between the decay width and the cross section. Even in a more general sense, the α s contributions that appear as corrections to the decay widths or cross sections expressed in terms of nonrelativistic wavefunctions should, in light cone perturbation theory, be thought of as perturbative corrections to the light cone wave function itself [39,92]. This can be understood in the sense that the degrees of freedom in the nonrelativistic wavefunction are constituent quarks as opposed to bare quarks in the light cone wave function, see the discussion in [39]. To obtain a consistent picture here, we will absorb the α s correction to the scalar part of the wavefunction φ(r), which is then transformed to the light cone wave function. We thus relate the value and derivative at the origin of φ(r) to the long distance matrix elements as The non-perturbative long distance matrix elements have been determined in Ref. [90] by considering simultaneously the J/ψ → e + e − and η c → 2γ decays. As a result of this analysis, the matrix elements for J/ψ read The analysis in Ref. [90] is done by using the charm quark mass m c,N R = 1.4 GeV. In general, the charm quark mass in NRQCD can differ from the charm quark mass used in the IPsat fits discussed in Sec. II. In our numerical analysis, we will use the NRQCD value for the charm quark mass in both the meson and photon wave functions when using the NRQCD results. Everywhere else in this work we use the charm mass m c = 1.3528 GeV obtained in the IPsat fit to the HERA structure function data. The uncertainties quoted above for the long distance matrix elements are not independent, and the correlation matrix is also provided in Ref. [90]. To implement these correlated uncertainties, we use a Monte Carlo method and sample parameter values from the Gaussian distribution taking into account the full covariance matrix. The uncertainty is then obtained by calculating the one standard deviation band with respect to the result obtained by using the best fit values.
To construct the meson wave function, we start from the meson rest frame where we can use the NRQCD matrix elements to constrain the wave function as discussed above. In the rest frame, we require that the quark spins are coupled into a triplet state, and the total spin and angular momentum to a J = 1 vector state. Thus we can in general write the spin-structure of the wave function in the following form: Here Y m L L are the spherical harmonics, ψ L is the radial wave function corresponding to the orbital angular momentum L and j 1 m j1 j 2 m j2 |Jm J are Clebsch-Gordan coefficients. In general, the conservation of spin-parity tells us that for J/ψ the orbital angular momentum can only take values L = S, D. Since J/ψ should be dominated by the S wave contribution, we will from now on consider the case where only the S wave component is non-zero. We note that in principle in the NRQCD approach one finds the D wave contribution to the vector meson wave function to be suppressed by v 2 compared to the S wave, and this is of the same order as the first relativistic correction included in terms of the wave function derivative above. However, the D wave contribution to the decay width is suppressed by an additional v 2 and as such the D-wave contribution is not constrained by the decay widths at this order. Thus it is most consistent to set it to zero. In this case the wave function simplifies to: where φ(r) is the scalar part of the wave function and related to the long distance matrix elements as shown in Eqs. (24) and (25). Using the 3-dimensional polarization vectors in Eq. (11) we can also write this as where in the case of transverse polarization and when the vector meson is longitudinally polarized. Here ξ + = (1, 0) and ξ − = (0, 1) are the two component spinors describing spin-up and spin-down states andξs = iσ 2 ξ * s is the antiquark spinor. The behavior of the quarkonium wavefunction at long distances is determined by non perturbative physics. This long distance physics affects short distances through the requirement of the normalization of the wave function. The NRQCD approach broadly speaking consists of parametrizing the nonperturbative long distance physics by measurable coefficients that serve as coefficients in the short distance expansion, which is used to calculate a physical process happening at short distance scales. In practice this amounts to expressing the wave function as a Taylor expansion around the origin: The linear term does not appear to ensure that the Laplacian of the wave function is finite at the origin. The coefficients can also be written as A = φ(0) and B = 1 6 ∇ 2 φ(0), and using equations Eqs. (24) and (25) we get the values The uncertainties in the long distance matrix elements are correlated as discussed above, and in our numerical calculations this correlated uncertainty is propagated to the coefficients A and B.
We then want to write our wave function ansatz (33) in light cone coordinates (k, z). We do this by first going to momentum space: where k = (k, k 3 ). We then want to change the longitudinal momentum variable from k 3 to the plus momentum fraction carried by the quark: z. Unfortunately there is no unique way to do this, due to the different nature of instant form and light cone quantization. In principle we would want to define z as the ratio of the quark k + to the meson P + = M V / √ 2, working in the rest frame of the meson. However, a quark inside a bound state described as a superposition of different k modes is not exactly on shell, its energy being affected by the binding potential. Thus we do not precisely know the k 0 required to calculate k + from k 3 . The rest frame wave function also includes values of k 3 that are very large, leading to values of k + that are larger than M V / √ 2. This is perfectly possible in instant form quantization with the time variable t. However, in light cone quantization k + is a conserved momentum variable, and has to satisfy 0 < k + < P + . The procedure that we adopt here is (similarly to e.g. [94]) to define the momentum fraction in practice as z = k + q /(k + q + k + q ) where k q and kq are the quark and antiquark momenta, with k + calculated as- In other words, we normalize by the total plus momentum of the quark-antiquark pair, instead of the meson plus momentum, and assume an on-shell dispersion relation. This choice has the advantage that it leads to 0 < z < 1 by construction. This leads us to the expression for the longitudinal momentum in the meson rest frame k 3 as where is the invariant mass of the quark-antiquark pair. We emphasize that since this choice is not unique, we might expect corrections or ambiguities proportional to powers of the difference between the meson mass and the quarkantiquark pair invariant mass M 2 V − M 2 to appear. Such corrections are, however, higher order corrections in the nonrelativistic limit and also numerically very small for J/ψ for the values of m c,N R and q 2 used here. We could also hope that since the invariant mass is a rotationally invariant quantity, these ambiguities would not lead to serious violations of rotational invariance (which expresses itself here as the equality of physical properties such as decay widths of transverse and longitudinal polarization states). We will see an example of such a correction explicitly in Appendix B.
To change the variables in our wave function, one needs to be careful with the delta functions and their derivatives. We therefore make the change by requiring that the overlap where ϕ is an arbitrary wave function, does not change under the change of variables. This requirement tells us that the scalar part φ( k) changes to Equation (40) is the scalar part of the NRQCD based vector meson wave function in the meson rest frame, expressed in momentum space. We note that this wave function is not normalizable due to the presence of the delta functions. However, as the NRQCD approach can only be used to constrain the coordinate space wave function and its derivative at the origin, we are forced to use the expansion of Eq. (33) which can not result in a normalizable wave function. However, for the purposes of this work this is not a problem, as the vector meson production is sensitive to the vector meson wave function overlap with the virtual photon wave function, and the photon wave function is heavily suppressed at long distances where the expansion (33) is not reliable.

B. Wave function on the light front
The NRQCD wave function obtained in the previous section is written in the vector meson rest frame in terms of the quark and antiquark spin states s,s. In order to calculate overlaps with the virtual photon wave function (17) and (18), we need to express it in terms of the light cone helicities h,h. The transformation between these two bases, usually expressed in terms of the 2-spinors, is known as the "Melosh rotation" [46,47].
The Dirac spinors that are used to factorize the nonrelativistic wavefunction into a a spin-and scalar part, are eigenstates of the spin-z operator in the zero transverse momentum limit. In terms of the two component spin vectors ξ defined above in Eqs. (31) and (32) they read The normalization factor N is determined from the conditionū s us = −v s vs = 2mδ s,s . Both the Dirac spinors in terms of the spin-z component u s and the helicity spinors u h (see Eqs. (7), (8)) are solutions to the Dirac equation, and as such can be obtained as linear combinations of each other. This mapping is the Melosh rotation R sh . It can be computed from the spinor inner products (see also Ref. [95]) as where k + = zq + and q + is the meson plus momentum, and s and h refer to the spin and light front helicity, respectively. The helicity spinors u h and v h can also be written in a similar form as the spinors in the spin basis, Eqs. (42) and (43), by introducing the two component helicity spinors χ h . To do this we write the helicity spinors (7), (8) in the form where N is again determined by the normalization requirement andχ h = iσ 2 χ * h . Using this form one can check that the Melosh rotation also connects the two component spin and helicity spinors as The coefficients R sh can also be expressed as a 2×2 matrix rotating the 2-spinors Here M is the invariant mass of the qq system from Eq. (38) and n = (0, 0, 1) is the unit vector in the longitudinal direction. In terms of this matrix the 2-spinors ξ s and χ h are related by Using Eq. (47) we can now express the NRQCD wave function in the light front helicity basis. We write where the scalar part is given in Eq. (40). The helicity structure U λ hh is obtained by applying the transform (47) in Eq. (31) and (32), i.e., After the Melosh rotation, we compute the Fourier transform to obtain the light front wave function in the mixed transverse coordinate -longitudinal momentum fraction space as The different helicity components of the final light front wave function resulting from this procedure are The first relativistic correction to the wave function, proportional to B or the wave function derivative, mixes the helicity and spin states. In particular, in the case of transverse polarization the h,h = ±∓ terms are nonvanishing when the relativistic correction is included. In general, we expect that if higher order corrections in v were included in the wave function parametrization, we would also find other components to be non-vanishing.
The Melosh rotation is crucial here, as it generates helicity structures that are not visible in the spin basis. This is in contrast to some early attempts to transform the wave functions obtained by solving the potential models to the light front as done e.g. in Ref. [78]. The role of the Melosh rotation in the context of vector meson light front wave functions and exclusive scattering was first emphasized in Ref. [47]. More recently it was applied to J/ψ production in the dipole picture in Ref. [96], and in [97] different quark-antiquark potentials were studied in this context. In the case of excited states such as ψ(2S) the role of the Melosh rotation is expected to be even more significant [98].
Let us in passing briefly compare our approach to the one in the recent work of Krelina et al in Ref. [96]. In our approach, we take the NRQCD wave function which only includes the S wave contribution (D wave part is suppressed by v 2 ). The quark spin dependence is now trivial, as the total angular momentum must be provided by the quark spins which gives us the structure of Eq. (30).
In Ref. [96], the authors assume, unlike we do here, that the spin structure of the vector meson wave function in the rest frame has the same form as the light cone helicity structure of the photon light cone wave function, Eq. (17) and (18). This structure is then supplemented by a wave function obtained from the potential model, and a Melosh rotation to the light front is applied at the end. Such a procedure leads to a large D-wave contribution in the wavefunction, which we do not have. We discuss the structure of the wavefunctions in terms of Sand D-waves in more detail in Appendix A.
To determine the role of the relativistic corrections in the vector meson wave function, we will also study for comparison the fully non-relativistic wave function where our starting point for the scalar part is Following the previous procedure, the final result for the light cone wave function can be read from Eq. (53) with the substitutions A = A , B = 0. One notices that this can now be written as In this extreme non-relativistic limit (k = 0, z = 1/2) the Melosh rotation simply corresponds to an identity matrix so that the spin and helicity bases are interchangeable here. The normalization A is obtained from the van Royen-Weisskopf equation for the leptonic width [99], which is also obtained from Eq. (22) by neglecting the relativistic correction proportional to q 2 J/ψ /m 2 c,N R , and the higher order QCD correction ∼ α s (note that parametrically α s ∼ v): By using the experimental value for leptonic width [100], we can calculate the coefficient A to be C. Overlap with photon Using the obtained J/ψ wave function on the light front, Eq. (53), we can directly compute overlaps with the virtual photon, Eqs. (17) and (18). In these overlaps, we also include the phase factor exp i z − 1 2 r · ∆ present in the vector meson production amplitude in Eq. (1). We also assume that the dipole amplitude does not depend on the orientation θ r of r as is the case in the IPsat parametrization, and integrate over θ r . The overlaps summed over the quark helicities read and where¯ 2 = Q 2 /4 + m 2 c,N R and ∆ = |∆| and r = |r|. In the case of transverse polarization, the result is identical in cases with λ = +1 and λ = −1. We will study these overlaps numerically in Sec. V A. We note that thanks to the delta function structure in z in our wave function (53), many phenomenological applications become numerically more straightforward as the z integral can be performed analytically.

IV. PHENOMENOLOGICAL WAVE FUNCTIONS
To provide a quantitative point of comparison for the effect of the relativistic corrections, we want to compare the light cone wave functions obtained in Sec. III to other parametrizations used in the literature. For this purpose, let us now discuss two specific alternative approaches used for phenomenological applications in the literature.

A. Boosted Gaussian
A commonly used phenomenological approach to construct the vector meson wave function is to assume that it has the same polarization and helicity structure as the virtual photon. This can be done by replacing the scalar part of the photon wave functions (17) and (18) by an unknown function as [48] with the explicit factor Q in the longitudinal wave function replaced by the meson mass as 2Q → M V . The scalar function φ(r, z) is then parametrized, and the parameters can be determined by requiring that the resulting wave function is normalized to unity and reproduces the experimental leptonic decay width. As we will discuss in more detail in Appendix B, this procedure does not correspond to the most general possible helicity structure. Nevertheless, our result at this order in the nonrelativistic expansion can in fact also be written in terms of the "scalar part of light cone wave functions." However, at higher orders in v different a different structure could appear.
In the Boosted Gaussian parametrization, the qq invariant mass distribution is assumed to be Gaussian, with the width of the distribution R and the normalization factors N T,L being free parameters. In mixed space, the parametrization reads In this work we use the parameters constrained in Ref. [87] by using the same charm quark mass m c = 1.3528 GeV as is used when fitting the IPsat dipole amplitude to the HERA data. The parameters are determined by requiring that the longitudinal polarization can be used to reprodcue the experimental decay width. The obtained parameters are R = 1.507 GeV −1 , N T = 0.589 and N L = 0.586 with M V = 3.097 GeV. The specific functional form and helicity structure of the Boosted Gaussian parametrization imply that in the vector meson rest frame there are both S and D wave contributions. This is demonstrated explicitly in Appendix A by performing a Melosh rotation from the light front back to the J/ψ rest frame. This feature is hard to describe in potential model calculations, and our NRQCD based wave function in particular has only the S wave component in the rest frame. The D-wave contribution in the Boosted Gaussian wavefunction is, however, quite small.

B. Basis Light-Front Quantization (BLFQ)
The second wave function we study here for comparisons is based on explicit calculations on the light front. In this approach, one constructs a light front Hamiltonian H eff , which consists of a one gluon exchange interaction, and a non-perturbative confining potential inspired by light-front holography. The formalism is developed in Refs. [101][102][103][104][105][106][107].
The quarkonium states are obtained by solving the eigenvalue problem As a solution, one obtains the invariant mass M 2 V spectrum and the light front wave functions in momentum space Here J, P, C and m J are the total angular momentum, parity, C-parity and the magnetic quantum number of the state, respectively. The free paramters, value of the coupling constant, strength of the confining potential, quark mass and the effective gluon mass, can be constrained by the charmonium and bottonium mass spectra [108,109]. In this work, we use the most up-to-date parametrizations from Ref. [109].
The obtained BLFQ wave functions have been applied in studies of the J/ψ production in the dipole picture at HERA [110], and in the context of exclusive J/ψ production in ultra peripheral heavy ion collisions at the LHC in Ref. [111]. Following the prescription used in Refs. [110,111], we consider the fitted quark mass m BLFQ c = 1.603 GeV in the "BLFQ wave function" to be an effective mass of the quarks in the confining potential, including some non-perturbative contributions. Consequently, when calculating the overlaps we use, as in [110,111], m c = 1.3528 GeV for the charm mass in the photon wave function, as constrained by the charm structure function data in the IPsat fit [87].

A. Photon overlap
The exclusive vector meson production cross section depends on the overlap between the cc component of the virtual photon wave function with the vector meson wave function, see Eq. (1). In Fig. 1 these overlaps for ∆ = 0 are shown as a function of the transverse size r = |r| of the intermediate dipole, using four vector meson wave functions: 1. NRQCD expansion, which is constructed by parametrizing the wave function and its derivative at the origin based on NRQCD matrix elements including corrections ∼ v 2 , and performing the Melosh rotation to the light front. This is our result from Sec. III.

2.
Delta, which is the fully non-relativistic limit (Eq. (55)) of the above wave function, without any information about the wave function derivative.

Boosted
Gaussian, the phenomenological parametrization discussed in Sec. IV A

BLFQ wave function based on Basis Light-Front
Quantization, discussed in Sec. IV B.
In Fig. 2 we show the same overlaps plotted as ratios to the fully nonrelativistic limit, i.e. the Delta parametrization. For the NRQCD expansion based wave function, we also show the model uncertainty related to the NRQCD matrix elements that control the value of the wave function and its derivative at the origin. The uncertainty band is in this case computed as discussed in Sec. III. The effect of the first relativistic correction can be determined by comparing the Delta and NRQCD expansion wave functions. At large dipoles the negative velocity suppressed ∼ r 2 contribution suppresses the vector meson wave function 1 compared to the fully non-relativistic form. This is especially visible at small Q 2 . At larger photon virtualities, the exponential suppression in the photon wave function becomes dominant before the relativistic −r 2 correction becomes numerically important. Thus, while the effect of the relativistic correction is dramatic in the ratio in Fig. 2, at large Q 2 it is insignificant for the actual overlap, as is seen in Fig. 1.
For small dipoles the wavefunctions are most strictly constrained by the quarkonium decay widths. The NRQCD parametrization does not, however, reduce exactly to the fully nonrelativistic Delta parametrization in the small r limit. This can be traced back to the fact that the gradient correction also affects the decay width, 1 The wave function would change sign at r 0 = 0.73 fm. As there should be no node in the J/ψ wave function, we set the wave function to zero at r > r 0 . We have checked that this cutoff has a negligible effect on our numerical results.
as seen in Eq. (22) (and from the fact that the constants A in (34) and A in (57) are different). A part of the 3-dimensional gradient correction becomes a correction to the functional form in z even at r = 0. This leads to the overlaps at small r being slightly different, even though the same decay width data is used to obtain the parameters of the rest frame wavefunctions.
Both the Boosted Gaussian and BLFQ wave functions are even more suppressed at large dipole sizes than the NRQCD parametrization. This is most clearly seen on the ratio plot, Fig. 2. This is a straightforward consequence of the fact that in these parametrizations the wavefunction normalization imposes an additional suppression at large r. For the Boosted Gaussian parametrization this additional suppression happens at such a large r that the overlap is already very small, and thus has a negligible effect on the overall overlap in Fig. 1. The Boosted Gaussian parametrization is very close to our NRQCD also for small dipoles. The BLFQ parametrization yields a a somewhat larger wave function overlap at small r than our NRQCD one, or the Boosted

Gaussian. 2
The suppression with respect to the nonrelativistic limit is larger for the longitudinal polarization state than for the transverse one. This can be understood as follows. The longitudinal virtual photon wave function depends on the quark momentum fraction as ∼ z(1 − z) (see Eq. (17)), and as such is peaked at z = 1/2. The z-structure of the fully non-relativistic wave function is δ(z − 1/2), and when the first relativistic corrections are included, the z = 1/2 region still dominates the overlap. 2 The difference between the BLFQ and Boosted Gaussian at small r is somewhat surprising given that the Boosted Gaussian and BLFQ parametrizations use the same decay width data to constrain the wavefunction. We do not know what is the main cause of this discrepancy. One possibility is that the different functional forms lead to a different relation between the decay constant on one hand and the photon wavefunction overlap on the other hand. Another possibility is that the difference is related to using different components of the electromagnetic 4-current J µ expectation value to calculate the decay constant from the wave function.
On the other hand, the transverse photon wave function is not peaked at z = 1/2, see Eq. (18). Thus, the suppression from the ∂ 2 z δ(z − 1/2) term in the relativistic correction is smaller for the transverse polarization.

B. J/ψ production
The total exclusive DIS J/ψ production cross section for a proton target at W = 90GeV is shown in Fig. 3, compared with the H1 [4] and ZEUS [8] data. The overall normalization of the cross section has a relatively large theoretical uncertainty. We note that the two corrections discussed in Sec. II, the real part and especially the skewedness correction are numerically significant, up to ∼ 50% (see e.g. Ref. [40]). As discussed in Sec. II, especially the skewedness correction is not very robust and its applicability in the dipole picture used here is not clear. In addition to the possibly problematic skewedness corrections, the fact that our NRQCD based wave functions are not normalized affects the absolute normalization of the vector meson production cross sections. Thus our fo- cus here is rather on the relative effects of different meson wave functions and the dependence on Q 2 . The vector meson cross section is dominated by dipole sizes of the order of 1/(Q 2 + M 2 V ) as can be seen 3 from Fig. 1. Consequently, it is more instructive to look at the dependence of the J/ψ cross section on Q 2 than the overall normalization. From Fig. 3 one sees that the fully nonrelativistic wave function results in a too steep Q 2 dependence compared to the HERA data. The first relativistic correction slows down the Q 2 evolution close to the photoproduction region and leads to a better agreement with the experimental data. This is a consequence of the basic behavior of the relativistic correction as a ∼ −r 2 modification that suppresses the vector meson wave function strongly at large dipoles. Thus the reduction from the relativistic correction is larger for smaller Q 2 . At large Q 2 the exponential suppression from the photon wave function starts to dominate at smaller dipole sizes, and the relativistic −r 2 correction becomes negligible. However, the relativistic contribution to the momentum fraction z structure is present at all Q 2 , and suppresses the longitudinal cross section more than the transverse one.
A similar trend in the Q 2 dependence is also visible with both the Boosted Gaussian and BLFQ wave functions. For the Boosted Gaussian case, the agreement with HERA data has been established numerous times in the previous literature, e.g. in Ref. [48]. The Q 2 dependence of the cross section is slightly weaker when the BLFQ wave function is used, but the difference is comparable to the experimental uncertainties. We note that in Ref. [110] the BLFQ wave function is found to result in a cross section underestimating the HERA data in the photoproduction region. In this work, compared to the setup used in Ref. [110], we use an updated BLFQ parametrization from Ref. [109] which was shown in Ref. [111] to result in a good description of the J/ψ production in ultra peripheral proton-proton collisions at the LHC, which in practice probe vector meson photoproduction [15,16].
To cancel normalization uncertainties, we next study cross section ratios. In Fig. 4 the longitudinal-totransverse ratio of the J/ψ production cross section is shown as a function of the photon virtuality. The results are compared with the H1 and ZEUS data from Refs. [4,8]. The first relativistic correction reduces the longitudinal cross section more than the transverse one. As discussed above, this is due to the fact that a part of the correction shifts the meson wave function away from the δ(z − 1/2), which is the structure preferred by longitudinal photons but not by transverse photons. This shows up as a decrease in the longitudinal to transverse ratio as a function of Q 2 . The effect is even stronger with the Boosted Gaussian and BLFQ wavefunctions.
Finally, we study vector meson production in the future Electron Ion Collider. As the diffractive cross section at leading order in perturbative QCD is approximatively proportional to the squared gluon density, exclusive vector meson production is a promising observable to look for saturation effects at the future Electron Ion Collider (see e.g. [112]).
The denominator corresponds to the so called impulse approximation, which is used to transform the photonproton cross section to the photon-nucleus case in the absence of nuclear effects, but taking into account the different form factors (transverse density profiles Fourier transformed to the momentum space). The A 4/3 scaling can be understood to originate from the fact that the coherent cross section at t = 0 scales as ∼ A 2 , and the width of the coherent spectra (location of the first diffractive minimum) is proportional to 1/R 2 A ∼ A −2/3 . The numerical factor c depends on the proton and nuclear form factors, and is found to be very close to c = 1 2 in Ref. [87]. In the absence of non-linear effects (or shadowing effects in the gluon distribution), with dipole amplitudes (19) and (21) that depend linearly on r 2 xg(x, µ 2 ), this ratio is exactly 1.
The obtained nuclear suppression factor is shown in Fig. 5 in the Q 2 range accessible at the Electron Ion Collider. We emphasize that all the nuclear modifications in this figure are calculated with exactly the same dipole cross sections, corresponding to the same nuclear shadowing (as measured e.g. by the nuclear suppression in F L or F 2 ). Thus the difference between the curves results purely from vector meson wave function effects. When using the NRQCD wave function with the relativistic corretion, the Boosted Gaussian wave function or the BLFQ wave function, the obtained nuclear suppression factors are practically identical. Even though large mass of the vector meson renders the scale in the process large, a moderate suppression ∼ 0.75 is found at small and moderate Q 2 . In the small Q 2 region the uncertainty obtained by varying the NRQCD matrix elements is large.
The fully non relativistic wave function results in a much stronger suppression at small Q 2 . This can be understood, as it was already seen in Fig. 1 that this wave function gives more weigh on larger dipoles compared to the other studied wave functions. As the larger dipoles are more sensitive to non-linear effects, a larger nuclear suppression in this case is anticipated. The first relativistic correction ∼ −r 2 suppresses the overlap at large dipole sizes, and consequently the nuclear suppression. At higher Q 2 the photon wave function again cuts out the large dipole part of the overlap in all cases, and as such the results obtained by applying the fully non-relativistic wave function do not differ from other wave functions any more. At asymptotically large Q 2 only small dipoles contribute and the dipole amplitudes can be linearized. Consequently, the suppression factor approaches unity at large Q 2 independently of the applied wave function.
The fact that the fully non-relativistic wave function results in a very different nuclear suppression demonstrates that the dependence on the meson wave function does not completely cancel in the nucleus-to-proton cross section ratios. Consequently, a realistic (and relativistic) description of the vector meson wave function is necessary for interpreting the measured nuclear suppression factors. This indicates that there is a large theoretical uncertaintly in using the fully nonrelativistic formula of Ryskin [3], not only for extracting absolute gluon distributions, but even for extracting nuclear modifications to the dipole cross section (or the gluon density) from cross section ratios.

VI. CONCLUSIONS
In this work we proposed a new parametrization for the heavy vector meson wave function based on NRQCD long-distance matrix elements. These matrix elements can be used to simultaneously constrain both the value and the derivative of the vector meson wave function at the origin using quarkonium decay data. This approach provides a systematic method to compute the vector meson wave function as an expansion in the strong coupling constant α s and the quark velocity v.
Compared to many phenomenological approaches used in the literature, our approach uses two independent constraints (the wave function value and its derivative). The obtained wave function is rotationally symmetric in the rest frame and contains only the S wave component. Consequently, we simultaneously obtain a consistent parametrization for both polarization states. This is unlike in some widely used phenomenological parametrizations where the virtual photon like helicity structure is assumed on the light front. Relating light cone wavefunctions to rest frame ones also provides a consistent way to discuss the effect of a potential D-wave contribution to the meson wavefunction. We do not see indications, neither theoretically nor phenomenologically, that a significant D-wave contribution would be required or favored for the J/ψ .
The first relativistic correction to the wave function, controlled by the wave function derivative at the origin, is found to have a sizeable effect on the cross section.
The negative ∼ −r 2 relativistic contribution in terms of the transverse size r suppresses the obtained wave function at larger dipole sizes. The momentum fraction part of the correction partially compensates for this effect for the transverse photon by shifting the wave function away from the fully non relativistic configuration where both quarks carry the same fraction of the longitudinal momentum, a configuration which is not preferred by the transverse photon.
A disadvantage in our approach is that it is not possible to obtain a wave function which is normalized to unity. In the NRQCD framework the value of the wavefunction at long distances is parametrized by a nonperturbative matrix element, whose effect is felt in the value of the wavefunction near the origin. This can lead to an overestimation of the cross section at Q 2 = 0, where one is most sensitive to the long distance behavior of the wave function. In practice, however, we obtain cross sections that are quite similar to what is given by e.g. the Boosted Gaussian parametrization. The wave function overlap with the photon is also smaller than with the BLFQ approach. Thus the lack of normalization in the wavefunction does not seem to be an important effect for J/ψ . The situation would be different for lighter vector mesons.
The structure of the wave function can be probed by studying cross sections (and cross section ratios) at different photon virtualities where the dipole sizes contributing to the cross section vary. The first relativistic correction is found to weaken the Q 2 dependence of the total J/ψ production cross section and the longitudinal-totransverse ratio. These effects are broadly similar to predictions obtained by the Boosted Gaussian parametrization, or by the BLFQ wave function that is based on an explicit calculation on the light front including confinement effects.
When comparing vector meson production off protons to heavy nuclei, we find that the wave function does not completely cancel in the nuclear suppression factor, which compares the γ * A cross section to the γ * p in the impulse approximation. This demonstrates that a realistic vector meson wave function is necessary to properly interpret the nuclear suppression results, and in particular a fully non-relativistic approach can not be reliably used to extract the non-linear effects on the nuclear structure.
In addition to the corrections in velocity, it would be important to include perturbative corrections in the strong coupling α s in the calculation of exclusive vector meson production. Indeed some recent advances [92,113] are gradually making it possible to do so in the dipole picture. However, a study of the phenomenological implications of these α s corrections remains to be done. In terms of understanding current and future experimental collider data, it would also be important to explore whether this approach can be extended to excited states such as the ψ(nS). wave function described in Sec. IV A. As shown in [114], in the rest frame, the squared J/ψ BLFQ wave function is dominated by the S wave component, the D wave contributing only a small fraction of the order of 0.1 . . . 4% (depending on the polarization). In heavier mesons, this contribution is even smaller. This is comparable to the Boosted Gaussian case discussed above.
Overall, based on neither the Boosted Gaussian nor the BLFQ parametrizations, we do not see any confirmation for the result of Ref. [96], where the D wave part of the J/ψ wave function was found to result in tens of percent contribution on the vector meson production cross section. Part of this discrepancy might be merely a question of terminology. In our discussion here, we have insisted that the terms S-wave and D-wave refer to the angular momentum components of the 3-dimensional wave function in the meson rest frame. Thus the mere presence, in the light cone wave function, of terms proportional to transverse momenta originating from the Melosh rotation cannot be taken as an indication of a D-wave component in the meson.
Let us now move to the case of a virtual photon. Since a spacelike virtual photon does not have a rest frame and is not a bound state, it is not customarily thought of in terms of an S -D wave decomposition. Now, however, we have an explicit light cone wave function for the photon just like for the meson, and we can use the same procedure to determine its S-and D-wave components in the meson rest frame. The resulting squared light front wave functions summed over quark helicities are shown in Fig. 7. The full photon wave function, written in Eqs. (17) and (18), is denoted by S + D, as it can be written as a sum of these two components. When compared to the full result, the squared D wave only contribution is found to be strongly suppressed. There is also a contribution originating from the overlap between the S and D wave contributions. This term would vanish if we integrated over all the angles. Here, we only integrate over the azimuthal direction of r. Integration over the momentum fraction z corresponds to the evaluation of the coordinate space wave function at x 3 = 0, and consequently one angular integral is not performed and the overlap does not vanish. The relative importance of different contributions is found to be approximatively independent of Q 2 . The S − D overlap contribution is numerically significant, which is reflected by the large difference between the full result and the S wave only contribution. This suggests that even though the D wave contribution is suppressed by the quark velocity, the charm quark mass is not large enough to render this contribution negligible. This is due to the fact that the photon wave function in the momentum space behaves as ∼ 1/k 2 where k is the quark transverse momentum, and this powerlike tail brings numerically large contributions from relatively large momenta. Additionally, the integration over the longitudinal momentum fraction z includes high momentum contributions, as the photon wave function has a support over a large range of z.
The contribution from the S − D overlap changes sign at large transverse separations in case of the longitudinal polarization. There is no node in the radial part of the wave function, but the spherical harmonic function describing the angular part of the D wave component changes sign, which explains the sign flip. In the S − D overlap mostly the helicity-+− and −+ components of the D wave contribute by coupling to the S wave. On the other hand, in the D wave squared wave function, one sums all helicity components. As the D wave component itself is a relativistic correction, none of the helicity structures dominates unlike in the S wave part. Moreover, only the +− and −+ helicity components change sign at large distances, and as the ++ and −− components do not vanish in this region, no node appears in the squared D wave result. In the case of transverse polarization with λ = ±1, the helicity component ±± in the D wave also changes the sign at large distances, but this effect is not easily visible in Fig. 7 as the other helicity components that do not change sign dominate.
Appendix B: Photon-like parametrizations of light cone wave functions As discussed in Sec. IV, an often used approach to parametrize vector meson wave functions is to start from the helicity structure of the virtual photon light cone wave functions (17), (18). One then replaces the Bessel function K 0 in the photon wave functions (17) and (18) by an unknown function as [48] with the explicit factor Q in the longitudinal wave function replaced by the meson mass as 2Q → M V . This leads, with our sign conventions, to the wave function being written as The scalar functions φ T,L (r, z) are then parametrized, and the parameters can be determined by requiring that the resulting wave function is normalized to unity and reproduces the experimental leptonic decay width. In terms of Lorentz-invariant form factors this means that the meson is assumed to have a nonzero Dirac form factor but a vanishing Pauli form factor, since this is the structure dictated by the gauge-boson-fermion vertex at leading order perturbation theory. The procedure therefore does not generate the most general possible helicity structure.
This photon-like parametrization approach starts from a spacelike photon, where the photon momentum breaks rotational symmetry that is manifested here as a symmetry between longitudinal and transverse meson polarization states. The common approach is to separately parametrize the longitudinal and transverse functions φ T,L (r, z). The helicity structure obtained by generalization from the photon wave function is of course consistent with rotational symmetry, since the decay of a timelike virtual photon is rotationally symmetric. Thus one could derive a constraint relating φ L (r, z) and φ T (r, z) by requiring the meson rest frame wavefunctions to be the same. To our knowledge this approach has not, however, been used in the literature. Using separate parametrizations for φ T,L (r, z) should be contrasted with the approach in this paper. Here, we maintain rotational invariance in the meson rest frame, in particular starting from the same decay constants calculated from the rest frame wave functions. Our procedure for going from the rest frame to the light cone wave function therefore simultaneously determines the wavefunction for both longitudinal and transverse polarization states.
One can take the parametrization (B2), (B3) in momentum space, perform the inverse Melosh rotation and separate the S-and D-wave components to get a rest frame 3-dimensional wavefunction. Assuming that the Fourier transforms of the scalar functions are rotationally invariant, i.e. φ T,L (k, z) depend only on k = | k| = (M/2) 2 − m 2 c , the result of this exercise in momentum space is: × N c π 2 (E + m c )(2E) 3/2 These expressions are written in terms of the energy of the quark in the meson rest frame E = k 2 + m 2 c = M/2, and k = | k|. Let us point out a few aspects of these expressions. Firstly, as discussed above, in the photonlike parametrization there is always a D-wave component in the meson wave function. It is, as expected, explicitly a relativistic correction, i.e. proportional to the squared 3-momentum of the quark. Secondly, the rest frame wave functions of the transverse and longitudinal polarizations are not the same, but differ by a factor where M is the invariant mass of the quark pair and M V the mass of the meson. As discussed earlier in Sec. III, the coordinate transformation from k 3 to z inevitably introduces ambiguities that are proportional to this difference, so this should not come as a surprise. In an NRQCD power counting, this difference is of the order of the binding energy of the meson, which is higher order than we are considering here.
In spite of this discussion, the wavefunction that we obtained in Sec. IV can in fact be written in a photon-like form in terms of scalar parts of light cone wave functions. In the notation of [48] these read We emphasize that we do not expect that writing down such a paramterization in terms of two scalar parts of a light cone wave function having the helicity structure of a photon would be possible at higher orders in the nonrelativistic expansion. As a side remark, we discussed above that the photonlike structure generically implies a non-zero D wave com-ponent, see Eqs. (B6) and (B7). On the other hand, our NRQCD based wave function by construction has no D component. However, the D-wave component resulting from inserting the scalar parts (B8) and (B9) into the formulae for the D-wave contribution, Eqs. (B6) and (B7), behaves as ∼ k 2 ∇ 2 k δ (3) ( k). Such a function actually yields zero when convoluted with any test function f ( k), since the angular integral picks out the = 2 component of f , which must vanish at k = 0. Thus the D-wave contribution corresponding to (B8) and (B9) is in fact zero in a distribution sense.
It is interesting to note, that when one calculates from these expressions the decay constants for the different polarization states using the light cone perturbation theory expressions (26) and (27) in Ref. [48], one obtains: The results are not exactly equal. However, as discussed above, if one approximates the meson mass by the quark pair invariant mass, as we did in transforming to the momentum fraction z, they do reduce to the same result. This can be seen explicitly by replacing M V in (B11) by M ≈ 2 m 2 c,N R + q 2 and using Eq. (35) to write, at lowest nontrivial order in the quark velocity, 2m c,N R /M V ≈ 1 + 3B/(m 2 c,N R A) (note that B < 0). In this approximation Eqs. (B10) and (B11) also give back the same decay width expression that we are using to determine the rest frame wavefunction. We reiterate that a difference such as this can be expected in our procedure. We are constructing our wavefunctions by requiring the decay widths calculated from the rest frame wave functions to have the correct value, and to be the same for the different polarization states. The coordinate transformation to light cone wave functions does not conserve these properties exactly, but only up to a given order in the nonrelativistic expansion.