Inclusive photoproduction of charmonia-bottomonia pairs

In this preprint we analyze the inclusive photoproduction of heavy charmonia-bottomonia pairs in the Color Glass Condensate framework and demonstrate that the cross-section of the process is sensitive to dipole and quadrupole forward scattering amplitudes (2- and 4-point correlators of Wilson lines). Using the phenomenological parametrizations of these amplitudes, we estimate numerically the production cross-sections in the kinematics of the forthcoming Electron Ion Collider and the ultraperipheral collisions at LHC. We found that the contribution controlled by the quadrupole amplitude is dominant, and for this reason, the suggested channel can be used as a gateway for studies of this nonperturbative object.


I. INTRODUCTION
Due to their heavy masses and small sizes, the quarkonia have been considered as important probes of the gluonic field of the target almost since their discovery [1,2].The modern NRQCD framework allows to evaluate systematically various perturbative corrections to processes which include quarkonia [3][4][5][6][7][8][9][10][11][12][13].This framework allows to get reasonable estimates for the cross-sections, although at present it includes some uncertainties in nonperturbative Long Distance Matrix Elements (LDMEs) of quarkonia states [10][11][12][13], as well as apparent mismatch of the LDMEs of different quarkonia which challenge the expectations based on heavy quark mass limit [13][14][15][16][17]. Furthermore, even in the heavy quark mass limit the production of single quarkonia provides only limited information about the gluonic field.For this reason, since the early days of QCD the theoretical efforts were dedicated to production mechanisms of multiple quarkonia in the final state [18][19][20][21], which could provide a significantly more exhaustive information about the gluonic content of the target.Recent experiments at LHC demonstrated a feasibility to study the production of quarkonia pairs experimentally.
In this paper we will focus on the associated production of charmonia and bottomonia, which has a very simple structure of the partonic amplitude.The theoretical studies of such processes (particularly, J/ψ +Υ hadroproduction) have been initiated in [22][23][24] and attracted a lof of theoretical attention as possible probes of the gluonic field of the target.In the collinear and k T factorization approach, such processes give a possibility to study the so-called double parton distribution functions (DPDFs), which encode the information about correlation of different constituents in the target.The recent experimental data [25,26] demonstrated that the double parton distributions might give significant contributions in this channel.However, due to a large number of different mechanisms in hadroproduction, these contributions have sophisticated structure, which complicates interpretation of experimental data in terms of partonic content of the target.At present the situation remains unclear: while inclusion of the DPDFs allows to explain the difference between the data and predictions based on single-parton scattering, the value of the so-called effective cross-section σ eff (parameter which controls the magnitude of the DPDFs) depends significantly on the channel used for its extraction [27][28][29].In view of these difficulties of hadroproduction channel, the photoproduction of the same states can serve as a simpler alternative for testing our understanding of the contributing mechanisms.Previously, the studies of the double quarkonia photoproduction mostly focused on exclusive channels [30][31][32][33][34][35][36][37][38][39].However, the inclusive production of such states also deserves interest as a potential gateway for studies of the gluonic distributions.Nowadays such processes may be studied in ultraperipheral pA collisions at the LHC, provided that the recoil nucleus is separated by large rapidity gap from the other hadrons.The use of heavy ions allows to boost significantly the flux of equivalent photons, and in this way facilitates measurement of tiny cross-sections.In the future, such processes may also be studied in electron-proton collisions at the Electron Ion Collider (EIC) [40,41], the Large Hadron electron Collider (LHeC) [42], and the Future Circular Collider (FCC-he) [43][44][45].
In the kinematics of all the above-mentioned experiments, due to enhanced role of multiparton distributions, which eventually lead to onset of saturation effects, it is appropriate to analyze this process in the frameworks with built-in saturation.In what follows we will use for our studies a Color Glass Condensate (CGC) framework [46][47][48][49][50][51][52][53][54][55][56], which naturally incorporates the saturation effects and provides a phenomenologically successful description of both hadronhadron and lepton-hadron collisions [57][58][59][60][61][62][63][64].The cross-sections of physical processes in this framework are expressed in terms of the forward multipole scattering amplitudes (n-point correlators of Wilson lines), which have probabilistic interpretation and present important physical characteristics of the target.Technically, these correlators generalize the multigluon distribution functions used in collinear and k T factorization picture, providing correct asymptotic behaviour in the high-energy limit.
The paper is structured as follows.Below in Section II we briefly describe the main components of the CGC framework and then present the theoretical results for the photoproduction of heavy quarkonia pairs in the CGC approach.In Section III we provide numerical estimates using the phenomenological parametrizations of dipole and quadrupole amplitudes.Finally, in Section IV we draw conclusions.

