BFKL Pomeron loop contribution in diffractive photoproduction and inclusive hadroproduction of $J/\psi$ and $\Upsilon$

We analyze contributions to the heavy vector meson production with large transverse momentum in proton--proton and diffractive photon--proton scattering driven by an exchange of two Balitsky--Fadin--Kuraev--Lipatov Pomerons in the squared amplitudes. The Pomerons couple to a single parton and form a Pomeron loop closed by the vector meson impact factors. For the photon--proton case the diffractive cut of the Pomeron loop contributes, and for the inclusive hadroproduction one finds the loop with two cut Pomerons. We compute both of these Pomeron loop contributions and study in detail their properties. The results are then used to calculate the cross sections for diffractive $J/\psi$ photoproduction with large transverse momentum at HERA and the correlated two Pomeron contribution for inclusive $J/\psi$ and $\Upsilon$ production cross sections at the LHC. Within a unified approach a good description of the photoproduction data is found, but correlated two Pomeron mechanism gives only a small contribution to hadroproduction of the vector mesons at the LHC.


Introduction
The heavy vector mesons with negative C-parity -charmonia and bottomonia -are classical probes of the QCD exchange at high energies. The signals of J/ψ, ψ and Υ mesons in their leptonic decay channels are clear and allow for accurate measurements of the corresponding differential cross sections. The underlying production dynamics is driven by gluonic degrees of freedom and their QCD evolution. The currently accepted pictures of the heavy vector meson production mechanisms in diffractive photon-hadron and inclusive hadron-hadron collisions are, however, quite different. The diffractive photoproduction data at high energies and large transverse momentum obtained by H1 [1] and ZEUS [2] collaborations at HERA have been successfully described [3][4][5][6][7] assuming an exchange of the non-forward Balitsky-Fadin-Kuraev-Lipatov (BFKL) Pomeron [8][9][10][11][12][13][14][15][16][17] in the diffractive amplitude. On the other hand, in the high energy inclusive hadroproduction of J/ψ and Υ, a good description of data from hadron colliders requires adopting the Color Octet Model (COM) [18][19][20][21][22][23][24][25], see also Refs. [26,27] for a review. In the present paper we shall investigate the diffractive photoproduction and a contribution to the inclusive hadroproduction of heavy vector mesons assuming the same underlying QCD dynamics of two BFKL Pomeron exchange.
The standard picture of diffractive photoproduction of vector mesons at HERA at large momentum transfer assumes exchange of gluonic hard color singlet across large rapidity distance between the incoming photon -outgoing meson vertex and the diffractive remnant of the proton. The kinematics of this process allows to apply the high energy limit in which the dominant contribution to the color singlet amplitude is given by the BFKL Pomeron [9][10][11][12][13][14]. By the BFKL Pomeron one understands the system of two Reggeized gluons in the t-channel interacting by exchange of usual gluons. The propagation of the Reggeized gluons and the effective interactions between them are derived in QCD in the high energy limit. In more detail, the exchange amplitude is described by the BFKL evolution equation that formally resums logarithmically enhanced perturbative corrections to all orders. In the BFKL approach one resums logarithms of a ratio of a large collision energy √ s and other, much smaller mass scales e.g. the meson mass or the momentum transfer. These logarithms are related to the rapidity distance Y between the projectile and the target in the high energy scattering process. So far the BFKL resummation in QCD was performed at the leading logarithmic (LL) [9][10][11][12][13][14] and next-to-leading logarithmic (NLL) approximation [15,16,[28][29][30][31][32]. In the LL approximation one resums terms δ (n) LL ∼ α n s Y n while δ (n) NLL ∼ α n+1 s Y n terms are resummed by the NLL BFKL evolution. The BFKL formalism assumes high energy (or k T ) factorization in which hard matrix elements are factorized in rapidity space from the BFKL evolution. In addition, matrix elements are off-shell, with initial quarks or gluons carrying non-zero transverse momenta unlike in the standard collinear approximation. Hence also the BFKL Pomeron may carry non-zero transverse momentum, and the corresponding amplitude is governed by the non-forward BFKL Pomeron. This formalism was applied [5,6,[33][34][35] some time ago to the data from HERA on J/ψ, ρ and φ diffractive mesons photoproduction [1,2,36] with large transverse momentum p T and was shown to describe the data well. In this paper we revisit the diffractive J/ψ photoproduction at HERA and use the established description of this process to estimate the BFKL Pomeron loop contribution to inclusive vector meson hadroproduction.
The COM of inclusive heavy meson hadroproduction assumes non-zero amplitudes for a change of the quantum numbers (in particular color and angular momentum) between the partonic phase and the meson [18][19][20]. More specifically, in the partonic subprocess heavy quark-antiquark pair QQ is produced with an arbitrary color and angular momentum quantum numbers, and the transition to the final state meson is governed by separate multiplicative coefficients for each set of the partonic quantum numbers. The values of these coefficients are obtained by fitting the predicted cross sections to experimental data. The theoretical basis for this mechanism comes from two complementary sources. First, within heavy quark effective theory one finds non-zero amplitudes of higher Fock states in the heavy vector meson wave function [18], for instance a state QQg with an additional gluon g. Obviously, the quantum numbers of QQ in this state do not match the quantum numbers of the meson. From another perspective, before the partonic QQ makes the meson a process of hadronization occurs, in which the quantum numbers may change. This is a picture of fragmentation of the primordial QQ state into the heavy meson [22]. The hadronization process does not have satisfactory perturbative description as it occurs at low (hadronic) scales and depends on non-perturbative properties of the QCD vacuum. In both scenarios the transition amplitudes of QQ to the meson cannot be derived from theory -only order-of-magnitude estimates can be obtained. Nevertheless the COM has a solid theoretical basis, and it is strongly supported by the successful fits of its predictions to the bulk of experimental data. However, since the parameters of the model are fitted, there is still room for other than COM possible mechanisms of the inclusive vector meson hadroproduction.
The classical alternative to the COM of the heavy meson production is the Color Singlet Mechanism (CSM) [37][38][39]. In fact, the CSM was considered to be the standard QCD prediction before it was contradicted by the Tevatron data [40,41]. In this approach one assumes the exact matching of the quantum numbers of the produced partonic QQ state and the final state meson. The main advantage of this mechanism is its completeness within perturbative QCD and no need for additional parameteres. For the C-odd mesons V at the leading twist the CSM is driven by the g + g → V + g partonic scattering. The predictions of the CSM, however, fail badly in describing the J/ψ hadroproduction at the Tevatron and the LHC, see e.g. [26,27]. For the total inclusive vector meson hadroproduction cross sections the CSM both at the leading order (LO) and at the next-to-leading order (NLO) are more than one order of magnitude below the data [26,27]. Moreover, the CSM predictions lead to the distribution in the meson transverse momentum p T which is much too soft, while the COM is able to describe well the p T dependence of the meson production cross section. The CSM and COM approaches were also extended from the collinear approximation to the k T -factorization framework [42][43][44][45][46][47][48]. It was found that also in the k T -factorization approach the CSM model is much below the data for the direct J/ψ production at large transverse momenta.
Beyond the leading twist approximation the CSM may be realized also with a fusion of three initial state gluons. At the partonic level the meson formation occurs by g+g+g → V diagrams with the coupling through the heavy quark loop. The Three Gluon Fusion (3GF) mechanism was considered in Ref. [49] as a contribution to J/ψ hadroproduction, and in Ref. [50] it was proposed as a possible leading contribution to heavy vector meson hadroproduction. Then it was further studied in Refs. [51][52][53]. Moreover, recently contributions of multiple gluon couplings in the J/ψ production were combined with the COM [54,55]. In this approach a good description of the J/ψ hadroproduction data was obtained, including the meson polarization. Although the uncorrelated 3GF mechanism is an implicit contribution of this framework, at larger p T the cross section is dominated by the COM contributions.
In the 3GF mechanism one of the gluons comes from one hadronic beam (the projectile), and two other ones from the other beam (the target). These two gluons in the t-channel can be taken either as completely independent (uncorrelated) or as coming from a single parton of the target (correlated). The two scenarios correspond to an uncorrelated double gluon distribution in the target, and to the correlated contribution in the double gluon distribution, respectively. A detailed study of the uncorrelated contribution to the 3GF mechanism showed that it may contribute to the total J/ψ hadroproduction cross section as a fraction of about 20 -25% [51]. The obtained p T -dependence of the meson production differential cross section dσ/dp T was found to be much steeper than the experimental data. Specifically, at larger values of p T the experimental data can be approximated with a power law: dσ/dp T ∼ 1/p n T with n 5, while with the uncorrelated 3GF mechanism n > 8 is obtained. This rather steep p T dependence is well understood as the uncorrelated 3GF contribution enters at a higher twist, and it is suppressed at large p T by an additional factor of (Λ h /p T ) 2 with respect to the leading twist, where Λ h is a small hadronic scale.
In this paper we analyze in detail the correlated 3GF mechanism. The two gluons that enter the vector meson production vertex from the target side are assumed to come from a splitting of a single parton: quark or gluon in the target. Since the parent parton is point-like, it does not introduce any additional hadronic scale and one expects a harder dependence of the meson p T -distribution than it was for the uncorrelated 3GF. At the lowest order this mechanism occurs through g + g → V + g (or g + q → V + q) partonic process with an exchange of two gluons between the g → V transition vertex at the side of the projectile and the g → g (or q → q) scattering at the target side. At the amplitude level the two gluon exchange in the t-channel carries the symmetric color octet. The important feature of the process at the lowest order is a flat dependence on the rapidity distance Y between the projectile gluon and the target gluon. Beyond the lowest order approximation the amplitude of this subprocess is modified by QCD radiative corrections. These corrections can be resummed by a QCD evolution equation. In the high energy limit corresponding to Y 1, the cross section of correlated 3GF cross section is driven by an exchange of four Reggeized gluons in the total color singlet state that interact with BFKL kernels. This system is described by Bartels-Kwieciński-Prasza lowicz (BKP) evolution equation [56][57][58]. It was shown in Ref. [59] that in the large N c limit the system of four Reggeized gluons in the color singlet channel may be approximated by two independent BFKL Pomerons. In consequence one expects a strong enhancement ∼ exp(2∆ P ) of the cross section at large Y by the double BFKL evolution with the BFKL Pomeron intercept ∆ P ∼ 0.3. In addition, the anomalous dimensions of the BFKL Green's functions could lead to a less steep p T dependence. So, despite this contribution enters at the O(α 5 s ) order and is a well defined part of the NNLO correction to the CSM contribution, it may be important due to strong effects of the BKP evolution. The first estimate of this contribution [50] suggested that it may reproduce the inclusive J/ψ hadroproduction cross section data from the Tevatron. In this paper we perform a detailed calculation of the correlated 3GF cross section to verify its importance.
In the analysis we shall use a connection between the correlated 3GF mechanism and the diffractive photoproduction of vector mesons at large p T . This connection originates from the kinematical identity of the impact factors corresponding to γ + (2g) → V and g + (2g) → V transitions. At the leading order the difference between these two impact factors comes only from the color factors and the gauge coupling constants. Hence, the cross sections for a vector meson production by gluon scattering off a partonic target and for the meson diffractive photoproduction are proportional to each other. This relation imposes important constraints on the correlated 3GF mechanism. When the BKP/BFKL evolution is taken into account however, the diffractive photoproduction and the 3GF mechanism are different due to the different color flow. The diffractive photoproduction cross section corresponds to the diffractive cut of two BFKL Pomeron exchange, and the correlated 3GF cross section to an exchange of two cut BFKL Pomerons. Thus independent evaluation is necessary for these two processes. In both cases, however, at the level of cross section one finds the topology of the BFKL Pomeron loop spanned between the meson impact factors at the projectile side and the parton impact factor at the target side. So, besides evaluating the correlated 3GF contribution with the BKP/BFKL evolution effects included, we shall also revisit the diffractive photoproduction case and analyse in detail the properties of the BFKL Pomeron loop in the t-channel.
The paper is organized as follows. In the next section we describe the theoretical framework for the diffractive heavy vector meson production in DIS, and the non-forward BFKL evolution. In Section 3 the two-Pomeron contribution to the vector meson hadroproduction is analyzed. The details of the color factors are discussed and the BFKL Pomeron in the conformal representation is described. In Section 4 the properties of the Pomeron loop at the parton level are analyzed, in particular the dependence on the transverse momentum. In Section 5 we present the comparison of our numerical calculations with the vector meson production in diffractive DIS as well as in hadroproduction, and we discuss the results in Section 6. Finally, in Section 7 we give summary and conclusions.

Diffractive heavy vector meson production
We begin with a short recollection of the perturbative QCD approach to the hard color singlet exchange in the diffractive heavy vector meson V photoproduction off a proton at large momentum transfer t. The process was investigated in detail at HERA in e ± p collisions with invariant c.m.s. energy √ S = 318 GeV. In the measurement of the process e ± p → e ± V X, a large rapidity gap devoid of particles is required. It separates the produced vector meson V and the dissociated proton remnant X. The e ± p cross section may be factorized into a universal flux function of quasi-real photons in the electron and the cross section for the diffractive photoproduction subprocess, with the rapidity gap. In this process the photon-proton invariant mass, √ s = √ zS, is assumed to be much larger than all the other scales present in the process, hence s |t| and s M 2 V , where M V is the meson mass. The applicability of the perturbative QCD is ensured by the conditions |t| Λ 2 QCD and M 2 V Λ 2 QCD . The diffractively produced C-odd state V is assumed to be a heavy vector quarkonium. In this work we focus on the J/ψ meson. Since in the available data the photon flux is strongly dominated by very low virtualities −q 2 = Q 2 we take the limit Q 2 → 0 in calculations of the QCD amplitudes.
Within the kinematic regime specified above, the color singlet exchange in the t-channel of the process (2.1) may be described in QCD as the perturbative Pomeron exchange (Fig. 1), that is governed by the non-forward BFKL equation [9-11, 13, 14]. Due to the large momentum transfer, the cross section can be factorized into a partonic cross section (dominated by the Pomeron exchange) and the collinear parton distribution functions (PDFs), that describe the structure of the proton target: whereŝ = xs is the photon-parton invariant mass squared, and f i is the PDF for the parton i which may be a quark (anti-quark) or a gluon. A natural choice for the factorization scale is µ ∼ |t|. Here x is the longitudinal fraction of the proton light cone momentum carried by the quark or gluon. In what follows, the transverse two-vectors in the light cone basis are denoted by the bold characters; for example, we denote the vector meson transverse momentum as p. In the high energy kinematics we have t −|p| 2 ≡ −p 2 . In Eq. (2.2) the coupling of the high-t Pomeron to the proton is assumed to occur through coupling to the individual partons. This approximation was studied in detail and motivated in [35,60,61].
In the high energy limit, the kinematic part of the BFKL Pomeron coupling to quarks and gluons is the same. The only difference between the quark and gluon partonic target comes from the color factors, where C γi , i = q, g, is the color factor and with A being the amplitude to produce the vector meson through a single Pomeron exchange (Fig. 2). The amplitude is dominated by the imaginary part, and the correction coming from the real part enters only at a subleading order in the logarithmic expansion and may be neglected. The imaginary part of the amplitude reads where k 1 and p − k 1 ≡ k 2 are transverse momenta of the exchanged gluons, and p is the transverse momentum carried by the Pomeron. Φ V , Φ q are the impact factors for the vector meson and for the quark, respectively. They are stripped off the color factors. In addition, the quark impact factor Φ q is evolved to rapidity y using the BFKL evolution. The exact form of the impact factors will be described below in this section. The rapidity evolution length y of the quark impact factor is defined by the relation where we set Λ = E T ≡ M 2 V + p 2 . This choice is different than the choice made in earlier studies [4][5][6][7], where Λ = M V was used. The latter value was selected as a result of fits driven mostly by the light vector meson high-t photoproduction data. The J/ψ data, however, were well described both with Λ = E T and Λ = M V . Hence in this paper we use Λ = E T , the value with a straightforward kinematic motivation. This choice was used e.g. in Ref. [3].
The small parameter s 0 is a phenomenological infrared cutoff that mimicks the effects of the color confinement; the results are, however, finite for s 0 → 0. We introduce this parameter following the approach proposed in Ref. [62]. The amplitude does not depend on the produced meson direction in the transverse plane, hence the phase space d 2 p in (2.4) may be trivially integrated over the azimuthal angle. We present the results in this form to keep the notation uniform with the more complicated case of the vector meson hadroproduction, which we shall describe in Section 3.
The diagrams contributing to the cross section for diffractive vector meson V photoproduction off a parton: a gluon (left) and a quark (right). The non-forward BFKL Green's function G is evolved along the rapidity gap between the vector meson and the proton remnants.
In the calculations we take the non-relativistic approximation for the meson wave function. With this assumption the lowest order photon to vector meson impact factor reads [33,63,64] where a, b are color indices of the exchanged gluons, and the kinematic part of the impact factor reads with e q being the charge of the quark in the meson in units of the elementary charge e, M V -the mass of the vector meson, and Γ V →ll its leptonic decay width. The photon to vector meson impact factor (2.8) is valid for the transverse polarizations of the photon and of the vector meson. The incoming photon is quasi-real hence its polarization is constrained to be transverse. The amplitude of the transverse photon to longitudinaly polarized vector meson was estimated in Ref. [33] to be small. Therefore the amplitude described by (2.8) provides the dominant contribution to diffractive vector meson photoproduction. The quark impact factor in the color singlet channel Φ 1,ab q (y, k 1 , p) is also factorized into the color part and the kinematic part, where the notation is the same as in the case of the photon-vector meson impact factor. The diffractive gluon impact factor Φ 1,ab g differs from the diffractive quark impact factor only by the color factor, The kinematic part of the quark impact factor Φ q (y, k 1 , p) in (2.5) is the solution of the non-forward BFKL equation with the initial condition given by the leading order quark impact factor [62] Φ q,0 (k 1 , p) = α s .
At a given momentum transfer p through the Pomeron, the evolved quark impact factor Φ q (y, k 1 , p) may be represented as convolution of the leading order impact factor Φ q,0 with the non-forward BFKL Green's function G y : In what follows, we shall use transverse momentum variables k, k that reflect the symmetry of the problem: k = (k 2 − k 1 )/2, k = (k 2 − k 1 )/2. The transverse momenta of gluons take the form 14) The explicit form of the leading logarithmic BFKL equation that defines the quark impact factor Φ q (y, k 1 , p) is the following The non-forward BFKL equation given in (2.15) is solved numerically, using an approximation introduced in [62]. The idea is to use the Fourier decomposition of the impact factor w.r.t. the angle φ k between k and p, where the Fourier coefficients for m > 0. We have checked that the full solution Φ q (y, k 1 , p) is well approximated by the leading component Φ q (y, k 2 , p 2 ), in accordance with results [5,62]. Since this approximation leads to much greater numerical efficiency, with a negligible effect on accuracy, we use it in the estimates of the cross section. The leading Fourier component with m = 0 obeys the equation, The independence of the leading order quark impact factor on the angles allows to set Finally, let us note that the cross section (2.2) can be written in terms of (2.4) as follows: (2.20) Using Eqs. (2.7), (2.10) and (2.11) it is straightforward to obtain the color factors for the diffractive production off the quark and gluon. The results read 3 Two Pomeron contribution to heavy vector meson hadroproduction

Direct approach
Let us now consider the heavy vector meson hadroproduction through the (cut) double Pomeron exchange. We focus on a hadronic analogue of the diffractive vector meson photoproduction -where the projectile is a gluon, and the two Pomerons couple to a single parton in the target. The general diagrams contributing to the partonic cross section, dσ 2−P , are depicted in Fig. 3. Note that at the lowest orders, the topologies of the considered partonic processes are the same as for the diffractive photoproduction. After inclusion of the evolution, however, the two processes correspond to different cuts through the two Pomerons; one finds the diffractive cut in the photoproduction, and the double cut Pomerons in the hadroproduction. The vector meson impact factor describes the fusion of three gluons into the meson. The coupling of both Pomerons to the single parton in the target leads to a correlation of the gluon distributions in the target, that enter the meson impact factor. This should be contrasted with another possible contribution, where the gluons in the target are uncorrelated, see [51].
In the calculations, the incoming partons from the projectile and target protons are treated as collinear and carry longitudinal momentum fractions x 1 and x 2 , whereas the gluon exchange in the t-channel is decribed in the k T -factorization framework, assuming the high energy limit. Hence the vector meson impact factor that enters the calculation describes the transition from the collinear projectile gluon to the vector meson by coupling of two t-channel gluons with non-zero transverse momenta. This impact factor may be obtained from the impact factor describing the fusion of three gluons into the vector meson that was derived in [64] by taking the collinear limit for the projectile gluon. Because of the odd C-parity of the meson, the impact factor is fully symmetric in color indices of the gluons, and the functional dependence on the external boson momenta is the same as in the diffractive impact factor for the exclusive vector meson photoproduction. Specifically, in the collinear limit for the projectile gluon, the three-gluon impact factor for inclusive vector meson hadroproduction Φ a 1 a 2 c where g s is the strong coupling constant (α s = g 2 s /4π), d a 1 a 2 c is the fully symmetric color tensor, a n with n = 1, 2, are the color indices of the t-channel gluons, and c is the color index of the projectile gluon.
The corresponding leading order impact factors of the target quark (i = q) and gluon (i = g) may be expressed in the following way: Figure 3. where matricesT bn R i are the generators of the color group in the color charge representation R i of the target particle. Note that after the projection on the color singlet channel performed by Tr( are recovered. Using the notation from the previous section we can express the cross section for the inclusive vector meson production in the considered two Pomeron mechanism as where S is the hadronic collision energy squared, the color factors C gq and C gg accommodate the color structure of the vector meson impact factor and the color projection onto the two-Pomeron state. Their calculation is straightforward, but requires certain assumptions on how the projection is made. We shall discuss this issue later in this section and give the values of the color factors. The diagrams in Fig. 3, stripped off the color factors, give the following expression for the partonic cross section for the vector meson production with the transverse momentum p: where we set the rapidity evolution length to y = log(x 1 x 2 S/(M 2 V +p 2 )), and the remaining notation is the same as in the previous section. In fact, the above equation describes the BFKL Pomeron loop with cut through both Pomerons, and the upper and lower Pomeron couplings given by the meson and quark four-gluon impact factors. The color part of the meson impact factor forbids the coupling of the gluon pair in the color singlet state to the meson, so the Pomerons are formed only between the gluons at the opposite sides of the unitarity cut. The transverse momentum of the Pomerons in the loop is ±q. In numerical estimates of (3.4) we shall use the approximation of the quark impact factors by the leading Fourier components, Φ q (y, cf. the discussion performed for the diffractive photoproduction. Note that we apply in Eq. (3.4) the vector meson impact factor Φ V corresponding to a transition of a transversely polarized gluon to a transversely polarized meson. The gluon is treated within the collinear approximation so it cannot carry the longitudinal polarization. As in the diffractive photoproduction case the transition of a transversely polarized gluon to a longitudinaly polarized meson is strongly suppressed [33].
Special attention should be paid to the target quark and gluon impact factors. At the leading order the parton four-gluon impact factor is a constant function of the gluon momenta, as it is for the two-gluon impact factor. This follows directly from the point-like nature of the partons. As a result, the four-gluon parton impact factor is proportional to the product of two-gluon impact factors. Next, in our approach we approximate the full Bartels-Kwieciński-Prasza lowicz (BKP) [56][57][58] evolution of the four gluon t-channel state in the amplitude squared by the independent evolution of two Pomerons. This approximation is valid in the large N c limit, as color reconnection between two Pomerons is suppressed by 1/N 2 c [59]. Hence, with the factorized form of the parton impact factor, and with the independent Pomeron evolutions, also the evolved four-gluon parton impact factor may be factorized (up to a constant factor) into a product of two-gluon evolved impact factors. Thus, in Eq. (3.4) both quark impact factors are evolved independently according to Eq. (2.19).
Let us now discuss the lowest order contribution to the correlated cross section (3.3), obtained by setting Φ q (y, k i , ±q) → Φ q,0 (k i , ±q) in Eq. (3.4). The corresponding diagrams are depicted in Fig. 4. We see that these diagrams are virtually the same as for the diffractive photoproduction at the lowest order, except for the vector meson vertex, which here contains an incoming gluon instead of a photon. Given the symmetry of the color part of the impact factor, the lowest order gq → V q and gg → V g cross sections may be obtained from the γq → V q and γg → V g cross sections by suitable modifications of the color factors and the coupling constants, while the momentum dependent part remains the same. The explicit result reads where This result is a useful benchmark for the cross section (3.4) with the BFKL evolution included.  Figure 4. The lowest order diagrams for the amplitude squared for the production of the heavy vector meson through the coupling of 1 + 2 gluons off a quark (left) and off a gluon (right).

The color factors
In terms of the BKP equation eigenstates, Eq. where a n , n = 1, . . . , 4, and c are color indices in the adjoint representation, a n describe the t-channel gluons in the natural order. The scalar product on the space of the color tensors may be defined as It is convenient to use the normalized color tensors The color tensorsP 0 ({a n }),P 2 ({a n }) andP d ({a n }) are normalized to one, and orthogonal up to 1/N 2 c corrections, P A |P B = δ AB + O(1/N 2 c ). Since the analysis is performed at the leading order in N c , the tensorsP A ({a i }) may be treated as an othonormal basis. Therefore the projectors may be defined on the color tensors, corresponding to the BKP eigenstates: (3.10) In our calculation we project the impact factors on the color state P 2 corresponding to the exchange of two cut Pomerons. The color tensors associated with the upper and lower impact factor are denoted by C α ({a n }) and C β ({b n }) correspondingly, and we define the color factor by where the color projector takes the following form The color tensors entering the cross sections are the following: for the color averaged gluon to meson transition amplitude squared, and for the color averaged scattering amplitude squared of the parton in color representation R i . The explicit form of color factors C αβ follows from Eqs. (3.11), (3.13), (3.14) and reads

Vector meson hadroproduction in conformal representation of the BFKL Pomeron
In what follows we shall recall the Lipatov solution to the non-forward leading logarithmic BFKL equation by means of the conformal eigenfunctions [14], and the application to diffractive scattering. Then we shall employ the formalism to describe the vector meson hadroproduction in the two Pomeron exchange mechanism. The solution to the nonforward BFKL evolution equation can be presented in the momentum or in the coordinate space. To this end, for a generic 2-dimensional transverse momentum k = (k x , k y ) we introduce the complexified momenta (k,k), where Similarly, for the transverse coordinate space, the 2-dimensional coordinates ρ = (ρ x , ρ y ) will be traded for the complex variables In what follows we set s 0 = 0 and use the original form of the non-forward BFKL equation. It was shown by Lev Lipatov [14,17] that the BFKL equation in the leading logarithmic approximation is invariant under the conformal transformations of the complexified transverse positions of the Reggeized gluons, for arbitrary complex parameters a, b, c, d. In analyses of high energy scattering amplitudes A(s, t) it is customary to use the Mellin moments ω conjugate to s [17], The solution for the gluon Green's function is then represented as , (3.20) where ω n (ν) is the LL BFKL eigenvalue and ψ are polygamma functions, the functions E n,ν are conformal eigenfunctions defined as where the powers h,h are the conformal weights The above form of the BFKL Green's function may be Fourier transformed to the transverse momentum representation, and inverse-Mellin-transformed in ω to the rapidity space are the conformal eigenfunctions in the momentum representation. The analytic form of the eigenfunctions E nν (k 1 , q) was derived in [65]. It is rather lengthy, so instead of listing it here we refer the reader to the original paper. In the limit y → 0, one finds (3.26) The dominant imaginary part of the amplitude for the BFKL Pomeron exchange between leading order impact factors Φ A,0 (k 1 , q) and Φ B,0 (k 1 , q) reads It is convenient to define the projection of the impact factors on the conformal eigenfunctions 28) and analogously for the index B. This leads to the following form of the Pomeron exchange amplitude The impact factors I q n,ν (q) and I V n,ν (q) were computed in the analytic form in Refs. [3,5,35,60]. In the calculation of vector meson hadroproduction we shall need the quark impact factor, which is treated within the Mueller-Tang scheme [35,60]. It takes the form where φ q is the polar angle of the vector q in the transverse plane. The formalism may be also applied to the two-Pomeron exchange process, assuming independent BFKL evolution of the Pomerons. Hence, we rewrite Eq. (3.4) as Using the representation of the BFKL Green's functions by the conformal eigenfunctions in the momentum space this may be rewritten as, A new non-trivial object that appears in this equation, I V ⊗V n 1 ,ν 1 ;n 2 ,ν 2 (p, q) is a projection of a pair of the vector meson impact factors (coming from the amplitude and its complex conjugate) on the conformal eigenfunctions: Note that in contrast to the diffractive scattering, in this expression the vector meson impact factors Φ V are coupled to two different conformal eigenfunctions each. In other words, the two meson impact factors and the two conformal eigenfunctions are all entangled and integral (3.33) cannot be factorized.
In order to cross check the numerical predictions for the vector meson hadroproduction by the two Pomeron exchange, obtained primarily by solving the Eq. (2.19) and Monte Carlo integration in Eq. (3.3), we calculated numerically integral (3.33). Then, the partonic cross section in the conformal representation was evaluated using Eq. (3.32). In this approach a good convergence of the emerging numerical sums and integrals was found for α s y 1, where the results were found to coincide with results obtained within the direct approach described in Section 3.1.

Properties of the Pomeron loop at the parton level
In this subsection we analyze the properties of the partonic cross section which includes Pomeron loop contribution (3.4). In all calculations we shall set the infrared cutoff s 0 to be equal to 0.5 GeV 2 .
In Fig. 5 the result of calculation of the partonic cross section from Eq. (3.4) is shown as a function of the transverse momentum p T of the produced vector meson for fixed values of rapidity. The normalization is given by Eq. (3.4), but the cross section is divided by α 5 s . The reason to divide out the strong coupling constant is that it introduces additional transverse momentum dependence, and we want to illustrate the transverse momentum dependence originating purely from the kinematical parts of the cross section and the BFKL Pomeron exchange. We observe that the cross section has a characteristic dip at values of transverse momentum around the vector meson mass p T ∼ M V . The dip appears because of a change of the sign of the partonic amplitude. Its appearance and evolution is well known from analyses of the diffractive vector meson production at large momentum transfer [5][6][7]. The dip appears at the lowest order for which the kinematical parts of the amplitudes for diffractive photoproduction and the triple gluon fusion hadroproduction are the same. For the diffractive photoproduction the dip moves to larger values of p T with increasing rapidity. In the hadroproduction where the BFKL evolution acts in a different way, the dip becomes shallower with increasing rapidity length of the BFKL evolution. This effect is especially visible at the highest plotted value of the LL evolution variablē α s y = 1. The fall-off with the transverse momentum is 1/p n T with power to be equal approximately n = 5 ÷ 6. In the next figure, Fig. 6 the result of the same calculation as a function of rapidity length of the evolved Pomerons is shown for fixed values of p T . The exponential growth with rapidity is clearly visible with the exponent being equal to approximately ∼ 0.5 which is consistent with the exchange of two LL Pomerons.
In order to better understand the properties of the Pomeron loop we also present calculations in which the integration of the momentum transfer q T is not performed. The results for  negativity of any observed differential cross sections. Of course, when integrated over the q T the cross section remains positive definite as it should be.

Numerical results
In this section we shall discuss the phenomenological implications of the correlated double-Pomeron exchange contribution to J/ψ and Υ hadroproduction. First, however, we will perform a numerical analysis of the diffractive photoproduction case described in Section 2 in the context of the most recent HERA data [1]. We treat this as a benchmark calculation and an opportunity to fix the (very few) theory parameters.
Our numerical calculations are performed in two ways. First, we use the Monte Carlo implementation of the cross section formulae which allows for a great flexibility in choice of distributions one can study. In particular it allows for usage of consistent observables with those measured. Second, we use direct implementations of the integrals, which serves as a cross-check and a useful tool for fast studies of smooth distributions.
The analysis of diffractive J/ψ photoproduction at √ S = 318 GeV [1] was performed by H1 collaboration and was found to be consistent with the previous measurements by ZEUS [2]. The latter was addressed in the literature at length in the BFKL context, see [4,5]. The newer data [1] cover larger γ-p c.m.s. energy range 50 < √ s < 200 GeV (vs 80 < √ s < 120 GeV at ZEUS) and much higher momentum transfers: 2 < |t| < 30 GeV 2 (vs |t| < 12 GeV 2 at ZEUS). We shall be interested in the cross section as a function of t.
Since we have the Monte Carlo implementation of the hadronic cross section, we can easily apply the kinematic cuts to be very close to the experimental ones. In particular, instead of calculating the cross section at fixed average γ-p c.m.s. energy as it was the case for previous studies [4,5], we can generate explicitly the equivalent photon flux [66] and calculate the weighted γ-p cross section. Above, m e is the electron mass and Q max ∼ M V . Following [1], we also restrict the mass of the diffractively produced state via the relation There are two main model parameters: the infrared cutoff s 0 and the slope of the trajectory in the evolution equation (2.19) given by the fixed strong coupling constant. The remaining freedom concerns: (i) the Λ parameter in Eq. (2.6) which we set to E T = p 2 T + M 2 V , (ii) the choice of the hard scale µ entering the strong coupling constant and the PDFs, (iii) the choice of the PDF set itself. Starting with the last point, we use the CT14nlo set [67]. The choice of the hard scale is, in principle, more subtle as α s may be running with a different scale in different components of the calculation. However, we simply set µ 2 = p 2 T + M 2 V as the scale for the coupling constant in the impact factors and vary the scale by a factor of two to estimate the theoretical uncertainty. In the LL BFKL kernel the fixed effective value of the strong coupling constantᾱ s = 0.14 is used. This value was adjusted to get a good description of the diffractive photoproduction data. It effectively parametrizes the effect of higher order corrections to the BFKL evolution. The results of calculations in this setup are presented in Fig. 9 together with H1 data. The sensitivity to the scale variation is visualized by the shaded boxes. We observe excellent agreement with data for the assumed parameters. In Fig. 10 we also show in one plot how different components of the model influence the results. In this study, the most drastic changes are generated either by decreasing the s 0 cutoff to 0.01 GeV 2 or by increasing the value of the intercept. Although in both cases the cross section gets increased significantly by a similar amount, the large |t| tail is different. The change of the Λ scale to the constant LL BFKL H1 data Figure 9. The differential cross section dσ/dt for the diffractive J/ψ production vs H1 data [1] for the coupling constant entering the LL BFKL trajectoryᾱ s = 0.14. CT14nlo, α s =0.14, s 0 =0.5, Λ=E T CT14lo, α s =0.14, s 0 =0.5, Λ=E T CT14nlo, α s =0.17, s 0 =0.5, Λ=E T CT14nlo, α s =0.14, s 0 =0.01, Λ=E T CT14nlo, α s =0.14, s 0 =0.5, Λ=M V Figure 10. Comparison of the differential cross sections dσ/dt for the diffractive J/ψ production for some configurations of the model parameters.
value equal to M V significantly hardens the spectrum. The change to the leading order PDFs has negligible effect.
Let us now turn to the double-Pomeron contribution to hadroproduction of J/ψ and Υ. The point we want to focus on is the p T -dependence of the vector meson in hadroproduction at large p T . This is because the various meson production mechanisms are characterized by distinct and robust resulting p T dependencies. This allows for the most stringent verification of the production mechanism and so far provides the strongest evidence for the dominance of the COM. Therefore we choose for the comparisons of the theoretical predictions with the experimental results, the dataset with the largest available p T range. These are measurements of prompt J/ψ and Υ production by the CMS collaboration in pp collisions at √ S = 13 TeV [68]. The kinematic setup for our calculation is adjusted to match this dataset. Accordingly, we set the c.m.s. proton-proton collision energy to 13 TeV and collect the Monte Carlo events in 2D bins of p T and rapidity Y of the vector meson reconstructed as Other parameters of the calculation are unchanged with respect to the diffractive photoproduction case (obviously, except the mass of the vector meson for the hadroproduction of Υ(1S) state). We make the slices on the 2D histogram in |Y | following the CMS setup.
Comparison of the two-Pomeron contribution to J/ψ hadroproduction and the respective data for four rapidity slices is shown in Fig. 11. The first striking observation is that the double-Pomeron contribution is very small in the kinematic region probed. On one hand, this is something one could somewhat naively expect from the Born level formula Eq. (3.5) due to the high power of strong coupling constant in the cross section, i.e. α 5 s . On the other hand, in principle, the Born cross section could gain strong enhancement due to the double BFKL evolution. It turns out, however, that the evolution length for the kinematics tested experimentally is not sufficiently long to compensate for the α 5 s suppression factor. To demonstrate this we calculate the distributions of the rapidity distance y, which actually sets the evolution length of the BFKL Green's functions. The results for the y-distribution in J/ψ proton-proton collisions at √ S = 13 TeV for typical meson p T = 30 GeV and rapidities Y = −0.6, −0.3, 0, 0.3, 0.6 are shown in Fig. 12. We note that, in this plot we only took one contribution to the cross section, i.e. the first term in Eq. (3.3). This was done in order to better illustrate the shift of the peak of the distributions with the changing rapidity of the vector meson, due to the change of the BFKL evolution. Of course in the full calculation the symmetry with respect to exchange of x 1 and x 2 is taken into account. The distributions peak around y ∼ 2, and decrease towards the kinematical cutoff at y ∼ 6 (depending on Y ). The mean evolution length in y is about 2 units, which is rather short and implies only moderate BFKL enhancement. A naive estimate using asymptotic enhancement factors of two the BFKL Pomeron exchange gives exp(2∆ P y) ∼ 4.5 for y = 2 and the LL BFKL intercept ∆ P 0.39 corresponding to the assumed effective valueᾱ s = 0.14 in the BFKL kernel. It is well known however, that preasymptotic effects lead to a decrease of the enhancement w.r.t. its asymptotic form.
The actual enhancement due to the evolution of the double BFKL ladder is shown in Fig. 14. We see that indeed the gain is roughly a factor of two. The moderate length y of  Figure 11. The double differential cross section dσ/dp T dY for the hadroproduction of the J/ψ, times the branching ratio for leptons B. Here Y is the rapidity of the J/ψ. The solid histograms represent the double-Pomeron contribution to the cross section calculated in this work. The error estimate (shaded boxes) comes from the scale variation. We also show the power-law behaviour C 1 /p 5 T and C 2 /p 7 T .
the BFKL evolution may be understood as result of a competition between x-dependencies of the PDFs that weight the partonic correlated 3GF cross section of the vector meson production off a quark or gluon, and the double BFKL evolution. Indeed, the PDFs grow towards small x as x −1−λ with λ ∼ 0.3, so they strongly push towards small x i values of the partons. The PDF factor ∼ (x 1 x 2 ) −1−λ squeezes the phase space for the partonic process and leads to only moderate average y.
We also study the hadronic cross section as a function of the transverse momentum q flowing in the double-Pomeron loop. In Fig. 13 we plot the contribution to the cross section, unintegrated in both p T and q T , similar to what has been done in Fig. 8 but at the hadronic level and within the kinematic window studied by the CMS experiment. We observe a non-negligible negative contribution to the integral coming from regions of non- zero momentum transfer flowing in the Pomeron loop. This is an exclusive property of the correlated two Pomerons and constitutes an additional source of the weaker cross section enhancement than naively expected. Note, that for any |p| slice the cross section will turn into a positive number upon integration over q.
The second observation regarding the comparison of the calculation with the experimental results in Fig. 11 is the power-like behavior. For the data we observe roughly 1/p 5 T fall-off, while the double-Pomeron contribution dies out faster, roughly like 1/p 7+ T with < 1. The obtained p T -dependence at the hadron level is significantly steeper than the parton level cross section Eq. (3.5) with a fixed coupling constant. This happens mostly because the strong dependence of the hadronic cross section to kinematic lower cutoffs x min i on the parton x values. One obtains x min i ∼ E T / √ S and these constraints lead to steeper dependence of the hadronic cross sections with E T . In addition α 5 s (E T ) in the impact factors leads to an additional suppression with increasing E T .
For completeness, we also perform the calculation for the Υ(1S) case. The result is presented in Fig. 15. For this case, all the previous observations apply, except the powerlike scaling of the correlated double-Pomeron contribution, which is now less accurate at lower values of p T . This is expected as in this region the effects of the second scale of the problem -the meson mass M V become relevant for Υ, that is much heavier than J/ψ.   Figure 14. The same theoretical distributions as in Fig. 11 but gathered on single plot and compared with the Born-level cross section.

Discussion
The analyses performed in this paper are based on the non-forward BFKL equation in the leading logarithmic approximation. The NLL corrections to the BFKL kernel are known to affect strongly the rapidity dependence of the amplitudes. We represent these effects by using an effective small value of the fixed coupling constant in the LL BFKL kernel. This approach is then confronted with the diffractive photoproduction data. The present paper is focused on estimating the correlated two Pomeron contribution to the inclusive  Figure 15. Same as in Fig. 11 but for hadroproduction of Υ(1S).
hadroproduction of heavy vector mesons. Therefore we use the diffractive photoproduction data to constrain the framework, and then apply it to the hadroproduction. We show that the J/ψ diffractive photoproduction data may be well described within the LL BFKL approach, so we do not find it necessary to consider the diffractive J/ψ photoproduction at large p T within the full next-to-leading logarithmic BFKL approximation. The tuned LL BFKL framework applied to inclusive hadroproduction gives a small contribution to the differential cross sections for J/ψ and Υ hadroproduction at large p T and does not explain the meson hadroproduction data. These main conclusions are not sensitive to details of the description. With this outcome it is not necessary to treat subleading effects, like the real parts of the BFKL exchanges or corrections due to off-diagonal effects -e.g. the Shuvaev factor [69]. In the present approach such correction factors may be efficiently absorbed into the effective value of α s driving the BFKL evolution. The strongest source of the theoretical uncertainty in the obtained results is the scale µ dependence of α s (µ) in the impact factors. The diffractive photoproduction cross section is proportional to α 4 s (µ), and the analyzed triple gluon fusion contribution to the hadroproduction is proportional to α 5 s (µ). Using the standard procedure of varying µ between the half and double of the natural process scale µ 0 = E T one finds rather large factor of about 2 (or more) of theoretical uncertainty due to the scale variation, see Fig. 9 for diffractive photoproduction, and similarly for the hadroproduction. These uncertainty factors are, however, correlated at different values of p T and they affect the overall normalization of the differential cross sections more than the shapes. Hence the agreement of the theoretical predictions and the data for the p T -shape of the diffractive photoproduction is rather robust. On the contrary -for the inclusive hadroproduction the correlated two Pomeron exchange leads to a much too steep p T dependence.
In a recent analysis [53] a good description of the direct J/ψ production in pp collisions at the LHC was obtained with the uncorrelated 3GF mechanism, without any other contributions included. Those results differ significantly from ones found in Ref. [51]. The   main reasons of the difference is using α s (m c ) in the impact factors in Ref. [53], while in Ref. [51] the running coupling constant α s (E T ) was used. This leads to a modification of the obtained transverse momentum distribution, because of a large span of the transverse momentum probed and the α 3 s dependence of the cross section. The choice of the scale µ = E T is also applied in the current paper, as it is related to the highest virtuality that occurs in the vector meson production vertex. In order to test the sensitivity of the results to the different choice of scale we also compute the cross section for the correlated 3GF mechanism of J/ψ hadroproduction in which the fixed α s (m c ) is used. The results are shown in Fig. 16 for the central rapidity bin |Y | < 0.3 . The effect of the change from α s (µ) to α s (m c ) is dramatic, as expected from the α 5 s dependence of the cross section. The obtained cross section is closer to the data at lower p T ∼ 20 GeV. The p T dependence is harder than it is for the running coupling, but still somewhat softer than the data.
Let us finally make connection of the computed amplitudes to the triple BFKL Pomeron vertex derived by Bartels [70,71]. In our approach two BFKL Pomerons are coupled to a single parton at the target side. The parton emerges from the target and undergoes the QCD evolution, which in the small x domain may be also described by the BFKL equation. For the gluon parton, this topology corresponds to some of the diagrams that enter the triple Pomeron vertex. In the computation we apply the Mueller-Tang approximation, that isolates the dominant contribution to the amplitude. A more complete treatment of this configuration should be possible using the results obtained for all cut triple Pomeron vertices [72].

Summary and conclusions
In this paper we have performed a detailed analysis of the contribution to the vector meson hadroproduction due to the correlated three gluon fusion mechanism (3GF). In this process, a vector meson, J/ψ or Υ in our analysis, is directly formed as a color singlet state from the interaction of three gluons. To be precise, in the process considered, one gluon originates from one hadron and two other gluons from another hadron. In the particular case of the correlated 3GF, we considered the contribution when two gluons came from the perturbative splitting of a parent quark or a gluon. The produced vector meson and the projectile quark and gluon were separated by a rapidity distance y which we modeled as an exchange of two cut non-forward BFKL Pomerons in the leading logarithmic approximation at high energy. Effectively this process contains a Pomeron loop contribution which is spanned between one collinear parton from one of the hadrons and the vector meson impact factor. For consistency and in order to fix the parameters as well as the normalization we have simultaneously calculated the diffractive photoproduction of vector mesons in deep inelastic scattering through the mechanism of exchange of the same BFKL Pomerons. The topology of this process, modulo color factors, is the same as the correlated 3GF vector meson hadroproduction. The important difference lies in the different cuts through the BFKL Pomerons: in the case of the diffractive DIS, one considers the diffractive cut, which leads to the rapidity gap between the vector meson and the target remnant, and in the case of the hadroproduction both BFKL Pomerons are cut. The appropriate color factors for the case of vector meson hadroproduction were obtained by the projection of the BKP states onto the color combination corresponding to two cut Pomerons and retaining terms leading in N c . The correlated Pomeron loop was numerically first analyzed at the level of the partonic cross section as a function of the transverse momentum of the produced vector meson and the rapidity. We found that, the cross section exhibits a dip in the transverse momentum, at the level of the Born cross section (i.e. without the BFKL evolution) which occurs at about the value corresponding to the vector meson mass. When the BFKL evolution is included, the dip is washed out and almost completely disappears for high values of rapidity. We checked that the rapidity dependence of the partonic cross section is consistent with the one expected from the two BFKL Pomerons, i.e. the rapidity dependence is exponential with the intercept that corresponds to twice the value of the intercept of the BFKL Pomeron. The momentum transfer dependence inside the Pomeron loop was also analyzed, and it turns out that the distribution in this variable has a sharp cutoff at the values which correspond to approximately half of the transverse momentum of the vector meson, i.e. q T ∼ p T /2, after which the distribution becomes slightly negative. This cutoff stems from the coupling of two Pomerons to the vector meson impact factor.
We have then evaluated the cross section for the diffractive photoproduction of J/ψ in DIS, for large values of momentum transfer. The non-forward BFKL evolution was solved numerically for fixed values of the strong coupling constant. Since the NLL corrections to the BFKL evolution are known to be large we have used the value of the strong coupling as a fitting parameter. The other freedom was in the choice of the infrared cutoff (mimicking the effective gluon mass due to the color confinement) inside the BFKL evolution and the choice of the scale in the running coupling in the impact factors. We also computed the cross section weighted by the equivalent photon flux. We obtained very good agreement with the experimental data from H1 collaboration as a function of the momentum transfer. Using the same parameters then, we have obtained essentially parameter free prediction for the 3GF correlated contribution to the vector meson hadroproduction. We made the comparison with the CMS data for J/ψ and Υ production at √ S = 13 TeV in bins of rapidity of the produced vector meson and as a function of the transverse momentum. We have found that this contribution is rather small, of the order of few percent at smallest values of the transverse momenta. The magnitude of the cross section is only slightly enhanced by the BFKL growth as compared with the Born cross section with just four gluon exchange, despite the fact that two Pomerons are exchanged. We further analyzed this, by evaluating the typical values of the rapidities probed in the CMS kinematics, and found that indeed the cross section is dominated by rather low rapidity intervals, and that the BFKL Pomeron exchange leads to about factor two enhancement over the Born cross section. The correlated 3GF contribution to the cross section falls off like 1/p 7 T , which is faster than the 1/p 5 T seen in the data, and is consistent with the evaluation of the Born cross section.
Several improvements could be made to the calculation presented in this paper. Nonforward BFKL evolution is known to NLL order, and this could be used to estimate better the value of the cross section. Combined with the NLL impact factors, this would be the first complete NLL evaluation of the correlated Pomeron loop contribution to this process.