NLO corrections to the deeply virtual meson production revisited: impact on the extraction of generalized parton distributions

,


Introduction
In the family of structure functions describing the quark-gluon structure of hadrons, generalized parton distributions (GPDs) [1,2,3] occupy a special place.Their understanding would provide insights into several tangible physical properties of hadrons, such as the composition of their spin [4], 3D tomography [5], or pressure [6].Simultaneously, the processes enabling their experimental determination are measurable in both current and forthcoming experimental facilities [7,8].Processes that in this sense are the subject of the most attention of the scientific community are deeply virtual Compton scattering (DVCS) and deeply virtual meson production (DVMP), for whose phenomenological status see [9,10,11] and references therein.Other processes have also been identified as possibly important sources of information about GPDs, such as timelike DVCS [12,13,14,15], double DVCS [16,17,13,18], or boson pair production [19,20,21,22,23], see also [24], but experimental measurements of these are so far scarce or nonexistent 1 .So in this paper, we focus on DVCS and DVMP, in particular, the production of longitudinally polarized vector mesons (DVV L P). Notably, transversely polarized vector meson production vanishes at leading-twist [31].Moreover, at experimentally accessible energies, pseudoscalar meson production appears to be predominantly influenced by higher-twist effects [32,33].
The significance and complementary nature of DVCS and DVMP in unveiling GPDs have long been recognized.This is evident in the ongoing efforts to enhance the precision of their theoretical 1 The list is not exhaustive and we refer to, for example [25], for the account on processes providing complementary information on GPDs from exclusive processes with a large momentum transfer to a nucleon (large t) such as nucleon form factors [26], wide-angle Compton scattering [27,28] and wide-angle meson production [29,30], as well as crossed processes.Their analysis is experimentally and theoretically more challenging.
descriptions by calculating higher-order perturbative QCD (pQCD) corrections to subprocess hardscattering amplitudes: NLO DVCS [34,35,36,37], NLO DVMP [38,39,40], and NNLO DVCS [41,42,43].Additionally, the power-corrections to DVCS have been determined [44,45,46,47,48], while the consistent inclusion of power-corrections in DVMP remains a challenging endeavor.The leading twist-2 contribution to DVMP accounts for the contribution of longitudinally polarized photons (γ * L N → MN ′ ).Currently, there remains a challenge in systematically distinguishing between experimental data associated with longitudinal and transverse contributions.Conversely, for DVCS the leading twist pertains to the contribution of transverse photons (γ * T N → γ N).Unfortunately, the application of the existing theoretical apparatus to phenomenology and to the extraction of GPDs is currently not at a satisfactory level.Analyses most often focus on only one of these processes, and even then most often they remain at the level of the leading order (LO) of the perturbation theory.Bearing in mind that in the most studied DVCS process the gluon contribution appears at NLO, it is clear that the resulting extractions of GPDs are not reliable.The researchers themselves are mostly aware of this; therefore, the largest number of previous studies stops at the determination of Compton form factors (CFF) that combine GPDs of all flavors and corresponding hard scattering amplitudes.
Notable earlier attempts to grasp GPDs by simultaneous analysis of DVMP and DVCS processes are: • Ref. [49], where it was shown that the GPDs obtained from an earlier fit to DVMP data [50,51,32] also relatively reasonably describe the observables of DVCS, in particular those from the HERA collider.The approach used relies on inclusion of transversal degrees of freedom in a way for which it is not clear how to naturally extend it to higher orders of perturbative QCD expansion.• Refs.[52,53,54], where it has been shown that a relatively acceptable simultaneous fit of GPDs to DVV 0 P and DVCS data obtained from the HERA collider at LO and NLO is possible.
The continuous increase in the volume and precision of available experimental measurements, along with the promising capabilities of future research facilities, necessitates the development of a robust theoretical framework.This framework should be on par with those established in related fields, such as the extraction of ordinary parton distribution functions (PDFs) and transverse momentum distributions (TMDs), extending at least to the NLO level, and potentially advancing to the NNLO level in the future.The foundation of such a framework, which is based on the conformal partial wave (CPaW) formalism [55,56] -specifically, the conformal partial wave expansion in conjunction with the Mellin-Barnes integral technique -was established and systematically organized for DVCS and DVMP in [57,58].In comparison to the traditional momentum fraction representation, it not only offers easier inclusion of GPD evolution to NLO order and beyond, but also opens the door to intriguing GPD modeling possibilities and paves the way for the development of stable and efficient computer codes capable of handling GPD evolution and fitting to both experimental and lattice data.
The comprehensive list of NLO contributions to deeply virtual production of vector (scalar) mesons in both the conformal moment and momentum fraction representations has been systematically organized in [58].Similarly, the corresponding expressions for pseudoscalar (axial-vector) meson production, were derived in [40], where additionally an omission in the NLO expressions for DVV 0 L P was recognized.Here we present the corrected NLO expressions for DVV 0 L P in the conformal moment representation.The main motivation of this paper is to continue the application of CPaW formalism to NLO analysis of DVMP and extraction of GPDs simultaneously from DVMP and DVCS data, for the first time systematically at the NLO level.The only predecessor is the unpublished work [53], with respect to which we bring here the following improvements: • The expressions for the NLO hard-scattering DVV 0 L P amplitude have been corrected in the meantime [40].
• Here we treat the measured cross sections as a reliable experimental result.In paper [53] the normalizations of the experimentally measured cross sections for DVV 0 L P were treated as fitting parameters.
• The GPD model we use here is of the same type as in [53], but simpler, with fewer free parameters.
In addition to the improved extraction of GPDs, the aim of this work is to study more carefully the impact of NLO corrections.We find that despite the very different roles of quarks and gluons in DVCS and DVV 0 L P, both processes at the NLO level give very similar GPDs, thus pointing out at the same time the universality of GPDs and the unquestionable need for a NLO theoretical approach.
In previous papers mostly dedicated to the calculation of NLO corrections, and not so much to applications, the impact of NLO corrections is usually studied in a simplified environment, using some fixed toy model, the same for LO and NLO.With such an approach it is possible only to get some idea about the size of the coefficients in a given order of pQCD expansion (see, for example, detailed analysis in [57,58]).Here, we want to answer a complementary, phenomenologically more important question: How does the shape of the GPD change when going from LO to NLO?In this sense, we exactly follow the proposition of the authors of [58] who consider their analysis of the impact of NLO corrections "not entirely realistic" and suggest that "appropriate method [. ..] would be the quantification of reparametrization effects that arise from LO and NLO fits to experimental data." To make our main message, that universal NLO GPDs can describe DVMP and DVCS processes at high energies (small Bjorken x B ) in the collinear pQCD approach, as clear as possible and free from technical complications, we assume that the only contribution to the amplitudes comes from the dominant sea quark and gluon GPDs H sea and H G , which is a good approximation in our chosen kinematics.Also, the only observables we consider are the cross sections for the processes γ * p → γ p (DVCS) and γ * p → ρ 0 p (DVMP).Other observables or lattice data, in other kinematic regimes, would require the addition of valence quarks, other GPDs, and the improvement of our GPD model, which we leave for future works.
In Section 2 we outline the main elements of collinear perturbative QCD framework and specifically of CPaW formalism.We summarize the LO and NLO contributions to subprocess hard-scattering amplitudes for DVV L P and DVCS, and explain the implementation of GPD evolution to NLO.In Section 3 we specify the GPD model we use, and in Section 4 the used experimental measurements are presented, with special emphasis on the problem of separation of the longitudinal and transverse components of the cross section for meson production.In Section 5, we present and discuss our results, while the concluding remarks are given in Section 6.

