Double heavy quarkonium production in diffractive processes at the Run 2 LHC energy

The double quarkonium production in single and double diffractive processes is investigated considering $pp$ collisions at the Run 2 LHC energy. Using the nonrelativistic QCD (NRQCD) factorization formalism for the quarkonium production and the Resolved Pomeron model to describe the diffractive processes, we estimate the rapidity and transverse momentum dependencies of the cross sections for the $J/\Psi J/\Psi$ and $\Upsilon \Upsilon$ production. The contributions of the color-singlet and color-octet channels are estimated and predictions for the total cross sections in the kinematical regions of the LHC experiments are also presented. Our results demonstrate that the double quarkonium production in diffractive processes is not negligible and that its study can be useful to test the underlying assumptions present in the description of the single and double diffractive processes.


I. INTRODUCTION
The study of hadronic collisions at the LHC provides a unique environment for precise measurements of poorly understood phenomena. In particular, the study of quarkonium production at the LHC is expected to provide important insight for improving the theoretical description of its mechanism of production, which represents one of the long-standing problems of Quantum Chromodynamics (QCD) [1]. During the last years, the study of the heavy quarkonium states in different reactions was a theme of intense interest, with distinct approaches proposed to describe the factorization of the small-and large-distance contributions present in the quarkonium production [2]. One of the possible frameworks is the Non-Relativistic QCD (NRQCD) formalism [3], which is a non-relativistic effective theory equivalent to QCD, which considers that the heavy quarkonium cross section receives contributions of the color singlet and color octet mechanisms. In this formalism, the cross section for the production of a heavy quarkonium state H factorizes as σ(ab → H + X) = n σ(ab → QQ[n] + X) O H [n] , where the coefficients σ(ab → QQ[n] + X) are perturbatively calculated short distance cross sections for the production of the heavy quark pair QQ in an intermediate Fock state n, which does not have to be color neutral. The O H [n] are nonperturbative long distance matrix elements, which describe the transition of the intermediate QQ in the physical state H via soft gluon radiation. Such quantities are determined from experimental data, with their values being, unfortunately, strongly dependent of the fit procedures [4][5][6]. Although the NRQCD formalism is able to describe a large set of experimental data, recent results indicate that the situation is far from being conclusive. Consequently, further tests for the color singlet and color octet mechanisms in NRQCD are still needed to clarify several aspects present in the heavy quarkonium production.
In the last years the production of quarkonium pairs became an alternative for the study of the underlying production mechanism (See e.g. Refs. [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21] ). The recent theoretical and experimental advances were strongly motivated by the large contribution of multiple parton interactions at the LHC energies [22]. With the increasing of the energy, the probability that one proton -proton collision leads to more than one scattering process also increases. As a consequence, a double quarkonium H 1 H 2 final state can be generated in a single partonic scattering (SPS) process, with σ(pp → H 1 H 2 X) ∝ σ(gg → H 1 H 2 ), as well as in a double parton scattering (DPS) process, where σ(pp → H 1 H 2 X) ∝ σ(gg → H 1 ) × σ(gg → H 2 ) [23][24][25][26][27][28], and in higher order multiple scattering processes [29,30]. Recent results for the double J/Ψ production indicate that the double scattering contribution is comparable to the single one. In contrast, the double Υ production is expected to be dominated by the SPS mechanism (See e.g. [31]). Considering the values predicted for the cross sections and the large luminosity present in the Run 2 of the LHC, a detailed experimental analysis of the double quarkonium production is expected to be feasible, which will allow to improve our understanding of the heavy quarkonium production.
The recent studies of the double heavy quarkonium production in single parton scattering processes are in general dedicated to the calculation of the production of this final state in inclusive reactions, where both initial protons dissociate in the interaction. A typical diagram is represented in Fig. 1(left panel). However, a double quarkonium can also be produced in diffractive interactions, where one (or both) of the protons remain intact and empty regions in pseudo-rapidity, called rapidity gaps, separate the intact very forward proton(s) from the H 1 H 2 state. Examples are the single and double diffractive processes, represented in the central and right panels of Fig. 1, respectively. Our motivation to study these processes is twofold. Firstly, the treatment of diffractive processes at the LHC have recently attracted much attention as a way of investigate the interplay of small-and large-distance dynamics within Quantum Chromodynamics [32]. Our goal is to complement previous studies [33,34] considering other final states and present a comprehensive analysis of the diffractive J/ΨJ/Ψ and ΥΥ production in pp collisions at the Run 2 LHC energy ( √ s = 13 TeV). Secondly, the LHCb Collaboration have measured during the Run 1 the double charmonium production [35] in exclusive reactions, which are characterized by two rapidity gaps and nothing else is produced except the leading protons and the central charmonium states. New data, with larger statistics and smaller uncertainties, are expected to be released soon. However, in order to have access to these exclusive processes, it is fundamental to have control of the background associated to the production of this same final state in double diffractive processes. The theoretical treatment of the exclusive double charmonium production was presented in Ref. [33]. In contrast, in this paper we will focus on the treatment of quarkonium pair production in double diffractive processes, which also exhibit two rapidity gaps in the final state. However, differently from the exclusive case, they contain soft particles accompanying the production of a hard diffractive object (double quarkonium), with the rapidity gaps becoming, in general, smaller than in the exclusive case. The presence of the soft particles should increase the number of tracks in the detector. Therefore, in principle, the double diffractive and exclusive contributions could be separated using the exclusivity criteria. If this separation is feasible, is also allows to study in more detail the modelling of the diffractive processes and the description of the quarkonium production. On the other hand, due to the large pile-up of events in each bunching crossing present in the Run 2, it is not clear if the separation of the diffractive and exclusive events will be possible by measuring the rapidity gaps and counting the number of tracks in the final state. In principle, the only possibility to detect events characterized by two rapidity gaps, as those associated to the double diffractive and exclusive double quarkonium production, is by tagging the intact hadrons in the final state. It implies that the key element to measure diffractive events at the LHC will be tagging the forward scattered incoming hadrons [36]. Moreover, if only one proton is tagged in the final state, the contribution of the double quarkonium production in single diffractive processes [ Fig. 1 (central panel)], should also be taken into account. Predictions for the single diffractive production of a quarkonium pair are presented here for the first time. Finally, it is important to emphasize that the DPS contribution for diffractive processes will be strongly suppressed, since the probability of two simultaneous interactions without the fragmentation of one (or both) of the incident proton is very small. The typical topology of the diffractive events (leading protons and rapidity gaps) is in general destroyed by multiple interactions. All these aspects motivate the analysis that will be performed in what follows. This paper is organized as follows. In the next Section we will review the formalism used to treat the double quarkonium production in single and double diffractive processes. The main aspects of the NRQCD formalism will be reviewed as well the basic assumptions of the Resolved Pomeron Model will be presented. The treatment of the gap survival probability will be also discussed. In Section III we will present our predictions for the rapidity and transverse momentum distributions for the single and double diffractive J/ΨJ/Ψ and ΥΥ production in pp collisions at the Run FIG. 2: Typical diagrams that contribute at leading order for the double quarkonium production by the gluon fusion subprocess.

