Measurement of inclusive J/$\psi$ pair production cross section in pp collisions at $\sqrt{s} = 13$ TeV

The production cross section of inclusive J/$\psi$ pairs in pp collisions at a centre-of-mass energy $\sqrt{s} = 13$ TeV is measured with ALICE. The measurement is performed for J/$\psi$ in the rapidity interval $2.50$. The production cross section of inclusive J/$\psi$ pairs is reported to be $10.3 \pm 2.3 {\rm (stat.)} \pm 1.3 {\rm (syst.)}$ nb in this kinematic interval. The contribution from non-prompt J/$\psi$ (i.e. originated from beauty-hadron decays) to the inclusive sample is evaluated. The effective double-parton scattering cross section is computed, neglecting the single-parton scattering contribution.


Introduction
In the quantum chromodynamics (QCD) parton model [1], hadrons are composed of elementary constituents, the partons.Due to the composite nature of hadrons, multiple parton hard scatterings can occur in a hadron-hadron collision.Thus, it is possible to have two or more hard parton interactions simultaneously.Multiple parton interactions (MPI) have been studied since the introduction of the parton model [2,3].Further studies included the generalization of the QCD evolution equations into multiparton distribution and fragmentation functions [4,5], and a discussion on the possible correlations in the colour and spin degrees of freedom [6].Double-parton scatterings (DPS) are the simplest case of MPI, and were found to play the most important role in processes with final states such as four jets, four leptons or n-jet + W/γ measurements [7][8][9][10][11][12][13][14][15].These studies were complemented by several measurements in hadron collisions at center-of-mass energies ( √ s) ranging from 63 GeV to 1.96 TeV [16][17][18][19][20][21][22].
At the LHC energies, the probability to have multiple parton interactions increases: as with the increase of collision energy, partons with smaller momentum fraction x are probed with larger fluxes.Recent measurements have shown the relevance of MPI at the LHC [23][24][25][26][27][28], and have contributed to stimulate recent progress in the theoretical understanding of MPI [29][30][31][32][33]. Nevertheless, a quantitative estimate of the DPS impact on observables remains challenging.Neglecting the parton correlations in the proton, the DPS contribution to a final state A + B can be evaluated as the product of the parton level cross sections ( σ ) divided by an effective cross section (σ eff ) [12,13,34] where the parameter m is a symmetry factor, m = 1 if A = B, and 2 otherwise.The effective cross section is a phenomenological parameter related to the transverse overlap function between the partons of the proton, and is thought to be universal.It was found to range between 2 and 25 mb [18, 20, 22-25, 27, 35-42].
In the quarkonium sector, quarkonium-pair production is a golden tool to probe the production mechanism of heavy quarkonia [53][54][55].The production mechanism of heavy quarkonia is not fully understood after more than forty years of study, and considered a long-standing puzzle of QCD.The colour-singlet model (CSM), which assumes the formation of an intermediate QQ state with the quantum numbers of the final state, underestimates the production cross section at high p T both at leading-order (LO) and next-to-leading-order (NLO) [56][57][58].The recent CSM next-to-next-to-leading-order NNLO ⋆ calculations have reduced the discrepancies [56,59].Non-relativistic QCD (NRQCD) calculations consider both colour-singlet (CS) and colour-octet (CO) states of the QQ pair [60], but fail to predict at the same time the production cross section and polarisation [61][62][63][64][65].The selection rules for pair production in the CS process of LO NRQCD forbid the feed-down from cascade decays of excited charge-conjugate-even states, e.g.χ c → J/ψ γ, whose contribution is significant in single quarkonium production, and makes the comparison of data to model calculations difficult.As a consequence, quarkonium-pair production provides stringent tests of model calculations.
In this letter we report the measurement of inclusive J/ψ pair production cross section in pp collisions at √ s = 13 TeV at large rapidity (2.5 < y < 4.0) with ALICE.Inclusive J/ψ results correspond to the sum Inclusive J/ψ pair production in pp collisions at √ s = 13 TeV ALICE Collaboration two contributions: the prompt contribution, originated from direct charm decays or decays of highermass excited states; and the non-prompt contribution, stemming from beauty decays.The results corroborate analogous measurements performed in a similar rapidity interval by LHCb [36].They constitute a probe to study quarkonium production mechanisms and the DPS contribution.