Perturbative QCD framework
In this paper, we use the standard collinear perturbative QCD approach [59,60] for both observed processes, DVCS and DVMP, specifically DVV L P, Fig. 1.This is worth emphasizing considering that the first attempts to apply this approach to DVMP [61,62] gave results for cross sections in discrepancy with measurements.The reason for this is, among other things, the effective ∼ 1/Q 4 scaling of measured cross sections, visible in Fig. 8, inconsistent with the canonical asymptotic ∼ 1/Q 6 scaling that follows from the collinear pQCD amplitude.In the currently most popular approach to the calculation of DVMP cross sections within GPD framework [50,51], this is fixed by introducing power corrections due to transverse momentum of partons.This approach has proved to be quite successful phenomenologically, but it is not clear how to apply it consistently at NLO.As we confirm in this paper, NLO corrections are significant, and the authors of the early studies of these corrections [63] concluded that for quantitative control in the high-energy regime, it is necessary to perform a resummation of large BFKL type log(1/x B ) logarithms [64].In the point of view to which we adhere [58], the change of effective scaling is achieved by log Q 2 perturbative logarithms.The collinear perturbative QCD approach relies on the property of factorization of the amplitude into the non-perturbative soft structure functions and the subprocess amplitude of hard scattering on the active parton [65,66], see Fig. 2.This approach provides the leading-twist contributions that describe the transversal DVCS (γ * T p → γ p) and longitudinal DVV L P (γ * L p → V L p) cross section components.In the following, we present the main points of this approach, taking the opportunity to collect and summarize the relevant expressions, with a particular focus on the CPaW formalism.Additional details can be found in [57,58,40].
The observables are expressed using process amplitudes termed Compton form factors (CFFs) for DVCS and transition form factors (TFFs) for DVMP.The nomenclature aligns with the nomenclature used for the contributing GPDs.Additionally, TFFs depend on meson distribution amplitudes (DAs), representing the meson's internal structure, making the analysis of the process more challenging, but also potentially more rewarding.
In the case of DVCS at twist-2, contributions arise from intrinsic parity even (vector) GPDs denoted as H a and E a , as well as intrinsic parity odd (axial-vector) GPDs denoted as H a and E a2 .Here, the variable a ranges over different quark flavors (u, d, s, . . . ) and gluon contributions (G).The gluon contributions to the subprocess hard-scattering amplitude of DVCS emerge at NLO.However, their impact on CFFs at LO arises through GPD evolution, attributed to the mixing of quark singlet and gluon GPDs.
In contrast, DVMP provides a natural distinction of GPDs of different parity.For the specific case of vector meson production (DVV L P), only the intrinsic parity even (vector) GPDs H a and E a contribute.On the other hand, the GPDs H q and E q come into play in the production of pseudoscalar mesons where there are no GPD gluon contributions 3 .Notably, in the production of neutral V L mesons, the contributions of gluon GPDs H G and E G are more significant, as, unlike DVCS, these gluon GPDs already contribute to the LO subprocess hard-scattering amplitude (see Fig. 2).
In the twist-2 approximation and for the small x B (regime of current interest in phenomenological 2 In this work we do not consider chiral-odd transversity GPDs.analysis) the differential cross section for DVCS (γ * T p → γ p) is given by where α em is the electromagnetic fine structure constant, M p is the proton mass, Mandelstam t = (p 1 − p 2 ) 2 , Q 2 = −q 2 1 is the virtuality of the initial photon, and x B = Q 2 /(2p 1 • q 1 ) is the Bjorken variable.Momenta p 1,2 and q 1,2 are specified on Fig. 1.In this approximation three, out of the four, CFFs contribute, H, E and H, and we neglect terms suppressed by further powers of small x B .For small x B arguments from Regge theory suggest H ∼ 1/x B and H ∼ 1/ √ x B , so it is justified to neglect the contribution of H.The situation with E is less clear, but we still rely on the suppression factor t/(4M 2 p ) and neglect E in the numerical treatment.
Similarly, at twist-2 the differential cross section for DVV L P (γ * L p → V L p) in the small x B regime simplifies to where H V L and E V L are TFFs (p → V L p) for the light vector meson V L .For the same reason as in the DVCS case in this work we ignore the contribution of E V L .