II. FORMALISM
In what follows we will present a brief review of the formalism needed to investigate the diffractive production of two spin-triplet S-wave quarkonia (J/Ψ and Υ) in pp collisions at the Run 2 LHC energy. Such final states provide a clean experimental signature when studied in their decay into muons. In order to include the color-singlet and color-octet channels for each quarkonium final state, we will assume the NRQCD formalism [3], which takes into account of both contributions in a systematic way. In this formalism the cross section can be factorized as the product of short-distance and long-distance parts. The short-distance coefficients are process dependent, but can be calculated using perturbative QCD. In contrast, the long-distance matrix elements (LDMEs) are expected to be process independent, but are associated to nonperturbative physics. One important characteristic of the NRQCD formalism, is that a heavy quark pair, created in a color-octet state in the perturbative part is also allowed to contribute for the quarkonium production, with the color being neutralized in the long-distance regime. The longdistance contribution can be expanded in terms of the relative velocity v of the heavy quarks in the rest frame of the heavy quarkonium. Moreover, in the NRQCD formalism, the Fock states are expressed in terms of the dynamical gluons and the heavy quark pair, which implies that the wavefunction of a heavy quarkonium H can be expressed as follows: As demonstrated in Refs. [12,13], the dominant double quarkonium Fock states are the following combinations: where the suppression present in the last term can be compensated by the enhancement of the large propagators present in the color octet channel for the production of a quarkonium pair with large transverse momentum. At the LHC energies, the double quarkonium production in inclusive and difffractive processes is dominated by the gluon fusion subprocess gg → H 1 H 2 . Some of the typical diagrams that contribute at leading order [O(α 4 s )] for the cross section are presented in Fig. 2. The first two diagrams are examples of the gg → QQ 1 ( 3 S 1 ) + QQ 1 ( 3 S 1 ) subprocess, where the heavy quark pairs are produced in a color-singlet state. In total, there are 31 diagrams that contribute for the color singlet channel. On the other hand, the color-octet channel, associated to the gg → QQ 8 ( 3 S 1 ) + QQ 8 ( 3 S 1 ) subprocess, contributes with 72 diagrams, with two of them being represented by the last two diagrams in Fig. 2.
In order to calculate the observables, the cross sections for the subprocesses, which are associated with the shortdistance physics, should be multiplied by the large-distance matrix elements O H 1 ( 3 S 1 ) and O H 8 ( 3 S 1 ) . These matrix elements account for the transition probability of the heavy quark pair QQ n ( 2s+1 L J ) in a given color state n, spin s, angular momentum L and total angular momentum J, to hadronize into quarkonium H. The color-singlet matrix elements O H 1 ( 3 S 1 ) are usually determined from the leptonic decay of H. On the other hand, the color-octet one, O H 8 ( 3 S 1 ) , are not known quite accurately and has been obtained by fitting the transverse momentum spectrum for quarkonium production in a global analysis of different sets of experimental data. As a consequence, the predictions for the color-octet contribution are, in general, strongly dependent on the underlying assumptions for the associated matrix elements. Such uncertainty has direct impact in the predictions of the observables that are dominated by the color-octet channel. As we will show below, that is not the case for the double quarkonium production at small transverse momentum, which determines the magnitude of the total cross section.
The leading order color-singlet and color-octet contributions for the differential cross sections of the gg → H 1 H 2 subprocess were calculated by several authors in the nonrelativistic approximation [7][8][9][10][11][12][13][14] and more recently taking into account the relativistic effects [15,17] and higher-order corrections [16,19,20]. As the impact of the relativistic effects is still a theme of debate [15,17] and the higher -order corrections are predicted to modify the transverse momentum distribution at large -p T [16,19,20], in our calculations for the diffractive production we will estimate the differential cross sections at leading order, disregarding these corrections. We will follow the notation from Refs. [10,13], where the expressions for the differential cross sections are explicitly presented. In the case of the color-singlet contribution, the differential cross section, already multiplied by the corresponding matrix element, is given by [10] dσ dt where |R H (0)| 2 is the square of the radial wave function of the quarkonium H at the origin. The variablesŝ,t andû are the usual Mandelstam variables for the partonic subprocess, and in the nonrelativistic approximation we assume that M = 2m Q . For the charm and bottom quark masses we will use m c = 1.5 GeV and m b = 4.7 GeV. The detailed expressions for the a jkl coefficients in Eq. (3) are given in Ref. [10]. Finally, the quarkonium radial function at the origin is related to the leptonic decay rate [37] as Γ(H → e + e − ) = . From recent PDG data [38] for Γ(J/Ψ → e + e − ) and Γ(Υ → e + e − ), we obtain |R J/Ψ (0)| 2 = 0.53 GeV 3 and |R Υ (0)| 2 = 4.6 GeV 3 for the J/Ψ and the Υ, respectively. For the color-octet contribution, the differential cross section, already multiplied by the corresponding matrix element, is given by [13] dσ dt where the a j coefficients can be found in Ref. [13]. As in Ref. [13] we will assume that O [40]. It is important to emphasize that more recent global analysis of the world data for charmonium production [4] indicate that smaller values for O J/ψ 8 ( 3 S 1 ) are prefered by the data. As consequence, our predictions for the associated color-octet contribution should be considered as an upper bound.
In order to estimate the double quarkonium production in single and diffractive interactions we will assume the validity of the factorization theorem for these processes. Consequently, the hadronic cross sections will be given by the convolution of the above differential cross sections with the inclusive and/or diffractive gluon distribution functions of the incident particles for the correspondent process. In the particular case of the single diffractive double quarkonium production, represented in Fig. 1 (central panel), the cross section can be expressed as follows with the final state being characterized by one intact hadron and one rapidity gap. In Eq. (5) we take into account that any of the two incident protons can remain intact in the interaction. Moreover, g p and g D are the inclusive and diffractive gluon distributions, probed at the scale µ 2 , which we assume to be equal to µ 2 = M 2 + p 2 T . On the other hand, the cross section for the double diffractive process will be given by This process is represented in Fig. 1 (right panel) and the final state will be characterized by two intact hadrons and two rapidity gaps. In our calculations we will use the CTEQ6L parametrization [41] for the inclusive gluon distribution. For the diffractive gluon distribution, it will be modelled using the Resolved Pomeron Model [42]. The basic idea is that the hard scattering resolves the quark and gluon content of the Pomeron [42] and it can be obtained analysing the experimental data from diffractive deep inelastic scattering (DDIS) at HERA, providing us with the diffractive distributions of singlet quarks and gluons in the Pomeron [43]. Consequently, the diffractive parton distributions are expressed in terms of parton distributions in the Pomeron and a Regge parametrization of the flux factor describing the Pomeron emission by the proton. The parton distributions have evolution given by the DGLAP evolution equations and should be determined from events with a rapidity gap or a intact hadron. The diffractive gluon distribution, g D p (x, Q 2 ), is defined as a convolution of the Pomeron flux emitted by the proton, f p I P (x I P ), and the gluon distribution in the Pomeron, g I P (β, Q 2 ), where β is the momentum fraction carried by the partons inside the Pomeron. The diffractive gluon distribution of the proton is then given by The Pomeron flux is given by   where t min , t max are kinematic boundaries. The Pomeron flux factor is motivated by Regge theory, where the Pomeron trajectory is assumed to be linear, α I P (t) = α I P (0) + α ′ I P t, and the parameters B I P , α ′ I P and their uncertainties are obtained from fits to H1 data [43]. The slope of the Pomeron flux is B I P = 5.5 −2.0 +0.7 GeV −2 , the Regge trajectory of the Pomeron is α P (t) = α P (0) + α ′ P t with α P (0) = 1.111 ± 0.007 and α ′ P = 0.06 +0.19 −0.06 GeV −2 . The t integration boundaries are t max = −m 2 p x 2 I P /(1−x I P ) (m p denotes the proton mass) and t min = −1 GeV 2 . Finally, the normalization factor A P = 1.7101 is chosen such that x I P × tmax tmin dt f I P /p (x I P , t) = 1 at x I P = 0.003. In our analysis we use the fit B obtained by the H1 Collaboration at DESY-HERA for the diffractive gluon distribution [43]. Finally, it is important to emphasize that the Resolved Pomeron model implies the presence of additional particles in the final state, associated to the remnants of the Pomeron, which dissociates in the interaction. As discussed in the Introduction, the presence of these particles can be used, in principle, to discriminate the diffractive from the exclusive double quarkonium production.
One important open question in the treatment of diffractive interactions in hadronic collisions is whether the cross sections for the associated single and double diffractive processes are somewhat modified by soft interactions which lead to an extra production of particles that destroy the rapidity gaps in the final state [44]. In the case of diffractive interactions in pp/pp collisions, the experimental results obtained at TEVATRON [45] and LHC [46,47] have demonstrated that one should take into account these additional absorption effects that imply the violation of the QCD hard scattering factorization theorem for diffraction [48]. One has that the soft rescattering corrections associated to reinteractions (often referred to as multiple scatterings) between spectator partons of the colliding hadrons produce additional final-state particles which fill the would-be rapidity gap and suppress the diffractive events. Consequently, in order to estimate the diffractive cross sections in hadronic collisions we need to take into account the probability that such emission does not occur. The fact that the diffractive factorization breaking is intimately related to soft multiple scattering in hadron-hadron collisions has motivated the modelling of these effects using a general purpose Monte Carlo [49][50][51]. A current shortcoming of these promissing approaches is that, due to the complexity of the diffractive interactions, their predictions are still strongly dependent on the treatment of the multiple interactions, the assumptions for the color flow along the rapidity gap as well as the modelling of possible proton excitations. Another possible approach to treat this problem is based on the assumption that the hard process occurs on a short enough timescale such that the physics that generate the additional particles can be factorized and accounted by an overall factor, denoted gap survival factor |S| 2 , multiplying the cross section calculated using the collinear factorization and the diffractive parton distributions extracted from HERA data. The modelling, magnitude and universality of this factor are still a theme of intense debate [52][53][54]. In general the values of |S| 2 depend on the energy, being typically of order 1 -5 % for LHC energies. Such approach have been largely used in the literature to estimate the hard diffractive processes at the LHC (See e.g. Refs. [55][56][57][58][59][60][61][62][63]) with reasonable success to describe the current data. However, as the effects that determine the gap survival factor have nonperturbative nature, they are difficult to treat and its magnitude is strongly model dependent (For recent reviews see Refs. [52,53]). In what follows we also follow this simplified approach, assuming |S| 2 = 0.02 for the double diffractive production and |S| 2 = 0.05 for the single diffractive one. Such values were estimated in Ref. [64] using a two-channel eikonal formalism that take into account the contributions of high-mass diffractive dissociation, possible nucleon excitations in the diffractive interaction as well as the contribution of pion loops for the bare Pomeron pole. It is important to emphasize that this choice is somewhat arbitrary, and mainly motivated by the possibility to compare our predictions with those obtained in other analysis. Recent studies from the CMS Collaboration [47] indicate that this factor can be larger than this value by a factor ≈ 4. Consequently, our results can be considered a lower bound for the diffractive contribution. However, it is important to emphasize that the uncertainty on |S| 2 only affect the normalization of the cross sections, with the shape of the distributions being a direct probe of the underlying assumption that the soft rescattering effects can be factorized of the hard process. In particular, if a different value for |S| 2 is constrained by the experimental data for e.g. diffractive heavy quark production, our predictions can be directly rescaled and compared with the diffractive double quarkonium data.

