Measurement of J/ψ -pair production in pp collisions at √ s = 13 TeV and study of gluon transverse-momentum dependent PDFs

,


Introduction
In the Standard Model (SM) of particle physics, the strong interaction is described by quantum chromodynamics (QCD).In the low-energy regime, the QCD coupling constant evolves to be so large that perturbation theory is not valid any more.The details of how the fundamental quarks and gluons are distributed inside hadrons and dynamically generate the hadron mass and spin are still largely unknown.As a building block of the physical world, the proton lies at the forefront of hadron structure studies.Moreover, the knowledge of the partonic structure of the proton is an essential input to the majority of measurements at hadron colliders, including searches for physics beyond the SM, among which the W -mass measurement is a typical example [1][2][3][4].
In the past, the internal structure of the proton was mainly studied in terms of one-dimensional parton distribution functions (PDFs) that parameterise the longitudinal momentum fraction (usually denoted x) distributions of partons inside the proton.Recently, significant progress has been made in constructing the theoretical framework for transverse-momentum dependent parton distribution functions (TMDs) [5][6][7], leading to a more comprehensive understanding of the proton structure.The quark TMDs have been studied through semi-inclusive deep-inelastic scattering and Drell-Yan measurements at HERMES, COMPASS and E866/NuSea experiments, and at a series of experiments at JLab [8][9][10][11][12][13].The gluon TMDs, however, are much less known, so probing gluon TMDs is one of the major objectives of future experimental facilities like the Electron Ion Collider (EIC) [14], LHCSpin [15] and LHC fixed-target experiments [16].In unpolarised protons, the gluon TMDs can be parameterised at leading twist using two TMDs [5]: the distribution of unpolarised gluons f g 1 (x, k 2 T , µ) and that of linearly polarised gluons h ⊥g 1 (x, k 2 T , µ), in which k T is the gluon transverse momentum and µ is the factorisation scale.In particular, the knowledge on the h ⊥g 1 (x, k 2 T , µ) function is still very limited.The production of J/ψ pairs in proton-proton (pp) collisions through single-parton scattering (SPS) has been proposed as the golden channel to probe gluon TMDs, in which the presence of linearly polarised gluons will lead to azimuthal asymmetries at the percent level [17,18].The transverse momentum (p T ) spectrum of J/ψ pairs also encodes information on f g 1 (x, k 2 T , µ).In fact, the differential production cross-section of J/ψ pairs as a function of p T measured by the LHCb experiment using the 2015 data was used to perform the first fit of f g 1 (x, k 2 T , µ) and obtain ⟨k 2 T ⟩ at an effective factorisation scale [17,19].Quarkonium production is also one of the best tools to study hadronisation.The nonrelativistic QCD (NRQCD) model provides the most successful description of quarkonium production so far, but it still can not describe coherently the production and polarisation of various quarkonium states measured in different collisions [20][21][22][23].The SPS production of J/ψ pairs (also referred to as di-J/ψ hereafter) can add valuable information to solve this puzzle [24][25][26][27].In addition to SPS, quarkonium pairs can be produced through double-parton scattering (DPS) [28].It is a process of great interest, which has been widely studied by many experiments via various reactions [29][30][31][32][33][34][35][36][37][38][39][40][41][42].It can be used, for instance, to reveal the profile and correlation of partons inside the proton, which are encoded in a characteristic parameter of DPS called effective cross-section, denoted as σ eff [43][44][45].The contribution from DPS can be estimated according to the formula [43][44][45] σ DPS di-J/ψ = 1 2 where σ J/ψ is the prompt J/ψ meson production cross-section, and the factor one-half accounts for the two identical particles in the final state.The effective cross-section σ eff characterises the transverse overlap area between the interacting partons.In this paper, the J/ψ-pair production cross-section in pp collisions at a centre-ofmass energy of √ s = 13 TeV is measured for J/ψ mesons with p T < 14 GeV/c and rapidity 2.0 < y < 4.5 using a subset of data collected by the LHCb experiment from 2016 to 2018 with specific trigger requirements, corresponding to an integrated luminosity of 4.2 fb −1 .The azimuthal asymmetry of J/ψ pairs is measured to probe the TMD function h ⊥g 1 (x, k 2 T , µ), presenting the first experimental measurement of linear polarisation of gluons inside unpolarised protons.The p T spectrum of the J/ψ pairs is measured in intervals of J/ψ-pair rapidity and mass, which will help to extract f g 1 (x, k 2 T , µ) with the TMD evolution effect considered [18].Updates of the differential production cross-sections given in the previous LHCb measurement using data with a luminosity of about 0.3 fb −1 at √ s = 13 TeV [19] are also provided, with the SPS and DPS contributions separated without dependence on any specific SPS production model.