Experimental apparatus and data sample
A description of the ALICE detector and its performances can be found in Refs.[66,67].At forward rapidity (2.5 < y < 4.0) the production of quarkonium states is measured in the muon spectrometer down to p T = 0 via their dimuon decay channel.The muon spectrometer of ALICE consists of a ten interaction length thick front absorber to filter muons, five tracking stations of two planes of cathode pad chambers each (MCH), a dipole magnet with a field integral of 3 Tm surrounding the third tracking station, a 1.2 m thick iron wall to absorb secondary hadrons escaping from the front absorber and low momentum muons coming mainly from π and K decays, and two trigger stations made of two planes of resistive plate chambers each (MTR) [68].The silicon pixel detector (SPD) and scintillator arrays (V0) are also used in this analysis.The V0 counters, two arrays of 32 scintillator tiles each, cover 2.8 ≤ η ≤ 5.1 (V0A) and −3.7 ≤ η ≤ −1.7 (V0C) and provide trigger information.The minimum bias (MB) trigger requirement consists of a logical AND of a signal in V0A and in V0C.The SPD, two cylindrical layers covering |η| ≤ 2.0 and |η| ≤ 1.4 for the inner and outer layers, respectively, is dedicated to the vertex reconstruction and allows estimating pile-up.The maximum interaction rate for the analysed data sample is 260 kHz, and the maximum pile-up probability is about 5 × 10 −3 , negligible for this measurement.
The J/ψ pair analysis is performed using data from pp collisions at √ s = 13 TeV collected from 2016 to 2018.The event sample was selected using the dimuon trigger condition, which is defined as the coincidence between the MB requirement and two opposite-charge sign track segments in the muon spectrometer trigger stations.Each track segment in the trigger stations is required to have a transverse momentum, evaluated online, larger than about 0.5 GeV/c.Only events passing a selection criterion to remove beam-background collisions contamination, based on the timing information from the V0 arrays, are considered in the analysis.
When multiple primary vertices are reconstructed by the SPD, the event is tagged as pile-up and removed from this analysis.In order to avoid acceptance biases on the reconstructed SPD tracklets, events with a displaced vertex with respect to centre of the SPD detector along the beam direction are discarded according to the requirement |v z | ≤ 10 cm.These selections allowed us to keep the pile-up below 0.3% for the analysed events, also for events with two muon pairs with an invariant mass above 2 GeV/c 2 .Considering the above selections, the total number of dimuon triggered events in the sample sums up to 587.4 × 10 6 events and corresponds to an integrated luminosity of 24.11 ± 0.01(stat.)± 0.80(syst.)pb −1 .