Factorization of form factors and CPaW formalism
As argued above, considering the specific energy regime of interest in this study, the differential cross sections are expressed solely in terms of the Compton form factor H for DVCS and the transition form factor H V L for DVV L P. In the following we illustrate the factorization of form factors using these two as examples.The given expressions remain valid across all energy regimes and are applicable to other form factors, in particular to E and E V L , taking care to use the subprocess hard-scattering amplitudes suitable for the intrinsic parity involved.
Both CFF H and TFF H V L , according to the factorization property, can be written as a convolution of the corresponding subprocess hard-scattering amplitudes T DVCS and T DVV L P and non-perturbative soft functions, GPD H and, in the case of DVV L P, the meson distribution amplitude φ V L : ⊗ and v ⊗ are defined in (9) below.The GPDs are functions of three variables: x -the parton's "average" longitudinal momentum fraction, ξ -the longitudinal momentum transfer (skewness parameter), and t -the momentum transfer squared, while their evolution with energy is encapsulated in the dependence on the factorization scale µ GPD .For the deeply virtual processes under consideration, skewness relates to x B through (5).The hard-scattering amplitude at twist-2 depends on x, ξ, Q 2 , renormalization and factorization scales and, in the DVMP case, v -the longitudinal momentum fraction of the leading Fock state parton in a meson 4 .In this collinear framework, the dependence of CFFs and TFFs on momentum transfer t solely arises from GPDs.
Formally, the all-order predictions of factorized expressions (3) and ( 4) are scale-independent with respect to renormalization (µ R ) and factorization (µ GPD , µ DA ) scales.However, finite-order predictions do exhibit residual scale dependence.Notably, the renormalization scale dependence is particularly pronounced in DVMP, as it comes into play already at LO in contrast to the DVCS where it appears at NLO (see Fig. 2).Consequently, the stabilization effect of higher-order corrections starts at NLO for DVMP and at NNLO for DVCS.In both processes, gluon contributions are proportional to α s and therefore sensitive to the choice of µ R .The scale settings are discussed in more detail in [57,58].It's worth noting that there are proposals in the literature regarding the 'optimal' scale settings and/or different resummation of factorization logarithms.However, these aspects are left for future investigation.For numerical evaluations in this study, we set the renormalization scale and both factorization scales to be equal to the photon virtuality, denoted as The evolution of GPDs and DAs with factorization scale is discussed below in Sec.2.3.
In (3) and ( 4) the superscript a denotes the flavor of the parton off which hard scattering occurs, 4 In the collinear approximation, the momentum fractions of the emitted and reabsorbed partons are u = (ξ + x)/(2ξ) and ū = 1 − u, respectively.Frequently, hard-scattering amplitudes are expressed in this notation [58].The appropriate "+iϵ" prescription in DVCS and DVMP kinematics is restored by employing ξ → ξ − iϵ. and the sum over active flavors is implied.Different flavor compositions of individual mesons make it, in principle, possible to separate the contribution of individual flavors to GPDs, which is one of the main reasons for studying the DVMP process.Flavor decompositions of CFFs have been systematically organized in [57] (Sec.3.1), while TFFs specifically for various mesons are systematized in [58] (Sec.2.2).In this work we are interested in neutral vector mesons.Taking into account the mixing of flavor singlet quark and gluon GPDs under evolution, it is convenient to decompose DVCS CFFs and DVV 0 L P TFFs into a flavor nonsinglet (NS) and singlet (S) part, where the singlet part is given as a sum of gluon (G) and flavor singlet quark (Σ) contributions: The factors Q 2 NS = 1/9 (1/6) and Q 2 S = 2/9 (5/18) depend on the number of active flavors depends on the meson type.In this work we use n f = 4, and For small x B , following Regge theory and by analogy with standard PDFs, the contribution of valence quarks is expected to be suppressed by ∼ √ x B , so we neglect it and assume that all of the amplitude comes from sea quarks and gluons.We further assume that the quark sea is flavor symmetric and given by the singlet combination and there is no nonsinglet contribution in our phenomenological analysis.The symmetry properties of the subprocess hard-scattering amplitudes and the meson DA project GPDs of well defined signature σ, as defined in [58], Eqs.(3.1), (3.9) and (3.10) 5 .In [58] the corresponding subprocess hard-scattering amplitudes were further simplified using the symmetry properties of quark (antisymmetric) and gluon (symmetric) GPDs and the vector meson DA (symmetric), and the use of GPDs of definite signature is understood.In [40] an analogous treatment of DV production of pseudoscalar mesons is described.Consequently, in the following the superscript a in (3) and ( 4), apart from information on flavor content, carries the information on signature.
The convolutions x ⊗ and v ⊗ in equations ( 3) and ( 4) entail integrations over the longitudinal momentum fractions of the partons, as expressed by: However, in this work, we apply a transformation from the longitudinal momentum space x into the space of conformal moments j (j-space) to all hard-scattering amplitudes, GPDs, and DAs.For singlet quark and gluon GPDs such a transformation is given by a convolution with Gegenbauer polynomials 5 For quark GPDs, σ = +1(−1) denotes antisymmetric (symmetric) behavior under x → −x, and the reverse holds for gluon GPDs.GPDs with well-defined symmetry under x → −x exchange are distinguished by superscripts (±).These superscripts indicate the charge parity of two partons, which, when multiplied by the signature, gives the intrinsic parity.Consequently, H q,G(+) and H q,G(+) contribute to DVCS, while H q,G(+) and H q(−) contribute to DV production of neutral vector and pseudoscalar mesons, respectively.The same holds for E and E.
for integer j, where the normalization coefficients are expressed via Euler's gamma function.For (the quark) meson DA the corresponding transformation is for an integer k also given by a convolution with Gegenbauer polynomials The normalizations of GPD and DA conformal moments differ.The former ensures that in forward kinematics (t = 0 and ξ = 0) conformal moments of GPDs coincide with corresponding PDF Mellin moments, while the latter is fixed by The symmetry properties of GPDs and the DA, whether symmetric or antisymmetric under x → −x and v → 1 − v, are transferred to even and odd parities of conformal moments.
As singlet quarks Σ and gluons G mix under evolution with scale, all relations for DVCS and DVV 0 L P are effectively matrices in (Σ, G) space.Thus the transformation of the hard-scattering amplitude into the space of conformal moments can be written as and in the case of DVV 0 L P additionally Applying these transformations, the convolution for H S becomes an infinite divergent sum and similarly for H S V 0 L , with additional summation over the meson DA conformal moments Σ ∞ k=0,even .In accordance with the CPaW formalism, we perform the resummation of the sum ( 14) through Mellin-Barnes integration in the complex plane.This yields the final expression for DVCS CFF and, completely analogously, for DVMP TFF where T and H are analytically continued in the complex j plane.The intersection c of the integration contour with the real axis should be to the right of the Regge poles and the poles of anomalous dimensions and hard-scattering amplitudes, and to the left of the poles of the tangent function.In this paper we use the value c = 0.35.The analogous expressions without gluon contributions are valid for H NS and from (7).As mentioned earlier, these contributions are negligible for the kinematics of interest and have not been considered in our numerical analysis.
We have elucidated here the transition from the usual momentum fraction representation to the conformal moment representation and the CPaW formalism, where the form factors ( 3) and ( 4) are expressed in the form [57,58] where Here one employs tan or cot in accordance with signature σ.The summation over k, i.e., DA conformal moments, can, in principle, also be replaced by a Mellin-Barnes integral when the analytical form is known.We will provide a summary of the relevant conformal moments (13) for DVCS and DVV L P subprocess hard-scattering amplitudes in the following section.

Hard-scattering amplitudes 2.2.1 DVCS
The DVCS hard-scattering amplitude contributing to the singlet intrinsic parity even (vector) CFF (15) and determined to NLO order in MS scheme reads in the conformal moment space [67, 57] with the LO part simply and the NLO part G C (1) j The parameter n f is the number of active quark flavors, (z) n is the Pochhammer symbol S 1 (z) is the harmonic number with Euler-Mascheroni constant γ E = 0.577 215 66, while ΣΣ γ (0) j and ΣG γ (0) j are elements of the LO anomalous dimensions matrix, given below in (28).
The hard-scattering amplitude contributing to H NS , which we do not include in the numerical evaluation, is, to this order, identical to Σ T DVCS j .The same subprocess hard-scattering amplitudes contribute to the CFF E, which is neglected in the current phenomenological analysis.Similarly, in this work we do not consider the parity-odd (axial-vector) hard-scattering amplitudes contributing to H and E TFFs.These amplitudes are listed in [57] (Section 4.2).