Detector and simulation
The LHCb detector [46,47] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b-or c-quarks.The detector includes a high-precision tracking system consisting of a siliconstrip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet.The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c.The minimum distance of a track to a primary pp collision vertex (PV), the impact parameter (IP), is measured with a resolution of (15 + 29/p T ) µm, where p T is in GeV/c.Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors.Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter.Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.
Simulated samples of J/ψ mesons are produced to study the expected behaviour of experimental signals and determine the detection efficiencies.The pp collisions are modelled using Pythia [48,49] with a specific LHCb configuration [50].In the Pythia model, J/ψ mesons are not polarised, and the leading order colour-singlet and colouroctet contributions [50,51] are included in prompt J/ψ meson production.Decays of unstable particles are described by EvtGen [52] with QED final-state radiation handled by Photos [53].The interactions of the generated particles with the detector are modelled using the Geant4 toolkit [54,55] as described in Ref. [56].

Candidate selection
The di-J/ψ candidates are reconstructed through the J/ψ → µ + µ − decays.The online event selection is performed by a trigger [57], which consists of a hardware stage (L0), based on information from the calorimeter and muon systems, followed by a two-step software stage (HLT1 and HLT2), which applies a full event reconstruction.At least one J/ψ meson is required to fulfil the selection criteria of the L0 and HLT1 triggers.The L0 trigger selects two muons with the product of their transverse momenta larger than 1.3 2 , 1.5 2 or 1.8 2 (GeV/c) 2 , depending on the data taking period.The HLT1 trigger requires two good-quality tracks with p T > 0.3 GeV/c and p > 6 GeV/c, that are loosely identified as muons and form a J/ψ candidate with an invariant mass m µ + µ − > 2.7 GeV/c 2 or 2.9 GeV/c 2 , depending on the data taking period.The HLT2 trigger requires both µ + µ − pairs to form a good vertex and have an invariant mass m µ + µ − within 120 MeV/c 2 of the world-average value of the J/ψ mass [58].Offline selections are applied to the di-J/ψ candidates to further reduce the combinatorial background.All the four muon tracks are required to have 1.9 < η < 4.9, p T > 0.65 GeV/c and p > 3 GeV/c.
In pp collisions, J/ψ mesons can be produced promptly at the PV, or in the decays of beauty hadrons.The nonprompt contributions with the J/ψ mesons originating from the decay vertices of beauty hadrons, which are typically separated from the PV, need to be subtracted in this analysis.The prompt and nonprompt J/ψ mesons can be distinguished by exploiting the pseudoproper time [59] where z J/ψ and z PV are the positions of the J/ψ meson decay vertex and its associated PV along the beam axis z, p z the component of the J/ψ momentum along the z-axis, and m J/ψ the world-average value of the J/ψ mass [58].The uncertainty σ tz is calculated by combining the uncertainties on z J/ψ and z PV since the uncertainty on p z is negligible in comparison.Candidates with both J/ψ mesons having −2 < t z < 10 ps and σ tz < 0.3 ps are retained.Finally, the four muon tracks are required to originate from the same PV, which reduces to a negligible level the number of candidates with the two J/ψ mesons originating from different pp interactions.
The di-J/ψ signals are extracted by performing a two-dimensional (2D) unbinned extended maximum likelihood fit to the distribution of the two J/ψ meson masses, . The two J/ψ → µ + µ − decays are labelled as J/ψ 1 → µ + 1 µ − 1 and J/ψ 2 → µ + 2 µ − 2 at random.There are four components of the 2D mass distribution: a signal di-J/ψ decay, a true J/ψ 1 → µ + 1 µ − 1 decay with a dimuon background µ + 2 µ − 2 , a dimuon background µ + 1 µ − 1 with a true J/ψ 2 → µ + 2 µ − 2 decay, and the association of two combinatorial dimuon backgrounds.The second and the third components have the equal fraction because the  mass distributions of the two J/ψ mesons are symmetric.For the mass distribution of each J/ψ meson, the signal component is described by the sum of a double-sided Crystal Ball (DSCB) function [61] with asymmetric tails and a Gaussian function with a common mean value but different widths.The tail parameters of the DSCB function, the fraction of the DSCB function and the ratio between the two widths are fixed from simulation.Only the common mean value and the width of the DSCB function are left as free parameters.The distribution of the combinatorial dimuon background component is modelled with an exponential function.Figure 1 shows the mass distributions for the two µ + µ − pairs and the projections of the fit result on m µ The subtraction of the nonprompt contributions is performed using a 2D unbinned extended maximum likelihood fit to the t z distribution of each J/ψ meson of the di-J/ψ signals with backgrounds subtracted using the sPlot method [62], taking m µ + 1 µ − 1 and m µ + 2 µ − 2 as the discriminating variables.The true t z distribution of prompt J/ψ mesons is expected to follow a Dirac delta function δ(t z ), while that of nonprompt J/ψ mesons should follow an exponential function.These are convolved with the sum of two Gaussian functions to model the detector resolution.The parameters of the resolution function are free parameters in the fit.A third component describes candidates with incorrectly associated PVs and is modelled using a binned histogram extracted from data by calculating t z with one J/ψ meson associated to the closest PV in the next event.Figure 2 shows the t z distributions for the two J/ψ mesons and the projections of the 2D fit result on t J/ψ 1 z and t . The yield of prompt di-J/ψ signal is (2.187 ± 0.020) × 10 4 , accounting for around 68% of the total di-J/ψ signal yield.
Since the kinematic distribution of the di-J/ψ signal is unknown a priori, the efficiency correction is performed on a per-event basis as where i is the event index, ω i is the sPlot weight corresponding to the component of prompt di-J/ψ signal, and ε tot i is the total efficiency for each candidate.Since no information related to the correlation of the two J/ψ mesons is used during the reconstruction and selection, the detection efficiency of the di-J/ψ candidate can be factorised into that of the two J/ψ mesons, which are determined from simulation as functions of the p T and y of the J/ψ mesons.The efficiency ε tot i is the product of the geometrical acceptance ε acc i , the reconstruction and selection efficiency ε rec&sel i , the particle identification (PID) efficiency ε PID i , and the trigger efficiency ε trig i , giving The efficiencies ε acc i , ε rec&sel i and ε PID i of a di-J/ψ candidate factorise as the product of that of the two J/ψ mesons, i.e.
while the trigger efficiency ε trig i of each di-J/ψ candidate factorises as where ε L0&HLT1 i is the L0 and HLT1 trigger efficiency.The offline selection criterion is tighter than the HLT2 requirements, making the HLT2 trigger fully efficient with respect to offline selected candidates.The track reconstruction efficiency, which is part of ε rec&sel i , and the PID efficiency ε PID i are calibrated using data-driven techniques to avoid known discrepancies between the simulation and data [63,64].With the detection efficiency corrected, the signal yield is determined to be N corr = (2.43 ± 0.04) × 10 5 , where the statistical uncertainty is verified by a bootstrapping approach [65].

