Exclusive photoproduction of $D$-meson pairs with large invariant mass

In this paper we analyze the exclusive photoproduction of the $D$-meson pairs with large invariant mass. We perform evaluations in the collinear factorization framework and in the leading order of the strong coupling $\alpha_{s}$, expressing the cross-section in terms of generalized parton distributions (GPDs) of different parton flavors in the proton. We focus on the photoproduction of the pseudoscalar-vector pairs, like e.g. $D^{\pm}D^{*\mp}$, $D^{0}\overline{D}^{*0}$, $D_{s}^{+}D_{s}^{*-}$, which gets the dominant contribution from the chiral even GPDs of the target, and estimate the cross-section in the kinematics of the future Electron Ion Collider (EIC). In all channels the amplitude of the process obtains comparable contributions from gluons and only one of the light quark flavors. This finding signals that the process potentially could be used to single out the contributions of the individual chiral even GPDs of light flavors. We found that the process is mostly sensitive to the behavior of GPDs in the so-called Efremov-Radyushkin-Brodsky-Lepage (ERBL) region. Numerically, the cross-section of the process is sufficiently large for experimental studies and thus can be used as a complementary probe for studies of the partons GPDs.


I. INTRODUCTION
In the last three decades the Generalized Parton Distributions (GPDs) of the nucleon turned into a standardized tool to encode information about the nonperturbative interactions of individual partons in the hadronic target [1][2][3][4][5][6], and for this reason have been in the center of theoretical and experimental studies.The GPDs allow to understand the contributions of different parton flavors to various observables which characterize the hadronic target.At present it is not possible to evaluate the GPDs directly from first principles, and for this reason studies of these objects rely on phenomenological extractions from experimental data or results of lattice simulations [7][8][9][10][11].However, the existing lattice studies, due to technical challenges, at present mostly focus on the special zero-skewedness (ξ = 0) limit and studies of some moments of GPDs, whereas phenomenological extractions suffer from various uncertainties, even for the cleanest and best understood channels [12].This motivates the search for new processes which could be used for extractions of the GPDs [13][14][15][16].
The expected high-luminosity experiments at the future Electron Ion Collider stimulated interest in various channels, which were previously disregarded due to the smallness of the cross-section.For studies of GPDs, a special interest is present in exclusive 2 → 3 processes, which have already been analyzed in the literature [17][18][19][20][21][22][23][24][25][26][27][28][29][30].Due to their different kinematic structure, these new probes could complement existing studies and provide new independent constraints on existing phenomenological models of GPDs.For amplitudes of such processes, the factorization theorem has been proven in the kinematics case when all the produced hadrons are well-separated kinematically, i.e. the pairwise invariant masses (≈relative velocities of the produced hadrons) are sufficiently large to avoid soft final-state interactions [31,32].Most of these studies focused on the production of pairs of light mesons and photons, and it has been discussed in detail how these novel channels could help to access new information about the GPDs [27,28,[33][34][35].
Potentially the production of the heavier meson pairs, like D-mesons and quarkonia, might be also used for the same purpose: in the kinematics where factorization theorems are applicable, the overall cross-section suppression by the heavy quark mass is comparable to the suppression by a large invariant mass.However, for these mesons the theoretical treatment should be adjusted, since the heavy quark masses cannot be disregarded, and should be considered as one of the heavy scales in the problem.This breaks the conventional twist suppression used for light quarks, and leads to new probes of the GPDs (compared to light meson production channels with the same quantum numbers).Furthermore, due to lack of substantial contributions from intrinsic heavy flavors in the proton, the heavy meson production channels might be used to disentangle the flavor structure of the GPDs, thus avoiding the usual superposition of GPDs of all light flavors.For heavy quarkonia pair production, this allows to single out the contribution of the gluon GPDs [36][37][38][39][40][41][42][43].For the D-mesons, the cross-sections get an additional contribution from one of the light flavors [13][14][15][16], which potentially allows to test individually the GPDs of light flavors.In this paper we will focus on the production of scalar-vector mesons pairs, like D + D * − , D 0 D * 0 and D + s D * − s which has not been discussed so far in the literature and might be used as a complementary probe of the target GPDs.This choice of quantum numbers also allows to avoid the contribution of the photon-photon fusion mechanism, which has been discussed previously in [44], and the contribution of the poorly known chiral odd transversity GPDs of the target, discussed in [27,28].We will analyze this process in the conventional collinear factorization approach, although maintaining the heavy quark mass as a hard scale, and will provide numerical predictions in the kinematics of low-and middle-energy electron-proton collisions at the forthcoming Electron Ion Collider (EIC) [45][46][47][48].
The paper is structured as follows.In the next Section II we introduce the framework and provide analytical expressions for the amplitude and cross-section of the process.In Section III we estimate numerically the crosssections, using publicly available parametrizations of the proton GPDs and D-meson distribution amplitudes.Finally, in Section IV we draw conclusions.

II. EXCLUSIVE PHOTOPRODUCTION OF MESON PAIRS
Below, in Section II A, we define the kinematic variables, discuss their typical ranges in EIC kinematics and introduce the light-cone decomposition for momenta of all particles, which will be used later.In the next Section II B, we evaluate analytically the amplitude of the process in the collinear factorization framework.

A. Kinematics of the process
For our evaluations we will use the photon-proton collision frame, in which the photon and proton move along the axis z, so the light-cone decomposition of their momenta are given by where the shorthand notation q stands for the momentum of the photon, Q2 = −q 2 is its virtuality, and P , P ′ , are the proton momenta before and after the collision.For the production of light mesons and meson pairs the analysis is frequently done in the so-called Bjorken kinematics [33-35, 49, 50], when the hard scale is set by the photon virtuality Q 2 , exceeding significantly the nucleon mass m 2 N , as well as all light quark masses.For processes involving heavy mesons this regime is hardly achievable experimentally, due to the rapid decrease of the flux of equivalent photons as a function of Q.Furthermore, due to the expected smallness of the cross-sections, most of the detected events will proceed via quasi-real photons with Q ≈ 0. For the sake of generality, till the end of this section we will keep Q ̸ = 0, assuming that the virtuality Q is bound by Λ 2 QCD ≪ Q 2 ≲ m 2 Q , although eventually we'll take Q = 0 in the numerical estimates.Since the spectrum of equivalent quasi-real photons emitted from the electron falls off rapidly as a function of the transverse photon momentum, at high energies the photon-proton frame should be close to the laboratory frame, in which case the electron-proton collision axis points in the direction of the axis ẑ.In the limit Q → 0, this frame, up to a trivial longitudinal boost1 , coincides with the frame used in earlier studies of exclusive photoproduction γp → γM p [17][18][19][20][21][22][23][24][25][26].In this frame, the polarization vectors of the virtual photons may be chosen as for the longitudinal and transverse polarizations respectively.The 4-momenta p 1 , p 2 of the produced heavy D-mesons can be parametrized in terms of the rapidities y a and transverse momenta p ⊥ a of these heavy mesons as where the positive rapidity is chosen in the photon direction.As we will see below, the cross-section falls rapidly as a function of transverse momenta; for this reason the dominant contribution to the cross-section in EIC kinematics comes from the region of relatively small momenta p ⊥ a .In these notations, the 4-vector ∆ of momentum transfer to the target is given by and the Mandelstam invariant t ≡ ∆ 2 is parametrized as The 4-momentum of the recoil proton should satisfy the onshellness condition (P + ∆) 2 = m 2 N , which provides an additional constraint The Equation ( 9) may be solved with respect to q − , yielding which, together with (1), allows to find the energy of the photon E γ ≈ q − /2 in terms of the kinematic variables y a , p ⊥ a of the produced D-mesons.In terms of these variables, the invariant energy W of the γp collision and the invariant mass M 12 of the produced heavy quarkonia pair can be rewritten as and respectively.In the high-energy limit q − , P + ≫ Q, M a ≫ {m N , |t|}, the results found earlier in this section may be simplified to In this kinematics the rapidities y 1 , y 2 and their difference can be rewritten in terms of the invariant Mandelstam variable(s) u 1 , u 2 , defined as In these notations, the conventional Bjorken variable x B can be represented as In the literature sometimes the variable x B is replaced by the so-called skewedness variable ξ, which is defined as [3] ξ This variable is directly related to the longitudinal light-cone momentum transfer ξ = −∆ + /2 P + = −∆ + / (2P + + ∆ + ).Using Eq. ( 16), it is possible to express ξ in terms of the rapidities y 1 , y 2 .As we can see from Figure 1, the typical values of ξ, x B are relatively small even for the lowest energy electron-proton beam at EIC; for this reason, we cannot disregard the gluon GPDs contributions.However, the kinematics of interest is still far from the saturation regime x B ≪ 1, and for this reason we may disregard saturation effects in our analysis.
The meson pair production in ep collisions is dominated by single photon exchange between the leptonic and hadronic parts; for this reason the cross-section of the process can be expressed as where dσ (T ) , dσ (L) are the photoproduction cross-sections of the transversely and longitudinally polarized virtual photons, y is the inelasticity (ratio of the energies of the virtual photon and electron), and dΩ h the phase volume of the produced heavy meson pair, which will be specified below.
The cross-section of the photoproduction process is related to the corresponding amplitude by where the δ-function in the right-hand side of ( 19) reflects the onshellness of the recoil proton.Using a light-cone decomposition (1-4), the argument of the δ-function can be rewritten as where ϕ is the azimuthal angle between the transverse momenta p ⊥ 1 , p ⊥ 2 , of the produced D-mesons.This allows to rewrite the δ-function in (20) as which permits to integrate out the dependence on ϕ.The condition |cos ϕ 0 | ≤ 1 at fixed invariant energies W of the photon-proton collision leads to a nontrivial constraint on the possible rapidities and transverse momenta of the produced D-mesons.In Figures (2,3) we illustrate the kinematically allowed domains for transverse momenta of the D-mesons, at several fixed rapidities of the D-mesons and photon-proton energies.From Figure (2) we may see that an increase of D-meson rapidities y 1 , y 2 leads to an increase of their longitudinal momenta, and due to energy conservation, this decreases the allowed transverse momenta of the produced mesons.These qualitative explanation also allows to understand Figure (3): increasing the rapidity or the transverse momentum of one of the D-mesons at fixed total invariant energy W inevitably decreases the allowed rapidities or transverse momenta of the other meson.The color coding in all plots reflects the value of cos ϕ 0 , fixed from (22).In the heavy quark mass limit, the difference of masses of the various D-mesons is suppressed as where Λ is some soft scale, and m Q is the mass of the heavy quark [51].For simplicity we disregard this difference altogether, assuming ξ -1.
-1.The complexity of the above-mentioned kinematic restrictions is a consequence of fixing the invariant energy W .In view of the symmetry of the final-state D-meson pair w.r.t.permutations of both mesons, in electroproduction experiments it might be easier not to impose conventional constraints on W and work with D-meson momenta as independent unconstrained variables.The energy W in each event can be rewritten in terms of these variables.The  δ-function in the right-hand side of ( 19) may be represented as The integration of (18), over all x B ∼ 1/ W 2 + Q 2 − m 2 N , allows to rewrite the electroproduction cross-section as dσ where the variables (y 1 , p 1⊥ , y 2 , p 2⊥ , ϕ) fully characterize the kinematics of the process, and dσ γp→M1M2p are the photoproduction cross-sections of longitudinal and transverse photons, for photons energy given by (10).

B. Amplitudes of the meson pair production process
As we have seen in the previous section, in the kinematics that we consider here, the typical values of x B , ξ are small, although still far from the saturation regime, and the dominant contribution to the cross-section comes from the region of small transverse momenta of D-mesons, p ⊥ ≲ M D .In this kinematic regime, it is convenient to use a collinear factorization framework for the evaluation of the amplitudes A (a) γp→M1M2p , and express the latter as a convolution of perturbative coefficient functions with distribution amplitudes of produced D-mesons and the GPDs of the target [1][2][3][4][5][6].The natural hard scales in this approach are the heavy quark mass m Q and the invariant mass M 12 ∼ m Q of the produced D-meson pair.In order to avoid nonperturbative final state interactions, we'll assume additionally that all hadrons in tthe final state are kinematically well-separated from each other, having sufficiently large (frame-invariant) relative velocities for each pair of final-state hadrons.For production of meson pairs at central rapidities, this constraint is satisfied almost everywhere, except in the small near-threshold region M 12 ∼ M 1 +M 2 .Since in the evaluation of the coefficient functions the transverse momenta of partons are disregarded, in order to guarantee the validity of ( 27), we will only consider the kinematics when the mesons are separated from each other at least by one unit of rapidity, ∆y ≳ 1.
In the collinear factorization picture all hadrons can be replaced with collinear partons, convoluted with nonperturbative distribution amplitudes which describe the momentum sharing between these hadrons, and contracted with appropriate spin projectors.For the proton target, in leading twist the Fock state is dominated by the two-partonic component, which is described by the quark and gluon GPDs.In what follows we disregard the contribution of the intrinsic heavy flavors in the proton, since they should be very small according to phenomenological estimates.Furthermore, we also do not take into account the contributions of the poorly known transversity GPDs H T , E T , HT , ẼT , since their contributions are accompanied by momentum transfer to the target ∆ µ and thus should be small for unpolarized observables in the small-∆ (small-t) kinematics which we study here (see [15,52] for more details).In the chiral even sector we should take into account the contributions of both quarks and gluons.The other (chiral-even) GPDs contribute to the amplitude of unpolarized process in the combination where the index a distinguishes the transverse and longitudinal polarizations of the photons, and, inspired by the previous studies of DVCS and DVMP [53,54], we introduced the double meson form factors Ha (ξ, ∆y, t) The variables z 1 , z 2 , are the light-cone fractions of the total momentum carried by the quarks in the D-mesons.In the evaluation of (29,30) we took into account that the final-state D-mesons are kinematically separated from each other; for this reason the Fock state of the final system is a direct product of Fock states of individual D-mesons, which in the heavy quark mass limit might be described by D-meson distribution amplitudes φ D (z).The detailed definitions of these distributions and discussion of their parametrizations may be found in Appendix A. The remaining partonic amplitudes a can be evaluated perturbatively, taking into account the diagrams shown in Figures 4, 5.The applicability of the perturbation theory is justified for this purpose in the heavy quark mass limit and large invariant mass M 12 of the produced D-mesons.The first diagram in the upper row of Figure 4 presents the dominant O (α s ) contribution; however, it is forbidden kinematically: the heavy quarks in the produced D-mesons carry significant fractions of the momentum, and thus are expected to have a large positive invariant mass (p c + p c) 2 ≳ 4m 2 c , whereas the photon in the electroproduction process has q 2 = −Q 2 ≲ 0. For the same reason, in the next-to-leading order there is no contributions which merely renormalize the propagators and vertices in this forbidden "leading order" diagram.Since we are mostly interested in the photoproduction regime, in what follows we will focus on the contribution of the transversely polarized photons; the contribution of the longitudinal photons is suppressed as ∼ Q/m Q and thus may be disregarded.The complete expressions for a and some technical details of their evaluation may be found in Appendix B. Other diagrams: representative (sub)leading ∼ O α2 s diagrams, which describe the (light) quark contributions to the D-meson pair production.The thin and thick lines correspond to light and heavy quarks respectively.Each diagram should be understood as a sum of diagrams with all possible permutations of the photon coupling to quark lines at fixed gluon vertices (so the second diagram in the upper row corresponds to 6 different diagrams, whereas each of the other diagrams should be understood as a sum of 7 different diagrams).We disregard the diagrams with heavy quark lines attached to the proton, since the intrinsic heavy flavors in the proton are negligibly small.
The quark and gluon GPDs of the target contribute in the amplitudes (29,30) integrated over the light-cone fraction x, and it is important to understand which region gives the dominant contribution in this convolution.The x-dependence of each individual Feynman diagram has a form of the rational function of the variable x, and for this reason the coefficient functions a can be represented as a sum of such contributions, where the functions P ℓ (x) , Q ℓ (x) are polynomials of the variable x.The polynomials Q ℓ (x) in the denominators can include up to n ℓ nodes x k in the region of integration, where n ℓ is the number of free propagators in the corresponding Feynman diagram 2 .The integral near the poles exists only in the principal value sense and should be evaluated using In Figure 6 we provide the density plots which show the dependence of the light quark coefficient function C q T (x, ξ, ∆y, z 1 , z 2 ) on some of its arguments.At fixed z 1 , z 2 , the poles show up as bright lines on a dark background.As could be seen from the analytic expressions in Appendix B, in the small-ξ limit all the poles scale as x (ℓ) k ∼ ξ, with proportionality coefficient which depends on ∆y, z 1 , z 2 ; for this reason, all the pole trajectories (bright lines in the density plot 6) are nearly straight.However, in the final result the coefficient functions contribute in the convolution with relatively broad D-meson distribution amplitudes, via the effective (integrated) coefficient function As we can see from the next Figure, 7, the function C q int does not have any singularities.The imaginary part of C q int appears due to a deformation of the integration contours in the integrals over z 1 , z 2 near the poles, which is carried out using the conventional ξ → ξ − i0 prescription3 .The integrated function C q int is mostly concentrated in the region |x| ≤ ξ, which suggests that the process is mainly sensitive to the behavior of GPDs in that domain (the so-called ERBL region).For the gluonic coefficient function, we observe a similar behavior: there are poles for fixed z 1 , z 2 , alghough they are smeared after convolution with the distribution amplitudes.

III. NUMERICAL RESULTS
For the sake of definiteness we will make predictions using the Kroll-Goloskokov parametrization of the GPDs [52,[55][56][57][58][59].This parametrization effectively incorporates the evolution of the generalized parton distributions, introducing a mild dependence of the model parameters on the factorization scale µ F .In Figure 8 we show how the cross-sections depend on the choice of this factorization scale.This dependence is mild at moderate energies, but becomes very pronounced at very high energies (small x B ).Such behavior is not surprising: it is known from studies of other channels [60][61][62][63] that this dependence exists, due to the omitted higher order corrections, which become especially important in the kinematics of small-x B . in what follows, for the sake of definiteness, we will choose the factorization scale In Figure 9 we can see the Q-dependence of the cross-section of the photoproduction subprocess 2 , since in the hard amplitudes the contribution of the O Q 2     terms is negligibly compared to the O M 2 12 contributions.However, for large 2 the virtuality Q 2 turns into the hard scale and leads to strong suppression of the cross-section.In electroproduction experiments the flux of equivalent photons decreases rapidly as a function of Q 2 , and thus the kinematics where Q-dependence becomes pronounced, is hardly achievable in the foreseeable future.For this reason in what follows we will focus only on the photoproduction regime Q ≈ 0.
In the Figures (10,11) we show the dependence of the cross-section (26) on the momentum transfer t to the target, and the related distributions of the produced D-mesons on transverse momenta and angle between them in the photon-proton frame.In the collinear factorization picture, the transverse momenta are disregarded in evaluation of the coefficient function, for this reason this dependence stems entirely from the t-dependence implemented in GPDs.The phenomenological analyses suggest that this dependence should exhibit a very pronounced (nearly exponential) suppression as a function of |t|.For this reason, the D-meson pairs are produced predominantly with small oppositely directed (ϕ ≈ π) transverse momenta, the so-called back-to-back kinematics which minimizes the momentum transfer t to the target.
In Figure 12 we show the dependence of the cross-section on the rapidities of the produced quarkonia.For the sake of definiteness, we assumed that the D-mesons are kinematically separated by a constant rapidity gap ∆y = 1 units.The growth of the cross-section with Y can be understood if we take into account that in the chosen kinematics, the increase of Y leads to an increase of the invariant energy W 2 , a corresponding decrease of x B , ξ and a growth of the gluon GPDs.For the same reason, the cross-section has a mild dependence on proton energy E p at constant average rapidity Y = (y 1 + y 2 )/2.The magnitude of the cross-section depends significantly on the quantum numbers (flavor content) of the produced D-mesons.It is instructive to understand the main sources of this dependence.As we discussed earlier in Section II A, our choice of kinematics corresponds to relatively small values of x B , ξ ≪ 1, and therefore we expect that all quark distributions are dominated by sea quarks, which have a mild dependence on flavor.The mass of the light quarks might be disregarded in the collinear approximation.However, some terms in the coefficient function, namely the diagrams which correspond to coupling of the photon to light quarks, include flavor-dependent electric charges e ℓ in prefactors.Their interference with diagrams which correspond to the coupling of the photon to heavy quarks can be constructive or destructive.Furthermore, there is additional flavor dependence introduced by the meson decay constant f D ≈ 209 MeV, f Ds ≈ 249 MeV, which can change the result for the crosssection by up to a factor of two, since the decay constant contributes as ∼ f 4 D .For the D + D * − mesons at small (Y, W ), the dominant contributions stem from the quark sector, and the contribution of the gluons becomes more pronounced at higher Y .For charged strange mesons D + s D * − s , the contribution of the quark sector is slightly smaller and is on par with the contribution of the gluons.Finally, for neutral D 0 D * 0 pairs, due to destructive interference of the contributions of quark and gluon GPDs, the cross-section is smaller than for the other mesons.In Figure 13 we compare side-by-side the cross-sections for different mesons, as well as show (in the upper horizontal axis) the Left plot: Dependence of the cross-section on the transverse momentum pT of produced D-mesons, at fixed azimuthal angle ϕ12 between the transverse momenta.Right plot: Dependence on the azimuthal angle ϕ12, at fixed transverse momentum pT of the produced D-mesons.The colored bands reflect the uncertainty due to choice of the factorization scale µF : the central line corresponds to factorization scale µF = 4 GeV ≈ 2mD, whereas the upper and lower limits were evaluated with µF = 8 GeV and µF = 2 GeV respectively.In the right plot, for better legibility, the cross-sections for Ep = 275 GeV and Ep = 41 GeV are multiplied by constant factors ×4 and ×1/4 respectively (shown near the right edge of each plot).A sharp peak at large angles ϕ12 ≈ π and relatively small transverse momenta corresponds to the so-called back-to-back kinematics, which minimizes the invariant momentum transfer |t|, as can be seen from Eq (7).
dependence on the invariant energy W γp .We may observe that the cross-sections have the same value for the same invariant energy W γp , independently of the proton energy E p or rapidity Y .
In Figure 14 we show the dependence of the cross-section on the rapidity difference ∆y between the two Dmesons at central rapidities.The cross-section decreases rapidly as function of ∆y, and similarly for the quark and gluon contributions.This behavior might be explained by an increase of the variables x B , ξ, and of the (minimal, longitudinal) momentum transfer |t min | at large ∆y and fixed Y ; this leads to a suppression of both quark and gluon contributions, due to t-dependence encoded in the partonic GPDs.In the limit ∆y → 0 the cross-sections remain finite, although numerically grow up to large values.Since our approach is justified only for the kinematically separated D-meson pairs, we do not consider that region.
Finally, in Figure 15 we provide predictions for the distribution of the produced D-meson pairs over their invariant mass M 12 .For all flavors of D-mesons and all energies W of γp pairs, the distributions have a very similar shape.Near the threshold, the cross-section grows due to increase of the phase volume.However, in this region the relative velocity v rel of the D-mesons is small, and thus potentially sizeable corrections might appear due to formation of the bound states (so-called tetraquarks).For higher values of M 12 ≳ 4 GeV the values of v rel become large enough to exclude such soft near-threshold effects, and the collinear theory becomes well-justified.The region of very large values of M 12 ≳ 6 GeV requires large rapidity difference ∆y or large transverse momenta of the heavy mesons.In this kinematics the momentum transfer t to the proton is large, so the cross-section becomes suppressed, in agreement with results shown in Figures (14,11).
We also analyzed a possibility to study the suggested channels in the kinematics of the proposed 22 GeV upgrade at JLab [64].However, we found that the cross-sections are extremely small in that kinematics and beyond the reach of experimental studies.This happens because in the kinematics of the proposed upgrade, the production would occur with relatively large values x B , ξ ∼ M 2  12 /W 2 ∼ 0.5−1, where the GPDs are strongly suppressed, both due to endpoint behavior of the underlying parton distributions ∼ (1 − x) n , n = 3 − 5, and additional suppression due to increase of the (longitudinal) momenta transfer to the target, |t| ≳ |t| min = m 2 N x 2 B /(1 − x B ).We also made similar estimates for photoproduction of B-meson pairs, yet found that even in the EIC kinematics the cross-sections are extremely small (sub-picobarn level).This could be understood from the structure of the coefficient functions, which scale as ∼ 1/m 5 Q in the heavy mass limit and thus are strongly suppressed.For this reason we do not include predictions neither for JLab 22 GeV kinematics, nor for B-meson pairs. .Dependence of the cross-section on the average rapidity Y = (y1 + y2) /2, at fixed rapidity difference y1 − y2 of the two mesons (positive rapidity is in direction of the photon/electron).The upper row corresponds to proton energy Ep ≈ 41 GeV, the middle row is for energy Ep ≈ 100 GeV, and the lower row is for Ep ≈ 275 GeV.The curves marked as "quarks" and "gluons" correspond to contributions of only quark or only gluon GPDs.The curve marked as "All" takes them all into account (as well as additional term due to interference of quark and gluon contributions in the cross-section).See the text for more detailed discussion.γp →D GeV) the relative velocity v rel of the produced D-mesons is small, and thus our approach might be not reliable (see the text for more explanation).

IV. SUMMARY AND CONCLUSIONS
In this paper, we analyzed the potential of exclusive photoproduction of D-meson pairs for studies of the GPDs of the target.We analyzed the photoproduction of pseudoscalar-vector pairs with opposite C-parities and different flavor content (D + D * − , D 0 D * 0 and D + s D * − s ), which get their dominant contribution from the chiral-even GPDs.We focused on the kinematics of large invariant masses of the produced meson pairs, moderate values of x B ∈ 10 −3 , 10 −1 and small photon virtuality Q ), achievable with low-and middle-energy ep beams at the Electron Ion Collider.We performed evaluations in the collinear factorization approach in leading order over the strong coupling α s (m Q ).In all channels the amplitude of the process obtains comparable contributions from one of the light quark flavors and gluons.This feature might present a special interest for phenomenological attempts to disentangle the flavor structure of the light quark GPDs.The sensitivity of the process to different GPD kinematics is controlled by the so-called coefficient functions (partonic level amplitudes), which have nontrivial behavior in the so-called ERBL region |x| < ξ, yet vanish rapidly outside of it.Due to convolution with sufficiently broad distribution amplitudes of D-mesons, the xdependence of the coefficient functions in the region |x| < ξ is relatively moderate, with a mild peak around x ≈ ξ.The inverse deconvolution apparently is not possible in view of complexity of the coefficient function, however we believe that experimental study of this process might present new constraints for phenomenological models of GPDs, especially in the ERBL region.The pronounced t-dependence of the cross-section, which stems from the phenomenological GPD models, implies that the D-meson pairs are produced predominantly in back-to-back kinematics, with relatively small and oppositely directed transverse momenta of the mesons.
The results of this study complement our previous analyses [42,43] of exclusive quarkonia pair production, which is sensitive only to gluonic GPDs.If compared in the same kinematics, the cross-sections for D-meson pairs are larger due to contribution of quark GPDs, and different structure of the gluonic coefficient functions due to replacement of one of the heavy quarks with nearly massless light quark.This fact should facilitate experimental studies of the suggested channel.The cross-section of the suggested process is comparable, by order of magnitude, to the crosssections of similar 2 → 3 processes (γ * p → γM p, M = π, ρ) suggested recently in the literature [17][18][19][20][21][22][23][24][25][26][27][28] for exhaustive studies of the light quark GPDs.This happens because the suppression due to heavy quark mass in D-meson pair production numerically is on par with a suppression by the additional fine-structure constant α em in processes which include additional emitted photon.For this reason both γ * p → γM p and D-meson production could be used as complementary tools for the study of both quark and gluon GPDs.