III. RESULTS
In what follows we will present our predictions for the single and double diffractive production of J/ΨJ/Ψ and ΥΥ in pp collisions at the Run 2 LHC energy ( √ s = 13 TeV). We will estimate the rapidity (y) and transverse momentum (p T ) distributions for these processes and a comparison with the inclusive predictions will also be presented. In the case of double diffractive processes, these distributions can be estimated directly from the differential cross section given by where x 1 min =x T e y 2−xT e −y , x 2 = x1xT e −y 2x1−xT e y ,x T = 2mT √ s and m T = M 2 + p 2 T . The Eq. (9) can be directly generalized for single diffractive and inclusive processes by the adequate replacing of the diffractive gluon distribution by the inclusive one.
In Fig. 3 we present our predictions for the rapidity distributions for the double J/Ψ (left panel) and double Υ (right panel) production in single and double diffractive interactions. The inclusive predictions are also presented for comparison. In our calculations we have included the color-singlet and color-octet mechanisms, which are denoted, respectively, by CSM and COM hereafter. The rapidity distributions are flat at the central rapidity region (y ≈ 0), with the inclusive prediction being a factor ≈ 30 (10 3 ) larger than SD (DD) one for the double J/Ψ production. On the other hand, for the double Υ production, we predict that the inclusive prediction is a factor ≈ 10 2 (10 4 ) larger than SD (DD) result. These differences are also observed in the predictions for the total cross sections, presented in Table  I for two rapidity ranges. In particular, we present our predictions for the kinematical range probed by the LHCb Collaboration. In this case, the single and double diffractive predictions are reduced by approximately one order of magnitude in comparison to the full rapidity range. Our predictions for the double J/Ψ production in the LHCb range for DD interactions are similar to those obtained in Ref. [33] for the exclusive production. As already emphasized in the Introduction, the topology of the final state of these two processes is different. While in central exclusive processes one has only the leading hadrons, two quarkonia and nothing else, in the double diffractive production one expect to have some extra particles coming from the Pomeron remnants. Consequently, in principle, the experimental separation between these two processes can be performed at smaller luminosities, as those presented in the LHCb detector. If it is feasible, a future analysis of the diffractive events can be useful to constrain the underlying assumptions present in description of these processes. In particular, it will allow to improve the treatment of the absorptive effects that lead to the breakdown of factorization in diffractive processes at hadronic colliders. On the other hand, if the separation of the exclusive and diffractive events is not possible due to large pile-up, our results indicate that the background associated to diffractive interactions cannot be disregarded in the selection of the exclusive events.
In Fig. 4 we present our predictions for the transverse momentum distributions for the double J/Ψ (left panels) and double Υ (right panels) production in single diffractive processes considering the full rapidity range (upper panels) and the LHCb range (lower panels). The contributions of the color-singlet and color-octet mechanisms are presented separately, as well as the sum of both for the single diffractive and inclusive production. The shape of the p T -distributions for single diffractive and inclusive processes are very similar, with the distribution vanishing for p T → 0, in agreement with the results obtained in Ref. [13]. For the full rapidity range (upper panels), one have that the contribution of the color-singlet mechanism is dominant at small values of p T . It implies that the magnitude of the total cross sections is determined by this contribution and are not affected by the current uncertainty in the determination of the color-octet matrix elements. On the other hand, our results indicate that the color-octet mechanism determines the behaviour of the distribution for p T ≥ 10 (20) GeV for the double J/Ψ (Υ) production. We obtain similar results for the LHCb range (lower panels), with the main difference being the normalization. The predictions for double diffractive processes are presented in Fig. 5. As in the SD case, the behaviour of the distribution at small -p T is determined by the color-singlet contribution and the color-octet one only contributes at large transverse momentum. Such result is expected, since the single and double diffractive processes are determined by the same differential cross section for the subprocesses.
As discussed before, in our calculations we are assuming that the contribution of the absorptive effects can be factorized and taken into account by the multiplicative factor |S| 2 , which is assumed to be a constant that is independent on the final state produced in the diffractive interaction and the kinematical range considered. Such strong assumption can be tested by a future analysis of the ratio between the cross sections for the double J/Ψ and double Υ production in single and double diffractive processes. If the assumption is correct, the transverse momentum and rapidity behaviours of the ratio will be independent of |S| 2 . Additionally, the impact of the next-to-leading order corrections and the dependence on the modelling of the inclusive and diffractive gluon distributions are expected to cancel in the ratio. Moreover, considering that the cross sections are dominated by the color-singlet contribution, which is reasonably well known, and that the inclusive and diffractive gluon distributions are also well determined in the kinematical range of interest (x ≈ 10 −3 ), a future analysis of the ratio can also be considered a direct probe of the framework considered to describe the double quarkonium production. In Fig. 6 we present our predictions for the transverse momentum (left panel) and rapidity (right panel) dependencies of the ratio. The predictions for the single and double diffractive production are similar. At large -p T (≥ 10 GeV) we predict that the ratio will be ≈ 1. Moreover, our results indicate that the ratio is almost rapidity independent. Similar results are obtained considering the LHCb kinematical range. If a different behaviour is observed in a future experimental analysis, it will give us some hint about the treatment and kinematical dependence of the soft gluon interactions in diffractive interactions. Some comments are in order before to summarize our main results in the next Section. In our calculations we have considered that the renormalization and factorization scales are equal to the transverse mass. As demonstrated in previous studies of the inclusive double quarkonium production [16,18,19], the predictions are strongly sensitive to this choice, since the cross section at leading order is proportional to α 4 s . As the dependence on the partonic cross sections is the same for inclusive and diffractive processes, we expect that a similar behaviour is also present in our predictions. In particular, a variation of the renormalization and factorization scales up and down by a factor 2 with respect to the central value used in our calculations is expected to modify our predictions by ≈ 70% [16]. Such uncertainty can be reduced in the future by precise measurements of the inclusive production. Another input parameter entering our calculations is the nonperturbative wave function of the quarkonium at the origin R H (0). As our predictions are dominated by the color singlet contributions, we have that if a different value for R H (0) is chosen, the new predictions can be obtained from the results presented here by multiplying our predictions by a factor of where |R H (0)| 2 new is the new value for the square of the radial wave function at the origin. Finally, in our analysis we have considered the CTEQ6L parametrization for the inclusive gluon distribution and the fit B of the H1 parametrization for the diffractive gluon distribution. As the main contribution for the cross sections come from values of x ≈ 10 −3 , where the inclusive gluon distribution is well determined, with different pametrizations predicting similar values, the impact of use a different model is small. In the diffractive case, we have verified that our predictions are modified by 9% if the fit A is used as input in the calculations. All these aspects imply an uncertainty of order of a factor 2 in our predictions, which can be strongly reduced when the ratio of cross sections is considered.

IV. SUMMARY
The experimental analysis of the double quarkonium production became a reality during the last years. These results have attracted a renewed attention to the theoretical description of this process, which is expected to provide important insights that will allow us to improve our understanding of the hadron structure and the production mechanism. Moreover, the recent installation of forward detectors can be considered the begin of a new era in the study of Diffractive Physics, whose description is still one of the great challenges of the QCD. Consequently, the study of the double heavy quarkonium in diffractive processes has the potentiality of provide results that can be used to improve the understanding of several aspects in the theory of strong interactions. That was the main motivation for the study performed in this paper, where we have presented a comprehensive analysis of the double heavy quarkonium production in single and double diffractive processes. We have provided predictions for the transverse momentum and rapidity distributions and for the total cross sections, derived using the NRQCD formalism for the quarkonium production and the Pomeron Resolved Model for the description of the diffractive interaction. We have demonstrated that the cross sections are dominated by the color-singlet mechanism, with the color-octet one being important only at large -p T , where the magnitude of the cross sections is strongly reduced. Moreover, our results indicate that the contribution of the diffractive production is non-negligible at the Run 2 LHC energy, which implies that a future experimental analysis is, in principle, feasible. We also have indicated that the analysis of the ratio between cross sections can be useful to probe the treatment of the absorptive corrections. Finally, our results indicate that a future analysis of the diffractive events can be useful to constrain the underlying assumptions present in the description of the double quarkonium production.