Systematic uncertainties
Systematic uncertainties are studied and summarised in Table 1.The uncertainty due to imperfect modelling of the J/ψ mass distribution is estimated by using an alternative function for the J/ψ signal component in the fit.A model derived from the simulation using kernel density estimation [66] is used instead.The relative variation of the extracted signal yield is 1.7%, which is taken as the systematic uncertainty.The uncertainty on the determination of the nonprompt contribution is evaluated by using an alternative variable for the discrimination.The log(χ 2 IP ) of the J/ψ mesons is used instead of t z , where χ 2 IP is defined as the difference in the vertex-fit χ 2 of the associated PV reconstructed with and without the J/ψ meson under consideration.The relative deviation from the nominal result, 2.4%, is taken as the systematic uncertainty.For the measurement of differential cross-sections, these two uncertainties are taken as common to all kinematic intervals.For a very small fraction of the di-J/ψ candidates, the J/ψ mesons may be associated to a wrong PV in two cases.In the first, the true PV is reconstructed but one of the two J/ψ candidates is associated to a wrong PV, while in the second the true PV is not reconstructed and both J/ψ mesons are associated to the same reconstructed PV in this event.For the first case, two J/ψ mesons are associated to different PVs and thus are rejected by the selection.The fraction of the first case is estimated using a simulated sample of Υ (1S) → J/ψJ/ψγ decays to be (0.65 ± 0.02)%.The second case is the wrong-PV component in the 2D t z fit, and its fraction is (0.16 ± 0.08)%.The total fraction of wrong-PV candidates adds up to 0.8%, which is taken as the systematic uncertainty.For the differential cross-sections, the fraction is studied in several intervals of the di-J/ψ rapidities using the same method.The maximum fraction, 1.5%, is conservatively taken as the uncertainty common to all kinematic intervals.
The uncertainty due to the finite sample size of the calibration samples used for efficiency determination is propagated to the final result using pseudoexperiments.It is determined to be 0.2% and varies up to 1.1% depending on the kinematic intervals for the differential cross-sections.The binning scheme that is used to determine the efficiencies could bias the signal yield N corr .For the PID efficiency, this effect is estimated by varying the binning schemes of the PID calibration sample.For the remaining efficiency terms, an alternative kernel density estimation approach [66] is used to determine the efficiencies as functions of (p J/ψ T , y J/ψ ).The relative difference in the cross-sections between the default and alternative approaches is 1.5% and quoted as the systematic uncertainty.For the measurement of differential cross-sections, it varies up to 6.6% depending on the kinematic intervals.The uncertainty on the track reconstruction efficiency consists of the statistical uncertainty due to the limited size of the calibration samples, which is estimated from pseudoexperiments to be 1.2% for di-J/ψ candidates, and the uncertainty due to the dependence of calibration factors on the event multiplicity, which is 0.8% per track.The two terms are added in quadrature to 3.4%.For the differential cross-sections, the first term varies up to 3.8%, while the second is considered to be common to all kinematic intervals.The trigger efficiency determined from the simulation is validated with the prompt J/ψ data, using a subset of events that fulfil the trigger requirement with the J/ψ signals excluded [57].The relative difference in the di-J/ψ cross-section calculated using the trigger efficiencies from the data and the simulation is 0.7%, and is taken as the systematic uncertainty.For the differential cross-sections, the uncertainty on the trigger efficiency varies up to 7.9% depending on the kinematic intervals.
The uncertainty on the J/ψ → µ + µ − branching fraction leads to an uncertainty of 1.1% on the di-J/ψ production cross-section.The relative uncertainty on the luminosity is determined to be 2.0%.With these uncertainties due to independent effects added in quadrature, the total systematic uncertainty on the di-J/ψ production cross-section is determined to be 5.4%.6 Production cross-sections The cross-section of the di-J/ψ production with both J/ψ mesons in the fiducial range p T < 14 GeV/c and 2.0 < y < 4.5 is measured to be where the first uncertainty is statistical and the second systematic, assuming negligible polarisation of the J/ψ mesons.The detection efficiency for the J/ψ mesons depends on their polarisation.In the helicity frame, there is a strong dependence on the polarisation parameter λ θ [67,68].For instance, when λ θ is assumed to be +0.2 (−0.2) for both J/ψ mesons, the di-J/ψ production cross-section changes by +6.2% (−6.3%) evaluated from simulation.
The differential di-J/ψ production cross-section as a function of a kinematic variable u is measured as where ∆N corr is the efficiency-corrected signal yield in an interval of the variable u, and ∆u is the interval width.In this analysis the cross-sections are reported as functions of: the absolute difference in rapidity between the two J/ψ mesons ∆y; the absolute difference in the azimuthal angle ϕ, defined in the laboratory frame, between the two J/ψ mesons ∆ϕ; the transverse momentum asymmetry A p T of the two J/ψ mesons, defined as the transverse momentum p di-J/ψ T , rapidity y di-J/ψ and invariant mass m di-J/ψ of the di-J/ψ signals; and the transverse momentum p J/ψ T and rapidity y J/ψ of either J/ψ meson.The binning scheme is chosen to have adequate and approximately even signal yield in each category.The differential cross-section as a function of p J/ψ T (y J/ψ ) is taken as the average of the two distributions of p J/ψ 1 T (y J/ψ 1 ) and p J/ψ 2 T (y J/ψ 2 ), taking advantage of the symmetry between the two J/ψ mesons.The 2D mass fit and the 2D pseudoproper time fit are performed independently for each kinematic interval to subtract the combinatorial backgrounds and the nonprompt contributions, respectively.The measured differential cross-sections of the di-J/ψ production are shown in Figure 3 as black data points (SPS+DPS), and summarised in Tables 2-9 in Appendix A.