II. THEORETICAL FRAMEWORK
Nowadays, the photoproduction processes may be studied in electron-proton, proton-proton and proton-nuclear collisions in ultraperipheral kinematics (UPC).The corresponding cross-sections of these processes are related to the photoproduction cross-section as where in (1) we use standard DIS notations, in which y is the inelasticity (fraction of electron energy which passes to the photon), q and P are the 4-momenta of the photon and proton, and y a , k ⊥ a , with a = 1, 2, are the rapidities and transverse momenta of the produced quarkonia with respect to electron-proton or hadron-hadron collision axis.The notation dσ T is used for the cross-section of the photoproduction subprocess, induced by a transversely polarized photon.The expression dn γ ω ≡ E γ , q ⊥ in (2) is the spectral density of the flux of equivalent photons with energy E γ and transverse momentum q ⊥ with respect to the nucleus, which was found explicitly in [65].The momenta p ⊥ a = k ⊥ a − q ⊥ are the transverse parts of the quarkonia momenta with respect to the produced photon.Due to nuclear form factor of the recoil nucleus, the spectral density dn γ ω, q ⊥ is strongly suppressed for the transverse momenta q ⊥ larger than the inverse nuclear radius R −1 A , so the average values of q ⊥ are quite small, , and the transverse momentum dependence of the cross-sections in the left-hand side of (2) almost coincides with the p ⊥ -dependence of the photoproduction cross-section in the right-hand side, and for this reason from now on we will tacitly assume that p ⊥ a ≈ k ⊥ a .For further evaluations we need to fix the reference frame and write out explicit light-cone momenta decompositions of the participating hadrons.In what follows we will use notations p 1 , p 2 for the 4-momenta of produced heavy quarkonia, and P for the momentum of the proton.We'll work in the reference frame in which the light-cone expansion of these 4-vectors is given by where m N is the mass of the nucleon, and M 1 , M 2 are the masses of produced quarkonia.In what follows we will be mostly interested in the high-energy collider kinematics q + , P − ≫ {Q, M a , m N }, when quarkonia are produced with relatively small transverse momenta.In this kinematics it is possible to use eikonal picture and assume that the plus-component q + of the photon light-cone momentum is shared only between the produced quarkonia, namely The invariant energy W of the γp collision and the invariant mass M 12 of the produced heavy quarkonia pair are given by and respectively.The photoproduction cross-section dσ T , which appears in (1,2), is the central quantity of interest for present study and will be evaluated in the Color Glass Condensate approach [46][47][48][49][50][51][52][53][54][55][56].In the following subsection II A we briefly remind the main assumptions of this theoretical framework, and in subsection II B we present final theoretical results of the evaluation.To avoid possible soft factorization breaking final state interactions, in what follows we will tacitly assume that the produced quarkonia are kinematically well-separated from each other, namely that the invariant relative velocity is sufficiently large, or equivalently formulated in terms of the invariant mass of hadron pairs, A. High energy scattering in CGC picture For heavy quarkonia production, the heavy quark mass m Q plays the role of the natural hard scale, which determines the interaction strength of heavy quarks with the gluonic field.In the heavy quark mass limit it is formally possible to develop a systematic expansion over the strong coupling α s (m Q ) ≪ 1.However, such expansion might be not very reliable in the small-x limit, when the gluon fields are enhanced due to saturation effects and reach values A a µ ∼ 1/α s .The interaction of the partons with the ultrarelativistic target (shockwave) in this kinematics is characterized by nearly instantaneous color exchange, which can be described in the eikonal approximation.In this picture the interaction of heavy quark with the target is described by a Wilson line where x ⊥ is the impact parameter of the parton, t a are the ordinary color group generators of pQCD in fundamental representation, and is the gluonic field in the target created by the color charge density ρ a .
According to classical CGC picture [47][48][49], for multiparton scattering processes each physical observable O should be averaged with a weight function W [ρ] which describes the probability of a given charge distribution ρ inside the target, namely where from now on we will use angular brackets ⟨...⟩ for such averaging.The analytic evaluation of the path integral Dρ is possible only for a few simple forms of W [ρ] (e.g. for gaussian).Fortunately, for many observables the averages ⟨O⟩ may be expressed in terms of a limited set of universal (process-independent) amplitudes, such as dipole or quadrupole forward scattering amplitudes, and thus avoiding explicit evaluation of the integral over all possible charge configurations.For example, the amplitude of the heavy quark pair or (single) quarkonia production from gluon may be described by the two diagrams shown in the Figure (see for details [66][67][68] ) and is given by where z represents the coordinate of the interaction point in the configuration space, and the corresponding parton fields are given by ū (p q , z) = 1 2 Diagrams which describe the inclusive heavy quark pair production in CGC picture in the leading order over αs.The subscript numbers 1,2 enumerate the heavy quark lines in our convention.The subscript letters z is the coordinate of the interaction vertex in the configuration space.The red block represents the interaction with the target (shockwave).
for the produced quarks, for produced antiquark, and for the initial-state gluon, where the vector n 2 in (18) is the unit light-cone vector pointing in the direction of the target momentum.The notation U ba is used for the components of the matrix (13) in adjoint representation.The expressions in square brackets in (16)(17)(18) merely reflect that corresponding U -matrices should be taken into account only if the colored parton exists at the moment of interaction with the shockwave of the target.In the limit U → 1 (no interaction) the amplitude (15) reduces to a mere perturbative spinors and gluon polarization vector.In what follows we will work in the mixed representation, describing the kinematics of each parton in terms of its light-cone momenta of the partons together with transverse coordinates in configuration space.Since in the eikonal approximation the interaction of any parton with the target is described by multiplicative factor U (X p ) in proper representation of the color group, evaluation of the amplitudes and cross-sections becomes straightforward.For a simple process shown in (1), the amplitude of the process is described by the S-matrix element [52,53,60] where the notation Y is used for the dipole rapidity.The dipole scattering amplitude N (x, r, b) can be related to S 2 (Y, x q , x q ) as where the variable r ≡ x q − x q is the transverse size of the dipole, and b ≡ α q x q + α q x q is the transverse position of the dipole's center of mass.In complete analogy we may introduce the multipole scattering S-matrix elements (correlator of Wilson lines) which characterize scattering amplitudes of 2n-particle ensemble of quarks.Each such correlator represents independent characteristics of the target and characterizes gluon distributions in the target.

