Open charm production in Double Parton Scattering processes in the forward kinematics

We calculate the rate of double open charm production in the forward kinematics studied recently in the LHCb experiment. We find that the mean field approximation for the double parton GPD (Generalized parton distributions), which neglects parton - parton correlations, underestimates the rate by a factor of two. The enhancement due to the perturbative QCD correlation \12 mechanism which explains the rate of double parton interactions at the central rapidities is found to explain 60 $\div$ 80 \% of the discrepancy. We argue that the nonperturbative fluctuations leading to non-factorized (correlated) contributions to the initial conditions for the DGLAP collinear evolution of the double parton GPD play an important role in this kinematics. Combined, the two correlation mechanisms provide a good description of the rate of double charm production reported by the LHCb. We also give predictions for the variation of the \effs (i.e. the ratio of double and square of single inclusive rates) in the discussed kinematics as a function of $p_t$. The account for two correlation mechanisms strongly reduces sensitivity of the results to the starting point of the QCD evolution.


I. INTRODUCTION.
It is widely realized now that hard Multiple Parton Interactions (MPI) play an important role in the description of inelastic proton-proton (pp) collisions at the LHC energies where MPIs occur with probability of the order one in typical inelastic collisions.
Also, in the past several years a number of Double Parton Scattering (DPS) measurements in different channels in the central rapidity kinematics were carried out [17][18][19][20][21][22], while many Monte Carlo (MC) event generators now incorporate MPIs.
The recent discovery by the LHCb of the double charm DPS production attracted a lot of attention since it expands the study of multiparton dynamics into a new kinematics region of large rapidities [23][24][25][26], and since the background from the leading twist processes is very strongly suppressed in this kinematics [27][28][29].
The LHCb data are available for the J/ψ DPS production: J/ψ − DD and for the DPS production of two DD pairs. According to the LHCb experiment results, the DPS rate in the studied kinematics, which is customarily parameterized by 1/σ eff is practically the same for all channels and σ eff ∼ 20 mb (see Fig. 10 in [24]). The observed universality of σ eff is consistent with expectations of the approximation outlined below. Here, as usual, σ eff is defined as where σ 1,2 are cross sections of elementary 2 → 2 processes and σ 4 is a cross section of a process pp → 1 + 2 final state. We will focus on production of two DD pairs since the data for this channel have the smallest errors [30]. Also, more complicated mechanisms than the gg → J/ψ + X process may contribute in the case of J/ψ production, i.e. ggg → J/ψ (see e.g. [31] for a recent discussion).
It was pointed out starting with [4,32,33] that the rate of DPS calculated under assumption that partons in nucleons are uncorrelated (and using information about the gluon GPDs available from the analysis [32,33] of the HERA data) is too low to explain the data. It was pointed out in [10,11,15,16,34] that correlations generated in the course of the DGLAP evolution -1 ⊗ 2 mechanism -explain the DPS rates in the central rapidity region [10,11,15,16,34] provided the starting scale for the QCD evolution -Q 2 0 = 0.5 ÷ 1GeV 2 is chosen. The remaining problem seems to be a strong enhancement of the processes involving J/Ψ production [35,36] at √ s = 2TeV which does not show up in the LHCb data.
In this letter we demonstrate that the new LHCb data [23][24][25][26] corresponding to the forward kinematics can be explained by taking into account two effects: buildup with increase of Q 2 of the perturbative correlations -the 1 ⊗ 2 mechanism, calculated using DGLAP formalism [4,7,9,10] and soft small x parton -parton correlations in the nucleon wave function which result in a nonfactorized contribution to the initial conditions of the double parton GPD which can be estimated using information on diffraction in lepton / hadron -nucleon scattering following the ideas first presented in [7].
The paper is organized as follows. In section two we describe the kinematics of the LHCb experiment. In the third chapter we present that the mean field approximation results for the rate of DD production and demonstrate that they are a factor of two lower than the data. In section four we present results for the 1 ⊗ 2 mechanism contribution (see Fig. 2) to the cross section, In section five we discuss the Reggeon model based estimate of the non-factorized contribution to the initial conditions at Q 2 0 ∼ 0.5 − 1GeV 2 , and its Q 2 evolution. In section six we present general formula for σ eff combining the mean field,1 ⊗ 2 and nonperturbative non-factorized contributions.
In section seven we demonstrate that the simultaneous account of all three DPS mechanisms leads to the σ ef f values consistent with the data. The results are summarized in section eight. So far LHCb experiment has presented results for σ eff integrated over a significant range of rapidities and transverse momenta. So in our analysis we will first perform calculations for the typical LHCb kinematics and later on present the results for the variation of σ eff within the LHCb kinematic range which turns out to be pretty weak.