Separation of DPS and SPS contributions
The distributions of the rapidity difference (∆y) between the two J/ψ mesons have different shapes for the SPS and DPS processes [31,39,69], so the DPS contribution can be extracted using the ∆y distribution with a data-driven template for the DPS process.The shape of the DPS component is obtained by combining two uncorrelated J/ψ mesons whose distributions follow the measured differential production cross-section of the single prompt J/ψ [68], assuming they are both uniformly distributed over the azimuthal angle ϕ.According to the NRQCD predictions [24][25][26][27], the SPS contribution to the di-J/ψ production in the range 1.8 < ∆y < 2.5 is negligible.The normalisation of the DPS contribution is thus determined in this range, while the remaining contribution is assigned to SPS.For the extraction of the DPS contribution, three sources of systematic uncertainties are considered: the uncertainty due to possible SPS remnant in the 1.8 < ∆y < 2.5 range, which is studied by varying the ∆y range for normalisation and estimated to be 3.5%; the uncertainty on the DPS template due to the binning scheme of the prompt J/ψ differential cross-sections, which is estimated to be 3.3%, determined through the relative difference between the results with and without interpolation across the intervals; the uncertainty propagated from the prompt J/ψ production cross-section, which is 1.8%.Consequently, the DPS cross-section of the di-J/ψ production with both J/ψ mesons in the fiducial range p T < 14 GeV/c and 2.0 < y < 4.5 is determined to be σ DPS di-J/ψ = 8.6 ± 1.2 (stat) ± 1.0 (syst) nb.
According to Eq. 1, the effective cross-section σ eff is measured to be where the prompt J/ψ production cross-section σ J/ψ in the range 0 < p T < 14 GeV/c and 2.0 < y < 4.5 is σ J/ψ = 15.03 ± 0.03 (stat) ± 0.94 (syst) µb [68], and the systematic uncertainties on σ DPS di-J/ψ and σ J/ψ are treated as uncorrelated.The σ eff result is compatible with existing measurements from different experiments in pp and pp collisions, as shown in Figure 4.With the DPS cross-section subtracted, the SPS cross-section of di-J/ψ production with both J/ψ mesons in the fiducial range p T < 14 GeV/c and 2.0 < y < 4.5 is determined to be σ SPS di-J/ψ = 7.9 ± 1.2 (stat) ± 1.1 (syst) nb.The differential di-J/ψ production cross-sections are shown in Figure 3 with the DPS and SPS contributions separated.The differential cross-sections for DPS are obtained by normalising the data-driven DPS template to the total DPS cross-section in the fiducial range, and the uncertainties are propagated from the total DPS cross-section.The extracted differential SPS cross-sections are listed in Tables 10-17 in Appendix A. By definition the SPS and DPS components are anti-correlated.In general, the difference  in distributions between DPS and SPS is due to the fact that the kinematics of two J/ψ mesons are uncorrelated for DPS while correlated for SPS. Figure 3(a) shows that the ∆y distribution for the DPS process is wider than that for the SPS process.As shown  in Figure 3(b), the ∆ϕ distribution for SPS peaks at ∆ϕ = 0 and π, while the DPS distribution is flat because two uncorrelated J/ψ mesons are assumed to be uniformly distributed over the angle ϕ in the data-driven template.The distributions of A p T and m di-J/ψ for DPS are both slightly wider than those for SPS, as shown in Figures 3(c) and 3(f).
sections are compared with the NLO* CS predictions, as shown in Fig. 5.The NLO* CS predictions include the uncertainties from the factorisation and renormalisation scales and subleading PDF uncertainties, correlated between the intervals.The measurements Figure 6: Differential cross-section of di-J/ψ production for SPS+DPS and SPS as a function of p di-J/ψ T for candidates with 9 < m di-J/ψ < 24 GeV/c 2 , compared with the PRA+NRQCD predictions for SPS [77].
are consistent with the NLO* CS calculations within the large theoretical uncertainties.The mass spectrum is also compared with the predictions combining parton Reggeization approach (PRA) [78] and NRQCD factorisation [77], which includes a subset of higherorder QCD corrections without ad-hoc kinematic cuts.Only the renormalisation and factorisation scale uncertainties are considered in the predictions.In the large m di-J/ψ region, where the NRQCD-based calculations are well justified, there is a good agreement with the data, while the same predictions exceed the SPS data at small m di-J/ψ , as shown in Fig. 5(f).The PRA+NRQCD calculations on the p di-J/ψ T spectrum are thus further compared with that of the SPS data with 9 < m di-J/ψ < 24 GeV/c 2 , where the PRA+NRACD approach is expected to be validated, as shown in Fig. 6.The predictions are consistent with the SPS cross-sections at small p di-J/ψ T , but they overestimate the SPS data at p di-J/ψ T larger than 3 GeV/c.