Dilute scattering limit
If the saturation effects are not very large, the interaction of heavy quarks with the target can be described perturbatively, making a Taylor expansion of the Wilson line (13), The dipole scattering amplitude (20) in this limit may be expressed as where we defined γ a as For further evaluations it is convenient to rewrite (24) in the form where we introduced a shorthand notation Γ (x) ≡ |γ a (x)| 2 , and r 12 , b 12 are the transverse size and impact parameter of the color dipole.For certain processes which involve partonic ensembles with net zero color charge, the contributions ∼ Γ (x i ) cancel, so the cross-sections eventually can be represented entirely as a linear superposition of the dipole amplitudes N (x, r, b).The multipole scattering amplitudes (21) in this limit may be rewritten as where c ijkℓ are some numerical coefficients.Using the identities (25), it is possible to express (26) as a linear superposition of the forward color dipole amplitudes.However, at higher orders in α s (m Q ) there are contributions which include multiple products of γ (x k ), and which cannot be reduced to a mere superposition of dipole amplitudes, demonstrating independence of the multipole contributions in general case.

B. Inclusive photoproduction of meson pairs
The inclusive heavy quarkonia pair production in general may proceed via a number of different mechanisms, and for this reason presents a sophisticated problem.In what follows we will focus on the production of charmonia-bottomonia pairs with opposite C-parity, in order to have vacuum quantum numbers in the t-channel.In the leading order of the perturbation theory at the amplitude level this process is described by a set of diagrams shown in the Figure 2.For comparison, the production of quarkonia pairs in general case requires inclusion of additional Feynman diagrams shown in the Figure 3, or emission of additional hard gluons, leading to significantly more complicated expressions for the cross-sections.
The evaluation of the diagrams shown in the Figure 2 is straightforward.The diagrams in the left and right rows might be related to each other by inversion of the quark line in the upper loop, and for this reason for quarkonia states with definite C-parity give the same (up to a sign) contributions.According to NRQCD, the quarkonia formation can proceed both from color singlet and color octet Q Q pairs, so the total cross-section may be represented as an incoherent sum dσ dΩ h = dσ (1)  dΩ h + dσ (8)  dΩ h , where we introduced a shorthand notation dϕ for the phase volume, ϕ is the azimuthal angle between the transverse momenta of the quarkonia, and dσ (1) , dσ (8) stand respectively for the color singlet and The dominant color singlet contribution comes from the diagrams in the last row (diagrams "5" and "6" respectively).The subscript numbers 1-4 enumerate the heavy quark lines in our convention.The subscript letters z1 ... z3 stand for the coordinates of the interaction vertices in the configuration space.The red block represents the interaction with the target (shockwave). .Additional diagrams which in general should be taken into account for inclusive quarkonia pair production in CGC picture in the leading order over αs, yet do not contribute to charmonia-bottomonia pairs.The subscript numbers 1-4 enumerate the heavy quark lines in our convention.The subscript letters z1...z3 stand for the coordinates of the interaction vertices in the configuration space.The red block represents the interaction with the target (shockwave).
octet contributions.The diagrams 1-4 lead to formation of QQ pair (partons "3", "4") in color octet state, whereas the remaining two diagrams lead to formation of both color singlet and color octet pairs.The evaluation of the color octet contribution, shown in diagrams 1-4 in the upper two rows may be interpreted as a production of the first quarkonium M 1 via γ * → M 1 g subprocess, with subsequent fragmentation of the emitted hard gluon into the second quarkonium M 2 .Indeed, we may observe that the gluon connecting quark lines of different heavy flavors is perturbative, and there is no instantaneous contributions to its propagator.For this reason, the contribution of diagrams 1-4 to the cross-section may be rewritten as dσ (8) where D g→M2 (z) is the fragmentation function of the gluon to the second quarkonium M 2 [69,70], and z is the fraction of the gluon momentum carried by the second quarkonium.In the leading order over α s (m Q ), the fragmentation function is given by where is the corresponding color octet matrix element for the quarkonium M 2 , however at higher orders in α s (m Q ), due to emission of soft gluons, the function D g→M2 acquires nonzero finite width, and its peak shifts towards smaller z < 1 [70].The evaluation of the cross-section dσ γ * →M1g may be found in [71] and largely repeats similar evaluation of the γ * → QQg subprocesses from [66,72].The final result of this evaluation reads as [66,72] where b 12 = (x 1 + x 2 ) /2 is the impact parameter of the quarkonium M 1 , r = x 1 − x 2 is the transverse distance between the quarks in the quarkonium, "primed" variables b ′ , z ′ , r ′ correspond to similarly defined variables in the conjugate amplitude, and the corresponding differential measure of integration dΠ i (...) is defined as The explicit expressions for the matrices Ξ [c] i, j and Γ λ i, j may be found in [71] and will be omitted here for brevity.According to modern estimates [10,12,69,[73][74][75], the color octet matrix elements are very small and constitute ≲ 10 −2 of the color singlet matrix elements.For this reason, a gluon fragmentation mechanism constitutes just a few percent correction and may be disregarded in the first approximation, especially in the kinematics of small transverse moments p ⊥ which we study in this paper.
In what follows we will focus on the contribution of the last two diagrams 5, 6 in the Figure 2 which represent the dominant color singlet contributions.In evaluations we may use the fact that the final-state quarks are nearly onshell, and in eikonal picture the free quark propagator becomes diagonal in transverse coordinates in configuration space, namely (see details in [76]).The corresponding cross-section in this limit may be rewritten as a convolution of wave functions and target-dependent multipole contributions, namely: , where h i and h i are the helicities of heavy quarks in the amplitude and its conjugate (summation is implied), ψ (γ)h1h2h3h4 QQ QQ (x 1 ; x 2 ; x 3 ; x 4 ; q) is the wave function which describes the 4-quark QQ QQ Fock state inside the quasireal photon (see details in Appendix A).We also use a shorthand notation for the conventional NRQCD projectors (in helicity basis) multiplied by the appropriate Long Distance Matrix Elements (LDMEs) of the corresponding meson M .In the heavy quark mass limit the largest LDMEs are associated with the 3 S 1 , 1 S 0 and 1 P 0 states, for which the corresponding projectors are given by [5,6,8,77] where P M is the 4-momentum of the corresponding quarkonium, m Q are the quarkonia masses, ⟨O M ⟩ are the corresponding long-distance matrix elements.We also introduced a shorthand notation for the density matrix which describes interaction of the 4-quark ensemble with the target , with subsequent formation of two color singlet dipoles.The angular brackets ⟨...⟩ denote averaging over the color sources, as defined in (14).
As we demonstrate in Appendix B, after convolution with color generators and color averaging, it is possible to rewrite (37) as which includes only dipole and quadrupole correlators.In the usual large-N c limit the averaging may be applied separately to different terms in the product in (38), so eventually the cross-section probes the quadrupole and dipole forward scattering amplitudes.We need to mention that the combination of quadrupole and dipoles, ), which appears bilinearly in (38), also controls the dominant color singlet contribution in inclusive hadroproduction of single quarkonia [78].In the latter channel the p T -dependence is predominanlty sensitive to the combination x 1 + x 2 − x 3 − x 4 (difference of impact parameters of the produced quarkonia in the amplitude and its conjugate), whereas dependence on other combinations of transverse coordinates is integrated out in the heavy quark mass limit [? ].The double quarkonia production cross-section (38) includes convolution with completely different dependence on coordinates, and for this reason provides independent observable for study of the quadrupole contributions.In the special limit when the distance between any quark-antiquark pair vanishes, the quadrupole amplitude reduces to a dipole amplitude, and for this reason (38) cancels exactly.Due to this, the integral (34) remains finite despite of possible singular behavior of the wave function ψ (γ) QQ QQ in this limit.While the differential cross-section (34) contains exhaustive information about the multipole elements of the target, due to smallness of the cross-sections, experimentally it may be interesting to study the p ⊥ -integrated cross-section, for which analytic integration over d 2 p ⊥ 1 and d 2 p ⊥ 2 yields dσ dy 1 dy 2 =4 (2π) Finally, we would like to discuss briefly evaluation of the diagrams shown in the Figure 3, which can contribute in general case of quarkonia pair production.The evaluation largely follows the same procedure, however the interaction of the shock wave with heavy quarks at early stages (diagrams in the left and middle columns of the Figure 3) does not allow to rewrite the cross-section as a convolution of ψ (γ) QQ QQ with universal target-dependent structure even in eikonal approximation: due to peculiar structure of the shock wave interaction with partons (16)(17)(18), there is no contributions with instantaneous virtual quarks and gluons which were taken into account in evaluation of ψ (γ) QQ QQ .Since in the final state we register only the color singlet heavy quarkonia (≈ take traces over the color indices of heavy quark lines, independently in the amplitude and its conjugate), the interaction with the target eventually can be expressed in terms of the dipole and quadrupole amplitudes (no sextupole and octupole contributions), however the structure of such expressions will be different for diagrams in each column of the Figure 3.A systematic evaluations for that case requires a dedicated study and will be presented elsewhere.