DVMP
Unlike DVCS, for DVMP the LO is already proportional to α s and in the singlet sector of interest, gluons contribute at LO (Fig. 2).The perturbative expansion of the DVV 0 L P hard-scattering amplitude contributing to the singlet intrinsic parity even (vector) TFF ( 16) takes the form [58,40] where factor 3 is a consequence of the normalization of the meson DA, and the LO term reads At the NLO order, the so-called pure singlet contribution pS T (1) appears for the first time Note that q T jk corresponds to the nonsinglet part that contributes to the, in our numerical analysis neglected, TFF H NS V 0 L from (7).It is convenient to separate amplitudes T jk according to contributions of individual classes of Feynman diagrams, i.e., to make a color decomposition (1) G T (1) where Conformal moments are given in [58]: the nonsinglet moments c in Eq. (4.53).However, as was stated in [40], the factorization of the gluon contribution employed in [58] (as well as in preceding work [39,68]), was a nonstandard one, i.e., differed from the one used in the case of DIS and DVCS.Hence, corrections were required for both the pure singlet and gluon contributions 6 , as elucidated in [40] (specifically in Eq. ( 20) and below).For convenience we list the corrected expressions for the pure singlet moment pS c (1) In the erratum [39,68], only the quark contributions were corrected.The mixing of quark and gluon contributions under evolution leads to corrections in both contributions: the correction of the NLO quark contribution depends on the LO gluon contribution, and vice versa.Thus, in the erratum [39], the correction to the gluon contribution was overlooked.In contrast, due to the vanishing LO quark contribution, the erratum [68] for the J/Ψ production is correct, as confirmed also by [69].and the gluon conformal moment 7 These replace pS c (1 jk and G c (1,F) jk in Eqs.(4.48) and (4.53b) in [58].Here ∆S 2 is defined as in [58], Eq. (4.13), while pS ∆c jk and G ∆c (1,F) jk are given in [58], Eqs.(4.48b) and (4.53f), respectively.The k are to be found in the next section.Finally, the same subprocess hard-scattering amplitudes contribute to the TFF E V 0 L , which is neglected in the current phenomenological analysis.The hard-scattering amplitudes for deeply virtual production of pseudoscalar mesons associated with H and E GPDs, are listed in [40].

Evolution of GPDs and DAs
In the context of the longitudinal momentum fraction space (often referred to as x-space), GPD evolution at LO has been implemented in computer code for a long time [70].However, this implementation has not yet been widely utilized for global fits.There is optimism that a new approach, as developed in [71], when combined with the PARTONS framework [72], may lead to advancements in this area.Additionally, there exists an implementation of NLO evolution in x-space [73], but this code has proven to be challenging to use effectively.One of the main advantages of working in the space of conformal moments (j-space) is the simplicity and numerical stability it offers for GPD evolution, which enables us to work at full NLO level.In conformal moment representation the evolution operator is diagonal at LO, and conformal moments evolve autonomously.Closed analytical expressions are available for both diagonal and non-diagonal terms appearing at NLO and beyond (for detailed account see [57]).The CPaW framework thus offers the potential for direct extension to the NNLO level.In the case of DVCS, preliminary analysis has been conducted in a specialized conformal scheme, leveraging DIS NNLO results [56,67,57].Additionally, there is already a significant number of components available in the MS scheme [41,42,43,74,75].

Evolution operators
The operator E(µ, µ 0 ; ξ) governing the perturbative evolution of the singlet intrinsic parity even (vector) GPDs of interest in this work, can be defined as In addition to the omission noted in [40], during the work on this article, in one of the terms in G c (1,F) jk in (24b) a typo in the sign was identified and corrected.
where at NLO accuracy and in the MS scheme, the operator is non-diagonal in both the two-dimensional flavor space (Σ, G) and in the infinite-dimensional space of conformal moments.We can write the perturbative expansion of this operator in the form Here, the summation is performed over the eigenstates a, b ∈ {+, −} of the LO evolution operator in (Σ, G) space, where the projectors P a j onto those states are Here γ (0) j is a 2 × 2 matrix of anomalous dimensions whose elements for the intrinsic parity even (vector) case read [76,77] ΣΣ γ while are its eigenvalues.The diagonal part of the NLO evolution operator is given as with β 1 specified below (41), and where the NLO anomalous dimensions γ j can be read from [78,79] with the convention change specified in [57] (Eq.( 134)).The non-diagonal part of the NLO evolution operator, which is the part that leads to mixing of conformal moments, can be written in the form where the matrices d jl and g jl can be read off [79] with the same convention change as applied to the diagonal anomalous dimension matrices.The expressions for the flavor nonsinglet sector can be derived from the ones provided above by reducing the matrix-valued quantities to scalar values associated with quark contributions, i.e., ΣΣ.Thus, one gets the nonsinglet evolution operator E(µ, µ 0 ; ξ) by employing the replacements outlined in [57], Eq. ( 144).We note that to NLO accuracy considered here The evolution of the DA is governed by the same evolution operators as for GPDs, with skewness ξ = 1 (ERBL evolution [60,59]), while properly accounting for the parity and charge quantum numbers.Therefore, for the vector mesons under consideration in this work, the evolution is governed by the nonsinglet intrinsic parity even evolution operator E. with The evolution of the intrinsic parity odd (axial) GPDs and of the pseudoscalar meson DA follow analogously, and the corresponding anomalous dimensions are specified in [57].

Evolution implementation
Finally, the evolution of GPDs and consequent convolution with the hard-scattering amplitude can be equivalently interpreted as the evolution of the hard-scattering amplitude 8 and a consequent convolution with GPDs T l (µ) This is practical for numerical evaluation, where evolved hard-scattering amplitudes can be stored in computer memory and called up during fitting of GPDs.Summation over conformal moment l in (36) can also be replaced by Mellin-Barnes integration, and after some reshuffling of the terms one gets with j ⊗ defined in (18a), and In this work we use c ′ = −0.25.The form specified in (38) is applicable to both intrinsic parity even and odd GPDs, regardless of their signatures.This is because the integral over l stems from the generic sum Σ ∞ l=0,even .Naturally, the hard-scattering amplitudes differ, and the operators 8 We note that the term 'evolution' in this context does not involve the resummation of terms are resummed.This should be viewed as an evaluation technique and the hard-scattering amplitude retains a residual dependence on the factorization scale µ

2
. For a more detailed discussion of factorization dependence, see [80].
include the corresponding vector or axial-vector anomalous dimensions, as elaborated in [57].The nonsinglet expression is analogous.
The distribution amplitude is at some input scale µ 0 often represented in terms of the finite number of conformal moments (11).The evolved DA (34) is then to LO also presented by a finite sum.However, NLO evolution introduces mixing between conformal moments and an infinite sum emerges.In practical applications, this infinite sum is 'truncated' assuming the suppression of higher conformal moments.When the analytical form of DA conformal moments is known, Mellin-Barnes integration can be employed, resulting in an expression similar to (38).For DVMP, it is advantageous to 'precalculate' both DA and GPD evolutions with the hard-scattering amplitude where l ⊗ corresponds to (38), while m ⊗ can be a summation over m or an analogue of (38).Often one takes µ 0 = µ ′ 0 .Finally, the CFF and TFF expressions suitable for practical use are given by ⊗ is given by (18a), while k ⊗ is a summation (18b) or (18a) as well.Moreover, expressions (40) allow us to quantify the influence of quarks and gluons present at initial scale µ 0 .Note that due to the fact that E (26) is nondiagonal in the flavor space (Σ, G), the components Σ T and G T represent the mixture of Σ T and G T .Consequently, in contrast to the gluon contribution at scale µ obtained using ) is different than zero even at LO.As is well-known in renormalization group theory, the dependence on the evolution scale is primarily manifested through the scale-dependent behavior of the strong coupling constant.In numerical evaluations, we utilize the equation where a s = α s /(4π), β 0 is given in (23), and . At LO we employ the analytical solution considering only the first term in (41), while at NLO we utilize the numerical fourth-order Runge-Kutta integration.At the initial scale µ 2 = 2.5 GeV 2 , we take the values 0.0606 (LO) or 0.0518 (NLO) for α s (µ)/(2π), and set n f = 4, consistent with the kinematical range of interest here.

Understanding contributions: simple lessons and insights
When examining the expressions presented in this section within the realm of phenomenology, it is crucial to note the following key points: • The gluon component of the NLO hard-scattering DVCS amplitude in ( 19) is negative.Consequently, gluonic contributions at NLO strongly suppress the quark contributions, as later illustrated in Fig. 11.The most dominant contribution in (19d) arises from the term: which in x-space corresponds to the most singular part of the amplitude proportional to The DIS structure function F 1 corresponds to the imaginary part of H in the forward limit 9 , and it's worth noting that precisely (42) does not contribute to its gluon Wilson coefficient ([57], Eq. ( 105)).• In DVMP, gluons appear already at LO (22b), and owing to their high density in the relevant kinematic regime of small x B , which is already established from the study of the DIS process, they are responsible for the dominant part of the DVMP amplitude.This is evident in Fig. 12. • The j = 0 pole, a well-known feature present in the LO gluon anomalous dimensions GΣ γ (0), j and GG γ (0), j (28), plays a crucial role in driving strong gluon evolution at small x B .This pole represents the j-space signature of the interplay between the Bjorken Q 2 → ∞ and the high-energy limit 1/x B → ∞.The impact of this gluon evolution is evident in Figs. 10 and 11, where the DIS structure function and DVCS CFF exhibit an increase with Q 2 .Additionally, Fig. 12 illustrates the flat behavior of the DVMP TFF in the relevant energy region.This flat behavior, as a consequence, facilitates the correct scaling behavior demonstrated in Fig. 8.
• The j = 0 pole is not only present in the intrinsic parity even evolution operator E (26) but also surfaces in the NLO DVMP contributions for the pure singlet moment pS c (1) jk (24a) and the gluon contribution G c (1,A) jk ( [58], Eq. (4.53a)).These contributions are parametrized using GΣ γ (0), j and GG γ (0), j , multiplying the factorization logarithm ln Q 2 /µ 2 GPD and constant terms.Similar contributions are expected for DVCS at NNLO.

