Exclusive photoproduction of heavy quarkonia pairs

In this paper we study the high energy exclusive photoproduction of heavy quarkonia pairs in the leading order of the strong coupling constant $\alpha_{s}$. In the suggested mechanism the quarkonia pairs are produced with opposite charge parities, and have predominantly oppositely directed transverse momenta. Using the Color Glass Condensate approach, we estimated numerically the production cross-sections in the kinematics of the forthcoming electron-proton colliders, as well as proton-ion colliders in ultraperipheral collisions. We found that the cross-sections are within the reach of planned experiments and can be measured with reasonable precision. The suggested mechanism has significantly larger cross-section than that of the same $C$-parity quarkonia pair production.


I. INTRODUCTION
The production of heavy quarkonia is frequently considered as a clean probe for the study of gluon dynamics in high-energy interactions, since in the limit of heavy quark mass m Q the running coupling becomes small, and it is possible to apply perturbative methods for the description of quark-gluon interactions. In many scattering problems the small size of the color singlet heavy quarkonium provides additional twist suppression [1,2], thus facilitating the applicability of perturbative treatments. The modern NRQCD framework allows to use quarkonia production as a powerful probe of strong interactions, systematically taking into account various perturbative corrections [3][4][5][6][7][8][9][10][11][12][13][14].
For precision studies of hadronic interactions, exclusive production presents a special interest in view of its simpler structure. However, up to now most of the experimental data on exclusive heavy quarkonia production were limited to channels with single quarkonia in the final state. This limitation was largely motivated by probable smaller crosssections of events with more than one quarkonia in the final state. Nevertheless, processes with two mesons in the final state present a lot of interest and have been the subject of studies since early days of QCD [15][16][17][18]. A recent discovery of all-heavy tetraquarks, which might be consider molecular states of two quarkonia, has significantly reinvigorated interest in the study of this channel [19][20][21][22][23][24][25][26][27][28][29].
In LHC kinematics most of the previous studies of exclusive double quarkonia production [30][31][32][33][34] focused on the so-called two-photon mechanism, γγ → M 1 M 2 , which gives the dominant contribution for the production of quarkonia pairs with the same C-parity in ultraperipheral collisions. Studies beyond the double photon fusion show that, in a TMD factorization approach, the exclusive double quarkonia production could allow to measure the currently unknown generalized transverse momentum distributions (GTMDs) of gluons [35]. However, in LHC kinematics the cross-section of this process can get sizable contributions from the so-called multiparton scattering diagrams. Such contributions depend on the poorly known multigluon distributions, leading to potential ambiguities in the theoretical interpretation of the data.
Electron-proton collisions have a significant advantage for studies of heavy quarkonia pair production, due to a smaller number of production mechanisms compared to hadron-hadron collisions. Moreover, precision studies of double quarkonia production in ep collisions could become possible after the launch of new high luminosity facilities, such as the forthcoming Electron Ion Collider (EIC) [36][37][38][39], the future Large Hadron electron Collider (LHeC) [40], the Future Circular Collider (FCC-he) [41][42][43] and the CEPC collider [44,45]. The main objective of this manuscript is the study of exclusive production of heavy quarkonia pairs, γp → M 1 M 2 p, in the kinematics of the above-mentioned electron-proton colliders. Potentially such production might be also probed in ultraperipheral heavy ion and protonion collisions. However, in these cases the analysis becomes more complicated in view of possible contributions of other mechanisms [30][31][32][33][34]. The large mass m q of the heavy flavors justifies the perturbative treatment in a wide kinematic range, without additional restrictions on the virtuality of the incoming photon Q 2 or the invariant mass of the produced quarkonia pair. In absence of imposed kinematic constraints, the dominant contribution to the crosssection will come from events induced by quasi-real photons with small Q 2 ≈ 0 and relatively small values of x B 1. In this kinematics it is appropriate to use the language of color dipole amplitudes and apply the color dipole (also known as Color Glass Condensate or CGC) framework [47][48][49][50][51][52][53][54][55]. At high energies the color dipoles are eigenstates of interaction, and thus can be used as universal elementary building blocks, automatically accumulating both the hard and soft fluctuations [56]. The light-cone color dipole framework has been developed and successfully applied to the phenomenological description of both hadron-hadron and lepton-hadron collisions [57][58][59][60][61][62][63][64], and for this reason we will use it for our estimates.
The paper is structured as follow. Below, in Section II, we evaluate theoretically the cross-section of exclusive photoproduction of heavy quarkonia pairs in the CGC approach. In Section III we present our numerical estimates, in the kinematics of the future ep colliders (EIC, LHeC and FCC-he) and ultraperipheral pA collisions at LHC. Finally, in Section IV we draw conclusions.