III. PHENOMENOLOGICAL ESTIMATES A. Parametrization of the dipole and quadrupole amplitudes
The numerical estimates depend crucially on the choice of the parametrization of the dipole and quadrupole amplitudes.For the dipole amplitude N , we will use the phenomenological impact parameter b-dependent "bCGC" parametrization [80,81] which has correct asymptotic behavior in the small-and large-size dipole limits, and has been fitted to reproduce various phenomenological data with reasonable precision.
The quadrupole amplitude in general is an independent nonperturbative object, not related to a dipole amplitude.However, as was discussed in [82][83][84], in the limit of large number of colors N c ≫ 1, it is possible to relate various multipole and multitrace matrix elements, and essentially express them all in terms of the dipole amplitude N .For heavy charmonia in QCD, the use of the large-N c limit is well-justified, since the discarded ∼ O (1/N c ) corrections parametrically are on par with strong coupling α s (m c ) ∼ 1/3.The derivation of these identities is cumbersome and may be found in the above-mentioned references [82,83].
The possibility to express the amplitude in terms of the color singlet dipoles may be also understood from the area enhancement argument introduced in [85][86][87]: The interaction of the color singlet multipole (e.g.quadrupole) with the target grows as a function of its size due to color transparency, for this reason we expect that in the phase space integration, the dominant contribution should come from the configurations with well-separated partons in the transverse space.On the other hand, color group averaging suppresses correlators in which the distance between partons exceeding the color correlation length set by saturation scale, namely s .However, this suppression does not work for configurations in which all partons are grouped into small-size color singlet pairs, separated by large distances between pairs [? ].After the integration over phase space, it becomes clear that such multi-pair configurations are enhanced compared to true or "genuine" multipoles (quadrupoles).The enhancement factor is given by ∼ S ⊥ Q 2 s n−1 , where S ⊥ is the transverse area of the target, and n ≳ 1 is the number of color singlet pairs.According to modern phenomenological parametrizations, the typical values of saturation scale Q s in the kinematics of interest are of order 1 GeV∼ (0.2 fm) −1 , for this reason the enhancement factor S ⊥ Q 2 s might be numerically large even for the proton.For this reason, any multipoint correlator may be approximated as a sum of pairwise products [85-87] where χ = {1, 2, ..., 2n} and Π(χ) is the set of partitions of χ with disjoint pairs.The result (40) significantly simplifies evaluation of the diagrams with multipole contributions and allows to apply Wick's theorem in order to introduce glasma graph approach.
As was suggested in [82,83], the result (40) for the the quadrupoles might be replaced with more accurate expression which takes into account the effects of JIMWLK evolution, as well as has more accurate behavior in a few physically important limits.Explicitly, the amplitude S 4 is given by where the function Γ Y (x i , x j ) is defined as The expression (41) includes denominators vanishing at the hypersurface however the function S 4 remains finite when approaching this hypersurface, with limiting value given by The combination S 4 (x 1 , x 2 , x 3 , x 4 ) − S 2 (x 1 , x 2 ) S 2 (x 3 , x 4 ) which appears in (38) may be rewritten as In the dilute limit, the expression for S 4 simplifies as [82,83] Substituting this result into the amplitude (37), it is possible to show that the nonzero contributions show up only when we take into account the nonlinear ∼ O N 2 terms in expansion and reads as The expressions in the second and the third lines of ( 47) may be rewritten using identity which clearly shows that (47) vanishes when any of the dipole sizes r 12 , r 34 , r ′ 12 , r ′ 34 goes to zero.Indeed, as was discussed in previous section, in this limit the quadrupole amplitude reduces to a dipole amplitude and exactly cancels in (38).In the small-r domain it is expected that the dipole amplitude N should have an asymptotic behavior where Q s is the saturation scale, and γ ≲ 1 is some numerical coefficient.For γ ≈ 1 (e.g."GBW" and "bSat" parametrizaitons [88]) it is possible to simplify drastically (47) and rewrite it as which demonstrates the dependence of the scattering amplitude N on sizes of individual dipoles in heavy quark mass limit.For γ ̸ = 1 a simple form ( 50) is not valid.However the antisymmetry of ( 48) with respect to permutations r 1 ↔ r 2 , r 3 ↔ r 4 indicates that the amplitude (47) should have a pronounced sensitivity to a relative orientation of vectors (r 12 , r 34 ) and (r ′ 12 , r ′ 34 ).