Study of gluon TMDs
The gluon TMD h ⊥g 1 (x, k 2 T , µ) inside unpolarised protons, representing the distribution of linearly polarised gluons, can be probed through the distribution of the azimuthal angle ϕ CS of either J/ψ meson in the Collins-Soper frame.This frame is the rest frame of the J/ψ pair with the polar axis (z-axis) bisecting the angle between the momentum of one proton and the reverse of the momentum of the other proton, the y-axis defined to be perpendicular to the plane spanned by the momenta of two protons, and the x-axis defined to complete a right-handed Cartesian coordinate system [79].The prediction of the differential di-J/ψ production cross-section as a function of ϕ CS through SPS is proportional to a + b × cos(2ϕ CS ) + c × cos(4ϕ CS ).The parameters a, b and c encode  information on the gluon TMDs as where F i ( ′ ) are hard-scattering coefficients, w i ( ′ ) are the TMD weights common to all gluon-fusion processes originating from unpolarised proton collisions, and C denotes the TMD convolutions [17,18].The calculation is valid in the TMD region with p di-J/ψ T < ⟨m di-J/ψ ⟩/2 [17,18].In this analysis, the ϕ CS distribution is measured in the TMD region p di-J/ψ T < 4.1 GeV/c, since the average value of m di-J/ψ in the whole fiducial range is ⟨m di-J/ψ ⟩ = 8.2 GeV/c 2 .The measured ϕ CS distributions with the SPS and DPS contributions separated are shown in Fig. 7(a).The expectation values ⟨cos 2ϕ CS ⟩ and ⟨cos 4ϕ CS ⟩ correspond to half of the ratio of the cos nϕ CS -modulations present in the TMD cross-section regarding its ϕ CS -independent component [18], i.e. ⟨cos 2ϕ CS ⟩ = b/2a and ⟨cos 4ϕ CS ⟩ = c/2a.They are calculated as where the index i denotes each interval, ∆ϕ CS i is the interval width and ϕ CS i is the interval centre