A. Kinematics of the process
We would like to start our discussion of the theoretical framework with a short description of the kinematics of the process. Our choice of the light cone decomposition of particles momenta is similar to that of earlier studies of pion pair [65][66][67][68][69] and single-meson production [70][71][72][73][74][75][76][77][78][79][80][81][82][83][84][85][86]. However, we should take into account that the mass of the quarkonium, in contrast to that of pion, is quite large, and thus cannot be disregarded as a kinematic higher twist correction. Besides, for photoproduction this mass can appear as one of the hard scales in the problem.
In what follows we will use the notations: q for the photon momentum, P and P for the momentum of the proton before and after the collision, and p 1 , p 2 for the 4-momenta of produced heavy quarkonia. For sake of generality we will assume temporarily that the photon can have a nonzero virtuality −q 2 = Q 2 , taking later that for photoproduction Q 2 = 0. We also will use the notation ∆ for the momentum transfer to the proton, ∆ = P − P , and the notation t for its square, t ≡ ∆ 2 . The light-cone expansion of the above-mentioned momenta in the lab frame is given by 1 where y a , p ⊥ a are the rapidity and transverse momentum of the quarkonium a, and M a is its mass. Using conservation of 4-momentum, we may obtain for the momentum transfer to the proton ∆ = P − P = q − p 1 − p 2 = (5) and for the variable t ≡ ∆ 2 After the interaction the 4-momentum of the proton is given by and the onshellness condition (P + ∆) 2 = m 2 N allows to get an additional constraint Solving Equation (8) with respect to q · P , we get In earlier theoretical studies [65-68, 80, 87] the evaluations were done in the so-called symmetric frame, in which the axis z is chosen in such a way that the vectors q andP ≡ P + P do not have transverse components. Besides, all evaluations were done in the Bjorken limit, assuming infinitely large Q 2 and negligibly small masses of the produced mesons (pions). In our studies we consider quasi-real photons, with Q 2 ≈ 0, and moreover the heavy mass of the quarkonia does not allow to drop certain "higher twist" terms. For this reason the kinematic expressions in the symmetric frame become quite complicated, and there is no advantage in its use for photoproduction.
which allows to express the energy of the photon E γ ≈ q + /2 in terms of the kinematic variables y a , p ⊥ a of the produced quarkonia. In the kinematics of all experiments which we consider below, the typical values q + , P − {Q, M a , m N , t}, and for this reason we may approximate (9) as From comparison of (3) and (8) we may see that at high energies the light-cone plus-component of the photon momentum q + is shared between the momenta of the produced quarkonia, whereas the momentum transfer to the proton (vector ∆) has a negligibly small plus-component, in agreement with the eikonal picture expectations. The expressions (9, 10) allow to express the Bjorken variable x B , which appears in the analysis of this process in Bjorken kinematics, using its conventional definition . As was discussed in [88][89][90][91], in phenomenological studies usually it is assumed that for heavy quarks all the gluon densities and forward dipole amplitudes should depend on the so-called "rescaling variable" which was introduced in [88] in order to improve description of the near-threshold heavy quarkonia production. While the color dipole framework is usually applied far from the near-threshold kinematics, the use of the variable x instead of x B for heavy quarks improves agreement of dipole approach predictions with experimental data. In Bjorken limit the variable x coincides with x B . For small Q 2 ≈ 0 (photoproduction regime) the variable x B vanishes, whereas x remains finite and is given by the approximate expression In this study we are interested in the production of both quarkonia at central rapidities (in lab frame) by high-energy photon-proton collision. In this kinematics the variable x is very small, which suggests that the amplitude of this process should be analyzed in frameworks with built-in saturation, such as color glass condensate (CGC). In contrast, in the Bjorken limit (Q 2 → ∞, Q 2 /2p · q = const) we observe that the variable x can be quite large, so it is more appropriate to analyze this kinematics using collinear or k T -factorization. The latter case requires a separate study and will be presented elsewhere.
In the photoproduction approximation the invariant energy of the γp collision can be written as whereas the invariant mass of the produced heavy quarkonia pair is given by In electron-proton collisions the cross-section of heavy meson pairs is dominated by a single-photon exchange between its leptonic and hadronic parts, and for this reason can be represented as where we use the standard DIS notation y for the elasticity (fraction of electron energy which passes to the photon, not to be confused with the rapidities y a of the produced quarkonia). The subscript letters L, T in the righthand side of (16) stand for the contributions of longitudinally and transversely polarized photons respectively. The structure of (16) suggests that the dominant contribution to the cross-section comes from the region of small Q 2 .
In this kinematics the contribution of d σ L is suppressed compared to the term d σ T . This expectation is partially corroborated by the experimental data from ZEUS [92] and H1 [93], which found that for single quarkonia production in the region Q 2 1 GeV 2 the longitudinal cross-section dσ L constitutes less than 10% of the transverse cross-section dσ T . For this reason in this paper we'll disregard the cross-section dσ L altogether, while the relevant cross-section dσ T is where A γ T p→M1M2p is the amplitude of the exclusive process, induced by a transversely polarized photon, and φ is the angle between the vectors p 1 and p 2 in the transverse plane. The δ-function in (17) reflects conservation of plus-component of momentum, discussed earlier in (8).
Similarly, for exclusive hadroproduction pA → pAM 1 M 2 in ultraperipheral kinematics we may obtain the crosssection using the equivalent photon (Weizsäcker-Williams) approximation, where dn γ (ω ≡ E γ , q ⊥ ) is the spectral density of the flux of photons created by the nucleus, q ⊥ is the transverse momentum of the photon with respect to the nucleus, and the energy E γ of the photon can be related to the kinematics of produced quarkonia using Eq. (9, 10). The explicit expression for dn γ (ω ≡ E γ , q ⊥ ) can be found in [94]. The momenta p * i = p ⊥ i − q ⊥ are the transverse parts of the quarkonia momenta with respect to the produced photon. Due to nuclear form factors the typical values of momenta q ⊥ are controlled by the nuclear radius R A and are quite small, For this reason, for very heavy ions (A 1) we may expect that the p T -dependence of the cross-sections in the left-hand side of (18) largely repeats the p T -dependence of the crosssection in the integrand in the right-hand side. For the special and experimentally important case of p ⊥ -integrated cross-section, the expression (18) simplifies and can be rewritten as where In the following subsection II B we evaluate the amplitude A γ T p→M1M2p which determines the cross-sections of photoproduction processes.
B. Amplitude of the process in the color dipole picture Since the formation time of rapidly moving heavy quarkonia significantly exceeds the size of the proton, the quarkonia formation occurs far outside the interaction region. For this reason the amplitudes of the quarkonia production processes can be represented as a convolution of the quarkonia wave functions with hard amplitudes, which characterize the production of the small pairs of nearly onshell heavy quarks in the gluonic field of the target. In what follows we will refer to these nearly onshell quarks as "produced" or "final state" quarks. For exclusive production the cross-section falls rapidly as function of transverse momenta p T of the produced quarkonia, and for this reason we expect that the quarkonia will be produced predominantly with small momenta. In this kinematical region it is possible to disregard completely the color octet contributions [8,9]. As was shown in [89][90][91], this assumption gives very good description of the exclusive production of single quarkonia.
The general rules for the evaluation of different hard amplitudes in terms of the color singlet forward dipole amplitude were introduced in [47,[49][50][51][52][53][54][55] and are briefly summarized in Appendix A. This approach is based on the high energy eikonal picture, and therefore the partons transverse coordinates and helicities remain essentially frozen during propagation in the gluonic field of the target. The hard scale, which controls the interaction of a heavy quark with the strong gluonic field, is its mass m Q , so in the heavy mass limit we may treat this interaction perturbatively. However, the interaction of gluons with each other, as well as with light quarks, remains strongly nonperturbative in the deeply saturated regime.
In the leading order over the strong coupling α s (m Q ), there are a few dozen Feynman diagrams which contribute to the exclusive photoproduction of meson pairs. In what follows, it is convenient to represent them as one of the two main classes shown schematically in Figure 1. For the sake of definiteness we'll call "type-A" all the diagrams in which quarkonia are formed from different heavy quark lines, as shown in the left panel of Figure 1. The opposite case, when quarkonia are formed from the same quark lines, as shown in the right panel of the Figure 1, will be referred to as "type-B" diagrams. This classification is convenient for a discussion of symmetries, as well as for analysis of quarkonia production with mixed flavors. For example, production of B + c B − c pairs clearly gets contributions only from type-A diagrams, whereas production of mixed flavor hidden-charm and hidden-bottom quarkonia (e.g. J/ψ + η b ) gets contributions only from type-B diagrams.
In configuration space the eikonal interactions with the target do not affect the impact parameters of the partons, so the interaction basically reduces to a mere multiplication of target-dependent factors, as discussed in Appendix A. This allows to express the amplitude of the whole process as a convolution of the 4-quark Fock component wave function ψ QQQQ of the photon with dipole amplitudes and wave functions of the produced quarkonia. The amplitude of the process γ * p → M 1 M 2 p can be represented as a sum where A 1 and A 2 stand for contributions of all type-A and type-B diagrams. Explicitly, these amplitudes are given by .
where in the expressions (22,23) we introduced a few shorthand notations, which characterize the pair of heavy partons i and j: the relative distance between them r ij = x i − x j , the light-cone fraction α ij = α i / (α i + α j ) carried by the quark in the pair (ij), and the transverse coordinate of its center of mass The notation˜ n in the first line of (22,23) implies summation over all possible attachments of t-channel gluons to the partons in the upper part of the diagram. For type-A diagrams the variables , n may take independently six different values, which correspond to connections to final quarks, virtual quark or virtual gluon. For type-B diagrams both produced quark pairs must be in a color singlet state, which translates into the additional constraint that , n should be connected to different quark loops (either upper or lower quark-antiquark pairs). The factors σ , σ n in the first line of (22,23) have the value +1 if the corresponding t-channel gluon is connected to a quark line or gluon, and -1 otherwise. On the other hand, the color factors c n depend on the topology of the diagram under consideration; more precisely, how the t-channel gluons are connected to the quark lines. For type-A diagrams, the color factor gluons are connected to the same quark line, or quark and antiquark lines of opposite color (e.g. quark-antiquark lines originating from colorless photon or leading to formation of colorless quarkonium). If the vertices of the t-channel gluons are separated by color changing vertex of a virtual gluon, then the color factor is given by For the diagrams with one 3-gluon vertex, when one of the t-channel gluons is attached to a virtual gluon, the corresponding color factor is c n = ±C 3 = ±N c /4, where the sign is positive for the diagram with attachment of the other t-channel gluon to the upper quark-antiquark pair (i.e. partons 1,2), and negative otherwise. Finally, for the diagram when both t-channel gluons are attached to a virtual (intermediate) gluon, the corresponding factor is c n = C 4 ≡ N c /2. For type-B diagrams, the corresponding color factor is c n = 1 where the summation is done over all final quarks j 1 , ...j n which stem from a given parton. The notations Ψ M1 , Ψ M2 are used for the wave functions of the final state quarkonia M 1 and M 2 (for a moment we disregard completely their spin indices), and ψ (γ) whereas for the type-B contribution it is given bỹ The variables Y ij in (22,23) stand for the lab-frame rapidity of quark-antiquark pair made of partons i, j. Explicitly Figure 2. Examples of higher order contributions, which become relevant for the exclusive production of quarkonia with the same C-parity. The left diagram corresponds to exchange of odderon (3-gluon ladder) in the t-channel, whereas the right diagram corresponds to photon exchange in the t-channel. In both plots it is implied summation over all possible attachments of t-channel gluons and photon (red) to black-colored partonic lines. As explained in the text, both types of contributions are suppressed compared to diagrams from Figure 1 and will be disregarded in what follows.
it is given by where α i and α j are light-cone fractions of the heavy quarks which form a given quarkonium. The dipole amplitude, which appears in (25,26), effectively takes into account a sum of different pomeron ladders [57,60], and for this reason it corresponds to exchange of vacuum quantum numbers in the t-channel. This fact imposes certain constraints on possible quantum numbers of heavy quarkonia produced via the γ + IP → M 1 M 2 subprocess. Since the C-parity of a photon is negative, the neutral quarkonia M 1 , M 2 must have opposite C-parities. This explicitly excludes production of quarkonia with the same quantum numbers (M 1 = M 2 ). For the case when quarkonia are charged (e.g. B + c B − c ), this implies that they necessarily must be produced with odd value of the mutual angular momentum L . Finally, we need to mention that at higher orders the interaction with the target should be supplemented by the exchange of C-odd three-gluon ladders (so-called odderons) in the t-channel [95] potentially giving contributions of odderon exchange, as shown in the right panel of Figure 2. Such interactions are suppressed at high energies, because the odderon has a smaller intercept than the pomeron. Besides, formally such contribution is also suppressed by O (α s (m Q )). Another possibility to produce C-even pair of quarkonia is via exchange of a (C-odd) photon, as shown in the right panel of Figure 2. Formally such contributions are suppressed by ∼ α em /α 2 s (m Q ), which is a small parameter for charm and bottom quarks, yet could get enhanced in the infinitely heavy quark mass limit m Q → ∞ due to suppression of α s (m Q ) in the denominator. Besides, this contribution can be enhanced in the very forward kinematics by the photon propagator ∼ 1/t, where t ≡ (p f − p i ) 2 is very small 2 . According to phenomenological analyses [30][31][32][33][34], the cross-sections of this mechanism numerically is much smaller than that of the mechanism suggested in this paper. For this reason in what follows we will focus on the production of opposite-parity quarkonia, and will disregard the contributions of t-channel odderons and photons altogether.

III. NUMERICAL RESULTS
The framework developed in the previous section is valid for heavy quarkonia of both c and b flavors. In what follows we will focus on the all-charm sector and present results for J/ψ + η c production, for which the cross-section is larger and thus is easier to study experimentally 3 .
2 Numerical estimates show that the invariant momentum transfer t for photoproduction of a pair of quarkonia M 1 , M 2 is restricted by where m N is the mass of the nucleon, W 2 ≡ sγp = (q + P ) 2 , and M 2 12 = p M 1 + p M 2 2 is the invariant mass of the quarkonia pair (clearly, M 12 ≥ M 1 + M 2 ). Already for EIC energies W ∼100 GeV, so we can see that it is possible to achieve the kinematics of very small t even for heavy quarkonia. 3 According to our estimates, for bottomonia the cross-sections are at least an order of magnitude smaller due to the heavier quark mass For the wave function of the J/ψ-mesons we will use a simple ansatz suggested in [96,97], where M is the helicity of J/ψ, r is the distance between the quark and antiquark, h,h are the helicities of the quark and antiquark, and f V , e V , ω are some numerical constants. This result can be trivially extended to the case of η cmeson, which differs from the J/ψ meson only by the orientation of the quark spins. Taking into account the structure of the Clebsch-Gordan coefficients for the 1/2 × 1/2 product, we may immediately write out the corresponding wave functions for η c , modifying the corresponding M = 0 component of the J/ψ wave function, Alternatively, the wave functions of quarkonia can be constructed using potential models or the well-known Brodsky-Huang-Lepage-Terentyev (BHLT) prescription [98][99][100] which allows to convert the rest frame wave function ψ RF into a light-cone wave function Ψ LC . It is known that in the small-r region, which is relevant for estimates, the wave functions of the S-wave heavy quarkonia in different schemes are quite close to each other [101][102][103][104], and for this reason in what follows we will use the ansatz of (28-31), in view of its simplicity. For our numerical evaluations we also need a parametrization of the dipole amplitude. In what follows we will we use the impact parameter (b) dependent "bCGC" parametrization of the dipole cross-section [89,105], We would like to start the presentation of numerical results from a discussion of the relative contribution of type-A and type-B diagrams introduced in the previous section. From the left panel of Figure 3 we can see that the dominant contribution comes from the type-A diagrams. Partially this enhancement can be explained by larger color factors in the large-N c limit. The interference of type-A and type-B contributions represents approximately a 10% correction and moreover, has a node, whose position depends on the produced quarkonia kinematics. As expected, the crosssection is suppressed as a function of p T (we considered p ⊥ J/ψ = p ⊥ η = p T for the sake of definiteness). In the right panel of the same Figure 1 we present the dependence of the yields on the azimuthal angle φ between the transverse momenta of the J/ψ and η c mesons. For definiteness, we assumed that the transverse momenta p ⊥ J/ψ , p ⊥ η of both quarkonia have equal absolute values. In order to make meaningful comparison of the cross-sections, which differ by orders of magnitude, we plotted the normalized ratio R(φ) = dσ (..., φ) /dy 1 dp 2 1 dy 2 dp 2 2 dφ dσ (..., φ = π) /ddy 1 dp 2 1 dy 2 dp 2 2 dφ , We can see that the ratio has a sharp peak in the back-to-back region (φ = π), which happens because in this kinematics the momentum transfer to the target |t| = ∆ 2 is minimal. In contrast, for the angle φ ≈ 0, which maximizes the variable |t| = ∆ 2 , the ratio has a pronounced dip. For p 1 = p 2 the dependence on φ is qualitatively similar, although the maximum and minimum are less pronounced.
In the left panel of the Figure 4 we analyze the p T -dependence, for the case when one of the quarkonia has a small transverse momentum p i ∼ 1 GeV. As expected, in this case the cross-section has a significantly milder suppression compared to the case when both quarkonia share the same transverse momentum. This result indicates that the EIC E p =275 GeV   (36), on the angle φ (difference between azimuthal angles of both quarkonia). The appearance of a sharp peak in back-to-back kinematics is explained in the text. For definiteness we considered the case when both quarkonia are produced at central rapidities (y1 = y2 = 0) in the lab frame; for other rapidities the φ-dependence has a similar shape.  quarkonia pair are predominantly produced with small transverse momenta p ⊥ 1 ∼ p ⊥ 2 1 GeV and opposite directions in the transverse plane ( φ ≡ φ 1 − φ 2 ≈ π). In the right panel of the same Figure 4, we show the p T -dependence of the cross-section in LHeC kinematics. While the absolute value increases in this case, we may observe that qualitatively the dependence on p T and angle φ remains the same.

EIC
In Figure 5 we analyze the dependence of the cross-section on rapidities of the quarkonia. In the left panel we consider the special case when both quarkonia are produced with the same transverse momenta p ⊥ 1 ∼ p ⊥ 2 ∼ 1 GeV and the same rapidities y 1 = y 2 in the lab frame. The variables y 1,2 in this case can be unambiguously related to the invariant photon-proton energy W γp ∼ √ s γp (shown in the upper horizontal axis), and as expected, the cross-section grows as a function of energy. In the right panel of the same Figure 5 we analyze the dependence of the cross-section on the rapidity difference ∆y between two heavy mesons. For the sake of definiteness we consider that both quarkonia have opposite rapidities in the lab frame, y 1 = −y 2 = ∆y/2. We observe that in this case the cross-section becomes suppressed as a function of ∆y, which illustrates the fact that the quarkonia are predominantly produced with the same rapidities.  (14). Right plot: The dependence on the rapidity difference between the produced quarkonia, y1 = −y2 = ∆y/2. For the sake of definiteness we assumed that both quarkonia are produced at central rapidities (y1 = y2 = 0) with transverse momenta p1 = p2 = 1 GeV in the lab frame.
Finally, in Figures 6, 7, 8 we show the results for the cross-section dσ γp→M1M2p /dy 1 dy 2 , which is integrated over transverse momenta p ⊥ i of both quarkonia. This observable can be the most promising for experimental studies, since it is easier to measure. We make the predictions in the kinematics of the ultraperipheral pA collisions at LHC, as well as future electron-hadron colliders. Largely, the dependence on y 1 , y 2 repeats similar dependence of the p T -unintegrated cross-sections. This happens because the p T -integrated cross-sections get its dominant contributions from the region of small p T m Q , where dependence on rapidity is mild. In Figures 6, 7 we have also shown the cross-sections of "master" processes ep → eM 1 M 2 p and Ap → AM 1 M 2 p. The expressions for these cross-sections differ from those of γp → M 1 M 2 p by a convolution with known kinematic factors, which correspond to fluxes of equivalent photons generated by the electron or heavy nucleus. These cross-sections have completely different behavior on the rapidity y 1 = y 2 of both quarkonia, which can be understood from (8)(9)(10)(11)(12)(13)(14)(15)(16). Indeed, mesons with higher lab-frame rapidities can be produced by photons of higher energy E γ , yet the flux of equivalent photons created by a charged electron or ion is suppressed and vanishes when the elasticity y = E γ /E e approaches unity. Finally, Figure 8 illustrates how the cross-section behaves as a function of y 1 , y 2 in general, when |y 1 | = |y 2 |. We can see that the cross-section has a typical ridge near y 1 ≈ y 2 , i.e. when quarkonia are produced with approximately the same rapidities.

IV. CONCLUSIONS
In this paper we studied in detail the exclusive photoproduction of heavy charmonia pairs. This process presents a lot of interest, both on its own, as a potential test of quarkonia production mechanisms in small-x kinematics, as well as a background to exotic hadron production. We analyzed in detail the leading order contributions and found that in this mechanism the quarkonia pairs are produced with opposite C-parities, relatively small opposite transverse momenta p T , and small separation in rapidity. This finding is explained by the fact that in the chosen kinematical region the momentum transfer to the recoil proton is minimal. As expected, the cross-section decreases rapidly as a function of p T , and grows as a function of photon-proton invariant energy (∼quarkonia rapidities), similar to single-photon production. However, the cross-section decreases as a function of the rapidity difference between the quarkonia. We estimated numerically the cross-section in the kinematics of ultraperipheral pA collisions at LHC, as well as in the kinematics of the future electron-proton colliders, and found that the cross-section is sufficiently large for experimental studies. Our evaluation is largely parameter-free and relies only on the choice of the parametrization for the dipole cross-section (32) and wave functions of quarkonia.
We need to mention that earlier studies of exclusive production focused on production of quarkonia pairs with the same quantum numbers (e.g. J/ψ J/ψ). In view of different quantum numbers, this process predominantly proceed via exchange of two photons at amplitude level, like e.g. via photon-photon fusion γγ → M 1 M 2 [30][31][32][33][34] or double photon scattering [46]. Due to extra virtual photon in the amplitude, the cross-sections of such processes  For UPC collisions the positive direction of rapidity is that of a heavy lead ion, and the cross-sections are given per nucleon. The solid curves correspond to cross-section of the γp → M1M2p subprocess, whereas dotted lines correspond to the cross-sections of the complete physically observable ep or Ap processes. We assume for definiteness that the rapidities of both quarkonia are equal to each other in the lab frame, y1 = y2 = y. The upper horizontal scale illustrates the corresponding value of invariant energy W ≡ √ sγp, as defined in (14).
are parametrically suppressed by ∼ α 2 em compared to the cross-section of opposite C-parity quarkonia, and thus numerically are significantly smaller. We hope that the process suggested in this paper will be included in the program of the future EIC collider, as well as ongoing studies at LHC in ultraperipheral kinematics.
Finally, we need to mention that it is quite straightforward to extend the framework developed in this manuscript to the case of all-heavy tetraquark production: for this it is only necessary that the product of final state quarkonia wave functions in (22,23) be replaced with the wave function of the tetraquark state. Estimates of the cross-sections for this case will be presented in a separate publication.  The dependence on rapidity difference for the pT -integrated cross-section, in the kinematics of ultraperipheral collisions, at LHC and future electron-proton colliders. The positive sign of rapidity is chosen in the direction of electron or emitted quasi-real photon. For UPC collisions the positive direction of rapidity is that of a heavy lead ion, and the cross-sections are given per nucleon. For the sake of definiteness we assume that in the lab frame the quarkonia have opposite rapidities, y1 = −y2 = ∆y/2. The upper horizontal scale illustrates the corresponding value of the invariant mass M12 ≡ p J/ψ + pη c 2 , as defined in (15). Dotted curves correspond to the cross-sections of the complete process (electron-proton or heavy ion-proton).

infrastructure of the NLHPC (ECM-02).
Appendix A: High energy scattering in the color dipole picture In this appendix for the sake of completeness we briefly remind the general procedure which allows to express different hard amplitudes in terms of the color singlet forward dipole scattering amplitude. While in the literature there are several equivalent formulations [47,[49][50][51][52][53][54][55], in what follows we will use the Iancu-Mueller approach [106].
The natural hard scale, which controls the interaction of a heavy quark with the gluonic field, is its mass m Q . In the heavy quark mass limit we may formally develop a systematic expansion over α s (m Q ) 1. Furthermore, for small color singlet dipoles there is an additional suppression by the dipole size, r ∼ 1/m Q , so the interaction of singlet dipoles with perturbative gluons is suppressed at least as ∼ α s (m Q ) /m Q . However, the interaction of gluons with each other, as well as with light quarks, remains strongly nonperturbative in the deeply saturated regime, so we expect that the dynamics of the dipole amplitudes should satisfy the nonlinear Balitsky-Kovchegov equation.
At very high energies the dynamics of partons can be described in the eikonal approximation. The transverse coordinates of the high energy partons remain essentially frozen during its propagation in the gluonic field dipole of the target. Similarly, due to eikonal interactions we may disregard completely the change of the quark helicities. In this picture the interaction of a dipole with the target is described by the S-matrix element [60,106] S(y, where we use the notation y = ln(1/x) for the dipole rapidity, x Q , xQ, are the transverse coordinates of the partons (quark or antiquark), and the factors V † (x Q ) and V (xQ) in (A1) are the Wilson lines, which describe the interaction of the partons with the color field of a hadron. They can be expressed as where A a µ is the gluonic field in a hadron. The impact parameter dependent dipole amplitude N (x, r, b) can be related to S y, x Q , xQ as where the variable r ≡ x Q − xQ is the transverse size of the dipole, b ≡ z x Q + (1 − z)xQ is the transverse position of the dipole center of mass, and z is the fraction of the light-cone momentum of a dipole which is carried by the quark Q. In view of the weakness of the interaction between heavy quarks and gluons, we can make an expansion of the exponent in (A2) over α s (m Q ). In this approximation the effective interaction of the quark or antiquark with the gluonic field of the proton can be described by the factor ±i t a γ a (x ⊥ ), where x ⊥ is the transverse coordinate of the quark, and t a are the ordinary color group generators of pQCD in the fundamental representation. Inspired by the color structure of the interaction, in what follows we will refer to these interactions as "exchanges of t-channel pomeron (gluons)", tacitly assuming that it can include cascades (showers) of particles. For the dipole scattering amplitude (A3), using (A1, A4), we obtain For further evaluations it is more convenient to rewrite this result in the form where we defined a shorthand notation ρ (x a ) ≡ |γ a (x)| 2 , and r 12 , b 12 are the distance and center-of-mass of the quark-antiquark pair located at points x 1 , x 2 . For many processes the contributions ∼ ρ (x i ) cancel, so the amplitude eventually can be represented as a linear superposition of the dipole amplitudes N (x, r, b). In what follows, we will see that the amplitude of the process considered in this manuscript can be represented as a bilinear combination of terms with structure ∼ [γ (x i ) − γ (x j )]. For this special case the substitution of (A6) allows to get a few important identities between bilinear expressions where r ij and b ij are the relative distance and center-of-mass of the quark-antiquark pair located at points x i , x j . For the impact parameter independent (b-integrated) cross-section the results (A5-A7) can be rewritten in a simpler form, The value of the constant term in the right-hand side of (A11) is related to the infrared behavior of the theory, and for the observables which we consider in this paper, it cancels exactly. In what follows we will apply this formalism to the evaluation of the exclusive dimeson production amplitudes.

Appendix B: Evaluation of the photon wave function
For evaluation of the photon wave function we follow the standard rules of the light-cone perturbaiton theory formulated in [16,107]. The result for theQQ component is well-known in the literature [96,108], yet below in Subsection B 1 we will briefly repeat its derivation in order to introduce notations. As we will see later in Subsection B 2, the wave function of theQQQQ-component can be expressed in terms of the wave function ofQQ-component. In Figure 9. Left plot: The leading order contributions to theQQ-component of the photon wave function ψ g→QQ . Right plot: the so-called gluon emission wave function, as defined in [109]. The momenta ki shown in the right-hand side are Fourier conjugates of the coordinates xi.
our evaluation we will focus on onshell transversely polarized photons, which give the dominant contribution, unless some specific cuts are imposed on its virtuality Q 2 . The momentum of the photon (1) introduced earlier simplifies in this case and has only light-cone component in the plus-axis direction, The polarization vector of the transversely polarized photon is given by where in (B2) we took into account that q ⊥ = 0. Before the interaction with the target, the photon might fluctuate into virtual quark-antiquark pairs, as well as gluons. In what follows we will use a convenient shorthand notation α i = k i /q + for the fraction of light-cone momentum of the photon carried by each parton, as well as k i⊥ for the transverse component of parton's momentum. In view of 4-momentum conservation we expect that α i , k i⊥ should satisfy an identity where summation is done over all partons. We may observe that the vector ε γ satisfies an identity and its scalar product with any 2-vector a yields ε γ · a = a x + iγa y √ 2 = |a| √ 2 e iγ arg(a) , arg(a) = arctan a y a x . (B6)

1.QQ component of the photon wave function
In this section for the sake of completeness we would like to remind the reader the main steps in the derivation of theQQ-component photon wave function [96,108] in the mixed (α, r) representation. In leading order the subprocess γ →QQ gets contributions only from the diagram shown in the left panel of the Figure 9. A bit later we will see that γ →QQ, as well as the closely related g →QQ subprocess, appears as constituent blocks in the more complicated 4-quark wave function. For this reason, in order to facilitate further discussion, temporarily in this section we will assume that the photon momentum q might has a nonzero transverse part q ⊥ , and will use notation z = k + 1 /q + for the fraction of light-cone momentum carried by the quark. In momentum space the evaluation is straightforward, using the rules from [16,107,109] and yields where λ is the helicity of the incoming photon, h,h, are the helicities of the produced quark and antiquark, c,c, are the color indices of Q andQ, respectively, and e q is the electric charge corresponding to a given heavy flavor. The momentum n is defined as n = k 1 − zq ⊥ = (1 − z)k 1 − zk 2 and physically has the meaning of the transverse part of the relative (internal) momentum of the QQ pair. The numerator of (B7) can be written out explicitly using the rules from [96,108], In configuration space the corresponding wave function can be found making a Fourier transformation over the transverse momenta, where the integral over k 2 was performed using the properties of the δ function, and before the integration over k 1 we shifted the integration variable as k 1 → n + zq. Explicitly, the integration over the variable d 2 n yields The structure of the Eq. (B10) clearly suggests that in a mixed representation the variable zr 1 +zr 2 plays the role of the dipole center of mass, whereas r 12 is its separation, in agreement with earlier findings from [110]. For the incoming offshell photon with virtuality −q 2 = Q 2 , straightforward integration yields in a similar fashion in the second line of (B10). The extension of this result for the production of a QQ pair by a gluon is straightforward and requires a simple replacement e q δ cc → g (t a ) cc . Finally, we would like to discuss briefly the so-called parton-level wave function of the gluon emission subprocess q → gq, as introduced in [109]. This object is useful for the analysis of different amplitudes, as we will see in the next section. In the leading order it gets contributions from the diagram shown in the right panel of the Figure 9. The evaluation of this object are quite similar to the derivation of (B7-B11). In momentum space we obtaiñ where λ is the helicity of the outgoing gluon; (h i , c i ) and (h f , c f ) are the helicities and color indices of the incident and final quark (before and after emission of a gluon); and similar to the previous case we have introduced the momentum n = k 1 − zq = (1 − z)k 1 − zk 2 which corresponds to the relative motion of the quark and gluon after emission of the latter. Using the rules from [96,108], we may rewrite the numerator as In configuration space the corresponding wave function is given by Figure 10. The leading order contribution to the wave function ψ where the integral over k 2 was performed using the properties of the wave function, and integration over the variable as k 1 = n + zq yields For the case of incoming offshell quark with virtuality Q 2 , straightforward generalization show that the second line of (B15) gets the form Similarly to the previous case, the structure of the Eq. (B16) clearly suggests that in a mixed representation the variable zr 1 +zr 2 plays the role of the dipole center of mass, whereas r 12 is its separation [110].

2.QQQQ component of the photon wave function
As was mentioned earlier in Section A, in the eikonal approximation the amplitude of the subprocess γ * →QQQQ in configuration space can be represented as a convolution of the wave function ψ (γ) QQQQ with linear combinations of dipole amplitudes (A7). In leading order over α s the amplitude of the process is given by the two diagrams shown in the Figure (10). It should be understood that these diagrams should be supplemented by all possible permutations of final state quarks. More precisely, for the production of different heavy flavors (e.g.ccbb) both diagrams should be supplemented by contributions with permuted pairs of momenta (k 1 , k 2 ) ↔ (k 3 , k 4 ). For the same-flavor quarkonia pairs (e.g.cccc) we should take into account contributions with independent permutations of the quarks and antiquarks, k 1 ↔ k 3 and k 2 ↔ k 4 . The evaluation of the corresponding process follows the standard lightcone rules formulated in [16,107]. We need to mention that some blocks, which will be needed for the construction of the amplitudem have already been evaluated in [109,111] (although in the chiral limit only). In these section we extend those studies and represent them in a form convenient for further analysis. According to the general light-cone rules [96,108], in the evaluation of the diagrams in Figure ( In leading order over α s , the amplitude of the process is given by the two diagrams shown in Figure (10) and depends on the momenta of the 4 quarks in the final state. In what follows we will use the standard notation α i = k i /q + for the fractions of photon momentum carried by each of these fermions, as well as k i⊥ for the transverse components of their momenta. We also will use a shorthand notation = k 3 + k 3 for the momentum of the virtual gluon connecting different quark lines. For the sake of generality we will assume that the produced quark-antiquark pairs have different flavors, and will use the notations m 1 for the current mass of the quark line connected to a photon, and m 2 for the current masses of the quark-antiquark pair produced from the virtual gluon.
Using the rules from [96,108], we may obtain for the corresponding amplitude of the subprocess where a i and c i are the helicity and color indices of the final state quarks, and D ij are the conventional light-cone denominators (with the first subscript index i = 1, 2 refers to the first and the second diagram in Figure 10 respectively, and the second index j = 1, 2, 3 numerates proper cuts shown with dashed vertical lines). Explicitly, these light-cone denominators are given by To simplify the structure of the expressions (B19-B23), we introduced a shorthand notationᾱ i ≡ 1 − α i , i = 1...4. The combination of momenta q 34 , defined in (B24), represents the relative motion momenta of quarks 3 and 4 (Fourier conjugate of a relative distance r 3 − r 4 ). Technically, the structure of the denominators, up to trivial redefinitions, agrees with the findings of [109]. The expressions in the numerator of (B18) can be written out explicitly using the light-cone algebra from [16,107,109], yielding for the amplitude where the momenta q i are defined as We may observe that the amplitude (B25) is antisymmetric with respect to permutation of the momenta and helicities of the first two quarks, (α 1 , k 1 , a 1 ) ↔ (α 2 , k 2 , a 2 ), and symmetric with respect to permutation of the the momenta and helicities of the 3rd and 4th quarks, (α 3 , k 3 , a 3 ) ↔ (α 4 , k 4 , a 4 ). This symmetry simply reflects that the amplitude (B25) was evaluated as a sum of the left and the right diagrams in Figure 10, which can be related by charge conjugation. This symmetry allows to simplify some evaluations.
For evaluations in the dipole framework we need to rewrite the amplitude in configuration space, making a Fourier transformation over the transverse components, In view of momentum conservation (B4), the wave function ψ (γ) QQQQ will be invariant with respect to global shifts i.e. should depend only on relative distances between quarks |x i − x j |. After straightforward evaluation of the integrals and algebraic simplifications it is possible to reduce (B27) to the form Figure 11. Graphical illustration of the transverse momentum dependence of the wave function ψ γ→QQQQ ({αi, ri}). The letters bij and b ijk stand for the center of mass position of the partons ij or ijk. See the text for more details.
The variable b j1...jn corresponds to the position of the center of mass of n partons j 1 , ...j n and was defined earlier in (24). The variables n i,j1...jn = (x i − b j1...jn ) / |x i − b j1...jn | are unit vectors pointing from quark i towards the center-of-mass of a system of quarks j 1 ...j n . It is not possible to do the remaining integrals over q 1 , k 2 analytically, nor present the wave function (B30) as a convolution of simpler "elementary" wave functions from Section B 1.. Technically, this happens because in the language of traditional Feynman diagrams the intermediate (virtual) partons are offshell, and the integration over q 1 , k 2 can be rewritten via integrals over virtualities of intermediate particles.
Nevertheless, the structure of the coordinate dependence of ψ γ→QQQQ ({α i , r i }) can still be understood using the simple rules suggested in Section B 1. Indeed, in the eikonal picture the transverse coordinates of all partons are frozen. The tree-like structure of the leading order diagrams 1, 2, in Fig. 10 and the iterative evaluation of the coordinate of the center of mass of two partons b ij = (α i r i + α j r j ) / (α i + α j ) allows to reconstruct the transverse coordinates of all intermediate partons, as shown in Figure 11. The variables r 1 − b 34 and r 2 − b 34 have the physical meaning of the relative distance between the recoil quark or antiquark and the emitted gluon. Similarly, the variables r 1 − b 234 and r 2 − b 134 can be interpreted as the size of theQQ pair produced right after splitting of the incident photon. These simple rules allow for the construction of the heavyQQQQ production amplitude in the gluonic field of the target.
The wave function ψ (γ) QQQQ ({α i , r i }) has a few singularities which require special attention in order to guarantee that the amplitudes of the physical processes remain finite. For the meson pair production, the choice of the quarkonia wave functions (28)(29)(30)(31), which vanish rapidly near the endpoints is sufficient in order to guarantee finiteness of the amplitudes (22)(23). In what follows we will refer to the diagrams in the first row as A1, B1, and the diagrams in the second row as A2, B2, respectively.

b. Instantaneous contributions
According to canonical rules of the standard light-cone perturbation theory [16,107], the evaluations from the previous section should be supplemented by the instantaneous contributions of virtual partons. The propagators of the instantaneous offshell quarks and gluons with momentum k are given by where n µ is the light-cone vector in minus-direction. The results for the instantaneous contributions of gluons are quite straightforward to get, essentially repeating the evaluations from the previous subsection. Since γ + γ + = 0, there is no diagrams with two instantaneous propagators (quark and gluon) connected to the same vertex. The numerators of amplitudes with instantaneous propagators have simple structure in view of identities [16,107,109] where the subscript indices q, g in the right-hand side denote the parton propagator, which should be taken instantaneous (q for quark, g for gluon), and Figure 13. Schematic illustration of the diagrams which contribute to a γ →QQQQ subprocess, via single-gluon exchange in t-channel. For the sake of simplicity we omitted a proton blob in the lower part. The square box with gluon connected in the middle stands for a coupling of a dipole (sum of the couplings to all partons which pass through the block, ∼ (±)γ(xi)ta). The center-of-mass bi 1 ...in of a system of partons i1...in is defined in (24). In all plots it is implied the inclusion of diagrams which can be obtained by inversion of heavy quark lines ("charge conjugation").
where the corresponding contributions A 1,2,3 are given by We may observe that all factors γ c (x i ) always appear in combination γ c (x i ) − γ c (x j ), which guarantees that in the heavy quark mass limit, when the distances between the quarks are small, the corresponding amplitude is suppressed at least as ∼ 1/m Q . The three-gluon coupling ∼ γ (x 34 ) always appears in combination [γ c (x 1 ) + γ c (x 2 ) − 2γ c (x 34 )] , in agreement with the earlier findings of [54].