B. Numerical estimates for cross-sections
For the sake of definiteness, we'll focus on the production of the J/ψ + η b and Υ(1S) + η c production, which are expected to have the largest cross-section among all the charmonia-bottomonia pairs.We disregard the color octet channels in view of the smallness of color octet LDMEs (see however discussion in [10][11][12]16] about their significant uncertainty).For the dominant color singlet LDMEs we will use the values [10,74,75] where the LDMEs of pseudoscalar mesons are fixed using the relations valid in heavy quark mass limit [3].As expected from potential models [89][90][91][92], these LDMEs approximately scale with heavy quark mass as ∼ (α s (m Q ) m Q ) 3 .In Figure 4 we illustrate the transverse momentum dependence of the cross-section for different quarkonia states.The significant difference of J/ψ η b and Υ(1S) η c cross-sections might be understood in the dilute limit: in this case the pseudoscalar quarkonium can originate only from the secondary (lower) quark loop in the Figure 2. The virtuality of the gluon which connects the two loops is largely controlled by the mass of this pseudoscalar quarkonium, thus yielding a relative suppression factor In Figure 5 we study the dependence of the cross-sections on the azimuthal angle ϕ between the transverse momenta of the two quarkonia states.In order to make meaningful comparison of the cross-sections, which differ by orders of magnitude, we plotted the normalized ratio of the cross-sections, where ϕ max is the angle which maximizes the numerators of (54).Since the quarkonia (dipole) sizes are small in the heavy quark mass limit, the dependence on transverse momenta in the small-p T kinematics is largely sensitive to the dependence on the (Fourier conjugate) impact parameters b 12 = (x 1 + x 2 )/2 and b 34 = (x 3 + x 4 )/2 in the implemented quadrupole amplitude.For the parametrization (41), the dependence on angle φ between the centers of mass of the dipoles, b 12 and b 34 , is largely controlled by prefactor The left column corresponds to cross-sections of the photon-proton subprocess, the right column includes predictions for the cross-section of the full process.For ultraperipheral collisions we assume proton-lead collisions with A = 82.For EIC we considered the electron-proton collisions only; for heavy nuclei these results should be understood as cross-sections per nucleon.
For the sake of definiteness we considered production at central rapidities (y1 = y2 = 0) in the lab frame.The shape of the pT -dependence has a very mild dependence on rapidity and azimuthal angle ϕ between the produced quarkonia.
where we introduced variables r ij = x i −x j for the distance between pairs of the heavy quarks.In the heavy quark mass limit, the dominant contribution comes from the region r 12 , r 34 ≪ b 12 , b 34 .For the angle φ ≈ 0 (collinearly directed impact parameters), the factor ( 55) is suppressed at b 12 ≈ b 34 : physically this corresponds to a tiny quadrupole passing at some distance from the nucleus.The cross-section increases homogeneously as a function of ϕ, which merely reflects growth of the quadrupole moment of the four-quark ensemble.The predicted pronounced ϕ-dependence suggests that predominantly the production occurs in the back-to-back kinematics, akin to exclusive photoproduction of the same pairs [37].In order to demonstrate that the expected ϕ-dependence is due to the implemented parametrization of the quadrupole amplitude, in the second row of the Figure 5 we compare predictions of the quadrupole parametrization (41) with predictions obtained with a trivial parametrization which corresponds to a free noninteracting quadrupole.For the latter parametrization, integration of the quadrupole terms in (34) yields δ-functions ∼ δ p ⊥ 1 , δ p ⊥ 2 , so for nonzero momenta p ⊥ 1 , p ⊥ 2 the quadrupole contribution effectively drops out.As we can see from the last row in the Figure 5, the ϕ-dependence in this case is negligibly small: it exists due to a mild dependence on orientation of the two dipoles r 12 , r 34 in the wave function ψ (γ) QQ QQ .This finding corroborates that the predicted ϕ-dependence is due to the chosen parametrization of the quadrupole Figure 5. Upper row: Dependence of the normalized ratio R(ϕ), defined in (54), on the angle ϕ (azimuthal angle between transverse momenta of quarkonia).The left plot corresponds to the kinematics of ultraperipheral collisions at LHC, the right plot is for the kinematics of the Electron Ion Collider.The cross-section has almost the same angular dependence for all rapidities and transverse momenta pT (see the text for more detailed explanation).Lower row: Comparison of angular dependence found with parametrization (41) (curve with label "Iancu et al.") and a trivial parametrization (56).
amplitude.Finally, in the Figures 6, 7 we show the dependence on quarkonia rapidities in LHC and EIC kinematics.The growth of the cross-section as a function of average rapidity Y = (y 1 +y 2 )/2 merely repeats the rapidity dependence (growth) of the dipole scattering amplitude, which contributes directly and indirectly via parametrization of quadrupole amplitude.In the second row of the Figures 6, 7 we show the dependence of the cross-section on the rapidity difference ∆y between the two heavy mesons.For simplicity we made plots for the quarkonia having opposite rapidities in the lab frame, y 1 = −y 2 = ∆y/2.The configurations with large rapidity difference correspond to highly asymmetric sharing of the photon momentum between the quarks, for which the wave function ψ (γ) QQ QQ is strongly suppressed.For this reason, the quarkonia are predominantly produced with close rapidities.Since the variable ∆y may be related to the invariant mass of the heavy quarkonia pairs (9), this behavior implies suppression of the cross-section for large invariant masses of charmonia-bottomonia pairs.