Model for GPDs and DA
To complete our theoretical framework, it is essential to specify a parameterization of GPDs that facilitates fitting to experimental measurements.We adopt a model wherein conformal moments of GPDs, as defined in (10), are characterized by an expansion in t-channel SO(3) partial waves [55,57,81,82].The foundational concepts for this approach were laid out in [55,57] and subsequently refined and applied, particularly in the context of small-x B phenomenology, in [81]  10 .In [82], the fundamental 9 This observation was exploited in [57] where NNLO DIS Wilson coefficients contributing to F 1 were utilized to access NNLO DVCS in the special conformal scheme.10 The model for GPD H conformal moments for small x B was utilized in [81], alongside the separate determination of Compton form factors in the valence region.It's worth noting that the KM-model nomenclature is used interchangeably for both results, even though they represent distinct approaches to DVCS phenomenology.steps connecting this j-space modeling with more conventional x-space GPDs were outlined.The former approach offers several advantages, including a clear connection with Mellin moments and lattice calculations, while the latter is more intuitively clear.It's crucial to emphasize that most of the hard-exclusive phenomenology and GPD modeling has traditionally been carried out in x-space, see recent examples in [83,84].Apart from the previously listed applications of the CPaW formalism to DIS, DVCS and DVMP at NLO and beyond, recently the conformal moment representation has been employed at LO [85,86], aiming to establish a parallel global analysis program encompassing PDFs, form factors, and DVCS measurements, with a foreseeable extension to lattice data.
We take generic expression [81] where J is the angular momentum in the t channel, and da J (ξ) are the corresponding Wigner functions.As in [81], we work in the approximation da J (ξ) ≈ 1, which is good enough for small-x B kinematics.In accordance with the fact that the forward limit of GPDs is equal to ordinary PDFs, e.g. for quark q, H q (x, 0, 0) = θ(x)q q (x) − θ(−x)q q (−x) , the amplitude H j,J=j+1 (t) of the leading partial wave is for t = 0 equal to the Mellin moment of the corresponding PDF In principle, the values q a (x) could be taken from the results of one of the collaborations specialized in the extraction of PDFs from experimental measurements.However, state-of-the-art PDF extraction usually employs variable flavor number procedures with sophisticated matching at heavy-quark production thresholds.Thus, naive usage of numerical values for q a (x) from those studies in our simplified fixedflavor-number framework would lead to inconsistencies.Therefore, as in [57,81], we determine the Mellin moments of PDFs ourselves, by fitting a simple ansatz to DIS measurements, where the counting rules determine β sea = 8 and β G = 6.This ansatz corresponds to a standard simple parameterization in x-space The normalization in (47) and ( 48) is chosen so that N a corresponds to the average longitudinal momentum fraction of parton a.These parameters are constrained by the summation rule and where, in accordance with results of DIS studies, we take N val = 0.4.Thus, we use only three parameters to fit small-x B DIS measurements: We complete the model of the leading partial wave by adding the dependence on t, and that, firstly, by completing the full Regge trajectory and, secondly, by adding the residual dependence on t via the dipole impact factor, so that the complete expression for the leading partial wave is with q a j given in (47).It turns out that the data we work with cannot distinguish the differences in tand j-dependencies of individual partial waves, so for subleading waves we assume that they are simply proportional to the leading wave ( 52), and at the end we take into consideration only two subleading waves, so our complete model is given by [57,81] where normalizations of the subleading partial waves s a 2 and s a 4 are additional parameters of the model.Thus, the total set of parameters, in addition to the three parameters from ( 50) is given by: These additional eight parameters are determined by fitting to measured DVCS and DVMP data.We also mention that the practically identical GPD model was also used in the first multichannel DVCS + DVMP fit in [53].The main difference is that the subleading partial waves were not taken as proportional to the leading one, but, for example, the Regge parameters α 0 and α ′ were different for each partial wave.In addition, the overall normalizations of the measured DVMP cross sections were considered as fitting parameters.A larger number of parameters makes the model in [53] more flexible, but the model used in this paper proved itself flexible enough for our needs.
The distribution amplitude φ ρ 0 (u) is relatively poorly known and for this study we decided to stick to its asymptotic form, so the sum over k in (16) has only one term and NLO evolution, being tiny for this form, is neglected.More advanced studies could also treat DA as an unknown function whose lower Gegenbauer moments could be simultaneously fitted to experimental measurements, or, at this moment probably more convenient, those lower moments could be taken from lattice QCD results [87].58) for the ρ 0 production obtained by separate fits of ansatz (60) to H1 measurements [88] (left), and to ZEUS measurements [89] (right) are plotted as the red solid line and the corresponding symmetric red uncertainty band.The function R(Q 2 ) determined by W -independent analysis [54] using the ansatz ( 59) is shown for comparison as the green dashed line and the corresponding symmetric green uncertainty band.