Analysis
J/ψ candidates are built from muon pairs of opposite-charge sign.Muons are identified by requiring that selected tracks in the MCH have a matching track segment in MTR.Only muon tracks within the detector acceptance are kept for analysis.Tracks are required to be within −4.0 < η µ < −2.5, and the radial distance from the beam axis at the end of the front absorber, R abs , is limited to 17.6 < R abs < 89.5 cm [69].J/ψ pair candidates are reconstructed from all combinations of double dimuon pairs (each dimuon consisting of an opposite-charge sign muon pair) per event.
The production cross section of inclusive J/ψ pairs is determined as Inclusive J/ψ pair production in pp collisions at where N is the signal estimate, ε is the acceptance-times-efficiency correction, B(J/ψ → µ + µ − ) = (5.961± 0.033)% is the branching fraction of J/ψ → µ + µ − [70], and L int is the integrated luminosity.
The J/ψ pair signal is evaluated from a fit to the 2-dimensional invariant mass distribution.A two-step procedure was chosen.The first step exploits the 1-dimensional distribution of all J/ψ candidates in the data sample analysed, to obtain a good description of the J/ψ line shape from data.A fit is performed with a superposition of J/ψ and ψ(2S) signal functions and a background function.The J/ψ mass, width, and normalisation are left free in the procedure.Instead, the ψ(2S) mass and width are bound to those of J/ψ as described in Ref. [71].The 2-dimensional invariant mass distribution of J/ψ pair candidates (m 1 (µ where N and R are the corresponding normalisation parameters.The ψ(2S) contribution is neglected in the 2-dimensional fit.The J/ψ pole mass and width determined from the first step are fixed in the second step of the fit, the rest of the fit parameters are left free.Different combinations of functional forms are used to determine the raw yield and its uncertainties.The signal S is modelled by a Crystal Ball function including a Gaussian core and two asymmetric power-law tails [72].The power-law tail parameters are obtained both from data or Monte Carlo and fixed in the fits [69].The background B contribution is described by either the sum of two exponentials, an exponential of a second order polynomial, or the ratio of a first-order to a second-order polynomials.The mass distribution is fit in two different mass intervals to test the results stability, i.e. [2.0, 4.5] and [2.2, 4.9] GeV/c 2 .As the candidates were assigned randomly, the fit function is symmetric under exchange of m 1 and m 2 .The projections of one of the fits on m 1 (µ + 1 µ − 1 ) and m 2 (µ + 2 µ − 2 ) are shown in Fig. 1.The J/ψ pair signal and statistical uncertainty are evaluated as the average of the values obtained in the twelve fit configurations.The systematic uncertainty is given by their standard deviation.The raw yield is N = 59.3 ± 13.5 (stat) ± 4.4 (syst).The acceptance, reconstruction and selection efficiency is evaluated assuming factorisation of the corrections of the J/ψ pair as ε(J/ψ J/ψ) = ε(J/ψ) × ε(J/ψ) .
The J/ψ ε, ε(J/ψ), is computed from Monte Carlo simulations as described in Ref. [69].An iterative procedure is used to generate input rapidity (y) and transverse momentum (p T ) distributions from data.
Inclusive J/ψ pair production in pp collisions at √ s = 13 TeV ALICE Collaboration The J/ψ are decayed into pairs of muons using EVTGEN [73] and PHOTOS [74].A GEANT3 [75] simulation is performed to transport the decay muons through the apparatus including a realistic description of the detector conditions during data taking.The validity of the factorisation approach for the efficiency calculation was tested.The invariant mass distribution was compared with the corresponding one after applying a 2-dimensional (y, p T ) acceptance-times-efficiency correction per J/ψ candidate.The shapes of the 2-dimensional invariant mass distribution, and their projections are not modified by the correction, confirming the validity of our assumption.In addition, a toy Monte Carlo was developed to study the possible influence of angular correlations between the two J/ψ of the pair.Two J/ψ were simulated per event, according to a (y, p T ) distribution extracted from single J/ψ measurements.To mimic possible correlations among the J/ψ, their rapidity difference was forced to follow either a triangular or a flat distribution.The average pair efficiency was computed for both cases.The resultant pair efficiency was found to be in agreement with the calculation from the factorisation approach for both cases.
Various sources of systematic uncertainties on the J/ψ pair production cross section are considered: (i) the signal extraction, (ii) the branching fraction uncertainty, (iii) the luminosity normalisation, and (iv) the acceptance-times-efficiency correction.Details on the signal extraction procedure were given previously in this Letter.The systematic uncertainty on the signal extraction, obtained as described above, amounts to 7.4%.The branching fraction uncertainty is 0.6% for single J/ψ [70], thus 1.1% for J/ψ pairs.The influence of the luminosity normalisation factor is evaluated by computing the equivalent number of minimum-bias events in the analysed dimuon sample with different methods as described in Ref. [76], which amounts to 2.9%.The uncertainty on the minimum bias cross section, evaluated in a van der Meer scan (1.6%), is also taken into account in the calculation [77].These two sources lead to a 3.3% systematic uncertainty for the luminosity.The systematic uncertainty on the acceptance-times-efficiency correction contains contributions from (i) the input p T and y distributions, (ii) the tracking efficiency in the MCH, (iii) the MTR trigger efficiency, and (iv) the matching of the reconstructed tracks in the MCH with the track segments in the MTR.The influence of the simulated J/ψ p T and y distributions is tested by comparing the corrected yield obtained via the iterative procedure, with the one obtained from an efficiency-corrected invariant mass distribution.For this exercise, a 2-dimensional ε(p T , y) correction is applied to each J/ψ candidate in order to build the efficiency-corrected invariant mass distribution, which was then fit to obtain the corresponding corrected yield.A 0.5% uncertainty is assigned to the MC input for J/ψ [69].The systematic uncertainties on the tracking efficiency in the MCH, the MTR trigger efficiency and the matching between the MCH and MTR are evaluated comparing single muon data and MC, as described in Ref. [78].The differences are then propagated to the dimuon case, being 4%, 2% and 1% respectively for the J/ψ [69].This results in a 4.6% acceptance-times-efficiency uncertainty for J/ψ, and is propagated to a 9.2% uncertainty for J/ψ pairs.The analysis requirement that all selected tracks in the MCH should match track segments in the MTR removes any possible dependence on which pair of tracks activated the trigger.Table 1 summarises the systematic uncertainties on the measurement of the J/ψ pair production cross section.
The ratio of the production cross section of the inclusive J/ψ pair to that of the inclusive J/ψ is considering dσ (J/ψ)/dy = (7533.3± 26.7 (stat.)± 491.6 (syst.))nb for p T > 0 and 2.5 < y < 4.0 [69], and assuming the systematic uncertainties to be uncorrelated.Likewise, the ratio can be calculated and interpreted as an effective cross section, according to Eq. ( 1).This interpretation assumes that all J/ψ pairs are produced via DPS processes.The relative contribution of SPS and DPS processes to J/ψ pair production is a topic of debate and intense studies, see e.g.Ref. [36].In addition, the understanding of this ratio gets challenged by the contribution of both the prompt and non-prompt components to the measured inclusive J/ψ cross section, where the non-prompt contribution originates from beauty-hadron decays.
The contamination from beauty-hadron decays to the J/ψ pair cross section is evaluated to assess the impact on the measurement according to The total beauty-hadron production cross section was measured to be σ total bb = 502 ± 16 (stat.)± 51 (syst.)+2 −3 (extr.)µb in Ref. [79].The branching ratio of a beauty hadron into a J/ψ is B(h b → J/ψ +X) = (1.16±0.10)%[70], and the acceptance correction factor α is estimated using PYTHIA 8.3 [80] simulations.Beauty hadrons are simulated according to three different configurations and forced to decay into J/ψ.The three configurations use the Monash 2013 tune for the calculation [81].Two of them also include a tuning of the parameters to get a good agreement with the NLO calculation by Mangano, Nason, and Ridolfi for the bb single and double differential distributions [82].The difference between the latter two is that one of them adds the ATLAS tune settings for multiple parton interactions [83].The α factor is obtained from the ratio of the J/ψ pair counts in the acceptance to the number of all J/ψ pairs in the simulation.The value of α = 0.044 +0.005 −0.007 is determined as the average of the factors obtained with all configurations, and the systematic uncertainty is conservatively set to the full spread of the values.This gives a non-prompt contribution of σ non-prompt (J/ψ J/ψ) = 2.97 ± 0.09 (stat.)+0.68 −0.76 (syst.)nb, and correspondingly the prompt J/ψ pair cross section is σ prompt (J/ψ J/ψ) = σ (J/ψ J/ψ) − σ non-prompt (J/ψ J/ψ) = 7.3 ± 1.7 (stat.)+1.9 −2.1 (syst.)nb.
A differential measurement of the prompt J/ψ pair production cross section and the corresponding ratios were previously reported by the LHCb collaboration in a slightly different kinematic interval, 2.0 < y < 4.5 and p T < 10 GeV/c [36,46].The results presented here are in agreement with the LHCb ones within uncertainties.
Despite the caveat caused by the calculation of this effective value considering both the SPS and DPS contributions to the production cross section, this value is consistent with the values obtained from quarkonium-pair production measurements, with σ eff values ranging from 2.2 to 12.5 mb [22,[35][36][37], and with the values obtained for quarkonium associated production at central rapidity (in the range 2.3-6.1 mb) [38,39,59].It is smaller than the values obtained for associated heavy-flavour production at large rapidity by LHCb (ranging from 12.8 to 18.0 mb) [23,27], or those from jet or electroweak associated production (whose values are between 12.0 and 21.3 mb) [18,20,24,25,[40][41][42].

Conclusion
The production cross section of J/ψ pairs at large rapidity in pp collisions at √ s = 13 TeV was studied by ALICE.The measurement exploits the full Run 2 data sample collected by ALICE.The production cross section of inclusive J/ψ pairs is reported to be 10.3 ± 2.3 (stat.)± 1.3 (syst.)nb, for J/ψ in the rapidity interval 2.5 < y < 4.0 and for p T > 0. The effective double-parton scattering cross section is evaluated neglecting the single-parton scattering contribution.The results are compatible with analogous measurements performed by the LHCb collaboration in a similar kinematic interval [36,46].The Run 3 data taking, with the upgraded ALICE detector and the larger accumulated luminosity [84], will allow us to perform this measurement with increased precision and separating the prompt and nonprompt contributions.This will also enable studying the kinematics of these events and probe model calculations. Inclusive

Figure 1 :
Figure 1: Projections of a fit to the 2-dimensional invariant mass distribution for inclusive J/ψ pair candidates for (a) m 1(µ + 1 µ − 1 ) and (b) m 2 (µ + 2 µ − 2 ).The (black) markers show data.The (black) solid line represents the total fit function.The (blue and green) dashed-dotted lines indicate the background contribution from a combination of a real J/ψ signal with a combinatorial candidate.The (yellow) dotted line represents muon pairs from combinatorial background.

Table 1 :
Sources of systematic uncertainty on the J/ψ pair production cross section measurement.
41 ± 0.04 (stat.)± 0.19(syst.)µb, J/ψ pair production in pp collisions at √ s = 13 TeV ALICE Collaboration Technologies, National Nuclear Research Center, Azerbaijan; Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Financiadora de Estudos e Projetos (Finep), Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and Universidade Federal do Rio Grande do Sul (UFRGS), Brazil; Bulgarian Ministry of Education and Science, within the National Roadmap for Research Infrastructures 2020¿2027 (object CERN), Bulgaria; Ministry of Education of China (MOEC) , Ministry of Science & Technology of China (MSTC) and National Natural Science Foundation of China (NSFC), China; Ministry of Science and Education and Croatian Science Foundation, Croatia; Centro de Aplicaciones Tecnológicas y Desarrollo Nuclear (CEADEN), Cubaenergía, Cuba; Ministry of Education, Youth and Sports of the Czech Republic, Czech Republic; The Danish Council for Independent Re-