IV. CONCLUSIONS
In this manuscript we analyzed in detail the inclusive photoproduction of heavy charmonia-bottomonia pairs.We focused on quarkonia pairs with opposite C-parities (Υ η c , J/ψ η b ), which do not require exchange of quantum numbers in the t-channel, and thus have the largest cross-sections.We found that the cross-section is sensitive to bilinear superposition of dipole and quadrupole contributions, which contribute in the same combination S 4 (x 1 , x 2 , x 3 , x 4 )− S 2 (x 1 , x 2 ) ⊗ S 2 (x 3 , x 4 ) as in inclusive single quarkonia production [78].Due to possibility to vary indepedently the  .The rapidity dependence of the cross-section in the kinematics of the ultraperipheral collisions at LHC.The plots in the upper row correspond to a configuration with equal rapidities of the produced quarkonia, y1 = y2, whereas the lower row corresponds to rapidities which differ by a sign in lab frame, y1 = −y2.In both rows the left plot corresponds to the cross-section of the photoproduction subprocess, and the right column shows predictions for the cross-section of the full process pA → A + M1M2X, assuming proton-lead collisions, as defined in (2).momenta of both quarkonia, the suggested process can allow to get better understanding of the quadrupole scattering amplitude S 4 (x 1 , x 2 , x 3 , x 4 ).We analyzed the role of the quadrupole scattering using parametrizations available from the literature [82,83] and found that numerically it gives the dominant contribution: its omission decreases the cross-section by a factor two-ten depending on the kinematics.The quadrupole term gives a pronounced dependence on the angle between transverse momenta of the quarkonia, which is almost negligible in its absence.
Numerically, the cross-sections of the suggested processes are small, but within the reach of the ultraperipheral collision experiments at LHC and at the future Electron-Ion Collider.The smallness of the cross-section happens due to heavy masses of charm and bottom quarks, which control the sizes of the produced cc and bb pairs.The suggested mechanism also contributes to charmonia-charmonia production, where its contribution is up to two orders of magnitude larger than for charmonia-bottomonia.However, as explained in Section II B, in that channel there are additional contributions which have a significantly more complicated structure and require a separate study.   .The rapidity dependence of the cross-section in the kinematics of the future Electron-Ion Collider.The left column corresponds to the cross-section of the photoproduction subprocess, and the right column shows predictions for the cross-section of the full electroproduction process eA → eAM1M2X.Due to heavy mass of the final state, the CGC approach is not applicable for y ≲ −1, and for y ≳ 1 the flux of the equivalent photons is suppressed due to leptonic factor (energy of the emitted virtual photon approaches from below the energy of the electron).The evaluation of the photon wave function follows the standard light-cone rules formulated in [19,93].The result for the QQ component is well-known in the literature [94,95] and is given by where λ, h, h are the helicities of the incoming photon and outgoing quark, ε µ (q) is the polarization vector of the photon, z is the fraction of the photon momentum carried by the quark, and r 12 is the transverse distance between quark and antiquark.The wave function of the QQ QQ-component can be expressed in terms of the wave function of QQ-component.In what follows we will focus on the photoproduction kinematics, assuming the transverse polarization of the incoming photons.We will use the reference frame in which the photon momentum has only plus-component, q ≈ q + , 0, 0 ⊥ .(A2) so the polarization vector of the incoming photon is given by In leading order over α s the wave function obtains contributions from the two diagrams shown in the Figure (8).In what follows 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.The evaluation of the diagrams follows the standard rules of the light-cone perturbation theory [19,93].The detailed evaluation of the QQ QQ component of the photon wave function ψ (γ) QQ QQ may be found in [36], and here for the sake of completeness we provide only the final result.This wave function can be represented as a sum where the first and the second terms correspond to contributions of non-instantaneous and instantaneous parts of propagators of all virtual particles, and for the sake of brevity we omitted color and helicity indices of heavy quarks (c i and a i respectively).The non-instantaneous contribution is given by the sum where where c a , c ′ a are the color indices in fundamental representation of the color group, and angular brackets stand for the color averaging Dµ [U x ].As explained in the text, this amplitude appears when we consider production of two color singlet quarkonia; for color octet channel the matrices U † , U are separated by additional color generators t a .Previously, similar correlator has been evaluated exactly in [86,87] in the context of studies of inclusive 3-quark production; for the 4-quark production the result was provided in implicit form (convoluted with hard amplitudes) in the large-N c limit, when all multipoles are expressed via dipole amplitudes.For the sake of completeness, below we provide a complete result for the amplitude.
Since the target is color singlet, it is expected that color averaging of the right-hand side should eventually include a sum over all possible color singlet irreducible representations which appear in direct product of U † U matrices, namely where summation is done over all possible permutations p of numbers (1,2,3,4), and we use a shorthand notation p k for the k-th element of permutation p. Overall, expression (B2) includes 4! = 24 unknown functions a p (...).In order to fix them, we may combine (B1,B2) and multiply both parts by where ℓ = (ℓ 1 , ℓ 2 , ℓ 3 , ℓ 4 ) is some arbitrary fixed permutaton.This yields a system of linear inhomogeneous equations where All the elements of the matrix M ℓ, p have very simple structure and are given just by N k c , where 1 ≤ k ≤ 4 is the number of disjoint cycles in permutation ℓ • p.For example, if ℓ 0 = (1, 2, 3, 4) and p 0 = (2, 1, 4, 3), then the corresponding element A contraction of δ-symbols in the right-hand side of (B5) allows to express all elements B ℓ in terms of color singlet multitrace operators.For example, for permutations p 0 = (2, 1, 4, 3), p 1 = (2, 3, 1, 4) and p 2 = (2, 3, 4, 1) we may get respectively B p 0 = tr U † (x 1 ) U (y 1 ) U † (x 2 ) U (y 2 ) tr U † (x 3 ) U (y 3 ) U † (x 4 ) U (y 4 ) (B7) B p 1 = tr U † (x 1 ) U (y 1 ) U † (x 2 ) U (y 2 ) U † (x 3 ) U (y 3 ) tr U † (x 4 ) U (y 4 ) , (B8) B p 2 = tr U † (x 1 ) U (y 1 ) U † (x 2 ) U (y 2 ) U † (x 3 ) U (y 3 ) U † (x 4 ) U (y 4 ) , (B9) which includes in the right-hand side quadrupole, sextupole and octupole contributions.The solution of nonlinear system (B3) is straightforward and was done with Mathematica (and FeynCalc package [96,97] for color group algebra).The result for the constants a p is very lengthy if written in terms of conventional notations S 2n (x 1 ), for this reason, for the sake of brevity we will introduce shorthand notations [i] = S 2 (x i , y i ) , [i, j] = S 4 x i , y i , x j , y j , [i, j, k] = S 6 x i , y i , x j , y j , x k , y k , (B10) [i, j, k, ℓ] = S 8 x i , y i , x j , y j , x k , y k , x ℓ , y ℓ , (B11) where each of the indices i, j, k, ℓ may take integer values between 1 and 4. In these notations the amplitude A is written as We need to mention that each term in the sum (B12) after contraction with partonic amplitudes can contribute with different factor N n c , n ∈ N, for this reason it might be incorrect to take the limit N c → ∞ directly in (B12).As can be seen from Figure 2 and (37), for diagrams 5,6 we will need only the amplitude convoluted with color group generators Using Fierz identity and performing convolutions with the help of FeynCalc, we may obtain a surprisingly simple result, )⟩ = (B15) = ⟨(S 4 (x 1 , y 1 , x 2 , y 2 ) − S 2 (x 1 , y 1 ) S 2 (x 2 , y 2 )) (S 4 (x 3 , y 3 , x 4 , y 4 ) − S 2 (x 3 , y 3 ) S 2 (x 4 , y 4 ))⟩ , namely all the contributions of sextupoles and octupoles cancel in the final result.This is an untrivial property of the photoproduction by color singlet photon: such cancellation does not occur if we consider hadroproduction or production of color octet QQ pairs.Note that the result of this appendix are exact, namely we haven't used the large-N c limit so far.

