title: Revisiting the production of $J/ψ$ pairs at the LHC
authors: A. A. Prokhorov A. V. Lipatov M. A. Malyshev S. P. Baranov
text: 
  
  : Semantic API Migrations
      Michael F.P. O'Boyle
      December 30, 2023
  =========================
  
  
  
  
  
  
  ^1Faculty of Physics, Lomonosov Moscow State University, 119991 Moscow, Russia
  
  ^2Joint Institute for Nuclear Research, 141980, Dubna, Moscow region, Russia
  
  ^3Skobeltsyn Institute of Nuclear Physics, Lomonosov Moscow State University, 119991, Moscow, Russia
  
  ^4P.N. Lebedev Institute of Physics, Moscow 119991, Russia
  
  
  
  
  
  
  
  
  
  
  
  Abstract 
  
  
  
  
  
  
  We consider the prompt double J/ψ production in pp collisions at the LHC
  in the framework of k_T-factorization QCD approach.
  Using the fragmentation mechanism, we evaluate the color octet 
  contributions to the production cross sections 
  taking into account the combinatorial effects of multiple gluon radiation
  in the initial state
  driven by the Ciafaloni-Catani-Fiorani-Marchesini 
  evolution equation.
  We demonstrate the importance of these contributions 
  in a certain kinematical region
  covered by the CMS and ATLAS measurements.
  On the other hand,
  the experimental data taken by the LHCb Collaboration
  at forward rapidities and 
  
  moderate transverse momenta
  
  can be described well by O(α_s^4) color singlet 
  terms and contributions from the double parton scattering mechanism.
  The extracted value of the effective cross section σ_ eff = 17.5 ± 4.1 mb 
  
  is compatible
  with many other estimations based on different final states.
  
  
  
  
  Keywords: charmonia, non-relativistic QCD, small-x, CCFM evolution, double parton scattering.
  
  
  
  
  
  § INTRODUCTION
   
  
  Prompt production of J/ψ meson 
  pairs at high energies is a very intriguing subject of studies<cit.>.
  It provides
  a unique laboratory to investigate the quarkonia production
  mechanisms predicted by the non-relativistic QCD (NRQCD) 
  factorization<cit.>,
  which is a rigorous framework for the description
  of heavy quarkonia production or decays.
  The NRQCD implies a separation of perturbatively 
  calculated short distance cross sections for
  the production of a heavy quark pair in an intermediate 
  Fock state ^2S+1L_J^(a) with spin S, 
  orbital angular momentum L, total angular momentum J and 
  color representation a from its subsequent non-perturbative transition 
  into a physical quarkonium via soft gluon radiation.
  The latter is described by the 
  long-distance non-perturbative matrix elements (LDMEs), 
  which obey certain hierarchy in powers of the relative heavy quark velocity v<cit.>.
  At the next-to-leading order (NLO), NRQCD can explain the LHC data on 
  the prompt J/ψ, ψ^' and χ_c transverse momentum distributions (see, for example,<cit.>).
  However, it has a long-standing challenge in the 
  J/ψ and ψ^' polarization and provides inadequate 
  description<cit.> of the η_c production data[One possible solution, which,
  however, implies certain modification 
  of the NRQCD rules, has been proposed<cit.> (see also<cit.>).] (see also discussions<cit.>).
  Studying the J/ψ meson pair production
  can shed light on the puzzling aspects above
  since cc̅ bound state formation takes place here twice.
  
  In the last few years, significant progress has been made
  in the NRQCD evaluations of prompt double J/ψ 
  production. 
  The complete leading-order (LO) calculations,
  including both the color singlet (CS) and color octet (CO) terms,
  were done<cit.>. The relativistic corrections to the J/ψ pair
  production are carried out<cit.>.
  The NLO contributions to the CS mechanism
  are known<cit.> and partial tree-level NLO^* contributions to
  the both CS and CO terms were calculated<cit.>.
  The latter were found to be essential for both low and 
  large transverse momenta, as compared to the LO results[At the moment, full NLO NRQCD predictions
  for double J/ψ production
  are not available yet.].
  However,
  being comparable with the LHCb measurements<cit.>, 
  all these evaluations  
  have sizeble discrepancies with the latest CMS<cit.> and ATLAS<cit.> data, 
  especially at large 
  transverse momentum p_T(J/ψ, J/ψ),
  invariant mass m(J/ψ, J/ψ) and rapidity 
  separation Δ y(J/ψ, J/ψ) of the J/ψ pairs.
  For example, the CMS data
  are underestimated by the NRQCD predictions with a factor of about 10<cit.>. The difference between the 
  theoretical calculations and more recent ATLAS data at large 
  p_T(J/ψ, J/ψ) or m(J/ψ, J/ψ)
  is typically smaller but still essential.
  It was argued<cit.> that new processes or mechanisms are needed to better
  describe the LHC data.
  
  At large invariant mass m(J/ψ, J/ψ) the processes with large angular separation
  between the J/ψ mesons could play a role.
  One of such processes are the gluon or quark fragmentation 
  shown in Fig. <ref>. The gluon fragmentation into ^3S_1^(8) intermediate state 
  scales as 1/p_T^4
  and govern the single J/ψ production at 
  high transverse momenta (see, for example,<cit.> and references therein).
  In the case of J/ψ pair production,
  such terms were found to be negligible since they 
  are suppressed by powers of QCD coupling 
  α_s<cit.>.
  However, at large p_T(J/ψ, J/ψ) or m(J/ψ, J/ψ)
  one can expect a sizeble combinatorial contribution
  to the fragmentation yield
  from the multiple gluon radiation originating
  during the QCD evolution of the initial gluon cascade.
  The latter determines the perturbative QCD corrections to the 
  production cross sections at 
  high energies, which can be effectively taken into account 
  using the Ciafaloni-Catani-Fiorani-Marchesini (CCFM)
  evolution equation<cit.>.
  Main goal of our study is to clarify this
  point and investigate the role of combinatorial cascade gluon
  fragmentation contributions to 
  the double J/ψ production
  in different kinematical 
  regimes at the LHC.
  
  
  
  Our other goal is connected
  with the investigation of additional 
  
  production mechanism, double parton scattering (DPS), which is
  widely discussed in the literature at present (see, for example,<cit.> and references therein).
  Apart from the single parton scattering (SPS), where J/ψ meson pair
  is produced in a single gluon-gluon collision, DPS events originate from 
  two independent parton interactions. 
  Studying the DPS mechanism is of great importance since 
  it can help in understanding various backgrounds in searches for new 
  physics at the collider experiments.
  Despite the relative low total production rate, the DPS mechanism is expected to be 
  important for double J/ψ production at forward rapidities<cit.>.
  Therefore, the latter can be used to determine the DPS key parameter, the effective 
  cross section σ_ eff, which is related to the transverse overlap function between partons in the proton and 
  supposed to be universal for all processes with different kinematics and energy scales.
  Most of the measured values of σ_ eff lie between 12 and 20 mb (see, for example,<cit.>).
  However, somewhat lower value σ_ eff = 8.8 - 12.5 mb
  was extracted from the latest LHCb data on J/ψ pair production
  within the NRQCD<cit.>.
  Moreover, the values σ_ eff = 8.2 ± 2.2 mb<cit.>, σ_ eff = 6.3 ± 1.9 mb<cit.>,
  σ_ eff = 4.8 ± 2.5 mb<cit.> and even σ_ eff = 2.2 ± 1.1 mb<cit.>, 
  σ_ eff = 2.2 - 6.6 mb<cit.> were obtained from 
  recent Tevatron and LHC experiments.
  Below we will try to extract the effective cross section σ_ eff from  
  combined analysis of the LHCb data<cit.> on double J/ψ production taken at √(s) = 7 and 13 TeV.
  
  To calculate the physical cross sections we use the k_T-factorization 
  approach<cit.>. 
  We see certain technical advantages in the fact that, even 
  with the LO amplitudes for hard subprocesses, one can include
  a large piece of higher-order pQCD corrections (NLO + NNLO + ...) taking them into 
  account in the form of CCFM-evolved Transverse Momentum Dependent (TMD)
  gluon densities in a proton[The description of this approach 
  can be found, for example, in review<cit.>.].
  In this way we 
  preserve consistency with our previous studies<cit.> and
  automatically incorporate the wanted effects of
  initial state gluon radiation.
  To reconstruct the CCFM evolution ladder,
  that is the key point of our consideration,
  we employ the TMD parton shower routine implemented into the Monte-Carlo 
  event generator cascade<cit.>.
  The k_T-factorization approach can be considered as a convenient 
  alternative to explicit high-order calculations
  in the collinear DGLAP-based scheme. The situation in J/ψ pair
  production is specific since calculating even the LO hard scattering amplitudes 
  is already complicated enough, so that extending to higher orders
  seems to be a rather cumbersome task. Thus, the k_T-factorization 
  remains the only way open to potentially important higher-order effects.
  To evaluate the DPS contributions to the 
  double J/ψ production we will use the results of 
  our previous analysis<cit.>.
  
  The outline of the paper is the following. In Section 2 we briefly
  describe the basic steps of our calculations. 
  In Section 3 we present the numerical results and discussion.
  Our conclusions are summarised in Section 4.
  
  
  
  § THE MODEL
   
  
  The neccessary starting point of our consideration is related with CS contribution
  to the double J/ψ production, that
  refers to O(α_s^4) gluon-gluon fusion subprocess
  
      g^*(k_1) + g^*(k_2) → cc̅[^3S_1^(1)](p_1) + cc̅[^3S_1^(1)](p_2),
  
  where the four-momenta of all particles are indicated in the parentheses.
  Some typical Feynman diagrams are depicted in Fig. <ref>.
  It is important that both 
  initial gluons are off mass shell.
  That means that they have
  non-zero transverse four-momenta k_1T^2 = - 𝐤_1T^2 ≠ 0 and
  k_2T^2 = - 𝐤_2T^2 ≠ 0 and an admixture of longitudinal 
  component in the polarization four-vectors (see<cit.> for more information).
  The corresponding off-shell (k_T-dependent) production amplitude contains widely used projection 
  operators for spin and color<cit.> 
  which guarantee the proper quantum numbers of final state charmonia.
  Below we apply the gauge invariant expression obtained earlier<cit.>.
  The derivation steps are explained in detail there.
  The respective cross section
  can be written as
  
      σ(pp → J/ψ J/ψ + X) = ∫1 16 π (x_1 x_2 s)^2 |A̅(g^* g^* → J/ψ J/ψ)|^2 ×
      × f_g(x_1,𝐤_1T^2,μ^2) f_g(x_2,𝐤_2T^2,μ^2) d𝐤_1T^2 d𝐤_2T^2 d𝐩_1T^2 dy_1 dy_2 dϕ_1 2πdϕ_2 2πdψ_1 2π,
  
  where ψ_1 is the azimuthal angle of outgoing J/ψ meson,
  ϕ_1 and ϕ_2 are the azimuthal angles of  
  initial gluons having the longitudinal momentum fractions x_1
  and x_2, y_1 and y_2 are the center of mass rapidities 
  of produced particles and 
  f_g(x,𝐤_T^2, μ^2) is the TMD 
  gluon density in a proton taken at the scale μ^2.
  
  
  
  In addition to the CS terms above, we have considered some 
  of CO contributions using the fragmentation approach.
  At high transverse momenta, p_T ≫ m_ψ,
  large logaritmic corrections proportional to 
  α_s^n ln^n p_T/m_ψ occur and, therefore, 
  description in terms of fragmentation 
  functions (FFs), evolving with the energy scale μ^2, appears to be appropriate.
  In general, the FF D_a^ H(z,μ^2) 
  describing the transition of parton a into the 
  charmonium state H 
  can be expressed as follows (see, for example,<cit.> and references therein):
  
      D_a^ H(z,μ^2) = ∑_nd_a^n(z,μ^2) ⟨ O^ H [n]⟩,
   
  where n labels the intermediate (CS or CO) state of charmed quark pair
  produced in the hard parton interaction
  and ⟨ O^ H [n]⟩ are the corresponding LDMEs.
  In the leading logarithmic approximation, 
  g^* → cc̅[^3S_1^(8)] transition is the only one
  giving a sizeble contribution 
  to S-wave charmonia production 
  at p_T ≫ m_ψ<cit.>, so that the 
  cross section of inclusive single J/ψ
  production in pp collisions could be approximately calculated as
  
      dσ(pp → J/ψ + X)  dp_T = ∫ dz dσ(pp → g^*)  dp_T^(g^*) d_g^[^3S_1^(8)](z,μ^2) ⟨ O^J/ψ [^3S_1^(8)] ⟩,
   
  where p = z p^(g^*) and p^(g^*)
  are the outgoing J/ψ meson and intermediate gluon momenta. One can easily obtain
  
      σ(pp → g^*) = ∫π x_1 x_2 s λ^1/2(m_ψ^2,k_1^2,k_2^2) |A̅(g^* g^* → g^*)|^2 ×
      × f_g(x_1,𝐤^2_1T,μ^2) f_g(x_2,𝐤^2_2T,μ^2) d𝐤^2_1T d𝐤^2_2T dy dϕ_1  2πdϕ_2  2π,
   
  where p^(g^*) = k_1 + k_2 and λ(m_ψ^2, k_1^2, k_2^2) is the known kinematical function<cit.>.
  Evaluation of the off-shell production amplitude 
  |A̅(g^* g^* → g^*)|^2 = (3/2) πα_s(μ^2) |𝐩_T^(g^*)|^2
  is an extremely straightforward and, in our opinion, needs no explanation.
  We only note that, according to the k_T-factorization prescription<cit.>,
  the summation over the polarizations of initial off-shell gluons is carried out with
  ∑ϵ^μϵ^* ν = 𝐤_T^μ𝐤_T^ν/𝐤_T^2.
  In the collinear limit 𝐤_T → 0 this expression converges to the ordinary one after averaging on the
  azimuthal angle.
  
  The key point of our consideration is
  that the gluon, produced in the hard scattering
  and fragmented into the J/ψ meson according to main formula (4), 
  is accompanied by a number of gluons radiated during the non-collinear QCD evolution,
  which also give rise to final J/ψ mesons.
  Thus, taking into account all their possible combinations into 
  the meson pairs, one 
  can calculate corresponding gluon fragmentation
  contribution to the double J/ψ
  production up to all orders in the pQCD expansion.
  At high energies, the QCD evolution of gluon cascade 
  can be described by the CCFM equation<cit.>, which smoothly 
  interpolates between the small-x BFKL gluon dynamics 
  and high-x DGLAP one, and, therefore, provides us with the suitable tool 
  for our phenomenological study.
  To reconstruct the CCFM evolution ladder,
  we generate a Les Houches Event file<cit.>
  in the numerical calculations according to (4) and (5) and then process the file with a 
  TMD shower tool implemented into the Monte-Carlo event generator cascade<cit.>.
  This approach gives us the possibility to take into account
  the contributions from initial state gluon emissions in a
  consistent way (see also<cit.>).
  
  Of course, the same scenario can be applied to fragmentation
  of charmed quark pairs into J/ψ mesons. 
  So, one can first simulate the perturbative production of c c̅ pair 
  in the off-shell gluon-gluon fusion and then reconstruct the 
  CCFM evolution ladder using the cascade tool.
  After that, one can easily produce J/ψ pairs by taking into account 
  all possible combinations of mesons
  originating from the charmed quarks and/or cascade gluon fragmentation. 
  Unlike the conventional (collinear) QCD factorization, where only fragmentation of both charmed quarks 
  into J/ψ mesons gives contribution, the model above 
  can lead to increase in the double J/ψ production 
  cross section due to additional combinatorial contributions from
  gluons and quarks.
  
  The charm and gluon FFs at the any scale μ^2, 
  D_c^J/ψ(z,μ^2) and D_g^J/ψ(z,μ^2), 
  can be obtained by solving 
  the LO DGLAP evolution equations:
  
      d/dlnμ^2([ D_c; D_g ]) = α_s(μ^2)/2π([ P_qq P_gq; P_qg P_gg ]) ⊗([ D_c; D_g ]),
   
  where P_ab are the standard LO DGLAP splitting functions. 
  The initial conditions for these FFs are calculated with<cit.>
  
      d_g^[^3S_1^(8)](z,μ_0^2)= α_s(μ_0^2)/24 m_c^3πδ(1 - z),
   
  
      d_c^[^3S_1^(1)](z,μ_0^2)= α_s^2(μ_0^2)/m_c^316z(1-z)^2  243(2-z)^6(5z^4-32z^3+72z^2-32z+16),
   
  where starting scale μ_0^2 = m^2_ψ. 
  As it was noted above, we keep only the leading 
  contributions to corresponding FFs (see, for example,<cit.> and references therein).
  According to the non-relativistic QCD approximation, we
  set the charmed quark mass to m_c = m_ψ/2
  and then solve the DGLAP equations (6) numerically.
  The obtained charm and gluon FFs, D_c^J/ψ(z,μ^2) and D_g^J/ψ(z,μ^2),
  are shown in Fig. <ref> as functions of z for several values of scale μ^2. 
  Using these FFs, we reproduce well the 
  results of calculations performed with the Monte Carlo 
  event generator pegasus<cit.> (see Fig. <ref>).
  
  Finally, we turn to the DPS contribution to the double J/ψ production.
  We apply a commonly used
  factorization formula (for details
  see reviews<cit.> and references therein):
  
      σ_ DPS(pp → J/ψ J/ψ + X) = 1 2σ^2(pp → J/ψ + X) σ_ eff,
  
  where factor 1/2 accounts for two identical particles in the final state. The 
  effective cross section σ_ eff can be considered as 
  a normalization constant which incorporates all "DPS unknowns" 
  in to a single phenomenological parameter.
  Derivation of the factorization formula (9) relies on the two approximations: first, 
  the double parton distribution function can be decomposed into 
  longitudinal and transverse components and, second, the longitudinal 
  component reduces to the diagonal product of two independent single parton 
  densities. The latter is generally acceptable for such collider experiments 
  where small-x values are probed.
  The typical values of the variable x
  in the considered process are of order
  x ∼ (2m_ψ^2 + p_T^2)^1/2/√(s)∼ 10^-3,
  that approximately corresponds to the kinematical region
  of CMS<cit.>, ATLAS<cit.> and even LHCb<cit.> measurements
  (due to relatively small invariant mass of produced J/ψ pair, see discussion below).
  Therefore, one can safely omit the kinematical constraint<cit.> often applied 
  at the edge of phase space[Phenomenological consequences of the kinematical constraint<cit.>
  at the large x were investigated<cit.>.].
  Detailed description of evaluation of inclusive cross section
  σ(pp → J/ψ + X) in the k_T-factorization approach supplemented with the 
  NRQCD formalism can be found<cit.>.
  
  In the numerical calculations below we will use TMD gluon density in a proton
  obtained<cit.> from the numerical solution of CCFM evolution equation (namely, A0 set), 
  where the input parameters have been fitted to the proton structure function F_2(x,Q^2).
  At present, the A0 gluon distribution function is widely used in the phenomenological 
  applications[A comprehensive collection of the TMD gluon 
  densities can be found in the tmdlib package<cit.>,
  which is a C++ library providing a framework and an interface to the different parametrizations.]
  (see, for example,<cit.>).
  The renormalization and factorization scales, μ_R and μ_F,
  were set to μ_R^2 = μ_F^2 = ŝ + 𝐐^2_T, where ŝ = (k_1 + k_2)^2 and 𝐐_T^2
  is the transverse momentum of initial off-shell gluon pair. 
  This choice is dictated mainly by the CCFM 
  evolution algorithm (see<cit.> for more information).
  As it is often done, 
  the fragmentation scale μ_ fr is choosen to be equal to μ_ fr = m_T, the 
  transverse mass of fragmenting parton.
  We use the one-loop formula for the QCD coupling α_s
  with n_f = 4 active quark flavors at Λ_ QCD^(4) = 250 MeV.
  Following<cit.>, we set the J/ψ meson mass m_ψ = 3.097 GeV.
  We take corresponding CS LDME from the known
  J/ψ→μ^+μ^- decay width: ⟨ O^J/ψ [^3S_1^(1)]⟩ = 1.16 GeV^3<cit.>.
  
  
  
  
  
  
  
  
  
  
  § NUMERICAL RESULTS AND DISCUSSION
   
  
  We are now in a position to present the results of our simulations.
  First we discuss the role of cascade gluon fragmentation
  
  in different kinematical regimes, which correspond to 
  the CMS, ATLAS and LHCb experiments.
  
  In Fig. <ref> we show the differential cross sections of double J/ψ 
  production calculated as a functions of J/ψ pair invariant mass m(J/ψ, J/ψ) and 
  difference in rapidity between the J/ψ mesons |Δ y(J/ψ, J/ψ)| 
  
  
  at √(s) = 13 TeV.
  We have required p_T(J/ψ) > 10 GeV for both produced mesons,
  that ensures the validity of the fragmentation approach used.
  Moreover, this restriction close to the CMS or ATLAS conditions.
  One can see that an accurate account of combinatorial
  contributions, originated from the cascade gluon fragmentation into
  the J/ψ mesons (labeled as "fragm. comb."), significantly (up to an order of magnitude) 
  increase the cross section compared to the 
  single gluon fragmentation, governed by the LO gluon-gluon
  fusion subprocess[Here we have applied the MMHT'2014 (LO) parton density set<cit.>.]
  (labeled as "fragm. coll.").
  For the latter, we reproduce the results<cit.>.
  To highlight the importance of the combinatorial gluon fragmentation, 
  we show the results obtained using the simplified selection of J/ψ pair in each event,
  where one of the J/ψ mesons is originated from the gluon produced in the 
  hard scattering subprocess and another one is produced from the leading cascade gluon (labeled as "fragm. lead.").
  This selection criterion almost corresponds to the 
  collinear limit, as it is clearly demonstrated in Fig. <ref>.
  Next, we find that the cascade gluon fragmentation plays a dominant 
  role at large invariant masses m(J/ψ, J/ψ) ≥ 25 GeV
  and |Δ y(J/ψ, J/ψ)| ≥ 1,
  where it greatly overestimates the CS contributions.
  Taking into account these combinatorial contributions
  results in the 
  drastical rise of the double J/ψ production cross sections 
  at large m(J/ψ, J/ψ), 
  where the strong discrepancy between 
  the NRQCD estimations (including both the CS and CO terms) 
  and experimental data, taken by the CMS and ATLAS Collaborations, is observed.
  Contrary, the combinatorial fragmentation effects should be significantly less pronounced
  at forward rapidities, which are covered by the LHCb measurements.
  To demonstrate it, we have repeated the calculations 
  under the requirements 4.5 < p_T(J/ψ) < 10 GeV and 2 < y(J/ψ) < 4.5.
  
  The upper limit of p_T(J/ψ) is set to be the same as in LHCb analyses<cit.> 
  while lower limit corresponds to the region, where the fragmentation approach is 
  valid.
  Our results for distributions in 
  m(J/ψ,J/ψ) and |Δ y(J/ψ,J/ψ)|
  are shown in Fig. <ref>.
  One can see that the cascade gluon fragmentation 
  gives only small contribution to the forward J/ψ pair 
  production and, in principle, can be safely neglected. 
  It can be easily understood 
  since at large rapidities (or, equivalently, at
  large momentum fraction x of one of the interacting gluons) the gluon 
  emissions in the initial state are insufficient. 
  
  Concerning the contributions from charm fragmentation,
  their role (compared to the LO predictions of conventional pQCD) is also a
  bit enhanced due to the multiple gluon emissions in
  the initial state. 
  We find that these processes amount several percent of 
  the J/ψ pair production cross section 
  (see Figs. <ref> and <ref>) and, of course, can be considered as additional non-leading
  terms[To generate cc̅ events 
  in the off-shell gluon-gluon fusion the Monte-Carlo event generator pegasus<cit.> has been used.].
  
  Thus, we have shown that taking into account the combinatorial contributions
  from the cascade gluon fragmentation
  could fill the gap between the NRQCD
  predictions and experimental data. 
  However, to perform the quantitative comparison with the available CMS<cit.> and ATLAS<cit.> measurements
  one has to include into the analysis a number 
  of other possible fragmentation channels
  playing role at low and moderate transverse momenta.
  Moreover, additional feeddown contributions 
  to the double J/ψ production 
  from the χ_c and ψ^' decays should be taken into account.
  An accurate theoretical description
  requires a rather long-time numerical calculations. 
  So, here we only claim the possible importance of the combinatorial
  fragmentation terms above and left their 
  further cumbersome analysis for a forthcoming dedicated study.
  
  Now we turn to available LHCb data
  collected at √(s) = 7 and 13 TeV<cit.>. 
  These data refer to p_T(J/ψ) < 10 GeV, m(J/ψ, J/ψ) < 15 GeV 
  and forward rapidity region, 2 < y(J/ψ) < 4.5.
  Since the combinatorial contributions from gluon and/or charmed 
  quark fragmentation are almost negligible there,
  only the CS terms and DPS production mechanism play the role.
  The latter give us the possibility to easily 
  extract the key parameter of DPS mechanism, the effective cross section σ_ eff, 
  from the LHCb measurements.
  The feeddown contributions from radiative χ_c and ψ^' decays 
  to the SPS cross section, which is governed by the subprocess (1),
  are also unimportant 
  at small transverse momenta 
  and invariant mass m(J/ψ, J/ψ), see discussions<cit.>.
  Thus, we neglect below all these terms for simplicity.
  To evaluate the DPS contribution
  to the J/ψ pair production
  we use the results of our previous studies
  and strictly follow the approach<cit.>
  for the inclusive cross section
  σ(pp → J/ψ + X), entering to the DPS factorization formula (9).
  So, the determination of σ_ eff can be performed 
  in a self-consistent way.
  
  The following kinematical variables have been investigated in the LHCb analyses<cit.>:
  transverse momentum p_T(J/ψ, J/ψ), rapidity y(J/ψ,J/ψ) and invariant mass of the J/ψ pair,
  transverse momentum and rapidity of J/ψ mesons, 
  differences in the azimuthal angle |Δϕ(J/ψ, J/ψ)|
  and rapidity |Δ y(J/ψ, J/ψ)|
  between the produced mesons and 
  transverse momentum asymmetry A_T, defined as
  
      A_T = | p_T(J/ψ_1) - p_T(J/ψ_2)  p_T(J/ψ_1) + p_T(J/ψ_2)|.
  
  The measurements have been performed for p_T(J/ψ,J/ψ) > 1 GeV, p_T(J/ψ,J/ψ) > 3 GeV
  and in the whole p_T(J/ψ,J/ψ) range.
  We consider σ_ eff as an independent parameter
  and perform a simultaneous fit to the LHCb data.
  The fitting procedure was separately done for each of the measured kinematical 
  distributions employing the fitting algorithm as implemented in the commonly used gnuplot package<cit.>.
  
  Not all of the existing data sets are equally informative for the σ_ eff extraction. 
  Using the data where the DPS contribution is smaller 
  than the uncertainty of the "main" contribution would only increase the total 
  error. So, our fit is based on the following distributions (all measured at
  √(s) = 7 TeV and 13 TeV):
  single J/ψ transverse momentum p_T(J/ψ);
  single J/ψ rapidity y(J/ψ);
  invariant mass m(J/ψ,J/ψ);
  transverse momentum of J/ψ pair;
  rapidity of J/ψ pair;
  transverse momentum asymmetry A_T;
  rapidity separation between the two J/ψ mesons |Δ y(J/ψ,J/ψ)|.
  For all observables except |Δ y(J/ψ,J/ψ)| we used the data 
  without cuts on p_T(J/ψ,J/ψ) and with p_T(J/ψ,J/ψ) > 1 GeV,
  while for the rapidity separation |Δ y(J/ψ,J/ψ)| we used the sets
  without cuts on p_T(J/ψ,J/ψ), with p_T(J/ψ,J/ψ) > 1 GeV, and with 
  p_T(J/ψ,J/ψ) > 3 GeV.
  
  The obtained mean-square average of the fitted values is σ_ eff = 17.5 ± 4.1 mb,
  where corresponding uncertainty is estimated in the conventional way using Student's 
  t-distribution at the confidence level P = 95%.
  
  Here we achieve a remarkable agreement 
  with the majority of other σ_ eff estimations  
  based on different final states, such as, for example, W + 2 jets<cit.>,
  2 γ + 2 jets<cit.>, γ + 3 jets<cit.>, 
  4 jets<cit.>, J/ψ + D^+, J/ψ + D^0, J/ψ + Λ_c^+<cit.>, 
  Υ(1S) + D^0<cit.>. Thus, our result  
  supports the expectation about the universality of this parameter for a wide
  range of processes with essentially different kinematics, energies and hard scales.
  The obtained value of σ_ eff significantly exceeds previous
  estimations based on the same final state, J/ψ + J/ψ,
  which are typically of about 2 - 5 mb<cit.>.
  Of course, the results<cit.> also contradict to the most of the measured σ_ eff
  values<cit.>.
  
  A comparison of our predictions with the LHCb experimental data is displayed
  in Figs. <ref> — <ref>. The theoretical uncertainty bands include both scale 
  uncertainties and uncertainties coming from the σ_ eff fitting procedure.
  First of them have been estimated in a usual way,
  by varying the μ_R scale
  around its default value by a factor of 2.
  This was accompanied with using the A0+ and A0- gluon densities instead of  
  default A0 distribution, in accordance with<cit.>.
  As one can see, we achieved a reasonably good
  agreement between the results of our calculations and 
  LHCb measurements, both for √(s) = 7 and 13 TeV.
  There is only exception in the threshold region, m(J/ψ, J/ψ) ≤ 9 GeV,
  where our predictions systematically overshoot the data.
  However, at such low m(J/ψ, J/ψ) 
  an accurate treatment of
  multiple soft gluon emissions, relativistic corrections and other
  nonperturbative effects becomes necessary
  to produce the theoretical estimations. All these issues are 
  out from our present consideration.
  Next, we find that neither the SPS terms, nor the DPS contributions alone
  are able to describe the LHCb data, but only their sum.
  In particular, the DPS contributions are essential to reproduce 
  the measured |Δ y (J/ψ, J/ψ)| distributions 
  at |Δ y (J/ψ, J/ψ)| ≥ 1 or 1.5,
  that confirms the previous expectations<cit.>.
  They are important to 
  describe also the normalization
  of J/ψ rapidity distributions and
  shape of transverse momentum asymmetry A_T at A_T ≤ 0.4, see Figs. <ref> — <ref>.
  
  The presented results, being considered altogether with the ones for inclusive
  single production of charmonia states<cit.>,
  can give a significant impact on the
  understanding of charmonia production 
  within the NRQCD framework and, in particular, on the further understanding of DPS mechanism.
  The most interesting outcome of our study is that
  the extremely low value of DPS effective cross section, σ_ eff∼ 2- 5 mb, obtained 
  in earlier analyses of double J/ψ production at the LHC, is not confirmed.
  
  
  
  § CONCLUSION
   
  
  We have considered the prompt production of J/ψ meson pairs in pp collisions at the LHC
  using the k_T-factorization approach of QCD.
  We employ the fragmentation mechanism to evaluate the color octet 
  contributions to the production cross sections and
  take into account the combinatorial effects of multiple gluon radiation
  in the initial state using the CCFM evolution equation.
  The latter could be essential
  in the kinematical region covered by the CMS and ATLAS measurements.
  On the other hand, we have demonstrated that
  the experimental data taken by the LHCb Collaboration
  at forward rapidities can be described well by the color singlet terms and contributions 
  from the double parton scattering mechanism.
  We determine the DPS effective cross section 
  σ_ eff = 17.5 ± 4.1 mb
  from the combined analysis of the LHCb data collected 
  at √(s) = 7 and 13 TeV.
  The extracted value is compatible with many other estimations 
  based on essentially different final states.
  The extremely low σ_ eff∼ 2 - 5 mb, obtained earlier
  from the double J/ψ production data, is not confirmed.
  
  
  
  § ACKNOWLEDGEMENTS
   
  
  The authors thank G.I. Lykasov for very useful discussions on the topic. 
  M.A.M. and A.A.P. were supported by grants of the foundation for the advancement 
  of theoretical physics and mathematics "Basis" 20-1-3-11-1 and 18-2-6-129-1, respectively.
  
  99
  
  1 CMS Collaboration, JHEP 09, 094 (2014). 
  2 ATLAS Collaboration, Eur. Phys. J. C 77, 76 (2017). 
  3 LHCb Collaboration, Phys. Lett. B 707, 52 (2012). 
  4 LHCb Collaboration, JHEP 06, 047 (2017). 
  5 G. Bodwin, E. Braaten, G. Lepage, Phys. Rev. D 51, 1125 (1995).
  6 P. Cho, A.K. Leibovich, Phys. Rev. D 53, 150 (1996); Phys. Rev. D 53, 6203 (1996).
  7 B. Gong, X.Q. Li, J.-X. Wang, Phys. Lett. B 673, 197 (2009).    
  8 Y.-Q. Ma, K. Wang, K.-T. Chao, Phys. Rev. Lett. 106, 042002 (2011). 
  9 M. Butenschön, B.A. Kniehl, Phys. Rev. Lett. 108, 172002 (2012). 
  10 K.-T. Chao, Y.-Q. Ma, H.-S. Shao, K. Wang, Y.-J. Zhang, Phys. Rev. Lett. 108, 242004 (2012). 
  11 B. Gong, L.-P. Wan, J.-X. Wang, H.-F. Zhang, Phys. Rev. Lett. 110, 042002 (2013). 
  12 Y.-Q. Ma, K. Wang, K.-T. Chao, H.-F. Zhang, Phys. Rev. D 83, 111503 (2011). 
  13 A.K. Likhoded, A.V. Luchinsky, S.V. Poslavsky, Phys. Rev. D 90, 074021 (2014). 
  14 H.-F. Zhang, L. Yu, S.-X. Zhang, L. Jia, Phys. Rev. D 93, 054033 (2016). 
  15 H. Han, Y.-Q. Ma, C. Meng, H.-S. Shao, K.-T. Chao, Phys. Rev. Lett. 114, 092005 (2015). 
  16 H.-F. Zhang, Z. Sun, W.-L. Sang, R. Li, Phys. Rev. Lett. 114, 092006 (2015). 
  17 M. Butenschön, Z. G. He, B.A. Kniehl, Phys. Rev. Lett. 114, 092004 (2015). 
  18 S.P. Baranov, Phys. Rev. D 93, 054037 (2016). 
  19 S.P. Baranov, A.V. Lipatov, Phys. Rev. D 100, 114021 (2019). 
  20 N.A. Abdulov, A.V. Lipatov, Eur. Phys. J. C 79, 830 (2019). 
  21 N.A. Abdulov, A.V. Lipatov, Eur. Phys. J. C 80, 5 (2020). 
  22 J.-P. Lansberg, H.-S. Shao, H.-F. Zhang, Phys. Lett. B 786 342 (2018).
  23 Y. Feng, J. He, J.-P. Lansberg, H.-S. Shao, A. Usachov, H.-F. Zhang, Nucl.Phys. B 945, 114662 (2019).
  24 J.-P. Lansberg, arXiv:1903.09185 [hep-ph].
  25 Z.-G. He and B.A. Kniehl, Phys. Rev. Lett. 115, 022002 (2015).
  26 Y.-J. Li, G.-Z. Xu, K.-Y. Liu, Y.-J. Zhang, JHEP 07, 051 (2013).
  27 L.-P. Sun, H. Han, K.-T. Chao, Phys. Rev. D 94, 074033 (2016).
  28 J.P. Lansberg, H.S. Shao, Phys. Rev. Lett. 111, 122001 (2013).
  29 M. Ciafaloni, Nucl. Phys. B 296, 49 (1988);
  
    S. Catani, F. Fiorani, G. Marchesini, Phys. Lett. B 234, 339 (1990);
  
    S. Catani, F. Fiorani, G. Marchesini, Nucl. Phys. B 336, 18 (1990);
  
    G. Marchesini, Nucl. Phys. B 445, 49 (1995).
  30 M. Diehl, D. Ostermeier, A. Schäfer, JHEP 12, 89 (2012).
  31 P. Bartalini et. al., arXiv:1111.0469 [hep-ph].
  32 H. Abramowicz et. al., arXiv:1306.5413 [hep-ph].
  33 S. Bansal et. al., arXiv:1410.6664 [hep-ph].
  34 R. Astalos et. al., arXiv:1506.05829 [hep-ph].
  35 S.P. Baranov, A.M. Snigirev, N.P. Zotov, Phys. Lett. B 705, 116 (2011).
  36 C.H. Kom, A. Kulesza, W.J. Stirling, Phys. Rev. Lett. 107, 082002 (2011).
  37 S.P. Baranov, A.M. Snigirev, N.P. Zotov, A. Szczurek, W. Schäfer, Phys. Rev. D 87, 034035 (2013).
  38 ATLAS Collaboration, JHEP 11, 110 (2016). 
  39 LHCb Collaboration, JHEP 07, 052 (2016)). 
  40 J.-P. Lansberg, H.-S. Shao, Phys. Lett. B 751, 479 (2015). 
  41 ATLAS Collaboration, Eur. Phys. J. C 77, 76 (2017). 
  42 D0 Collaboration, Phys. Rev. D 90, 111101R (2014). 
  43 D0 Collaboration, Phys. Rev. Lett. 116, 082002 (2016). 
  44 CMS Collaboration, JHEP 05, 013 (2017). 
  45 S. Catani, M. Ciafaloni, F. Hautmann, Nucl. Phys. B 366, 135 (1991);
  
    J.C. Collins, R.K. Ellis, Nucl. Phys. B 360, 3 (1991).
  46 L.V. Gribov, E.M. Levin, M.G. Ryskin, Phys. Rep. 100, 1 (1983);
  
    E.M. Levin, M.G. Ryskin, Yu.M. Shabelsky, A.G. Shuvaev, Sov. J. Nucl. Phys. 53, 657 (1991).
  47 R. Angeles-Martinez et al., Acta Phys. Polon. B 46, 2501 (2015).
  48 H. Jung, S.P. Baranov, M. Deak, A. Grebenyuk, F. Hautmann, M. Hentschinski, A. Knutsson,  
    M. Kramer, K. Kutak, A.V. Lipatov, N.P. Zotov, Eur. Phys. J. C 70, 1237 (2010).
  49 C.-H. Chang, Nucl. Phys. B 172, 425 (1980);
  
    E.L. Berger, D.L. Jones, Phys. Rev. D 23, 1521 (1981);
  
    R. Baier, R. Rückl, Phys. Lett. B 102, 364 (1981);
  
    S.S. Gershtein, A.K. Likhoded, S.R. Slabospitsky, Sov. J. Nucl. Phys. 34, 128 (1981).
  50 S.P. Baranov, Phys. Rev. D 84, 054012 (2011).
  51 Y.-Q. Ma, J.-W. Qiu, H. Zhang, Phys. Rev. D 89, 094029 (2014).
  52 E. Bycling, K. Kajantie, Particle Kinematics, John Wiley and Sons (1973).
  53 J. Alwall, A. Ballestrero, P. Bartalini, S. Belov, E. Boos, A. Buckley, J.M. Butterworth, 
    L. Dudko, S. Frixione, L. Garren, S. Gieseke, A. Gusev, I. Hinchliffe, J. Huston, B. Kersevan, 
    F. Krauss, N. Lavesson, L. Lönnblad, E. Maina, F. Maltoni, M.L. Mangano, F. Moortgat, S. Mrenna, 
    C.G. Papadopoulos, R. Pittau, P. Richardson, M.H. Seymour, A. Sherstnev, T. Sjöstrand, P. Skands, 
    S.R. Slabospitsky, Z. Wcas, B.R. Webber, M. Worek, D. Zeppenfeld, Comput. Phys. Commun. 176, 300 (2007).
  54 A.V. Lipatov, M.A. Malyshev, H. Jung, Phys. Rev. D 100, 034028 (2019).
  55 A.V. Lipatov, S.P. Baranov, M.A. Malyshev, Eur. Phys. J. C 80, 330 (2020).
  56 V.L. Korotkikh, A.M. Snigirev, Phys. Lett. B 594, 171 (2004).
  57 A.M. Snigirev, Phys. Rev. D 83, 034028 (2011).
  58 S.P. Baranov, A.V. Lipatov, M.A. Malyshev, A.M. Snigirev, N.P. Zotov, Phys. Rev. D 93, 094013 (2016).
  59 H. Jung, arXiv:hep-ph/0411287.
  60 http://tmd.hepforge.org
  61 PDG Collaboration, Phys. Rev. D 98, 030001 (2018).
  62 L.A. Harland-Lang, A.D. Martin, P. Motylinski, R.S. Thorne, Eur. Phys. J. C 75, 204 (2015).
  63 S.P. Baranov, A.H. Rezaeian, Phys. Rev. D 93, 114011 (2016). 
  64 
  65 CMS Collaboration, JHEP 03, 32 (2014). 
  66 ATLAS Collaboration, New J. Phys. 15, 033038 (2013). 
  67 D0 Collaboration, Phys. Rev. D 93, 052008 (2016). 
  68 D0 Collaboration, Phys. Rev. D 89, 072006 (2014). 
  69 LHCb Collaboration, JHEP 06, 141 (2012). 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  