Experimental data and L/T separation
Since we focus on the high energy regime (small x B ) and large values of virtuality Q 2 , the only experimental data that we use are those of the H1 and ZEUS collaboration realized using the HERA collider.In particular, we use the results of • H1 measurements of F 2 (x B , Q 2 ) DIS structure function from [90]  11  • H1 and ZEUS measurements of the DVCS cross section σ γ ≡ σ(γ * p → γ p) from [92,93,94,95], • H1 and ZEUS measurements of the DVMP cross section for production of the ρ 0 meson σ(γ * p → ρ 0 p) from [89,88].
In order to be in the regime where perturbative QCD can be used with confidence in the leading twist, for DVCS data, we impose a cut Q 2 > 5 GeV 2 , and for DVMP, which is more problematic due to two hadrons in the final state, we use a data cut Q 2 > 10 GeV 2 .For the possible importance of higher-twist contributions see recent results [96] at lower Q 2 .We could also straightforwardly include measurements of the production of ϕ mesons.Neglecting the contribution of valence quarks leads to a simple SU (3) flavor relation between the cross sections for the production of ρ 0 and ϕ mesons: where the meson decay constants are f ρ 0 = 0.209 GeV and f ϕ = 0.221 GeV.However, deviations from this relation are visible in the experimental data in Fig. 9, where the cross sections for ϕ, scaled according to (57), are consistently smaller than those for ρ 0 .To avoid tensions that are entirely an artifact of our simplified model, in this paper we do not use the data for ϕ (nor for ω), but only those for ρ 0 that have better statistics.Neglected influence of valence quark contributions results in a small reparametrization of the sea-parton GPDs.One of the main goals of this study is to establish the usability 11 This data has been superseded by the combined H1 and ZEUS data [91].However, these more recent measurements are reported in terms of cross-sections that combine all DIS structure functions F i , whereas older measurements report F 2 so we can use them directly within our framework which is limited to describing only the F 2 part of the DIS cross section.
of the collinear twist-2 approach to describe the longitudinal DVMP process 12 , and that can already be judged on the basis of the ρ0 meson production alone.
In contrast to the DVCS situation where the dominance of the twist-2 amplitude is a safer assumption and the measured data can be directly compared with the theoretical prediction, it is known that the measured DVMP cross section has a large higher-twist contribution from the exchange of transversely polarized virtual photon.Consequently, for comparison with the twist-2 theory presented in Sec. 2 it is necessary to separate the contributions of longitudinal and transverse photons (L/T separation) and compare the theory with measurements of the longitudinal cross section σ ρ 0 L = σ(γ * L p → ρ 0 p).Direct observation of the polarization of the virtual photon is of course not possible, so L/T separation of the DVMP process is done indirectly, typically by measuring the polarization of the meson in the final state and using the assumption of s-channel helicity conservation (SCHC), which implies that the polarization of the photon follows that of the meson, σ The SCHC assumption can also be experimentally verified by measuring and comparing individual spin-density matrices elements (SDME), which additionally enables the improvement of the L/T separation procedure itself.Using such an approach, both HERA collaborations determined the ratio and studied its dependence on the kinematic variables of the process.For example, the results in [88] show • lack of apparent dependence of R on W within uncertainties, and • certain indications of dependence of R on t, visible for higher values of L are given in a binned form only in Q 2 , and for a quality extraction of GPDs it is important to know the dependence on other kinematic variables, W and t.Therefore, in this work we mainly rely on measurements of the quantity R and, using it together with measurements of the total cross section σ ρ 0 , we determine the values of σ ρ 0 L which we use in our fits.Since R is mostly not given for the same kinematic points as σ ρ 0 , measured values for R cannot be used directly, but we interpolate them by fitting, thus obtaining R as a continuous function of kinematic variables.
In [52], motivated by the fact that only the dependence on Q 2 is unquestionable, an ansatz for the function R is proposed of the form where m ρ 0 = 0.776 GeV, and p and a are the fitting parameters.In [52,53,54] these parameters were determined by fitting to the H1 and ZEUS values for R. In the present work, we have reexamined the dependence of R on kinematic variables W and t.First, as our interest is focused on data having large Q 2 , we observe that precisely at these values the uncertainty of R is greatest, and thus it is the least 12 From this point forward, for the sake of readability, we will refer to the process as DVMP, although the technically accurate term is DVρ clear whether or not there is a dependence on W (cf. Fig. 39 in [88]).Also, precisely at higher values of Q 2 both H1 and ZEUS measurements, which are at different W , show more significant divergence in values for R, see Fig. 3. Therefore, for the purposes of this work, we made a new determination of the function R, where we extended the ansatz (59) by a factor that also describes the dependence on where b is the new parameter.For the numerical efficiency of the fitting procedure, W is not squared in the denominator, which makes the parameter b dimensional.Also, because of tension between the measurements, we determined the function R separately for H1 data, and for ZEUS data: Thus obtained functions R(W, Q 2 ) H1 and R(W, Q 2 ) ZEUS are used in this work to determine σ ρ 0 L from H1 and ZEUS measurements for σ ρ 0 , respectively.They are shown in Fig. 3.The question of dependence of R on the variable t also arises.The situation is again problematic for large values of Q 2 , even more than in the case of the dependence on W . Precisely at higher values of Q 2 a dependence on t begins to appear, and additionally, measurements of the production of light vector mesons at slightly higher x B measured by the COMPASS collaboration (cf.Fig. 11b in [97]) unambiguously indicate that R also depends on t.On the other hand, we cannot determine this functional dependence because measurements of R(t) are available only for smaller values of Q 2 , below the limit Q 2 = 10 GeV 2 that we set in this paper.In this unsettled situation, we decided to take a conservative stance and we do not use DVMP differential cross section measurements dσ ρ 0 L /dt, but only measurements of the cross section integrated over t.
As discussed above, one of the important features of the DVMP cross section is its scaling with Q 2 .Being able to correctly reproduce the scaling is an important test of a theoretical model.The scaling of the total cross section σ ρ 0 was usually considered in the literature, but here we consider the scaling of its longitudinal component σ We will see later that NLO models (unlike LO) have no problem reproducing this value.However, what is usually considered "the scaling" is the scaling with Q 2 for fixed x B .Four measurements of the H1 collaboration at x B = 0.0018 are shown in Fig. 8 and fitting the power function Such a considerable deviation of scaling from the asymptotic value w = 6 is considered to be difficult to describe within the collinear perturbative QCD approach.However, we find that our LO and NLO models are capable of reproducing this effective behavior in the measured data region.[90], while the blue squares are combined measurements of the H1 and ZEUS collaborations [91].