Figure 2 .
Figure 2. Diagrams which describe the inclusive charmonia-bottomonia pairs production in CGC picture in the leading order over αs.The diagrams in the right column differ from diagrams in the left column just by inversion of the quark line in the upper loop (charge conjugation).The dominant color singlet contribution comes from the diagrams in the last row (diagrams "5" and "6" respectively).The subscript numbers 1-4 enumerate the heavy quark lines in our convention.The subscript letters z1 ... z3 stand for the coordinates of the interaction vertices in the configuration space.The red block represents the interaction with the target (shockwave).

Figure 3
Figure 3.Additional diagrams which in general should be taken into account for inclusive quarkonia pair production in CGC picture in the leading order over αs, yet do not contribute to charmonia-bottomonia pairs.The subscript numbers 1-4 enumerate the heavy quark lines in our convention.The subscript letters z1...z3 stand for the coordinates of the interaction vertices in the configuration space.The red block represents the interaction with the target (shockwave).

Figure 4 .
Figure 4.The cross-sections of inclusive production of different quarkonia pairs with opposite C-parities.For the sake of definiteness we considered that both quarkonia are produced with the same absolute value of the transverse momenta p ⊥ 1 = p ⊥ 2 = pT .Upper and lower rows correspond to kinematics of ultraperipheral collisions and the Electron-Ion Collider, respectively.The left column corresponds to cross-sections of the photon-proton subprocess, the right column includes predictions for the cross-section of the full process.For ultraperipheral collisions we assume proton-lead collisions with A = 82.For EIC we considered the electron-proton collisions only; for heavy nuclei these results should be understood as cross-sections per nucleon.For the sake of definiteness we considered production at central rapidities (y1 = y2 = 0) in the lab frame.The shape of the pT -dependence has a very mild dependence on rapidity and azimuthal angle ϕ between the produced quarkonia.