ACKNOWLDGEMENTS
We thank our colleagues at UTFSM university for encouraging discussions.This research was partially supported by Proyecto ANID PIA/APOYO AFB220004 (Chile) and Fondecyt (Chile) grants 1220242 and 1230391.

Appendix A: Distribution amplitudes of the D-mesons
In this appendix, for the sake of completeness, we provide definitions and parametrizations of the D-meson distribution amplitudes.For light mesons the distribution amplitudes are conventionally characterized (ordered) by the twist of the corresponding quark-antiquark operator [65][66][67][68][69].For spinless pseudoscalar meson M at leading twist, there is only one distribution amplitude, defined as where p is the momentum of the meson (we assume that the meson moves in the plus-direction), z is the fraction of the momentum carried by the quark, and L is the standard path-ordered gauge link.Similarly, for the vector (spin-1) mesons there are 2 independent leading-twist distributions ⊥ defined as where ε M (p) is the (transverse) polarization vector of the meson, which satisfies n • ε J/ψ = p • ε * J/ψ, µ (p) = 0.The definitions (A1-A4) may sometimes include distribution amplitudes normalized to unity, taking out the explicit meson decay constants f M in the right hand side.To avoid ambiguity, we will use the notation φ(z) (with corresponding other indices) for such normalized DAs.These findings imply that during perturbative evaluation of the coefficient function, the meson formation from perturbative quarks might be described using the effective vertices Sometimes this scheme is also extended to heavy D-mesons.However, for processes involving heavy quarks, the twist-based classification of distribution amplitudes requires adjustments, since the heavy quark mass limit breaks twist-based suppression of higher twist contributions in physical amplitudes.Instead of this, the distribution amplitudes are ordered by the inverse powers of the heavy quark mass m Q , using the Heavy Quark Effective Theory (HQET) framework [51].For the pseudoscalar and vector D-mesons containing c-quark, this implies that we should replace (A5) with [44,51,70,71] where v µ = p µ /M D is the 4-vector D-meson velocity, ε µ (p) is the 4-vector polarization for the vector mesons, f D is the corresponding decay constant.For the charge-conjugate D-mesons (containing c quark), the definitions (A6) . Typical shape of the D-meson distribution amplitude φD(z).For definiteness this plot was done using "model 2" from [76]; other models described in [71,74,75] have very similar shapes.
should be modified as v µ → −v µ .Due to heavy quark spin-flavor symmetry [51,72,73] it is that the functional form φ D should be approximately the same for pseudoscalar and vector mesons.While there is a plethora of different phenomenological and model-based parametrizations [71,[74][75][76], they all suggest that φ D (z) should have a broad peak near z ≈ 2/3 and vanish at the extremes, as shown in Figure 16.This shape suggests that up to 1/3 of the momentum of the D-meson is carried by the light-quark.

Appendix B: Evaluation of the coefficient functions
The coefficient functions (partonic amplitudes) may be evaluated using standard light-cone rules, which may be found in [1,3,49,[77][78][79].In the leading order, this requires the evaluation of all the diagrams shown in Figures 4,5.This evaluations drastically simplify in the collinear factorization approach, assuming that in our hierarchy of scales the mass of the heavy quark m Q and photon virtuality Q are large parameters, Q ∼ m Q ∼ W ≡ √ s γp , and omitting the mass of the proton m N and all the transverse momenta in the coefficient functions evaluation.Below, in subsections B 1, B 2 we discuss some technical details and provide final results for the coefficient functions of the light quarks and gluons.

Light quark contribution
The leading twist chiral even quark GPDs H q , E q , Hq , Ẽq , which are expected to give the dominant contributions, are conventionally defined from quark-antiquark correlators [3] In what follows we work in the light-cone gauge, so the gauge link becomes trivial, L (ζ 1 , ζ 2 ) = 1.The skewedness variable ξ = (P ′+ − P + ) / (P ′+ − P + ) ca be related to the rapidities of the heavy mesons using (17).For the sake of brevity we do not show explicitly the dependence on the factorization scale µ, which is described by the conventional DGLAP evolution equations.The factorization theorems allow to rewrite the amplitudes of physical processes as convolution of these GPDs with cross sections of partonic processes γq → M 1 M 2 q, where the quarks before and after interactions have light-cone momenta k + 1,2 ∼ (x ± ξ) P + .The transverse part of the quark momentum ∆ ⊥ is small compared to typical hard scales (Q, m Q , M 12 ) and thus may be disregarded in the evaluation of the coefficient function.The evaluation of the amplitude requires calculation of all the diagrams shown in Figure 4, together with all possible permutations of photon vertices, assuming that gluon vertices are fixed.The explicit expressions for the quark spinors and projectors onto F q and F q may be found in [3]; the spinor algebra was done using FeynCalc package for Mathematica [80,81].While the total number of diagrams is huge, in the collinear factorization picture many diagrams vanish.This happens because the transverse momenta of partons are disregarded in the collinear approach, so all the light quark propagators are given by the gluon propagators in light-cone gauge satisfy the light-cone gauge condition n µ G µν = 0, and for the γ + , γ − matrices in the numerator we have standard relations of Dirac algebra For example, this leads to a cancellation of all diagrams with 3-gluon vertices, whose representative is shown in the second diagram in the first row of Figure 4: the gluon vertex itself is proportional to the difference of gluon momenta ℓ i − ℓ j , which in the collinear limit has only light cone ± components.The gluon propagator for collinear gluons is proportional to the transverse part of the metric tensor, g ⊥ µν = g µν − n , where n (±) µ are the light-cone vectors in the ± direction.After contraction of such propagators with the ± components of ℓ µ i − ℓ ν j , we can see that the 3-gluon contribution vanishes (beyond collinear approximation, giving corrections ∼ O t/M 2 to the coefficient function).The evaluation of the remaining diagrams is straightforward and was done using the FeynCalc package for Mathematica [80,81].In the evaluation of the coefficient function, we replaced the final-state D-mesons with a system of collinear heavy-light quark-antiquarks (qQ and Qq respectively), adding appropriate spinor projectors, as required by the definitions A6 (see Appendix A for details).To avoid confusion, we will use the notations z 1 , z 2 for the light-cone momenta fractions of quarks and z1 ≡ 1 − z 1 , z2 ≡ 1 − z 2 for antiquarks in the first and the second mesons.In what follows we focus only on the Q ≈ 0 limit and provide the results for the transversely polarized photons (extensions for the case Q ̸ = 0 is straightforward, although tthe final expressions become too lengthy for publication).The final result for the light quark coefficient functions are where T and ε D * are the polarization vectors of the incident photon and produced vector meson, e ℓ , e H , are the charges of the light and heavy quarks (in units of proton charge), and the constants a k , b k , c k , d k , ãk , bk , ck , dk are given by 2 = e 2∆y 1 − e ∆y 2m 5  Q z 1 (e ∆y z1 + z 2 ) (1 + e ∆y z 2 ) (z 1 + e ∆y z 2 ) × (B10) e 2∆y e ∆y − 1 2m 5 Q z 1 (1 + e ∆y z1 ) (e ∆y z1 + z 2 ) (z 1 + e ∆y z 2 ) × (B11) 2m  We may observe that, as a function of variable x, each of the contributions has at most one or two poles, whose position depends on the skewedness ξ, rapidity difference ∆y, and the light-cone fractions z 1 , z 2 carried by quarks in the corresponding D-mesons.The endpoint singularities at z 1 = 0, z 1 = 1, z 2 = 0 and z 2 = 1 do not present any difficulties, neither for the factorization property at conceptual level, nor for numerical integration, since in the amplitude these coefficient functions contribute multiplied by the distribution amplitudes of the D-mesons, which vanish rapidly at these points.

Gluonic contribution
The definition of the leading twist gluon GPDs is very similar to the quark GPDs definitions and differ only by the structure of the operators inserted between the initial and final states, namely [3,60] where L is the standard path-ordered gauge link.In what follows we'll use the standard light-cone gauge A + = 0, in which L = 1, and the remaining two-gluon operators in (B21, B22) may be rewritten in the form The integration over z in (B21, B22) effectively corresponds to a transition to the momentum space, where the derivatives ∂ + z1 , ∂ + z2 convert into the multiplicative factors k + 1,2 ∼ (x ± ξ) P + , and thus (B21, B22) may be rewritten in a joint form as [60] The factors ±i0 in the denominator of (B26) reflect the conventional ξ → ξ − i0 prescription used to define the deformation of the integration contour near the poles of the amplitude.The structure of the numerator of Eq. (B26) implies that the coefficient functions C a and Ca may be extracted from the corresponding partonic amplitude of γg → D 1 D 2 g, contracting the Lorentz indices of the gluons with g ⊥ µν and ε ⊥ µν respectively.In the leading order over α s , the partonic amplitude γg → D 1 D 2 g gets contribution from the Feynman diagrams shown in the Figure 5.Each diagram in that Figure should be understood as a sum of two contributions, with photon coupled to heavy or light quark line, and complemented with diagrams which differ by permutation of the t-channel gluons, as demonstrated in Figure 17.The latter permutation leads to contributions which differ only by the sign in front of the light-cone variable x and interchange of the Lorentz indices µ ↔ ν.As explained earlier, the coefficient functions C a , Ca , may be found contracting the free Lorentz indices µ, ν with symmetric g ⊥ µν or antisymmetric ε ⊥ µν ; for this reason C a and Ca will be even or odd functions of the variable x respectively.In the collinear picture we assume that the quarks and antiquarks in the produced D-mesons carry the light-cone momenta z a p a and (1 − z a ) p a respectively, so all the parton momenta in leading-order diagrams of Figure 5 may be fixed by energy-momentum conservation and represented as linear combinations of the momenta of the D-mesons and t-channel gluons.Using (B26), we may rewrite the coefficient functions of the transversely polarized photons as where the constant κ was defined earlier in (B7), and the factors (x ± ξ ∓ i0) −1 in (B27, B28), come from (B26).
Since the gluon GPDs H g , E g , are even functions of the variable x, and Hg , Ẽg are odd functions [3], in tthe convolution over x both terms in the numerators of (B27, B28) give equal nonzero contributions.Numerically the dominant contribution to unpolarized cross-section comes from the GPD H g , whereas the contribution of Hg is negligibly small.The evaluation of the diagrams from Figure 5 was done using the FeynCalc package for Mathematica [80,81] (see the beginning of the previous Appendix B 1 for some technical details); the contributions to C a and Ca were extracted, contracting the free Lorentz indices µ, ν of the partonic amplitude with g ⊥ µν or ε ⊥ µν , respectively.Explicit expressions for the functions C T , CT are given by where e ℓ , e H , are the charges of the light and heavy quark (in units of the proton charge), and the expressions in front of them correspond to a sum of the diagrams in which photon is connected to the light or heavy quark lines respectively.Explicitly, these contributions are given by a 1 = 2e 2∆y e 2∆y − 1 −e ∆y (ξ + 1)(x − ξ) + (ξ + 1)z  Similar to the case of quark coefficient functions, all the contributions can have up to three non-coinciding poles as a function of variable x.The position of the poles depends on the kinematics of the produced D-mesons (variables ∆y, ξ), as well as the light-cone fractions z 1 , z 2 carried by the quarks inside the D-mesons.As we mentioned in the main text, all the convolutions (integrals) which include these coefficient functions, should be interpreted in the principal value sense, using the standard ξ → ξ − i0 prescription [60] for contour deformation near the poles.

Figure 1 .
Figure 1.The contour plots illustrate the dependence of the skewedness variable ξ = − P + f − P + i / P + f + P + i on the rest-frame rapidities y1, y2, of the D-mesons, for different fixed photon virtualities Q and proton energies Ep.Each dashed line corresponds to a line ξ = const in the y1, y2 plane, with value of constant ξ shown as a label on the contour line.For simplicity, we disregard the transverse momenta of the produced D-mesons.The upper and lower rows differ by the choice of the proton energy Ep (41 and 100 GeV respectively).

Figure 2 .
Figure 2. (Color online) The colored bands show the kinematically permitted regions for D-meson pair production, for fixed rapidities y1, y2, at fixed photon energy Eγ, virtuality Q and proton energy Ep.The color of each pixel reflects the cosine of the angle ϕ (azimuthal angle between the transverse momenta of the D-mesons), fixed from (22).The variable Y = (y1 + y2)/2 is the average rapidity of the D-mesons.The left and right columns differ by the choice of the rapidity difference (∆y = 0 and ∆y = 1 respectively).The upper and lower rows differ by the choice of the photon energy Eγ in the lab frame (3 and 5 GeV respectively).The borders of the kinematic domains have very mild dependence on Q 2 (up to 2-3 units of M 2 D ) and mild sensitivity to the proton energy Ep.

Figure 3 .
Figure 3. (Color online) Kinematic constraints on the rapidity and transverse momenta of one of the D-mesons, when the momentum of the other (its rapidity y2, transverse momentum p 2⊥ ) is fixed.The photon virtuality Q and energy Eγ, as well as proton energy Ep are fixed to values shown in the upper part of the Figure.The color of each pixel reflects the value of the angle ϕ (azimuthal angle between the transverse momenta of the D-mesons), fixed from (22).The left and right columns differ by the choice of the rapidity interval for y2.

Figure 4 .
Figure 4. First diagram in the upper row: the leading-order ∼ O (αs)-diagram, which describes the (light) quark contributions to the D-meson pair production.As explained in the text, a diagram of this type does not contribute in the collinear approximation.Other diagrams: representative (sub)leading ∼ O α 2s diagrams, which describe the (light) quark contributions to the D-meson pair production.The thin and thick lines correspond to light and heavy quarks respectively.Each diagram should be understood as a sum of diagrams with all possible permutations of the photon coupling to quark lines at fixed gluon vertices (so the second diagram in the upper row corresponds to 6 different diagrams, whereas each of the other diagrams should be understood as a sum of 7 different diagrams).We disregard the diagrams with heavy quark lines attached to the proton, since the intrinsic heavy flavors in the proton are negligibly small.

Figure 5 .
Figure 5.The leading order diagrams, which describe the contribution of the gluons to D-meson pair production.Each diagram should be understood as a sum of two diagrams with complementary assignment of quark flavors (heavy-light and light-heavy) .

Figure 6 .
Figure 6.Density plot which illustrates the light quark coefficient function C q T (in relative units) as a function of the variables x and skewedness ξ, at fixed rapidity difference ∆y, between the heavy D-mesons and the fixed quark light-cone fractions z1, z2 in D-mesons.The first and the second row differ by the choice of the values of z1, z2; the left and the right columns differ by the values of ∆y.For the sake of definiteness we consider D + D * − in all plots).Thick white lines effectively demonstrate the position of the poles x ℓ k of the coefficient function (31).For reference, we also added red dashed lines x = ±ξ, which separate DGLAP and ERBL regions.

Figure 7 .
Figure 7. Left: Density plot which shows the integrated coefficient function C q int (33), as a function of the variables x and skewedness ξ, at fixed rapidity difference ∆y between the heavy D-mesons.Right: The real, imaginary and absolute parts of the coefficient function C q int as a function of x/ξ, at fixed ξ, ∆y.For the sake of definiteness we consider D + D * − in all plots.Both plots illustrate that the singularities were smeared after convolution with the final-state D-meson distribution amplitudes.Similar behavior is observed for other choices of kinematics and other mesons.

4 Figure 8 .
Figure 8.The factorization scale dependence of the cross-section, for different D-meson pairs.For the sake of definiteness we consider the production at central rapidities (Y = (y1 + y2)/2 ≈ 0), for different energies Ep of the proton beam.All these values in the photoproduction regime correspond to small xB ≲ 5 × 10 −2 ≪ 1.All frame-dependent variables are given in the reference frame described in Section II A.

4 Figure 9 . 4 Figure 10 .
Figure 9. Dependence of the photoproduction cross-section(26) on the virtuality Q of the photon.In the left and right plots we compare predictions for different rapidities y J/ψ , yη c and different proton energies Ep.Both plots clearly illustrates the transition from photoproduction to Bjorken regime in the region Q ∼ 1 − 2 M J/ψ .In both plots the photon energy is evaluated from(1,10).All frame-dependent variables are given in the reference frame described in Section II A.
Figure 11.Left plot: Dependence of the cross-section on the transverse momentum pT of produced D-mesons, at fixed azimuthal angle ϕ12 between the transverse momenta.Right plot: Dependence on the azimuthal angle ϕ12, at fixed transverse momentum pT of the produced D-mesons.The colored bands reflect the uncertainty due to choice of the factorization scale µF : the central line corresponds to factorization scale µF = 4 GeV ≈ 2mD, whereas the upper and lower limits were evaluated with µF = 8 GeV and µF = 2 GeV respectively.In the right plot, for better legibility, the cross-sections for Ep = 275 GeV and Ep = 41 GeV are multiplied by constant factors ×4 and ×1/4 respectively (shown near the right edge of each plot).A sharp peak at large angles ϕ12 ≈ π and relatively small transverse momenta corresponds to the so-called back-to-back kinematics, which minimizes the invariant momentum transfer |t|, as can be seen from Eq(7).

Figure 13 .
Figure 13.Side-by-side comparison of the cross-sections, for different mesons and different proton energies Ep.The lower horizontal axis shows the average rapidity Y = (y1 + y2) /2, whereas in the upper axis we show the corresponding values of invariant energy W ≡ √ sγp.A comparison of different plots shows that the cross-sections for the same W , but different proton energy Ep and Y coincide with each other.

Figure 14 .
Figure 14.Dependence of the D + D * − production cross-section on the rapidity difference ∆y, at fixed average rapidity Y = (y1 + y2) /2 and different proton energies Ep.For other mesons we observe a similar behavior.

Figure 15 .
Figure 15.The distribution of the produced D-meson pairs over their invariant mass M12, for several fixed invariant energies W of the γp collision (left plot) and for different flavors at fixed W (right plot).In the near-threshold region M12 ≈ 2MD (M12 ≲ 4 GeV) the relative velocity v rel of the produced D-mesons is small, and thus our approach might be not reliable (see the text for more explanation).

Figure 17 .
Figure 17.Schematic illustration of the diagrams with direct and permuted t-channel gluons, which are related to each other by an inversion of sign in front of light-cone fraction x ↔ −x, and permutation of the Lorentz indices µ ↔ ν.