Fits
In order to understand the influence of NLO corrections of different processes, we use six different GPD fits, listed in Tab. 1.All of them were obtained by fitting the ansatz described in Sec. 3, at LO or NLO of perturbative QCD.Every model is first fitted to the DIS F 2 (x B , Q 2 ) data of the H1 collaboration [90].After that, we fix the parameters (50) that affect the form of GPDs in the forward limit, and the remaining parameters (55) of the three LO and three NLO models are fitted to only DVCS data, only DVMP data and to the total DVCS+DVMP set, as specified in Tab. 1.For the sake of brevity, in the names of all models we suppressed the DIS label, so, for example, model NLO-DVCS-DVMP should have been called NLO-DIS-DVCS-DVMP.
The initial values of the parameters are listed in the first row of Tab. 2. They were chosen very generically, using nevertheless to some extent the knowledge of the corresponding Regge trajectories, so initial intercepts are set to 1, and initial slopes to 0.15 for both quarks and gluons.Also, the first subleading quark partial wave normalization is set to a negative initial value s sea 2 = −0.2based on earlier LO DVCS studies [81] which preferred it so.Some parameters are restricted as stated in the second row of Tab. 2. The squares of masses m 2 a and Regge slopes α ′a are for physical reasons limited to positive values of the order of one.Negative values of α ′a enable solutions with discontinuous behavior of the form factor as a function of t, which should be avoided.Partial-wave normalizations of the quarks are constrained to achieve a natural hierarchy 1 ≫ s sea 2 ≫ s  [92,93] and ZEUS [94,95] measurements of the dependence of the DVCS cross section on t, W and Q 2 compared to LO (thin) and NLO (thick) models which are, besides DIS, fitted only to these DVCS measurements (blue dashed), and also additionally to DVMP measurements (red solid).The three H1 lines on left panels correspond, from top down, to Q 2 = 8, 15.5, and 25 GeV 2 , respectively.
the gluon sector prevented a good fit, so gluon partial waves are left more flexible than quark ones.Future studies should show whether the final relatively large values of formally subleading gluon partial waves s G 2 and s G 4 , which were also observed in earlier studies that also included data from experiments with a fixed target [98], are a signal of a problem with our model.The resulting values of these and other parameters for the best of the six models, NLO-DVCS-DVMP, are displayed in the third row of Tab. 2, and the uncertainties are given in the fourth row, where everything is the result of a standard least squares fit with the software routine MINUIT [99,100].The most correlated pairs of parameters for that model are • α ′sea and m 2 sea with correlation 0.937 which shows that it is not possible to clearly distinguish the t dependence associated with x ("shrinkage" effect, controlled by parameter α ′ ) from the residual    Regarding the goodness of the six fits, Tab. 3 shows the obtained values of χ 2 /n pts where we observe as follows.
• Models describe data well only if they are fitted to them (see Tab. 1).For other data not included in the fitting, the obtained values χ 2 /n pts are in the range 10 to 30 000, and these cases are marked as ≫ 1 in Tab. 3. • Only the NLO model NLO-DVCS-DVMP describes all the data well and therefore we consider it the best model presented in this paper.• LO models can describe DIS and DVCS data well, as shown already in [57], but they cannot describe DVMP well, with χ 2 /n pts > 3 (although the total χ 2 /n pts = 1.4 of the model LO-DVCS-DVMP for the total data set looks satisfactory).
In Tab. 3, we give the values for χ 2 divided by the number of data points and not, as is common, by the number of degrees of freedom.This is because for particular subsets of the total dataset, it is not possible to distinguish the corresponding number of degrees of freedom.Let us give therefore for the We repeated the entire analysis on the data where L/T separation was performed using the universal function R(W, Q 2 ), i.e., not using separate R(W, Q 2 ) H1 and R(W, Q 2 ) ZEUS .In that case, the conclusions are unchanged, the plots are barely distinguishable, and the values χ 2 /n d.o.f. also change only a little.For example, the value χ 2 /n d.o.f. of the best model NLO-DVCS-DVMP becomes 1.3 instead of 1.2.
As stated in Sec. 3, instead of setting a forward limit of our GPD model to be equal to PDF functions taken from one of the modern PDF sets, we also fit models to DIS data, i.e, F 2 structure function whose form in terms of the Mellin-Barnes integral is given in [57], Eq. (198).The fit result is shown in Fig. 4 and in the first line of Tab. 3, where it can be seen that the description of these data is perfect.This may seem a relatively trivial test of our model, but given the subtleties when choosing the number of quark flavors and crossing the corresponding thresholds, we believe that it has not been given enough attention in other attempts to extract GPDs in the literature.In Fig. 4, we also present more precise combined H1 and ZEUS data, which we have not utilized in our fits for reasons outlined in footnote 11.In the illustrated region, where due to the small value of electron inelasticity the F 2 structure function dominates the DIS cross-section, even this highly accurate data is reasonably described by our models.As can be seen in the second row of Tab. 3 and in Fig. 5, the description of DVCS data is good in both LO and NLO.There are certain problems visible when describing ZEUS data in the third panel of Fig. 5, but it is for a low value Q 2 = 3.2 GeV 2 and these data were not even included in the fit.
In Figs. 6 and 7, and in the third row of Tab. 3 results are given related to one of the central questions that we deal with in this work: Is it possible to describe the cross sections for DVMP in the canonical collinear GPD approach?It can be seen immediately, from the problem in the description of the Q 2 dependence for higher values of Q 2 > 30 GeV 2 , especially for the H1 data in Fig. 6 on the left, and from the problem in the description of the dependence on W , especially for the ZEUS data in Fig. 7 on the right, and from large values χ 2 /n pts > 3, that at the LO level the description is not good.However, at NLO, the situation improves noticeably.True, the right panels of Figs. 6 and 7 suggest that the W dependence of the NLO model is not steep enough (interestingly, for the LO models it is too steep), but the description is mostly satisfactory, the behavior for large values of Q 2 is correct and χ 2 /n pts values are acceptable.It is worth focusing additionally on the problem of the Q 2 dependence which was often highlighted as the main problem of a purely collinear approach, because in that approach, the canonical scaling of the longitudinal cross section at a fixed value x B is of the form ∼ 1/Q 6 , cf. ( 2) and ( 4), whereas experimental data on the total cross section reportedly suggest a milder dependence ∼ 1/Q 4 .There, of course, the first problem is to distinguish the contributions of the purely longitudinal cross section, which we dealt with in section 4. If we accept the L/T separation function proposed there, in fact the longitudinal component of the cross section also shows the behavior σ however, it should be noted that experimental data for fixed x B are relatively scarce.They are shown in Fig. 8 where we see that within the experimental errors both LO-DVMP and NLO-DVMP show satisfactory effective Q 2 scaling in the measured range, which is a consequence of perturbative Q 2 logarithms.Only for much larger values of Q 2 the models obey the canonical ∼ 1/Q 6 behavior.We can better assess the quality of the description of the Q 2 dependence of the cross section if we study it for fixed W , for which there is more abundant experimental information.As discussed in section 4, the longitudinal cross section for fixed W shows the effective behavior σ ρ 0 L ∝ Q −5 , so in Fig. 9 we show the rescaled Q 5 σ ρ 0 L data compared to four of our models fitted to DVMP measurements.It can be seen that only NLO

Quark vs gluons at LO and NLO
If we grant that the description of the experimental data presented in the previous section is satisfying, let us see what we can learn about the quark-gluon structure of the proton responsible for that description.As the DIS and DVCS processes are related in the sense that they use a photon probe that couples to the same flavor singlet combination of quarks and gluons, separation of the quark and gluon structure using only those two processes is necessarily indirect.Adding the DVMP process to the analysis brings a new, meson probe, with, in principle, a different flavor structure, and at the same time a more direct approach to the gluonic GPD that already contributes to the LO amplitude, see Fig. 2. Therefore, DVMP rightly stands out as the key to separating individual flavors and gluonic GPD in the proton.
However, the influence of the gluon GPD on the DVCS cross section should not be underestimated either.This influence is reflected through two effects.The first of them is strong evolution of gluons at high energies, which is already present in DIS observables and is visible in Fig. 10, where the contribution of gluons to the total DIS F 2 (x B , Q 2 ) at x B = 0.001 is of the same sign and becomes equal to the quark one already for Q 2 ∼ 50 GeV 2 .Another effect is the significant contribution of gluons to the total NLO non-forward DVCS amplitude of hard scattering.It is negative, i.e. of the opposite sign to the quark contribution, and despite the formal α s /(2π) suppression is equally large, which was already pointed out earlier [101,102].This can be clearly seen in the second panel of the Fig. 11, where already at the initial scale µ 2 0 = 4 GeV 2 a large negative contribution of gluons would completely cancel out the quark LO contribution, so the NLO quark contribution must be twice as large as at LO so that the final CFF Im H could be in accordance with the experiment.This second effect is not present in DIS, which can be clearly seen in Fig. 10 where the displacement of gluons on the right NLO panel is positive and very small, so that the LO and NLO cases are barely distinguishable.It should be noted that what corresponds to the DVCS CFF function H in deep-inelastic scattering is not the structure function F 2 but F 1 , as discussed in Sec.2.4.Unlike F 2 , where the NLO correction is positive and small, see Fig. 10, for F 1 , the relative correction is negative, therefore in sign agreement with DVCS.It is also larger than for F 2 , but still turns out to be about four times smaller than for the CFF H.It is also important to highlight that the quark and gluon contributions in Figs.10-12 were separated according to (40) i.e., they reveal the influence of quarks and gluons present at the initial scale µ 0 .Thus although there is no true gluon contribution for DVCS at LO, the influence of initial gluons on the quark GPD, which they exhibit through evolution, and thus on the CFF, is revealed.As can be seen in Fig. 12, gluons equally dominate the DVMP amplitude at both LO and NLO and the main NLO effect is a relative suppression of the real part of the amplitude.Some models of the longitudinal production of light vector mesons have predicted that the gluon GPD alone is almost sufficient to describe the amplitude of this process at high energy and our results approximately confirm this.
Such different behavior of GPDs and corresponding PDFs at LO and NLO can be succinctly described using the so-called skewness ratio of GPDs and PDFs: which for small values of x practically does not depend on x.To derive x-space GPDs from those modeled in j-space, one must perform an inversion of the Mellin-Barnes integral.This process is discussed in detail in [55,82].An example of a GPD on the cross-over line is provided in [81], specifically in Eq. (43).In papers [103,104] it was proposed that a GPD, up to the dependence on t, be completely determined by the value of the corresponding PDF and the conformal values of the skewness ratio: Obviously, the real situation is more complex and this quantity varies a lot from LO to NLO.An additional variation is brought by a completely different interplay of quarks and gluons in the DVMP amplitude.We can gain additional insight by observing the skewness ratio (66) of all six models, shown in Fig. 13.The thin lines show the three LO models and it can be seen that the quark and also the gluon GPDs fitted to DVMP, DVCS or combined DVMP+DVCS data are drastically different.However, almost universal GPD functions emerge at NLO for the combined DVCS+DVMP fit and for fits where DVCS or DVMP data are turned off (thick lines).The NLO skewness ratios are consistently higher than conformal ones (67) and are asymptotically slightly above the value of 2 for quarks and about 1.5 for gluons.We see that the gluon GPD from the combined NLO DVCS+DVMP fit is closer to that from the DVMP fit, which is probably a reflection of the fact that DVMP has more influence on the extraction of gluonic GPD than DVCS.
Considering the diverse roles of quarks and gluons, successfully achieving a comprehensive simultaneous description of all three observed processes using unique universal GPD functions, exemplified by

