QCD analysis of near-threshold photon-proton production of heavy quarkonium

The near-threshold photo or electroproduction of heavy vector quarkonium off the proton is studied in quantum chromodynamics. Similar to the high-energy limit, the production amplitude can be factorized in terms of gluonic generalized parton distributions and the quarkonium distribution amplitude. At the threshold, the threshold kinematics has a large skewness parameter ξ , leading to the dominance of the spin-2 contribution over higher-spin twist-2 operators. Thus, threshold production data are useful to extract the gluonic gravitational form factors, allowing studying the gluonic contributions to the quantum anomalous energy, mass radius, spin, and mechanical pressure in the proton. We use the recent GlueX data on the J= ψ photoproduction to illustrate the potential physics impact from the high-precision data from the future Jefferson Laboratory 12 GeV and Electron-Ion Collider physics program. DOI:


I. INTRODUCTION
Recently, there has been rising interest in measuring the photoproduction of J=ψ particles near the threshold. Dedicated experiments with such a purpose are running at Jefferson Laboratory, and similar experiments are planned at the future Electron-Ion Collider [1]. It has been proposed in Refs. [2,3] that the photoproduction of J=ψ can be used to measure the gluon matrix element hPjF 2 jPi in the nucleon and provide crucial information about the trace anomaly contribution to the nucleon mass and mass radius [2,[4][5][6][7].
In the literature, different methods have been adopted to analyze the process. In the vector dominance model [2,3,7], the vector-meson photoproduction is related to the forward meson-nucleon scattering where a direct operator product expansion (OPE) in terms of gluonic matrix elements is applicable [8][9][10][11]. Recently, this process has also been approached using dispersive analysis [12,13] and holographic QCD [14][15][16][17]. However, a comprehensive understanding of the process in perturbative QCD in the threshold region is still lacking except a few early attempts [18,19]. This is sharply contrary to the large Q 2 and large W diffractive region where the process has attracted attention since the mid-1990s [20][21][22][23][24] with well-established allorder factorization [25] in terms of generalized parton distributions (GPDs) [26][27][28], and the Q 2 ¼ 0, large W, and small t region where the factorization has been explicitly shown at next to leading order [29]. In this paper, we study the near-threshold photoproduction of the heavy vector-meson with mass M V in the heavy-quark mass limit M V → ∞. In this limit, the process is dominated by the direct photon coupling with the heavy quarks, and the heavy-quark production through gluonic subprocesses including possible intrinsic heavy flavor is suppressed by α s ðM V Þ → 0. We show that the QCD factorization in terms of gluon GPDs in Ref. [29] remains valid in the threshold region as well. Different from the proposals in Refs. [2,3] and recent calculations in a different limit [30,31], the leading contribution comes from the tensor part of gluonic energy momentum tensor (EMT) and high-dimensional twist-2 gluonic operators, due to the emergent light-cone structure in the large M V limit. The near-threshold region is characterized by large skewness parameter ξ ∼ 1 regardless of the Q 2 [31], and the gluon EMT dominates over highdimensional operators as well as three-gluon exchanges [18,30]. Therefore, the process can be used to probe the gluonic gravitational form factors of the proton, which provide important information about the gluon contributions to the proton's mass and spin as well as pressure structures [4,5,[32][33][34].
The organization of the paper is as follows. In Sec. II, we introduce the near-threshold kinematics of the process, paying attention to the emergent light-cone structure in the initial state. In Sec. III, we perform the analysis of the two-gluon exchange diagram and express the amplitude in terms of gluon GPDs. In Sec. IV, we Taylor expand the amplitude in terms of moments of the GPD and show that the cross section can be expressed in terms of the gluonic gravitational form factors of the proton for the large skewness parameter ξ and the decay constant of the heavy meson. In Sec. V, assuming the dominance of the spin-2 matrix element, we compare our prediction to experiment data of J=ψ production using the results of the gravitational form factors extracted from lattice calculation in Ref. [35] and fit a parametrization of the form factors with the GlueX [36] data, especially the M A and Cð0Þ, which are important to the proton mass radius. With those fitted form factors, we predict the cross section of ϒ production near threshold. In Sec. VI, we study the polarization effects and make predictions for polarization-dependent cross sections which are important for disentangling various form factors. Finally, we make several comments and conclude in Sec. VII.

II. NEAR-THRESHOLD KINEMATICS
We first investigate the near-threshold kinematics of the process. Without loss of generality, we work in the c.m. frame as shown in Fig. 1, though the final result is frame independent. The four-momenta of the incoming photon, outgoing vector meson, incoming proton, and outgoing proton are denoted by q, K, P, and P 0 , respectively. We restrict ourselves to the real photon and near-threshold region, although similar analysis can be extended to a finite virtual photon mass Q 2 . We choose the incoming proton to move in the þz direction, as in the laboratory frame. In c.m. frame, the magnitudes of the three-momenta can be expressed in Lorentz scalars as where we label the nucleon mass M N , the c.m. energy squared W 2 ≡ ðP þ qÞ 2 ≥ ðM N þ M V Þ 2 , the average proton four-momentumP ≡ ðP 0 þ PÞ=2, four-momentum transfer Δ ¼ P 0 − P, and associated invariant t ≡ Δ 2 . In the heavy-quark limit, the vector-boson mass M V ≫ M N . As a result, W ≫ M N , and the incoming proton travels almost along the þ light-cone direction in terms of the light-front coordinate , where x μ is a four-coordinate. One has with p ¼P þ ffiffi 2 p ð1; 1; 0 ⊥ Þ and n ¼ 1 ffiffi 2 pP þ ð1; −1; 0 ⊥ Þ being two opposite light-cone unit vectors. The four-momenta of the final-state proton and vector meson are where ⃗ e can be in any spatial direction. In this paper, we mainly focus on the threshold region defined by the condition that the velocity of the final-state proton β ≡ j ⃗ Kj=P 00 is of order 1. This condition implies that the velocity of the heavy meson is of order Oð M N M V Þ and Right at the threshold, the invariant momentum transfer t equals which is of order M V M N in the heavy-quark limit. As W increases, the allowed region of −t forms a band ½jtj min ðWÞ; jtj max ðWÞ between the backward case with −t ¼ jtj min and forward case with −t ¼ jtj max . Near the threshold, both of them are of order OðM N M V Þ, much larger than M 2 N . In the standard notation of GPDs [27], this implies that the skewness is close to 1, which has also been observed for the large Q 2 region [31]. In Figs. 2 and 3, we show the ξ value on the ðW; −tÞ plane for J=ψ and ϒ production in the kinematically allowed region. This condition ξ → 1 near threshold will allow us to study the form factors of the gluon EMT as we discuss later.

III. LEADING-ORDER QCD FACTORIZATION
We now proceed to consider the amplitude of the process written as where J μ is the electromagnetic current operator, ε μ is the polarization vector of the photon normalized as ε 2 ¼ −1, and jKðε V Þi denotes the vector-meson state with momentum K and polarization vector ε V normalized as ε 2 V ¼ −1 and satisfying K · ε V ¼ 0. All states are normalized covariantly. In the heavy-quark limit, the leading-order contribution to M is given by two-gluon exchange diagrams as illustrated in Fig. 4. To calculate these contributions, we perform the following approximations which are justified in leading order in Oð M N M V Þ and Oðα S Þ. First, since the momentum transfer is mostly in the lightcone direction p, only the þ components of the loop momenta l and −Δ − l for the gluons are kept in the heavyquark loop and can be expressed as where −1 < x < 1 is the momentum fraction for l. In the threshold region, ξ approaches 1 in the heavy-quark limit, though we will keep the ξ dependence in general. Second, at the large M V limit, we perform the nonrelativistic approximation to the heavy-meson wave function, which amounts to the substitution Here, v V is the four-velocity of the vector meson, and the scalar functionφðkÞ corresponds to the distribution amplitude of the vector meson. The relative momentum between the quark and antiquark is of order Oðα S M V Þ for hydrogenlike systems, and to the leading order in α s , we can further approximate the wave function as where the velocity v V of the final-state heavy-meson can be approximated as the static one is proportional to the nonrelativistic hydrogenlike wave function ψ NR at the origin. The above meson wave function can be studied more systematically with nonrelativistic QCD, which shows that the leading correction to the wave function is indeed suppressed by order α S ðM V Þ [37][38][39].
Finally, we notice that for collinear gluons in light-cone gauge A þ ¼ 0, the A − component is suppressed by the large boost factor. Therefore, to the leading-order approximation, it is sufficient to keep the transverse components of the gluon gauge potential only.
Given these approximations, we can evaluate the Feynman diagrams in light-cone gauge straightforwardly. The result for the leading amplitude Mðε V ; εÞ reads Here, the function Gðt; ξÞ implicitly depends on the polarization of the initial and final proton which is suppressed and can be expressed in terms of the gluon GPDs as where the hard kernel Aðx; ξÞ reads The standard gluon GPD F g is defined as [40] F g ðx; ξ; tÞ where the states are normalized as hPjPi ¼ 2E P ð2πÞ 3 δ ð3Þ ð ⃗ 0Þ and the renormalization scale has been omitted. Similar results as in Eqs. (13) and (14) have been derived in the literature both in the high-energy limit [41,42] and in the heavy-quark limit [29], but not in the threshold region as we are interested in here.
In Eq. (14), GPDs are evaluated at a generic ξ. Near the threshold and in the heavy-quark limit, t approaches infinity, and ξ is close to 1, as shown in Eq. (7). Therefore, in principle, one should set ξ ¼ 1 in our equations above. However, there are two reasons for us to keep the generic ξ dependency.
First, the expression with generic ξ agrees with the leading-order result in the kinematic region at large W and small jtj [29,41,42]. What we have shown is that in the heavy-quark limit, the validity region of the leading-order factorization formula (13), (14) in terms of GPDs can be smoothly extended to the threshold region along the jtj min line. It shows that the large M V is sufficient to generate light-cone structure even in the threshold region. Therefore, our result can be viewed as a generalization of Ref. [29].
Second, although in the heavy meson limit ξ approaches 1 near threshold, in reality, especially the production of J=ψ, the kinematic finite meson mass corrections are important. By using the physical ξ, one can expect to take part of these kinematic corrections into account.
Given the amplitude above, we can calculate the cross section, and the result is with e Q the charge of the quark in the unit of proton charge, and the photon and meson polarization are summed over. The kinematic prefactor is the same as that in Refs. [14,15], and in the second line it is shown that the cross section of near-threshold meson production is related to the nonrelativistic wave function at origin and the gluon GPDs. Again, the renormalization scales in GPDs and wave function have been omitted and shall be on the order of M V .

IV. EXPANSION IN MOMENTS OF GPD
As shown in Eq. (14), the entire gluon GPD enters in the amplitude for the threshold region heavy quarkonium production because of the light-cone dominance. However, because ξ ∼ 1, one can Taylor expand the Aðx; ξÞ and express the Gðt; ξÞ in terms of the even moments of the GPD, The moments can be related to the matrix element of the gluonic twist-2 operator O μμ 1 …μ n−2 ν g [40] hP 0 jO μμ 1 …μ n−2 ν g jPi ¼ SūðP 0 Þγ μ uðPÞ X i;even By parametrizing the GPD, the moments for the scalar functions H and E can then be written as [40] where C-terms represent the highest power in ξ. In particular, for n ¼ 0, we have for the leading moments of the GPDs. As we will argue, the higher moment contributions shall be suppressed in the threshold region. If one only keeps the leading moment n ¼ 0, the Gðt; ξÞ becomes in terms of the matrix elements of gluon EMT in the proton. We then can derive a formula for the cross section in terms the gravitational form factors defined through following parametrization of the matrix element of gluon EMT [27], where A ¼ A 2;0 , B ¼ B 2;0 , and C ¼ C 2;0 are the same form factors as in Eq. (23). It leads to the following form of jGðt; ξÞj 2 after summing/averaging over the final and initial proton spin, where H 2 ≡ H 2 ðt; ξÞ and E 2 ≡ E 2 ðt; ξÞ are defined in Eq. (23). Combining with Eq. (17), the cross section of heavy vector-meson photoproduction can be expressed in terms of those gravitational form factors. This result agrees with the holographic QCD predictions [16,17] that the leading contribution to the cross section is due to exchange of 2 þþ excitations, or the spin-2 twist-2 operators, instead of the 0 þþ F 2 operators suggested in Refs. [3,7]. However, in the holographic approach, the twist-2 part of the gravitational form factor is dual to the graviton exchange and is free from the C g contribution. This differs from the generic QCD parametrization in Eq. (24). Because of the EMT conservation, the quantum anomalous energy F 2 form factor can be related to the twist-2 ones here. If the further limit ξ → 0 is taken in Eq. (26), only the A form factors are leading, and all the results agree. However, this is inconsistent with our approximation here.
Here, we return to the validity of the moment expansion in Eq. (18) and the leading moment approximation in Eq. (24). Notice that for generic ξ < 1, the GðtÞ has imaginary part ImGðt; ξÞ ∝ F g ðξ; ξ; tÞ þ F g ð−ξ; ξ; tÞ; which implies that the moment expansion in Eq. (18) has some limitations. However, as ξ gets closer to 1 in the heavyquark limit, the above imaginary part vanishes in power of 1 − ξ. The renormalization group evolution will also help to improve the convergence of the expansion at large renormalization scale μ ∼ M V . Indeed, it has been shown that the asymptomatic form of gluon GPD reads [43,44] F asym which vanishes quadratically at x ¼ ξ. With this form of asymptotic behavior, one can show that the expansion in Eq. (18) is convergent. To summarize, Eq. (18) should be a good approximation for Gðt; ξÞ in the threshold region where ξ is close to 1.
Regarding the validity of the leading moment approximation in Eq. (24), one can read from Eq. (18) that, due to the ξ 2n in the denominator, Eq. (18) is simultaneously a moment expansion and a ξ −1 expansion. The closer the ξ to 1, the larger the contribution from the leading term n ¼ 0.
Using the asymptotic form in Eq. (28), in the exact μ ∼ M V → ∞ limit, the ratio between the contribution from higher (n ≥ 1) moments and the leading moment is and the total higher-moment contribution is 25%. Notice that in the large Q 2 region, a similar estimation in favor of the leading moment dominance was first performed in Ref. [31]. Realistically, when μ ∼ M V is not very large and ξ is not exactly at 1, we can use the ratio between the second and the leading terms to estimate the effect of the high-order terms. Using the explicit expressions for GPD moments in Eq. (21), the ratio between second and leading moments can be calculated as for the moments of H, and B g 4;0 þ 4ξ 2 B g 4;2 − 16ξ 4 C g for the moments of E. For a quick estimation, one can keep only the A 0 form factors which equal to moments of gluon parton distribution functions (PDFs) at t ¼ 0. Away from t ¼ 0, one can use the following dipole model: Here, f g is the gluon PDF at zero momentum transfer, and m A and m A4 are unknown dipole masses. Neglecting the difference in dipole masses and using the recent CTEQ global analysis [45], at renormalization scale μ ¼ 1.3 GeV, one has For J=ψ production, one has ξ ¼ 0.6 right at the threshold, and the ratio is around 0.1, consistent with the dominance of the leading moment. Away from the threshold, the ratio increases as ξ drops, and the dominance of the leading moment becomes less pronounced. As one increases the renormalization scale, more momentum fraction is carried by small x gluons, and the ratio becomes smaller. However, since the form factor C is also important for the amplitude, a more realistic estimation needs input from lattice calculations.

V. GRAVITATIONAL FORM FACTORS FROM THRESHOLD DATA
As we have discussed in the previous section, threshold quarkonium photoproduction in the heavy-quark limit is dominated by the nucleon's gravitational form factors in the leading moment approximation. In this section, we study the phenomenology of determining these gravitational form factors from realistic heavy quarkonium production data, neglecting various higher-order corrections which we will take into account in the future publication. We consider the J=ψ photoproduction where there are more data available near the threshold and predict the near-threshold ϒ photoproduction cross section where the heavy-quark expansion works better.

A. J=ψ photoproduction
To start with, we consider the J=ψ photoproduction total cross sections from SLAC [46], Cornell [47], and the most recent GlueX experiments [36]. Given that the GlueX data seem to deviate from the SLAC and Cornell data, we will focus on comparing with the GlueX result, which includes more data as well as some measurements of differential cross section.
As for Gðt; ξÞ, the input of gravitational form factors is required. At large momentum transfer, it has been argued that the form factors decay polynomially [51][52][53] based on power-counting methods. Here, we use the gravitational form factors from lattice calculations, which model those form factors with dipole expansion In a recent calculation, those parameters are found approximately, m A ¼ 1.13 GeV, m C ¼ 0.48 GeV, A g ð0Þ ¼ 0.58, and C g ð0Þ ¼ −1.0 [35], while the B g ðtÞ form factor is numerically small. These inputs without any fitting parameters yield the curve shown in Fig. 5.
On the other hand, our theoretical formulas allow us to extract the gravitational form factors from the J=ψ photoproduction data. Since the data are quite limited at this stage, We choose to fit the GlueX total cross section and differential cross section combined. If we have both A g ð0Þ and m C fixed to be the lattice values in order to avoid overfitting, we get m A ¼ 1.64 AE 0.11 GeV, C g ð0Þ ¼ −0.84 AE 0.82. In Figs. 6 and 7, we compare the fit results with the data from GlueX. The uncertainties of m A and C g ð0Þ indicate that the data are more sensitive to m A rather than C g ð0Þ in the region of GlueX data. One can of course explore other fitting scenarios as well.

B. ϒ photoproduction
The above result can be used to predict the ϒ photoproduction rate near threshold with M ϒ ¼ 9.46 GeV and the wave function at origin for ϒ [49,50], Considering the effect of running coupling constant, α S ¼ 0.2 is used for the ϒ production. Then, we have the predicted total cross section as shown in Fig. 8 and the predicted differential cross section in Fig. 9. The large mass of ϒ implies that our calculation in the heavy meson limit works better for ϒ production, as there are fewer theoretical uncertainties from higher-order correction. However, the production rate is suppressed by the heavy meson mass, and thus the cross section is much lower than J=ψ.
FIG. 6. Fit total cross section for J=ψ production compared with the total cross section measured at GlueX [36]. The 95% confidence band is shown as the shaded region hereafter.
FIG. 7. Fit differential cross section for J=ψ production compared with the differential cross section at W ¼ 4.58 GeV measured at GlueX [36].
C. Quantum anomalous energy, mass radius, and pressure distribution The gravitational form factors measured from J=ψ photoproduction can be used to study mass, spin, and pressure properties of the nucleon [4,7,27,32,34].
There are two methods to access the quantum anomalous energy through the matrix elements of the EMT. The first is through the matrix element of F 2 suggested in Ref. [2]. In our factorization, this is a subleading contribution. Alternatively, one can also access the same quantity by considering the trace of the full EMT, where A q ð0Þ and A g ð0Þ are traditionally interpreted as the momentum fractions carried by quarks and gluons, respectively. The above relation comes from the conservation of the EMT. Thus, fitting A g ð0Þ through heavy-quarkonium production data, one obtains the gluonic contribution to the quantum anomalous energy, and the result can be compared with the vector-dominance model analysis [2,54]. Of course, it also provides an alternative determination of the gluonic momentum fraction in the nucleon. It shall be clear that our formula for the cross section works for the large t. To extrapolate to t ¼ 0, one has to learn the t dependence through lattice QCD calculations. The scalar and mass radii of the proton are defined as [34] The two different mass radii are related with the Cð0Þ term as The gravitational form factors here are for the total EMT, so one has AðtÞ ¼ A q ðtÞ þ A g ðtÞ and CðtÞ ¼ C q ðtÞ þ C g ðtÞ.
For the quark contribution, we use the lattice data C uþd ð0Þ ¼ −0.267 or C uþd ð0Þ ¼ −0.421 depending on the extrapolation method [55]. As for A uþd ðtÞ, one could also fit the form factor with the above dipole expansion as for which one has the dipole mass m Di;uþd ≈ 1.8 GeV and A uþd ð0Þ ≈ 0.5 [56]. Then, the quark form factors combined with the above gluon form factors from fitting give and then we have One should be aware that those results are associated with large uncertainties resulting from the lack of precision for the fitted value C g ð0Þ ¼ −0.84 AE 0.82. The gravitational form factor C also provides a direct access to the pressure and shear force distributions of the nucleon [32,[57][58][59]. Given our fit value C g ð0Þ ¼ −0.84 AE 0.82 and combining with the lattice data m C ¼ 0.48 GeV in the dipole assumption (38), the gluon contribution to the FIG. 9. Solid line shows the predicted differential cross section at W ¼ 11.5 GeV for ϒ photoproduction near threshold. pressure distribution is shown as Fig. 10. Here, we have neglected theC contribution.
The above results serve as an example of measuring the gravitational form factors from the photoproduction of heavy vector mesons. Although the results rely on the input from lattice calculations as well as the models for the form factors, it can be a useful tool when more data are available and higher-order corrections are taken into account. One can also make a full GPD-based analysis directly using Eq. (14).

VI. CROSS SECTION WITH POLARIZATION
Up to now, we have only considered the unpolarized cross section. In this section, we study the polarization effect. By measuring the polarization of the photon and the vector meson, we are able to verify the specific structure ε Ã V · ε of the factorization formula (13) which will serve as a crucial test of the factorization formula itself. Since Eq. (13) also predicts that the matrix element Gðt; ξÞ contains all the proton polarization dependence, more information on gravitational form factors can be extracted by measuring the polarization of the proton. Notice that the measurement of the polarization of the final-state proton can be nontrivial [60,61]. In this section, we shall study the polarized cross section from two different point of view, the vector-meson polarization and the proton polarization.
Polarization observables allow us to separate out the contributions from different form factors. In particular, even though the B q;g ðtÞ is small, the polarized cross sections can be sensitive to the B g ðtÞ form factor which is related to the gluon angular momentum [27]. If given sufficient data with high precision, it might be possible to shed some light on the gluon angular momentum contribution to the proton spin.

A. Vector meson polarization
The total cross section is proportional to jMðε V ; εÞj 2 , and thus one has from Eq. (13) In general, one could define the photon polarization tensor and meson polarization as ρ μν ðqÞ and ρ Vμν ðKÞ so dσ dt ∝ ρ μν ðqÞρ μν V ðKÞ. For a physical photon or meson, the polarization tensor is simply while in the case in which they are virtual particle coupled with leptonic current, the polarization tensor can be written with those leptonic current as ρ μν ðqÞ ¼ūðl 1 ; S 1 Þγ μ uðl 2 ; S 2 Þūðl 2 ; S 2 Þγ ν uðl 1 ; S 1 Þ; ð53Þ with l i and S i the momenta and the polarization vectors of the lepton and antilepton pair. They satisfy S i · l i ¼ 0, S 2 i ¼ −1, and similarly for the vector meson. If we consider the photoproduction process when the photon is on shell and unpolarized, the polarization tensor can be written as whereq i is the unit vector in the direction of photon threemomentum. As for the vector meson, assuming the final leptons are unpolarized and averaging over the lepton polarization, we have ρ Vμν ðKÞ ¼ − 1 2 K 2 g μν þ l 1;μ l 2;ν þ l 2;μ l 1;ν ; ð55Þ with l 1 , l 2 the momenta of the two leptons satisfying l 2 þ l 1 ¼ K, and thus we have The last termq · ⃗ l 1q · ⃗ l 2 indicates that the cross section depends on the angle between the momentum of outgoing lepton and the momentum of initial photon and can be examined by experiments.

B. Proton polarization
Here, we consider the polarization of the initial and final proton. For a demonstration, let us consider the asymmetric difference δjGðt; ξ; S 0 ; SÞj 2 as δjGðt;ξ;S 0 ;SÞj 2 ≡jGðt;ξ;S 0 ;SÞj 2 −jGðt;ξ;S 0 ;−SÞj 2 ; which measures the difference in the cross section when the spin S μ of the initial proton has been flipped. Notice the polarization vector satisfies S · P ¼ 0 and S 2 ¼ −1 and similarly for S 0 . Then, one has FIG. 10. The gluon contribution to the pressure distribution inside the proton from the combined lattice calculation [59] and our fit to GlueX data.
δjGðt; ξ; S 0 ; SÞj 2 ¼ where E 2 ≡ E 2 ðt; ξÞ and H 2 ≡ H 2 ðt; ξÞ are leading moments of GPDs defined in Eq. (23). When the initial proton is transversely polarized, one can show that In particular, if both S and S 0 are perpendicular the scattering plane, then only the first term in Eq. (59) survives, and the cross section is nonzero only when S, S 0 are parallel. For generic transverse spin S, S 0 , the second term is also nonzero. By measuring the cross section at different S, S 0 , one is allowed to extract all three form factors.

VII. CONCLUSION
In this paper, we have made a heavy-quark expansion for the vector-meson production at the threshold region in QCD. The result shows that the cross section can be used to measure the form factors of the gluonic EMT. Before ending the paper, we would like to make several comments.
First, our calculation is performed only at leading order in Oð M N M V Þ and Oðα S Þ. For application to J=ψ production, since the J=ψ mass is not heavy enough, one expects large mass correction and high-twist effect. A comprehensive study of the mass correction is definitely very important but is beyond the scope of the current paper. Beyond leading order in Oðα S Þ, there are quantum corrections both from the emission of virtual gluon and from the internal motion of the vector meson [29]. For production of J=ψ, the order Oðα s Þ corrections can be significant [29]. A calculation of these effects will be left to a future work.
Second, we should emphasize again that our result should be viewed as a generalization of the leading-order factorization in Ref. [29]. The kinematic region considered therein is large W and small momentum transfer t, corresponding to the large W part of the t min line. The heavy-quark limit is performed after taking the large W limit. In this work, we have extended the factorization along the t min line into the entire threshold region where W is of the same order of M V and is a single ultraviolet scale. While the factorization seems to work in leading and next leading orders, the validity of the factorization theorem to all orders in perturbation theory remains to be established.
Third, although in this paper we considered only the small Q 2 ∼ 0 region, at large Q 2 , the factorization formula in terms of GPD should remain valid. In fact, as Q 2 increase, the incoming proton becomes faster, and the momentum transfer becomes larger. The light-cone structure shall be more relevant. This picture differs qualitatively from the Euclidean OPE approach in Ref. [30].
Forth, we comment on the scalar contributions. In Eq. (13), only the F þi F þi component has been kept due to the standard power-counting rule of collinear gluon. If one keeps all the terms, then the F þ− F þ− contribution, corresponding to the F 2 contribution, is proportional to 1 − ξ 2 and is suppressed in the heavy-quark limit. A detailed study of F 2 and the three-gluon contribution is left to future work.
Finally, one should also notice that, although in the current work only two-gluon exchange has been considered, starting from twist-3, there can be three-gluon contributions as well. It has been argued in Ref. [18] that in the near-threshold region the three gluon contributions might be important. To further clarify the importance of the three and more gluon contributions requires a thorough power-counting analysis in the threshold region. A thorough study of three-gluon contribution will be left to a future work.
In conclusion, we have shown that the factorization formula in terms of gluon GPD for the photo or elector production of heavy vector meson in Ref. [29] remains to be valid in the near-threshold region. At a large vectormeson mass, the skewness ξ is close to 1, and the amplitude is dominated by the leading moment of the gluon GPD, the gravitational form factors of the proton. This allows us to extract alternatively the gluonic contribution to the momentum and anomalous quantum energy. It also allows us to extract the gluonic C form factor (or D terms) and determine the mass radius and pressure distribution of the proton. We have applied our formalism to the case of unpolarized J=ψ production with GlueX data, and the result is encouraging.