Figure 6
Figure 6.The rapidity dependence of the cross-section in the kinematics of the ultraperipheral collisions at LHC.The plots in the upper row correspond to a configuration with equal rapidities of the produced quarkonia, y1 = y2, whereas the lower row corresponds to rapidities which differ by a sign in lab frame, y1 = −y2.In both rows the left plot corresponds to the cross-section of the photoproduction subprocess, and the right column shows predictions for the cross-section of the full process pA → A + M1M2X, assuming proton-lead collisions, as defined in(2).

Figure 7
Figure 7.The rapidity dependence of the cross-section in the kinematics of the future Electron-Ion Collider.The left column corresponds to the cross-section of the photoproduction subprocess, and the right column shows predictions for the cross-section of the full electroproduction process eA → eAM1M2X.Due to heavy mass of the final state, the CGC approach is not applicable for y ≲ −1, and for y ≳ 1 the flux of the equivalent photons is suppressed due to leptonic factor (energy of the emitted virtual photon approaches from below the energy of the electron).

3 qℓFigure 8 .
Figure 8.The leading order contribution to the wave function ψ (γ) QQ QQ defined in the text.The momenta ki shown in the right-hand side are Fourier conjugates of the coordinates xi.It is implied that both diagrams should be supplemented by all possible permutations of final state quarks (see the text for more details).