fb SPS
Figure 8: Normalised p T spectrum of di-J/ψ production in different y di-J/ψ intervals, compared with TMD predictions [18] in the TMD region p di-J/ψ T < ⟨m di-J/ψ ⟩/2.The average values of the p di-J/ψ T distributions in three y di-J/ψ intervals are presented at the top of the figure .fixed to the values calculated by Eqs. 13 and 14, and the normalisation is fixed to that of the SPS measurement.The results are consistent with zero, but the presence of an azimuthal asymmetry at a few percent level is allowed.The prediction of ⟨cos 2ϕ CS ⟩ varies from 0.009 to 0.016 due to nonperturbative uncertainties [18], also consistent with the measured result given the large uncertainty so far.
The p T spectrum of the di-J/ψ signals from SPS can also be used to probe the gluon TMDs, especially f g 1 (x, k 2 T , µ) [17,18].It was pointed out in Ref. [18] that the variation of the momentum fractions of the two interacting gluons, x 1,2 = m di-J/ψ e ±y di-J/ψ / √ s, do not have significant impact on the shape of the p di-J/ψ T spectrum.The p di-J/ψ T spectrum is thus measured in three different intervals of y di-J/ψ for the SPS process, and the cross-section results are listed in Tables 18 and 19 in Appendix A for SPS+DPS and SPS separately.The distributions are normalised for comparison in Fig. 8.They are consistent with each other within the uncertainties.The average values of the p di-J/ψ T distributions in three y di-J/ψ intervals are also presented at the top of Fig. 8, and show no significant variations.The TMD predictions [18], which are only applicable in the TMD region p di-J/ψ T < ⟨m di-J/ψ ⟩/2, are also shown in Fig. 8, and peak at higher p di-J/ψ T than the measured distributions.
In addition, the study of the dependence of TMDs on the renormalisation and rapidity scales, requires a measurement of the p T spectrum at different m di-J/ψ [18].The differential cross-sections dσ/dp di-J/ψ T in the three intervals 6 < m di-J/ψ < 7 GeV/c 2 , 7 < m di-J/ψ < 9 GeV/c 2 and 9 < m di-J/ψ < 24 GeV/c 2 , are listed in Tables 20 and 21 in Appendix A for SPS+DPS and SPS separately.The normalised p T spectra of the di-J/ψ production for SPS in different m di-J/ψ intervals with the expected values of ⟨m di-J/ψ ⟩ = 6.6, 7.9 and 11.0 GeV/c 2 , respectively, are compared in Figure 9, with the TMD predictions [18]