Conclusions
In this work we have revisited the NLO corrections to DVMP.The main ingredients and advantages of the CPaW formalism [55,56,57,58] have been provided along with revised expressions for deeply virtual production of longitudinally polarized vector mesons.The impact of NLO corrections has been analyzed through global fits to DIS, DVCS and DVV 0 L P data, in particular, ρ 0 meson data, at small-x B .We conclude that a reliable and consistent description of the longitudinal component of the DVMP cross sections at high energies and for Q 2 > 10 GeV 2 is possible in the standard collinear pQCD approach, only if NLO corrections are included.A simultaneous description of DIS, DVCS and DVMP processes becomes feasible at the NLO level, thus revealing the proton structure described by universal GPD functions.

Figure 1 :
Figure 1: Kinematics of the DVCS and DVMP processes.

Figure 2 :
Figure 2: Generic diagrams contributing to the factorized DVCS (top) and DVMP (bottom) process at LO.In the DVMP case, the gluon GPD at input scale contributes already at LO (lower right).

Figure 3 :
Figure 3: R(W, Q 2 ) functions(58) for the ρ 0 production obtained by separate fits of ansatz(60) toH1 measurements  [88]  (left), and to ZEUS measurements[89] (right) are plotted as the red solid line and the corresponding symmetric red uncertainty band.The function R(Q 2 ) determined by W -independent analysis[54] using the ansatz (59) is shown for comparison as the green dashed line and the corresponding symmetric green uncertainty band.

Furthermore
, in[88], the H1 collaboration specifically determined values for σ

ρ 0 L . A fit of the function σ ρ 0 L
∝ Q −w to the data from the left panel of Fig.6gives w = 5.1 ± 0.1 (fixed W ) .

Figure 4 :
Figure 4: Description of the DIS structure function F 2 (x B , Q 2 ) at LO (red, thin) and NLO (red, thick) by models obtained by a multichannel fit to DIS, DVCS and DVMP data.The black points are measurements of the H1 collaboration[90], while the blue squares are combined measurements of the H1 and ZEUS collaborations[91].

sea 4 .Figure 5 :
Figure 5: Description of H1[92, 93]  and ZEUS[94,95] measurements of the dependence of the DVCS cross section on t, W and Q 2 compared to LO (thin) and NLO (thick) models which are, besides DIS, fitted only to these DVCS measurements (blue dashed), and also additionally to DVMP measurements (red solid).The three H1 lines on left panels correspond, from top down, to Q 2 = 8, 15.5, and 25 GeV 2 , respectively.

Figure 6 :
Figure 6: Description of H1 DVMP measurements [88] by models fitted to H1 and ZEUS DVMP data (green, dot-dashed lines) and additionally to DVCS data (red, solid lines).It can be seen that, unlike the NLO models (thick lines), the LO models (thin lines) have problems describing the data at larger values of Q 2 .

Figure 7 :G 4
Figure 7: Description of ZEUS DVMP measurements[89].Meaning of lines is the same as in Fig.6.

Figure 8 :
Figure 8: Longitudinal cross sections for production of ρ 0 mesons [88], at fixed x B = 0.0018, together with LO-DVMP (thin green dot-dashed line) and NLO-DVMP (thick green dot-dashed line) models fitted to DVMP data.Both models are able to reasonably reproduce the effective experimental Q 4 scaling in data region.

Figure 9 : 0 L ∼ 1 /Q 5 ,
Figure 9: Longitudinal cross sections for production of ρ 0 mesons (black circles) [88], plotted according to the approximate effective experimental fixed-W scaling σ ρ 0L ∼ 1/Q 5 , together with LO (thin lines) and NLO (thick lines) models fitted to DVMP only (green dot-dashed) and DVMP+DVCS data (red solid lines).One observes that only NLO models are able to reproduce the experimental large Q 2 scaling.We also display the SU(3) flavor rescaled ϕ production data (blue triangles) which was not used in this work.

Figure 10 :
Figure 10: Contributions of the quark (blue dot-dashed) and gluon (red dashed) PDF to the total DIS structure function F 2 (black solid) at LO (left) and NLO (right), for x B = 0.001.

Figure 11 :
Figure 11: Contribution quark (blue dot-dashed) and gluon (red dashed) GPD to the total DVCS CFF (black solid) in models LO-DVCS-DVMP (left) and NLO-DVCS-DVMP (right), in dependence on Q 2 for x B = 0.001 and t = 0. Imaginary (top) and real (bottom) part of the leading CFF H are displayed.

Figure 12 :
Figure 12: Contributions of the quark (blue dot-dashed) and gluon (red dashed) GPD to the total DVMP TFF (black solid) in models LO-DVCS-DVMP (left) and NLO-DVCS-DVMP (right) in dependence on Q 2 for x B = 0.001 and t = 0. Imaginary (top) and real (bottom) part of the leading TFF H ρ 0 L are displayed.The steep behavior below Q 2 = 10 GeV 2 is a reflection of strong evolution effects.It lies below our data range and is likely outside the scope of validity for our twist-2 collinear approach.

Figure 13 :
Figure13: Skewness ratio(66) for the quark (left) and gluon (right) GPD H, for six models from Tab. 1. Depending on the processes used, the LO models (thin lines) show large mutual differences of the resulting GPDs.The dependence of NLO GPDs (thick lines) on the process is significantly smaller.

Table 1 :
List of models used in this work, pQCD order used, and datasets to which the model was fitted.References to used DIS, DVCS and DVMP experimental data are given in Tab.3.

Table 2 :
The first two rows give the initial values of the parameters and limits on their ranges, if set.The final fitted values and their uncertainties (one standard deviation) of the best NLO-DVCS-DVMP model are given in the last two rows.The corresponding values of χ 2 /n pts are given in the last column of Tab. 3. The parameters correspond to the scale µ 2 0 = 4 GeV 2 .

Table 3 :
Values of χ 2 /n pts for each LO or NLO model (columns) for the total DIS + DVCS + DVMP dataset and for subsets corresponding to different processes (rows).(The values denoted by ≫ 1 are greater than 10.) 2/n d.o.f. for models fitted to the maximal dataset DIS+DVCS+DVMP: