Double-quarkonium production at a fixed-target experiment at the LHC (AFTER@LHC)

We present predictions for double-quarkonium production in the kinematical region relevant for the proposed fixed-target experiment using the LHC beams (dubbed as AFTER@LHC). These include all spin-triplet S -wave charmonium and bottomonium pairs, i.e. Psi(n_1S) + Psi(n_2S), Psi(n_1S) + Upsilon(m_1S) and Upsilon(m_1S) + Upsilon(m_2S ) with n_1,n_2 = 1,2 and m_1,m_2 = 1,2,3. We calculate the contributions from double-parton scatterings and single-parton scatterings. With an integrated luminosity of 20 fb-1 to be collected at AFTER@LHC, we find that the yields for double-charmonium production are large enough for differential distribution measurements. We discuss some differential distributions for J/Psi + J/Psi production, which can help to study the physics of double-parton and single-parton scatterings in a new energy range and which might also be sensitive to double intrinsic c-bar(c) coalescence at large negative Feynman x.


Introduction
Heavy-quarkonium production is typically a multi-scale process, which involves both shortand long-distance aspects of the strong interaction. This particularity makes heavy-quarkonium production an ideal probe to study Quantum Chromodynamics (QCD) in its perturbative and nonperturbative regimes simultaneously. Studies have extensively been performed at collider and fixed-target energies in proton-proton, proton-nucleus and nucleus-nucleus collisions (see reviews e.g. Refs. [1,2,3]). The associated production of heavy quarkonium is a very interesting process not only because it provides a way to pin down the heavy-quarkonium production mechanism but also because it can help to understand a new dynamics of hadron collisions appearing at high energies, where multiple-parton scatterings (MPS) happen simultaneously, among which the most likely is of course two short-distance interactions from a single hadron-hadron collision -doubleparton scattering (DPS). The relevant DPS analyses with heavy quarkonium are J/ψ + W [4], J/ψ + Z [5], J/ψ+charm [6] and J/ψ + J/ψ [7] production.
In particular, double-quarkonium production is of specific interest. It provides an original tool to study quarkonium production from the conventional single-parton scatterings (SPSs). The relevant earlier studies in the literature can be found in Refs. [8,9,10,11,12,13,14,15,16,17,18,19]. Moreover, it was claimed in Refs. [20,21,22,23,24,25,18,19] that DPS contributions should be significant source of J/ψ + J/ψ, especially at high energies where there is a high parton flux. In addition, the spin-triplet S -waves (e.g. J/ψ, Υ) provide clean signatures with their small background when they are studied in their decay into muon pairs. Hence, they are easy to trigger on experimentally, in contrast to hadronic jets and open charm mesons productions, which require either good calorimetry or good particle identification.
A first comprehensive comparison between experiments [26,7,27] and theory for J/ψ-pair production at the Tevatron and the LHC was performed in Ref. [18], where we pointed out that this observable could be used to probe different mechanisms in different kinematical regions. We noted that the direct DPS measurement by D0 collaboration [7] -looking at the rapidity-difference spectrum-is consistent with the J/ψ-pair measurement by the CMS collaboration [27] and, as we will discuss later on, compatible with rather larger DPS rates. On the other hand, as we advocated in [16], one cannot draw a definite conclusion on the presence of DPS in the LHCb data [26] with relatively low statistics.
In this context, we found it important to study the potentialities offered by the use of the 7 TeV proton LHC beams used in the fixed target mode to study quarkonium pair production. Its multi-TeV beams indeed allow one to study p + p, p + d and p + A collisions at a centre-of-mass energy √ s NN 115 GeV as well as Pb + p and Pb + A collisions at √ s NN 72 GeV, with the high precision typical of the fixed-target mode. It was indeed advocated in [3,28] that such a facility, referred to as AFTER@LHC, would become a quarkonium, prompt photon and heavyflavour observatory thanks to its large expected luminosity (for recent phenomenological studies, see [29,30,31,32,33,34,35,36,37,38]). A first feasibility study for quarkonium production was presented in [39] and demonstrated that a LHCb-like detector would perform extremely well in the fixed-target mode. Similar performances are expected for quarkonium-pair production. Integrated luminosities as large as 20 fb −1 [3] can be delivered during a one-year run of p + H collisions with a bent crystal to extract the beam [40]. The LHC beam can also go through an internal-gas-target system 1 . Conservatively sticking to gas pressures already reachable now, yearly integrated luminosities reach 100 pb −1 . With a designed target cell similar to that of HER-MES [44], a few fb −1 yr −1 are probably also reachable. We have reported in Tab. 1 the instantaneous and yearly integrated luminosities expected with the proton beams on various target species of various thicknesses, for both options.
The structure of this paper is as follows. In section 2, we detail and justify our methodology to compute both DPS and SPS contributions to quarkonium-pair production. Section 3 contains a general discussion of the interest to look at DPS vs SPS contributions at different energies. This prepares the discussion of our results at √ s = 115 GeV relevant for AFTER@LHC in Section 4. Section 5 gathers our conclusions.

Target
Thickness perfect gas 100 10 −9 10 100 Table 1: Expected luminosities obtained for a 7 TeV proton beam extracted by means of a bent crystal or obtained with an internal gas target with a pressure similar to that of SMOG@LHCb [42].

Methodology
In this section, we will explain the main ingredients used to compute the rates of doublequarkonium production at AFTER@LHC, which closely follows from our previous work in Ref. [18].

Double-parton scatterings
The description of such a mechanism is usually done by assuming that DPS can be factorised into two single-parton scatterings (SPS); this can be seen as a first rough approximation which can however be justified by the fact that possible unfactorisable corrections due to parton correlations should be small at small x. For the double-quarkonium production, the general formalism with a factorisation assumption is (see e.g. Ref. [24]) where δ Q 1 Q 2 is the Kronecker delta function, Γ i j (x 1 , x 2 , b 1 , b 2 ) is the generalised double distributions with the longitudinal fractions x 1 ,x 2 and the transverse impact parameters b 1 and b 2 . A further factorisation assumption is to decompose Γ i j (x 1 , x 2 , b 1 , b 2 ) into a longitudinal part and a transverse part where D i j (x 1 , x 2 ) is the double-parton distribution functions (dPDF) [45]. Moreover, by ignoring the correlations between partons produced from each hadrons, one can use where f i (x 1 ) and f j (x 2 ) are the normal single PDFs. This yields to If one further ignores the parton flavour dependence in T i, j,k,l (b) and defines the overlapping function one reaches the so-called "pocket formula" where σ Q 1 and σ Q 1 are the cross sections for respectively single Q 1 and Q 2 production and σ eff is a parameter to characterise an effective spatial area of the parton-parton interactions via Hence, it is only related to the initial states and should be independent of the final state. However, the validation of its universality (process independence as well as energy independence) and the hold of the factorisation in Eq.(6) should be cross checked case by case. Thanks to its larger luminosity and probably wide rapidity coverage, AFTER@LHC provides a unique opportunity to probe DPS and to extract σ eff from double-quarkonium final states.
To perform our predictions, we will use σ eff = 5.0 ± 2.75 mb, which was determined from J/ψ-pair production data at the Tevatron by D0 collaboration [7]. The reason for such a choice is that all of the double-quarkonium production processes share the same gluon-gluon initial states and the typical x are not that much different. It guarantees that we only need to assume the energy independent of σ eff . However, we do not claim that this value is the only one possible; we only take it as our reference number. If one wants to use another value of σ eff , one can just simply perform a rescaling of our presented numbers.  Table 2: χ 2 results after a combined fit of d 2 σ/dP T dy to (a) the ψ(nS ) PHENIX data [46] by fixing n = 2 and P T = 4.5 GeV and (b) the Υ(nS ) data CDF [47], ATLAS [48], CMS [49] and LHCb [50,51] data by fixing n = 2 and P T = 13.5 GeV.
Due to the fact that the description of single heavy quarkonium production at hadron colliders in the whole kinematical region is still a challenge to theorists, using ab initio theoretical computation of σ Q would significantly inflate theoretical uncertainties. Instead, we will work in a data-driven way to determine σ Q .
Our procedure is as follows. We start from the cross section σ Q i which can be written as  Figure 1: Comparison with PHENIX measurement [46] for J/ψ (a,b) and ψ(2S ) (c) production and with LHCb measurement [50] for where f a , f b are the parton distribution functions (PDF) of the initial partons a and b, dLIPS Q+X is the Lorentz-invariant phase space measure for pp → Q + X andŝ is the partonic centre-of-mass energy (i.e.ŝ = x 1 x 2 s). For single quarkonium production in p + p collisions at √ s = 115 GeV, the gluon-gluon initial state is dominant. The initial colour and helicity averaged amplitude square for gg → Q + X can be expressed in the form of a crystal ball function [20] |A where K = λ 2 κŝ/M 2 Q . The parameters κ,λ,n and P T can be determined by fitting the (differential) cross sections to the experimental data. The dedicated codes to performing the fit and to compute DPS contribution to double-quarkonium production have been implemented in HELAC-ONIA [52].
With the same conditions as in Ref. [20], we reproduced their results. However, after a combined fit of the charmonium data measured at the Tevatron and the LHC, the same parameter set cannot perfectly reproduce the low-energy data measured by PHENIX collaboration [46] at RHIC. Since the collision energy of RHIC √ s = 200 GeV is very close to the centre-of-mass energy of the fixed-target experiment at the LHC (AFTER@LHC) √ s = 115 GeV, we used the PHENIX data alone to determine the parameters in Eq. (9). A combined fit of d 2 σ/dP T dy to the PHENIX data [46] for J/ψ and ψ(2S ) production gives the χ 2 results presented in Tab. 2a having fixed n = 2 and P T = 4.5 GeV. We also show the comparisons of the P T spectra in Fig. 1a There is no such measurement of Υ at RHIC. We therefore performed a combined fit of d 2 σ/dP T dy to CDF [47], ATLAS [48], CMS [49] and LHCb [50,51] data. The results for Υ are presented in Tab. 2b having fixed n = 2 and P T = 13.5 GeV. For illustration, the comparisons between the fit and the LHCb data [50] are shown in Fig. 1d-f. However, we should recall that the results from the Tevatron/LHC fit may underestimate the RHIC total Υ production cross section by a factor of 2. Hence, we present the conservative predictions for the corresponding DPS yields of Υ at AFTER@LHC. All of the above fits are performed with MSTW2008NLO PDF set [53] available in LHAPDF5 [54] and the factorisation scale µ F = M 2 Q + P 2 T . The physical mass M Q for quarkonium is taken from PDG data [55] as well as the branching ratios.
Once a fit is done, |A gg→Q+X | 2 is fixed allowing us to evaluate σ(pp → Q+X) (or its differential counterparts in any variable) which can then be injected into the "pocket formula" Eq. (6) in order to predict the DPS yield. Since we do not apply any muon cuts, we do not need to make any assumptions regarding the polarisation of the production quarkonia.

Single-parton scatterings 2.2.1. Double-charmonium and double-bottomonium production
The SPS contribution for J/ψ-pair production have systematically been investigated in our previous works [16,18]. We have shown that the leading order (LO) in the strong coupling constant α s calculation is enough to account for the low P T data as well as the P T -integrated cross section, the bulk of the events lying at low P T . However, if one goes to mid P T (e.g. P T > 5 GeV), O(α 5 s ) contribution start to be dominant. As a consequence , the yield and the polarisation changes significantly compared to a LO calculation. Since we are only interested in the data in the low P T regime and since our aim is essentially to assess the feasibility of measuring quarkonium pair production with AFTER@LHC, we only perform a LO calculation. Colour-octet contributions in this regime are also negligible: it is suppressed by powers v without any kinematical enhancement 6 at variance with the single-quarkonium-production case. Yet, the feed-down contributions from higher excited spin-triplet S -wave quarkonium might be substantial as already shown for the J/ψpair production in Ref. [18]. We will take these into account. The branching ratios that will be used in this context are taken from PDG [55] and we have listed them in Tab. 3 for completeness.
(a) Decay within a family decay channel branching ratio (%) The general formula for the amplitude of the production of a pair of colour-singlet (CS) S -wave quarkonia Q 1 and Q 2 with as initial partons a and b is where we denote the momenta of quarkonia Q 1 and Q 2 as P 1 and P 2 respectively and their polarisations as λ 1,2 , N(λ 1,2 |s 1,3 , s 2,4 ) are the two spin projectors and R 1,2 (0) are the radial wave functions at the origin in the configuration space for both quarkonia. In the above equation, we have defined the heavy-quark momenta to be q 1,2,3,4 such that P 1,2 = q 1,3 + q 2,4 and p 1,2 = (q 1,3 − q 2,4 )/2. s 1,2,3,4 are then the heavy-quark spin components and δ c i c j / √ N c is the colour projector. The spin-triplet projector N(λ|s i , s j ) has, in the non-relativistic limit, v → 0, the following expression All these computations can be performed automatically in HELAC-ONIA [52] framework based on recursion relations. The radial wave functions at the origin R(0) are taken from Ref. [56], which were derived in the QCD-motivated Buchmüller-Tye potential [57]. We also listed their values in Tab. 4.