II. KINEMATICS OF THE LHCB STUDY OF THE DOUBLE CHARM PRODUCTION
In the case of production of two D mesons, the main mechanism is production of two pairs of DD mesons in two hard process (DPS) (see Fig.1) with two D (D s ) mesons originating from two DD pairs. D-mesons are observed in the rapidity interval y = 2 ÷ 5, the average rapidity interval between D andD meson is of the order ∆y = 0.5. A cutoff of p t ≥ 3GeV was introduced in the DPS analysis leading to the average transverse momenta of the D mesons of the order of p t ∼ 4 GeV.
Hence D-mesons are created in the interaction of two gluons with virtualities Q 2 ∼ (2p 2 t + m 2 c ) ∼ 34 GeV 2 . The factor of two takes into account the fragmentation of c → D in which D-mesons carry, in average, ∼ 0.75 fraction of the jet momentum [30] (see Fig.1).
The invariant mass squared of the created D meson pair is with a small uncertainty for the channels with the highest statistics.
The important advantage of these processes as compared to the processes experimentally studied before is that in this kinematics the SPS production of D-meson pairs is practically negligible [27][28][29] and the dominant process is the DPS production of cc pairs by gluons, thus permitting to use the methods developed in [4,7,9,10].
Similar calculations can be carried out for double bb pair production and bbcc pair production.
The only difference is that the corresponding transverse scale for b-pairs is Q 2 = m 2 b + 1.5p 2 t ∼ 50 GeV 2 where we shall take p t ∼ 4 GeV below as characteristic momenta. The corresponding invariant mass squared is of order 200 GeV 2 and x 1 ∼ 0.003, x 3 ∼ 0.014.

III. MEAN FIELD APPROXIMATION ESTIMATE OF e↵
Recall that in the mean field approach (see Fig.2 left ) double parton GPDs, describing the DPS, are where the one particle GPDs 1 D are known from the analyses [31,37] of exclusive J/ photoproduction at HERA. They are parametrized as Here D(x, Q 2 ) is the conventional gluon PDF of the nucleon, and F 2g ( , x) is the two gluon nucleon form factor. The e↵ective cross section e↵ is then given by We shall use exponential parametrization [37] F 2g ( , x) = exp( B g (x) 2 /2), where

III. MEAN FIELD APPROXIMATION ESTIMATE OF σ eff
Recall that in the mean field approach (see Fig.2 left ) double parton GPDs, describing the DPS, are where the one particle GPDs 1 D are known from the analyses [32,37] of exclusive J/Ψ photoproduction at HERA. They are parametrized as Here D(x, Q 2 ) is the conventional gluon PDF of the nucleon, and F 2g (∆, x) is the two gluon nucleon form factor. The effective cross section σ eff is then given by We shall use exponential parametrization [37] where gives a very similar numerical result for σ eff in our kinematics, decreasing σ eff by 4-5% which is well within the uncertainties of the current knowledge of the t-dependence of the gluon GPD in the studied x, Q 2 range ).
Integrating over ∆ 2 , we obtain for σ eff in the mean field approximation: where x i are the longitudinal momentum fractions of the four partons involved in the 2 ⊗ 2 mechanism. Hence we find for the mean field value of σ ef f in the LHCb kinematics which, as we already mentioned, is a factor of two larger than the the value reported by the LHCb.