fb SPS
Figure 9: Normalised p T spectrum of di-J/ψ production in three m di-J/ψ intervals with ⟨m di-J/ψ ⟩ =6.6, 7.9 and 11.0 GeV/c 2 , compared with TMD predictions [18] in the TMD region p di-J/ψ T < ⟨m di-J/ψ ⟩/2.The average values of the p di-J/ψ T distributions in three m di-J/ψ intervals are presented at the top of the figure .overlaid in the TMD region.According to the prediction, the p T spectrum would broaden as m di-J/ψ increases [18], but no obvious broadening of the p T spectrum can be seen in the TMD region due to the large uncertainties.The average values of the p di-J/ψ T distributions in three m di-J/ψ intervals are also presented at the top of Fig. 9, and slightly increase with mass.

Conclusion
The J/ψ-pair production cross-section in pp collisions at √ s = 13 TeV is measured to be 16.36 ± 0.28 (stat) ± 0.88 (syst) nb using a data sample corresponding to an integrated luminosity of 4.2 fb −1 collected by the LHCb experiment, with both J/ψ mesons in the range of p T < 14 GeV/c and 2.0 < y < 4.5.The contributions from DPS and SPS are separated based on distinctive ∆y dependences of their corresponding crosssections.The effective cross-section characterising the DPS process is determined to be σ eff = 13.1 ± 1.8 (stat) ± 2.3 (syst) mb, and is consistent with most of the existing measurements.The differential cross-sections in SPS are consistent with the NLO* CS predictions which are plagued by large theoretical uncertainties.The cross-sections predicted by PRA+NRQCD overshoot the SPS data at small m di-J/ψ and agree with them at large m di-J/ψ .
The gluon TMDs are probed via the ϕ CS distribution and the p di-J/ψ T spectrum from the SPS process.The extracted values of ⟨cos 2ϕ CS ⟩ and ⟨cos 4ϕ CS ⟩ are consistent with zero, but the presence of an azimuthal asymmetry at a few percent level is allowed.The p di-J/ψ T spectra are consistent with each other in different y di-J/ψ intervals, and peak at lower p di-J/ψ T values than the theoretical predictions.No significant broadening of the p T spectrum with increasing m di-J/ψ is seen in the TMD region, but the average value of the p di-J/ψ T distribution increases slightly with m di-J/ψ .The results provide important experimental inputs to study gluon TMDs and to improve theoretical models.The ongoing data taking by the LHCb experiment will enable more precise measurements of production cross-sections, and allow to test the predictions of gluon TMDs with higher precision.