Charmonium-bottomonium pair production
The simultaneous production of a charmonium and a bottomonium has been studied in Refs. [13,19]. Its CSM contributions are expected to be suppressed because the LO contributions in CS mechanism (CSM) are O(α 6 s ), i.e. α 2 s suppressed compared to double-charmonium and doublebottomonium production. Hence, it is expected to be a golden channel to probe colour-octet mechanism (COM) at the LHC [13]. However, such a statement is valid only if one can clearly separate DPS and SPS events experimentally since the DPS contributions would be dominant. For a thorough discussion, the reader is guided to [19]. In contrast, colour octet (CO) contributions can appear at O(α 4 s ), which however are suppressed by the small size of the CO long distance matrix elements (LDMEs). If one follows the arguments of Ref. [13], one is entitled to consider only the cc( 1 ) channels. This approximation is however based on the validity of the velocity scaling rules of the LMDEs which may not be that reliable. Since we have also followed this approximation, our SPS results for charmoniumbottomonium production should only be considered as a rough estimate of the yields.
The formula of the amplitude is similar to that of CS state production with the following replacements for CO in Eq.(11) where T a c i c j is the Gell-Mann matrix and O i ( 1 ) should be determined from experimental data. We took the values from a LO fit to the large P T data as described in Refs. [58,59,60] as our reference numbers; these are consistent with those used in Refs. [13,19]. The values of O(  Finally, we describe our parameters for our SPS calculations. In the non-relativistic limit, the mass of heavy quarkonium can be expressed as the sum of the corresponding heavy-quark-pair masses. In our case, we have where m Q = m c for charmonium and m Q = m b for bottomonium. The masses of charm quark and bottom quark are taken as m c = 1.5 ± 0.1 GeV and m b = 4.75 ± 0.25 GeV. The factorisation scale µ F and the renormalisation scale µ R are taken as µ F = µ R ∈ [ 1 2 µ 0 , 2µ 0 ] with µ 0 = (M Q 1 + M Q 2 ) 2 + P 2 T . The advantage of using µ 0 = (M Q 1 + M Q 2 ) 2 + P 2 T is that we are 8 able to recover the correct mass threshold M Q 1 + M Q 2 in the low P T regime. Finally, the PDF set for the SPS calculation is CTEQ6L1 [61] with the one-loop renormalisation group running of α s .

Energy dependence of the ratio DPS over SPS
Due to the very large integrated luminosity of AFTER@LHC (up to 20 fb −1 per year) compared to the experiments performed at RHIC, the measurement of double-quarkonium production at AFTER@LHC will provide a unique test of the interplay between the DPS and SPS production mechanisms in a new energy range. The energy dependence of σ eff will be explored at a wide energy range when combined with the LHC collider and Tevatron data. Due to the double enhancement of the initial gluon-gluon luminosity with the energy, √ s, DPS contributions are expected to be more and more important with respect to the SPS ones at larger √ s. This can be observed on Fig. 2. The fact that we have used a crystal ball fit on the partonic amplitude (gg → QX) to the Tevatron and LHC data to perform the DPS predictions allows us to predict the non-trivial energy (and rapidity) dependence of the DPS predictions.
One however sees on Fig. 2 that a change of σ eff from 15 mb -which seems to be the favoured value for jet-related observables-to 5mb -which is the value extracted by D0 from the J/ψ + J/ψ data [7]-results in a significant change in the point where both contributions are equal. In the former case, it occurs very close to the energy of AFTER@LHC, in the latter case, it occurs between the Tevatron and the LHC energies. All this clearly motivates for measurement and σ eff extractions at low energies.

Predictions at AFTER@LHC
We are now in the position to present our numerical results at √ s = 115 GeV in p+ p collisions. The total cross section we obtained are given in Tab. 5, 6 and 7. The results have been multiplied by the branching ratios into a muon pair and they are all in unit of fb. In general, we have The DPS contributions decrease quickly when the mass threshold M Q 1 + M Q 2 increases because of its square dependence of the initial-state parton luminosity. With the nominal integrated luminosity of 20 fb −1 proposed to be collected at AFTER@LHC, we find that the measurement double-bottomonium production is out of reach 2 and one may be able to record a few J/ψ + Υ(1S ) events, which are DPS dominated. However, one should always keep in mind that σ SPS for ψ + Υ production strongly depends on the CO LDMEs O( ) . An upper limit determination would already be insightful. The quoted theoretical uncertainties include σ eff = 5 ± 2.75 mb for the DPS yields and scale uncertainties as well as heavy quark mass uncertainties for the SPS yields, which were already discussed in Sec.2.
As regards double-charmonium production, about 10 thousand events could be collected per year -which is more than what has so far been collected by LHCb and CMS. In the analysis  of the differential distributions, we therefore only focus on these and, in particular, on J/ψ-pair production. In Fig. 3, we show three interesting distributions without kinematical cuts. Along the lines of [39], we also used the LHCb kinematical acceptance, i.e. the rapidity of J/ψ restricted to be in the interval of [2,5]. The corresponding distributions are shown in Fig. 4.
As we already discussed in Ref. [18], at LO, a 2 → 2 kinematics for SPS would result in the transverse momentum of the J/ψ-pair P ψψ T to be zero. However, it can be smeared by the intrinsic k T of the initial partons. We accounted for such an effect with a Gaussian distribution with k T = 2 GeV as also done in Refs. [16,18]. However, due to the relative smaller yields at AFTER@LHC than at CMS, one can only access to P ψψ T < 10 GeV regime, which can be clearly seen from Fig. 3a and Fig. 4a. In such a kinematical region, the k T smearing effect makes the spectrum of SPS as broad as the DPS one. This difference mostly depends on the value of k T which is essentially empirical.
The absolute rapidity difference between the J/ψ pair is expected to be a good observable to discriminate the DPS and SPS contributions. This was first pointed out in Ref. [20] and this was used later on by D0 collaboration [7] to extract σ eff from double-J/ψ production at the Tevatron. The DPS events should have a broader distribution in ∆y than the SPS ones, because two (relatively) independent hard interactions happen simultaneously in DPS while the two J/ψ from SPS are more correlated. The situation still does not change at AFTER@LHC with or without cut as Fig. 3b and Fig. 4b show. In the latter case, the restriction to negative rapidities in the centre-ofmass obviously reduce the ∆y range (see Fig. 4b). Starting from ∆y = 2, the DPS events dominate the SPS events. A ratio DPS/SPS of 10 is obtained for ∆y > 2. The distribution of the invariant mass for the J/ψ pair M ψψ reflects a similar information as the ∆y distribution. Hence, it follows that the M ψψ spectra of DPS are also broader than those of SPS, which can be seen on Fig. 3c and Fig. 4c.
Finally, we present the cross section as a function of the total rapidity of the J/ψ pair, Y ψψ , and of the sub-leading P T between the J/ψ pair in Fig. 3d and Fig. 4d. One sees that the sub-leading P T spectrum may be measured up to 6 GeV with AFTER@LHC. As regards the rapidity distribution, its maximum is obviously located at Y cms = 0, that is Y = 4.8 in the lab frame. One sees that one can expect some counts down to Y ψψ 2.5 where x F 2M ψψ √ s sinh(Y ψψ − 4.8) −0.5. This is precisely the kinematical region where double intrinsic cc coalescence contributes on average [10]. Any modulation in the pair-rapidity distribution would sign the presence of such a contribution.

Conclusion
We have discussed double-quarkonium production in proton-proton collisions at a fixed-target experiment using the LHC proton beams, AFTER@LHC. These processes have lately attracted much attention, both in the theorist and experimentalist communities. They are expected to be good observables to further constrain the various models describing heavy-quarkonium production. Double-quarkonium production also provides a good opportunity to study DPS since the yields of single quarkonium production is large and their decay to four muons is a clean signal at a hadron colliders. AFTER@LHC provides very appealing opportunities to study these observables with a LHCb-like detector and in new energy region.
In this paper, we have studied both DPS and SPS contributions for double-quarkonium production. These processes include ψ(n 1 S ) + ψ(n 2 S ), ψ(n 1 S ) + Υ(m 1 S ) and Υ(m 1 S ) + Υ(m 2 S ) with n 1 , n 2 = 1, 2 and m 1 , m 2 = 1, 2, 3. DPS contributions are estimated in a data-driven way, while SPS ones are calculated in LO non-relativistic QCD [62]. From our calculations, we find that ten thousand of double-charmonium events can indeed be measured at AFTER@LHC with the yearly integrated luminosity of 20 fb −1 . In the most backward region, a careful analysis of the rapidity distribution could also uncover double intrinsic cc coalescence contributions. In general, future measurements on double-charmonium production can provide extremely valuable information on QCD, in particular important tests on the factorisation formula for DPS and the energy  (in)dependence of σ eff .