IV. 3 TO 4 MECHANISM
The mechanism for the enhancement of the rate of DPS (increase of 1/σ ef f ) as compared to its mean field value was suggested in [4,7,9,10], where it was shown that taking into account the pQCD DGLAP ladder splits leads to a decrease of σ ef f -the 1 ⊗ 2 mechanism, see the right hand side of Fig. 2.
We calculate R by solving by iterations the evolution equation for 2 GP D [7,10,11]. The cross section due to the 1 ⊗ 2 mechanism is calculated as [7] 1 Note here that the 1 ⊗ 1mechanism contribution must be excluded [7]. Here 2 D 1 corresponds to 1 ⊗ 2 mechanism, while 2 D 2 to 2 ⊗ 2 contribution (with generic initial conditions -either factorized, or including non-factorized terms).
The distribution 2 D 1 corresponding to Fig. 2 (the right hand side) is obtained by solving by iterations of the evolution equation for 2 GP D [7,9,11]. The mean field distribution 2 D 2 gets corrections from the QCD evolution due to 1 ⊗ 2 mechanism. The DPS effective cross section is then parametrized as In [7,9] it was assumed that the factorized form given by eq. 3 is valid at the starting point of the evolution, Q 2 0 , which is essentially the parameter separating soft and hard dynamics. The enhancement coefficient increases with decrease of Q 2 0 . Numerical results for the enhancement coefficient R pQCD for charm pair production are given in Figs. 6,7 below. The direct calculations of σ eff pQCD = σ M F /(1 + R pQCD ) shows that the pQCD correlation leads to a decrease of σ eff to 24-28 mb for the double charm production, slightly larger reduction for charm + bottom, and a more significant reduction for double bottom.
Thus the 1 ⊗ 2 mechanism significantly improves the agreement with experimental data, but its relative contribution is smaller than in the central rapidity range covered by the ATLAS and CMS experiments.
There is an additional contribution to the DPS at small x which is absent in the case of processes involving x i ≥ 0.01 (production of jets, etc at the central rapidities). This contribution was first discussed in [9]. It results in a non-factorized contribution to 2 GPD at the initial scale Q 2 0 that separates soft and hard physics and which we consider as the starting scale for the DGLAP evolution. In the previous sections we assumed that at this scale 2 GPD factorizes into the product of two 1 GPDs. It is natural to expect that transition from soft to hard QCD regime is smooth and occurs at scales Q 2 ∼ 0.5 − 1 GeV 2 . In this case one expects that at such a scale the single parton distributions at small x below 10 −3 are given by the soft Pomeron exchange. In this picture the two soft partons may originate from two independent "multiperipheral ladders" represented by cut Pomerons, see Fig. 3. The soft Pomeron amplitude is practically pure imaginary [38] see also [39] We can write where factors x i /x take into account a smaller rapidity intervals occupied by the ladders in the case of transition to inelastic diffractive states. The factor corresponds to the cut Pomeron that splits into two Pomerons in diagram 4. It is equal to the product of the triple Pomeron vertex and square of proton -Pomeron residues, cf. [38,40]. Here we use α IP (0) = 1.1 corresponding to a soft effective Pomeron [41].
If x 1 = x 2 the right hand side of Eq.12 is equal to that is to the ratio of inelastic and elastic diffraction in DIS for the invariant γp energy s = m 2 0 /x. Using the triple reggeon parametrization of the cross section we can determine normalization of the three pomeron vertex C 3IP in eq. 13 from the HERA data [42,43] for the ratio of inelastic and elastic diffraction at t = 0 in the processes of vector meson production: The constant C 3IP is roughly the same for diffractive production of light mesons and J/ψ in a wide range of Q 2 , thus confirming the hypothesis of a smooth transition between soft and hard regimes. It is determined from the condition ρ(x 1 , x 1 , Q 2 0 ) = ω, where x 1 ∼ 0.001, that corrsponds to HERA data in [42,43]. Note here that to have a smooth connection with the low Q 2 gluon density model of GRV we take the x-dependence of gluon density at small x from this model. This may corresponds to relatively hard effective Pomeron in the lower legs though a priori density of partons in the Pomeron may grow more rapidly at small x than the overall Pomeron dominated amplitude.
In the Reggeon calculus [38] the effective triple Pomeron coupling is expected to decrease slowly with energy due to screening corrections somewhat reducing the rate of the increase of ω expected in the unscreened triple Pomeron model.
In any case, our procedure involves normalizing parameters of the model for x ∼ 10 −3 and studying a relatively narrow x range 10 −4 < x < 10 −2 . As a result our results are not sensitive to the variation of the Pomeron intercept between the soft and hard values.
Accordingly, for the parton density in the ladder we use: where the small x intercept of the parton density λ is taken from the GRV parametrization [44] for the nucleon gluon pdf at Q 2 0 at small x. Numerically λ(0.5 GeV 2 ) ∼ 0.27, λ(1.0 GeV 2 ) ∼ 0.31 Using eqs. 15, 12, 13 and the above values of λ(Q 2 0 ) we obtain C 3IP = 0.125 ± 0.025 GeV −2 for Q 2 0 = 0.5 GeV 2 , and C 3IP = 0.14 ± 0.025 GeV −2 for Q 2 0 = 1 GeV 2 . As a result we can estimate 2 D(x 1 , x 2 , Q 2 0 ) nf as where we introduced an additional factor of a = 0.1 in the limit of integration over x (or, equivalently, the limit of integration over diffraction masses M 2 ) to take into account that the Pomeron exchanges should occupy at least two units in rapidity, i.e. x > max(x 1 , x 2 )/0.1. The dependence on rapidity gap cutoff is weak, of order 10 %, and is present in all inelastic diffraction calculations [40]. The constant c 3IP = m 2 0 C 3IP , where m 2 0 = m 2 N = 1 GeV 2 is the low limit of integration over diffraction masses.
Consider now the t = −∆ 2 dependence of the above expressions. Strictly speaking all eqs. 11,12,17 have to include the explicit dependence on t. Here we shall however assume the factorization of the t-dependence, that reveals itself in the form where the function U does not depend on t and all t-dependence is given by the form factor F (t), for which we will use the exponential parametrization. Then we can use eqs. 11,12,17 at t = 0 (with corresponding functions, given by these equations, and the t-dependence given by the exponential form factors F (t)). Note that these form factors depend on x and the resolution scale only weakly and the scale dependence can be neglected while performing integration in eqs. 12,17.
Such factorization is known to work well for a pure diffraction case (diagram 4 for x 1 = x 2 , and we expect it to work in the general case as well.
The t-dependence of elastic diffraction is given by Thus the t dependence of the factorized contribution to 2 D f is given by where F 2g is the two gluon nucleon form factor.
The t-dependence of the non-factorized term eq. 17 is given by the t-dependence of the inelastic diffraction: exp((B in (x 1 ) + B in (x 2 ))t/2.).
Using the exponential parameterization exp(B in t) for the t-dependence of the square of the inelastic vertex pM X IP , the experimentally measured ratio of the slopes B in /B el 0.28 [42] translates into the absolute value B in = 1.4 ÷ 1.7 GeV 2 . A much weaker t-dependence of inelastic diffractive residue as compared to the elastic vertex is observed also for reaction pp → p + M X , see e.g. [45].
In the language of the Reggeon calculus this is a consequence of the well known observation that the t-dependence of three Pomeron vertex is much weaker than of the square of the ppIP vertex, see e.g. [40].
The evolution of the initial conditions, eq. 17, is given then by where G(x 1 /z 1 , Q 2 1 , Q 2 0 ) is the conventional DGLAP gluon-gluon kernel [46] describing evolution from Q 2 0 to Q 2 1 , Q 2 2 . In our calculations we neglect initial sea quark densities in the Pomeron at scale Q 2 0 (obviously Pomeron does not get contribution from the valence quarks). Let us define the quantity K (generalizing ρ from eqs. 11,12 to arbitrary Q 2 1 , Q 2 2 ): .
The nominator of this quantity is given by integral 21, while the denominator is a product of the conventional PDFs. We carried the numerical calculation of K for Q 2 0 = 0.5 GeV 2 and Q 2 0 = 1.0 GeV 2 . The typical results are presented in Fig. 5. (the corresponding x i are taken in accordance with analysis of section 2, and the calculations are carried out at t=0.).
One can see that K grows with increase of Q 2 0 and that the QCD evolution leads to the suppression of the nonperturbative contribution. We perform calculation neglecting the PPR (Pomeron-Pomeron-Reggeon) contribution. Inclusion of this term would increase the result by ∼ 10%. Overall we estimate the errors in the K-factor due to uncertainties in the input parameters are ∼ 25-35%.
The characteristic feature of K-factor is its increase as one considers more forward kinematics for charm production. Moreover, if we start from smaller x 1 , x 2 the rate of decrease of K with the increase of transverse momenta decreases. We illustrate these features in Fig. 5, where we consider K for the charm production kinematics described in section 2 (in Fig.5 Q 2 0 = 0.5 GeV 2 , the behaviour for Q 2 0 = 1 GeV 2 is similar. Upper curve is K for the 2 GPD with small x ∼ 10 −4 gluons and lower one -for larger x ∼ 10 −2 . We can see that main non-factorizable contribution originates from a smaller x gluon pair. The same is true for production of ccbb and bbbb. One can see from Fig.5 that K(x, Q 2 ) decreases strongly with increase of Q 2 . This reflects the increase of typical x at Q 2 0 scale contributing to K(x, Q 2 ) with increase of Q 2 and a fast decrease of K(x, Q 2 0 ) with increase of x (remember that K(x ≥ 0.05 − 10 −1 , Q 2 0 ) ≈ 0 and grows strongly with decrease of x less than 10 −2 . We can now write the general expression for σ eff taking into account non-factorized contribution to the initial conditions, the 1 ⊗ 2 mechanism and the mean field contribution.
Here B i ≡ B(x i ), and Also is the ratio of 2 GPD obtained from non-factorized and factorized terms at the scale Q 2 1 , Q 2 2 . After carrying out integration over ∆ 2 we obtain the expression for σ eff in terms of R pQCD , K, B el and B in . For simplicity we will write it only for the case of kinematics under considerations where the K term enters only for the partons with smaller x's.

VII. σ eff FOR PRODUCTION OF THE HEAVY QUARK PAIRS
.
We can now return to the analysis of the process of production of two charmed pairs. We consider the symmetric kinematics, i.e. x 1 ∼ x 2 ; x 3 ∼ x 4 .
In this case we can neglect terms proportional to K 34 since it corresponds to a negligible Regge mechanism contributions at x 3 , x 4 ∼ 0.01 ÷ 0.1, and in particular neglect S 12 -K 34 interference terms). Then we have Carrying out the integration we obtain for the full rescaling of σ eff including all three mechanisms discussed above: where R pQCD is the cross section enhancement due to 1 ⊗ 2 mechanism, i.e. proportional to is the enhancement due to nonperturbative correlations and interference of nonperturbative and perturbative contributions.
Note that the main sources of large R tot are the presence of the pQCD enhancement -1 ⊗ 2 for two partons with larger x and nonperturbative enhancement for smaller x's. The latter enhancement is amplified by the fact that the only ∆ 2 dependence in this case due to exp(−B in ∆ 2 ), whose slope is almost three times smaller than that of the mean field term, leading to the major enhancement of the corresponding contribution, compensating relatively small K. (The smallness of K is connected with a rapid decrease of the effect of non-perturbative correlations with the increase of Q 2 .) Thus the enhancement we obtain is essentially due to asymmetric (between upper and lower parts of diagram Fig. 2) kinematics of two pairs of x's.
altogether the Regge type contribution to R is ∼ 0.3,R pQCD ∼ 0.7 For Q 2 0 = 1GeV 2 we find the Regge contribution to R to be larger-∼ 0.4, while R pQCD ∼ 0.4, As a result for both choices of the initial conditions we obtain R ∼ 1.8 − 2., leading to σ eff ∼ 20 − 22mb (29) Note that numerically variation of the values of R pQCD , with a choice of the starting point of the Q 2 evolution is practically completely compensated by the variation of the soft non-factorizable contribution.
Note that eq.29 does not include additional uncertainties in the Reggeon calculation. For example, uncertainty in the ratio of inelastic and elastic diffraction of order 25% will lead to 19 − 23 mb in eq.29 and so on. There is a similar uncertainty due to the input t-dependence of the gluon GPDs.
The same calculation for the production of two bottom and two charm pairs in the LHCb kinematics [23,24] also gives R ∼ 1.9 − 2. In this case σ eff mean field ∼ 38mb, and we find σ eff ∼ 19 mb, in a good agreement with the LHCb data.
We show different contributions to σ eff enhancement as a function of the transverse momentum of D meson p t for 3.5 TeV and 6.5 TeV runs in figures 6 and 7. We see that the R pQCD slowly decreases with energy, but this is compensated with increase of R sof t , whose relative contribution also increases with the increase of energy.
The corresponding σ eff for two LHC runs are depicted in Figs. 8. We see that the σ eff increases by less than 1 mb for small p t when we move from 3.5 to 6.5 TeV, i.e. it effectively remains constant with the increase of energy, due to increase of soft correlations contribution compensating the decrease of pQCD contribution and increase of mean field σ M F eff . In fact of course such small changes are beyond the accuracy of our model, and we can only conclude that σ eff are approximately constant in this interval of energies for given transverse momenta p t .
We obtain very similar results for the production of two pairs of bb (Fig. 9). Note that in our approach the same σ eff are expected for production of two Υ and Υbb, cf discussion in section 2 of the case of charm production.

VIII. CONCLUSIONS
We have demonstrated that the rate of DPS of the production of two pairs of D mesons in the pp collisions in the forward kinematics studied by the LHCb can be explained by taking into account two types of correlations in the nucleon double GPD -the pQCD mechanism of [4,7,9,10] which allowed previously to describe the rate of DPS at the central rapidities and new nonperturbative correlation mechanism specific for small x which is related to the phenomenon of the inelastic diffraction.
Account for two correlation mechanisms significantly reduces sensitivity of the results to the starting point of the QCD evolution, both for forward and for central kinematics.  Fig.10 in [24]). We obtain similar results for other charm DPS production processes (2 J/P si, J/Ψ and DD pair), and this is indeed observed in experiment [23] in the forward kinematics (within experimental accuracy). For the DPS production of the bottom quarks we find (see Fig.9) σ eff ∼ 21-23 mb which is nearly a factor of two smaller than the mean field estimate of σ eff =38 mb. Thus we observe that combining pQCD correlation mechanism and the Regge inspired model for initial conditions we find approximately constant σ eff of order 20-22 mb for the LHCb kinematics.
Our calculations of σ eff were performed both for the 3.5 × 3.5 TeV 6.5 × 6.5 TeV runs. (The corresponding differences with 4 and 7 TeV runs respectively are negligible). We obtain practically the same values of σ eff since the decrease of R pQCD is compensated by increase if R soft . The actual difference is of order 1mb, slightly increasing to 2 mb (σ eff slightly decreases with increase of energy, but this change may be artifact of our model assumptions, i.e. it is obviously beyond the accuracy of our model).
Clearly, the role of soft correlations increases with the decrease of typical Bjorken x in the process. The same is true for transverse scale where the soft correlations start to be relevant, we see that it increases with energy. On the other hand the changes in the scale of pQCD and soft correlations tend to compensate each other with the increase of energy. This means that from the theoretical point of view it will be extremely helpful to carry the measurement of σ eff for new 6.5 TeV run at LHCb, as well as to measure dependence of σ eff on rapidity of the forward quark pair.
Obviously the calculations of soft correlations presented in this paper can be considered only as a semi-quantitative estimate. In particular this is connected with large uncertainty in the parameters of the model (see section 4) that are known with the accuracy of 25−30%, leading to corresponding inaccuracy in R sof t . Additional inaccuracy (although significantly reduced) is due to the choice of Q 2 0 scale. Nevertheless, our estimate clearly reveals importance of soft correlations in forward kinematics and increase of their contribution with the energy of the collision.
Finally , let us note that our model can be used also for the central kinematics, where in particular it can be applied to calculate σ eff in the underlying event (UE). Preliminary results show that it will not influence significantly the MC simulations of UE given in [15,16], although it may lead to stabilization of σ eff in the region of small p t characteristic for UE. The detailed results for the central kinematics, as well as comparison of our predictions for σ eff with that of [47], that developed a different model also based on observations of [9] and applied it to the central kinematics, will be given elsewhere [48].