Figure 2 :
Figure 2: Distributions of (a) t J/ψ 1 z and (b) t J/ψ 2 z for di-J/ψ signals with backgrounds subtracted.The projections of the fit result are overlaid.

Figure 3 :
Figure 3: Differential cross-section of di-J/ψ production for SPS+DPS, DPS and SPS as a function of (a) ∆y, (b) ∆ϕ, (c) A p T , (d) p di-J/ψ T , (e) y di-J/ψ , (f) m di-J/ψ , (g) p J/ψ T and (h) y J/ψ .The error bars represent the statistical and systematic uncertainties added in quadrature.The purple dashed lines in (a) and (c) indicate the baseline of zero.

Figure 4 :
Figure 4: Summary of σ eff measurements in pp or pp collisions [29-42].The legend entries marked with an asterisk are taken from a third-party calculation based on the original experimental result [28, 70-73].

Figure 5 :
Figure 5: Differential cross-section of di-J/ψ production for SPS+DPS and SPS as a function of (a) ∆y, (b) ∆ϕ, (c) A p T , (d) p di-J/ψ T

Figure 7 :
Figure 7: Distribution of ϕ CS (a) with the SPS and DPS contributions separated in the TMD region p di-J/ψ T < 4.1 GeV/c and (b) for SPS with the function described in the text overlaid.The systematic uncertainties correlated between intervals are excluded from the error bars.

Table 1 :
Summary of the systematic uncertainties on the measurement of the di-J/ψ production cross-section.The total systematic uncertainty is a quadratic sum of these uncertainties.

Table 2 :
Differential cross-sections dσ/d∆y of di-J/ψ production.The first uncertainties are statistical, and the second systematic.

Table 3 :
Differential cross-sections dσ/d∆ϕ of di-J/ψ production.The first uncertainties are statistical, and the second systematic.

Table 4 :
Differential cross-sections dσ/dA p T of di-J/ψ production.The first uncertainties are statistical, and the second systematic.

Table 6 :
Differential cross-sections dσ/dy di-J/ψ of di-J/ψ production.The first uncertainties are statistical, and the second systematic.

Table 7 :
Differential cross-sections dσ/dm di-J/ψ of di-J/ψ production.The first uncertainties are statistical, and the second systematic.

Table 9 :
Differential cross-sections dσ/dy J/ψ of di-J/ψ production.The first uncertainties are statistical, and the second systematic.

Table 10 :
Differential cross-sections dσ/d∆y of di-J/ψ production in SPS.The first uncertainties are statistical, and the second systematic.

Table 11 :
Differential cross-sections dσ/d∆ϕ of di-J/ψ production in SPS.The first uncertainties are statistical, and the second systematic.

Table 12 :
Differential cross-sections dσ/dA p T of di-J/ψ production in SPS.The first uncertainties are statistical, and the second systematic.

Table 13 :
Differential cross-sections dσ/dp di-J/ψ T of di-J/ψ production in SPS.The first uncertainties are statistical, and the second systematic.

Table 14 :
Differential cross-sections dσ/dy di-J/ψ of di-J/ψ production in SPS.The first uncertainties are statistical, and the second systematic.

Table 15 :
Differential cross-sections dσ/dm di-J/ψ of di-J/ψ production in SPS.The first uncertainties are statistical, and the second systematic.

Table 17 :
Differential cross-sections dσ/dy J/ψ of di-J/ψ production in SPS.The first uncertainties are statistical, and the second systematic.