Measurement of associated J/ψ - ψ (2 S ) production cross-section in pp collisions at √ s = 13 TeV

The cross-section of associated J/ψ - ψ (2 S ) production in proton-proton collisions at a centre-of-mass energy of √ s = 13 TeV is measured using a data sample corresponding to an integrated luminosity of 4 . 2 fb − 1 , collected by the LHCb experiment. The measurement is performed for both J/ψ and ψ (2 S ) mesons having transverse momentum p T < 14 GeV/ c and rapidity 2 . 0 < y < 4 . 5, assuming negligible polarisation of the J/ψ and ψ (2 S ) mesons. The production cross-section is measured to be 4 . 5 ± 0 . 7 ± 0 . 3 nb, where the first uncertainty is statistical and the second systematic. The differential cross-sections are measured as functions of several kinematic variables of the J/ψ - ψ (2 S ) candidates. The results are combined with a measurement of J/ψ - J/ψ production, giving a cross-section ratio between J/ψ - ψ (2 S ) and J/ψ - J/ψ production of 0 . 274 ± 0 . 044 ± 0 . 008, where the first uncertainty is statistical and the second systematic.


Introduction
Quantum chromodynamics (QCD) is the fundamental theory of the strong interaction between quarks and gluons.One of the most important properties of QCD is that the coupling constant increases with decreasing energy.A perturbative treatment can be applied at large momentum transfer, while the non-perturbative effects at small momentum transfer make theoretical calculations difficult and usually need to be described by effective theories or phenomenological models with input from experiment.The study of heavy quarkonium production in proton-proton (pp) collisions can provide important information to improve QCD calculations in the non-perturbative regime.The process involves production of a QQ pair, where Q denotes a beauty or charm quark, followed by its hadronisation into a heavy quarkonium state.In the nonrelativistic QCD (NRQCD) approach [1][2][3], both colour-singlet and colour-octet states of the intermediate QQ pair are taken into account, and the transition probabilities from the QQ states to the heavy quarkonium are described by non-perturbative long-distance matrix elements.Although the NRQCD approach describes many aspects of quarkonium production successfully, no satisfactory solution has been achieved thus far to simultaneously describe differential cross-section and polarisation data over the whole kinematic region [4,5].This puzzle can be probed via the production of heavy quarkonium pairs through the single-parton scattering (SPS) process [6][7][8][9], which can provide new tests of the NRQCD approach.
Besides SPS, heavy quarkonium pairs can be produced through double-parton scattering (DPS) [10], which makes the study more complicated.The DPS process is of great relevance since it provides valuable information on the profile and correlations of partons inside the proton.A thorough understanding of DPS is crucial in some searches for physics beyond the standard model because several important sources of background, e.g.Z + bb, W ± W ± , have significant DPS contributions.In DPS, the heavy quarkonium pair production is usually assumed to come from two independent partonic interactions, and the cross-section can be estimated according to [11][12][13] where δ Q 1 Q 2 is one if the two quarkonium states Q 1 and Q 2 are the same and zero if they are different, σ Q 1 and σ Q 2 are the production cross-sections of single quarkonium Q 1 and Q 2 respectively, and σ eff is an effective cross-section characterising an effective transverse overlap area between the two pairs of interacting partons.Experimentally, it is difficult to distinguish the SPS and DPS contributions.The measurements of quarkonium pair production so far mostly separate them according to the kinematic relation between the two quarkonium decays, which is model dependent.This paper focuses on charmonium pair production.The production of J/ψ mesons associated with the χ c state, whose charge parity is even, through SPS is forbidden at leading order in NRQCD, resulting in a large difference in the fraction of J/ψ pairs produced from feed-down of χ c and ψ(2S) mesons between SPS and DPS.The fraction of feed-down from ψ(2S) production is expected to be large (around 46%) in SPS, while it is only about 20% in DPS [10].The fraction of feed-down from χ c production is suppressed in SPS, while it is around 50% in DPS [10].Thus the feed-down fraction of J/ψ pairs from J/ψ-ψ(2S) or J/ψ-χ c is considered as another important test to pin down the DPS or SPS dominance.
The first evidence for J/ψ-pair production in hadron collisions was reported by the NA3 collaboration in the 1980s [14,15].In recent years, the J/ψ-pair production crosssection in pp collisions was measured by the LHCb collaboration at √ s = 7 TeV [16] and 13 TeV [17], by the CMS collaboration in pp collisions at √ s = 7 TeV [18] and by the ATLAS collaboration at √ s = 8 TeV [19].It was also measured in pp collisions by the D0 collaboration at √ s = 1.96TeV [20].However, the production of charmonium pairs involving other states, such as J/ψ-ψ(2S) and J/ψ-χ c has not been measured before.
In this paper, the measurement of J/ψ-ψ(2S) production cross-sections in pp collisions at √ s = 13 TeV is presented, 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 production cross-sections are measured for J/ψ and ψ(2S) mesons having transverse momentum p T < 14 GeV/c and rapidity 2.0 < y < 4.5, assuming negligible polarisation of the J/ψ and ψ(2S) mesons, since all the LHC measurements so far indicate that the polarisation of quarkonia is small [21][22][23][24][25].The differential cross-sections of J/ψ-ψ(2S) production as functions of several kinematic variables are also measured, as well as the cross-section ratio between J/ψ-ψ(2S) and J/ψ-J/ψ production.

Detector and dataset
The LHCb detector [26,27] 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 pp interaction point (i.e.primary 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.
The J/ψ-ψ(2S) candidates are reconstructed using final states with two µ + µ − pairs.The online event selection is performed by a trigger [28], which consists of a hardware stage (L0), based on information from the calorimeter and muon systems, followed by two stages of software trigger, HLT1 and HLT2, at which full event reconstruction is performed.It is required that at least one µ + µ − pair fulfils the selection criteria of the L0 and HLT1 triggers.The L0 trigger selects two muons with the product of their transverse momentum p T 1 × p T 2 > 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, which are loosely identified as muons, to form a J/ψ or ψ(2S) candidate with high 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 J/ψ and ψ(2S) candidates to be reconstructed with good vertex-fit quality, and the invariant mass m µ + µ − is required to be within 120 MeV/c 2 of the known J/ψ or ψ(2S) masses [29].
Simulated samples are produced to study the behaviour of signal candidates and determine the detection efficiencies.The pp collisions are modelled using Pythia [30,31] with a specific LHCb configuration [32].Decays of unstable particles are described by EvtGen [33] with QED final-state radiation handled by Photos [34].The interactions of the generated particles with the detector are modelled using the Geant4 toolkit [35,36] as described in Ref. [37].Due to limited knowledge of the associated J/ψ-ψ(2S) production mechanism, the candidate selection is designed not to utilise information based on the correlation between the J/ψ and ψ(2S) mesons.Therefore, the detection efficiency can be factorised into that of J/ψ and ψ(2S).The efficiency factorisation is validated using a simulated sample of Υ → J/ψJ/ψγ decays.The J/ψ and ψ(2S) efficiencies are evaluated using simulated samples of single-J/ψ and single-ψ(2S) candidates.In the Pythia model, J/ψ and ψ(2S) mesons are generated with no polarisation, and the leading order coloursinglet and colour-octet contributions [32,38] are considered in prompt J/ψ and ψ(2S) production.

Candidate selection
Further selections are performed offline to J/ψ-ψ(2S) candidates to reduce the combinatorial background.Four muon tracks with good track-fit quality are selected and required to have 1.9 < η < 4.9, p T > 0.65 GeV/c and p > 3 GeV/c, and to form two good-quality di-muon vertices.The muon identification requirement is tightened to further suppress the contribution of mis-identified tracks.The background from fake tracks is reduced by a neural-network based algorithm [39].All events are required to have at least one reconstructed PV.For candidates with multiple PVs in the event, the one with the smallest χ 2 IP is taken as the associated PV, where χ 2 IP is defined as the difference in the vertex-fit χ 2 of a given PV reconstructed with and without the candidate under consideration.The four muon tracks are required to originate from the same PV.This reduces the number of pile-up candidates, i.e.J/ψ and ψ(2S) candidates originating from two independent pp interactions, to a negligible level.A selection is performed on the pseudoproper time t z [40] of the J/ψ and ψ(2S) candidates, defined as where z ψ and z PV are the coordinates of the J/ψ or ψ(2S) decay vertex and the PV along the beam axis z, p z is the projection of the J/ψ or ψ(2S) momentum along the z axis, and m ψ is the known J/ψ or ψ(2S) mass [29].The t z uncertainty σ tz is calculated by combining the uncertainties from z ψ and z PV since the uncertainty on p z is negligible in comparison.Candidates with both J/ψ and ψ(2S) mesons having −2 < t z < 10 ps and σ tz < 0.3 ps are selected.After the selection, the events with multiple candidates, in particular those in which the four muons can be combined in two different ways to form a J/ψ-ψ(2S) candidate, account for 1.2% of the total.According to simulation, the candidates with wrong combination of muons follow the background distribution, so their inclusion do not introduce any bias to the signal yield determination.

Cross-section determination
The cross-section of associated J/ψ-ψ(2S) production is given by where N corr is the efficiency-corrected yield, L is the integrated luminosity, and B is the appropriate branching fraction.The integrated luminosity is L = 4.18 ± 0.08 fb −1 measured using the van der Meer scan method [41].The branching fraction of the J/ψ decay is B(J/ψ → µ + µ − ) = (5.961± 0.033)% [29].Under the assumption of lepton universality in electromagnetic decays, the branching fraction of the ψ(2S) decay is taken to be B(ψ(2S) → e + e − ) = (7.93 ± 0.17) × 10 −3 [29], taking advantage of its much smaller uncertainty compared to the muonic decay.
The J/ψ-ψ(2S) signals are extracted by performing an unbinned extended maximum likelihood fit to the two-dimensional (2D) mass distribution of J/ψ and ψ(2S) candidates, . The 2D mass distribution function can be factorised into the product of the two mass dimensions.Each one-dimensional mass distribution consists of the signal and combinatorial background components.The signal component is described by a sum of a double-sided Crystal Ball (DSCB) function [42] and a Gaussian function with a common mean value but different widths.All the tail parameters, 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 shape parameters.The combinatorial background distribution is modelled with an exponential function.
Figure 1 shows the projections of the 2D mass distribution, together with the fit result, of the m µ + 1 µ − 1 and m µ + 2 µ − 2 invariant mass combinations.In pp collisions, J/ψ or ψ(2S) mesons can be produced either directly from hard collisions of partons, through the feed-down of excited charmonium states, or via decays of beauty hadrons.The J/ψ or ψ(2S) mesons from the first two sources originate from the PV and are called prompt mesons, while those from the last source originate from decay vertices of beauty hadrons, which are typically separated from the PV, and are called nonprompt mesons.The nonprompt contributions with the J/ψ or ψ(2S) meson originating from b-hadron decays are subtracted by performing a maximum-likelihood fit to the 2D t z distribution of J/ψ-ψ(2S) candidates with backgrounds subtracted.The background subtraction is implemented with the sPlot method [43] 2 as the discriminating variables.The t z distribution of prompt J/ψ or ψ(2S) mesons is described by the Dirac delta function δ(t z ), while that of nonprompt mesons by an exponential function.These distributions are convolved with a Gaussian detector resolution function.For a very small fraction of candidates, the true PV is not reconstructed and both J/ψ   and ψ(2S) candidates are associated to the nearest reconstructed PV.The positions of the true PV and the wrongly associated PV are not correlated, which results in a broad distribution of t z .This component can be modelled from data using event mixing, i.e. calculating t z with the J/ψ or ψ(2S) candidate associated to the closest PV in the next event of the sample.Figure 2 shows the projections of the 2D t z distribution with backgrounds subtracted, together with the fit result, on t J/ψ z and t Since the kinematics of J/ψ-ψ(2S) decay are not known a priori, an efficiency correction is performed event-by-event as where the index i denotes each event in data, w i is the sWeight [43] corresponding to the component of prompt J/ψ-ψ(2S) signal candidates, and ε tot i the detection efficiency.The efficiency ε tot i is evaluated as the product of four components: the geometrical acceptance ε acc i , the efficiency of reconstruction and selection ε rec&sel i , the PID efficiency ε PID i , and the trigger efficiency ε trig i , giving For a J/ψ-ψ(2S) candidate, each efficiency factor can be parameterised as a function of the p T and y of the J/ψ and ψ(2S) mesons.As long as the real p T and y values of both J/ψ and ψ(2S) mesons are taken, it makes no difference for the efficiencies whether the candidate is produced by SPS or DPS.The possible kinematic correlation between two mesons, reflected in p T and y, carries information on the production mechanism.On this basis, the efficiencies of each J/ψ-ψ(2S) candidate can be factorised into those of the J/ψ and ψ(2S) mesons separately, since no information related to the correlation between them is used during the reconstruction and selection process.The efficiencies ε acc i , ε rec&sel i and ε PID i of each J/ψ-ψ(2S) candidate are factorised as the product of efficiencies of the J/ψ and ψ(2S) mesons as The trigger efficiency ε trig i of each J/ψ-ψ(2S) candidate is factorised as where ε L0&HLT1 i is the L0-and-HLT1 trigger efficiency.The online HLT2 reconstruction algorithm is equivalent to the offline algorithm and the offline selections are tighter than the HLT2 requirements, making the HLT2 trigger fully efficient with respect to offline selected candidates.All the efficiency terms are estimated in (p T , y) intervals of the J/ψ (ψ(2S)) meson using the simulated single-J/ψ (single-ψ(2S)) samples.The track reconstruction efficiency, which is part of ε rec&sel i , and the PID efficiency, ε PID i , are calibrated using data to eliminate differences between the simulation and data [44,45].The efficiency-corrected signal yield is determined to be N corr = (8.9± 1.4) × 10 3 , where the statistical uncertainty is verified by the bootstrapping approach [46].

Systematic uncertainties
Table 1 summaries the systematic uncertainties on the cross-section measurement.Due to the small yields of the J/ψ-ψ(2S) data, a larger pseudo J/ψ-ψ(2S) sample is produced by combining the single-J/ψ and single-ψ(2S) candidates in data randomly to reduce the statistical fluctuations in the evaluation of systematic uncertainties.The same treatment of the selection, signal extraction and the efficiency correction is applied to this pseudo sample as for the J/ψ-ψ(2S) data sample.
An uncertainty is attributed to the imperfect modelling of the J/ψ (ψ(2S)) mass distribution.As an alternative to the sum of a DSCB function and a Gaussian function, the invariant-mass distribution of J/ψ (ψ(2S)) signals is described by a model derived from simulation using the kernel density estimation approach [47].To account for the resolution difference between data and simulation, the alternative model is convolved with a Gaussian function with a mean value of zero and the width varied freely.The relative difference of the extracted signal yields using the default and alternative models is 1.7%, which is taken as the systematic uncertainty.For the measurement of differential cross-sections, this uncertainty is taken to be common to all kinematic intervals.An uncertainty is attributed to the subtraction of nonprompt contributions.The quantity log(χ 2 IP ) of the J/ψ (ψ(2S)) meson is used for the discrimination of nonprompt components instead of the pseudoproper time t z .The log(χ 2 IP ) distribution of the prompt J/ψ (ψ(2S)) signals is described by the sum of two Bukin functions [48] with the same peak position and asymmetry, while the nonprompt component is described by one Bukin function.The shape parameters of the Bukin functions other than the peak position and width are fixed from simulation.The relative difference of the yield of prompt J/ψ-ψ(2S) signal candidates with respect to the default result, 1.9%, is taken as the systematic uncertainty.For the measurement of differential cross-sections, this uncertainty is considered common to all kinematic intervals.
A small fraction of J/ψ or ψ(2S) candidates in the J/ψ-ψ(2S) data may be associated to a wrong PV, which can be divided into two cases.The first is that the true PV is reconstructed but one of the J/ψ and ψ(2S) candidates is associated to a wrong PV accidentally.In this case, J/ψ and ψ(2S) mesons are associated to different PVs and thus are rejected by the selection.The fraction is estimated using the simulated sample of Υ → J/ψJ/ψγ decays to be (0.65 ± 0.02)%.The second is that the true PV is not reconstructed and both J/ψ and ψ(2S) candidates are associated to the nearest reconstructed PV in this event.This is a component in the 2D t z fit, with the fraction determined to be (0.8 ± 1.5)%.The total fraction of wrong PV events adds up to 1.5%, which is taken as a systematic uncertainty.For the measurement of differential crosssections, the fraction is studied in several intervals of J/ψ-ψ(2S) rapidities using the same method.The maximum, 2.2%, is conservatively taken as the uncertainty common to all kinematic intervals.
The limited size of the calibration samples used to determine the efficiencies can introduce a systematic uncertainty.This is studied by using pseudoexperiments, in which the efficiency in each (p T , y) interval is sampled by Gaussian distributions with the corresponding efficiency as the mean value and the statistical uncertainty of the efficiency as the width.The width of the cross-section distribution obtained from pseudoexperiments is determined to be 0.7%, and taken as the systematic uncertainty.It varies up to 3.0% depending on the kinematic intervals for the measurement of 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 studied by varying the binning schemes of the calibration sample, and the relative difference in the cross-sections is taken as an uncertainty.For the other efficiencies, ε acc , ε rec&sel and ε trig , the kernel density estimation approach [47] is used as an alternative, in which the efficiencies are determined as smoothed functions of (p T , y).The relative difference in the cross-sections between the default and alternative approaches is quoted as the systematic uncertainty separately for each efficiency factor.The uncertainty on the total efficiency is taken as a quadratic sum of these uncertainties on the four efficiency factors, and its value is 2.8%.For the differential cross-sections, it varies up to 5.5% depending on the kinematic intervals.
The systematic uncertainty from the track reconstruction efficiency has two sources.One is propagated from the statistical uncertainties of the correction factors due to the limited size of the calibration sample.It is evaluated to be 0.9% using pseudoexperiments.The correction factors also depend on the event multiplicity, which is not perfectly matched between the simulation and calibration samples.This introduces a systematic uncertainty of 0.8% per track, adding up to 3.2% for the J/ψ-ψ(2S) signals containing four muon tracks.The track detection efficiency uncertainty is in total 3.3% with the two sources added in quadrature.For the differential cross-sections, the first term varies up to 4.2%, while the second is considered to be common to all kinematic intervals.The systematic uncertainty associated with the trigger efficiency is evaluated by comparing the efficiency determined from the simulated single-J/ψ sample to that from the single-J/ψ data.The trigger efficiency is determined using a subset of events that fulfil the trigger requirement with the J/ψ signals excluded [49].Due to the small yields of the ψ(2S) sample, the ratio of trigger efficiencies between data and simulation determined using the J/ψ events is also applied to ψ(2S) mesons.The relative difference in the J/ψ-ψ(2S) cross-section between using the trigger efficiency obtained from data and simulation is 0.7%, and is taken as a systematic uncertainty.For the differential cross-sections, the quoted uncertainty varies up to 7.5% depending on the kinematic intervals.
The uncertainties on the branching fractions lead to an uncertainty of 2.2% on the measured cross-section.The luminosity is determined with a relative uncertainty of 2.0% using the methods described in Ref. [41].With these uncertainties due to independent effects added in quadrature, the total systematic uncertainty on the J/ψ-ψ(2S) production cross-section is determined to be 5.9%.

Differential cross-sections
The differential cross-section of J/ψ-ψ(2S) production is measured as a function of the absolute value of the rapidity difference between the J/ψ and ψ(2S) mesons ∆y, the absolute value of the difference in the azimuthal angle ϕ between the J/ψ and ψ(2S) mesons ∆ϕ, the transverse momentum p J/ψ-ψ(2S) T , rapidity y J/ψ-ψ(2S) and invariant mass m J/ψ-ψ(2S) of the J/ψ-ψ(2S) candidates.The binning scheme is chosen to have adequate and approximately equal signal yields in each bin.Due to the limited signal yields, the differential cross-sections of J/ψ-ψ(2S) production are calculated using the sWeights based on the global mass fit and global t z fit rather than from the fit in each kinematic interval.The potential correlation between these kinematic variables and the discriminating variables used in the sPlot method may bias the signal yields.Correction factors are thus obtained using a pseudo J/ψ-ψ(2S) sample produced by combining single-J/ψ and single-ψ(2S) candidates in data randomly.They are calculated as the ratio of yields obtained from the fit in each kinematic interval and from the global sWeights.The uncertainties on the correction factors come from the limited size of the pseudosample, resulting in another systematic uncertainty on the differential cross-sections of J/ψ-ψ(2S) production.It varies up to 10.7% depending on the kinematic intervals.The other systematic uncertainties on differential cross-sections are presented in Section 5.The measured differential cross-sections of J/ψ-ψ(2S) production are shown in Figure 3, and listed in Tables 2-6 in Appendix A.
In the NRQCD approach, explicit calculations [7] indicate the colour-octet contributions are much smaller than colour-singlet contributions.The incomplete (not including loop diagrams) next-to-leading order colour-singlet (NLO* CS) predictions [54] of J/ψψ(2S) production in the SPS process can be obtained from HELAC-Onia [55,56], an automatic matrix element generator for heavy quarkonium physics.The NLO* CS predictions are compared with the measured J/ψ-ψ(2S) production cross-sections as shown in Figure 3.The theoretical uncertainties include the uncertainty due to factorisation and renormalisation scales, which is dominant, and the PDF uncertainty.The data and theory  , (d) y J/ψ-ψ(2S) and (e) m J/ψ-ψ(2S) , compared with the NLO* CS predictions for SPS [54][55][56].The m J/ψ-ψ(2S) spectrum is also compared with the PRA+NRQCD predictions for SPS [57].The data contain contributions from both SPS and DPS.
are consistent within the large theoretical uncertainties, albeit the DPS contribution is not subtracted from the measurements.
As shown in Figure 3(e), the m J/ψ-ψ(2S) spectrum is also compared with the SPS predictions combining the parton Reggeization approach (PRA) [58] and the NRQCD factorisation approach [57], which includes a subset of higher-order QCD corrections in an infrared-safe way without ad-hoc kinematic cuts.Only the scale uncertainties are considered in the predictions.The PRA+NRQCD predictions for SPS are larger than the SPS+DPS data at small m J/ψ-ψ(2S) , and consistent with them at large m J/ψ-ψ(2S) .6.2 Cross-section ratio between J/ψ-ψ(2S) and J/ψ-J/ψ production The cross-section ratio between J/ψ-ψ(2S) and J/ψ-J/ψ production is measured to be σ J/ψ-ψ(2S) σ J/ψ-J/ψ = 0.274 ± 0.044 (stat) ± 0.008 (syst), in which the J/ψ-J/ψ production measurement is performed under a similar strategy as this analysis, as described in Ref. [59].The systematic uncertainties due to the signal mass model, nonprompt contribution, wrong PV association, efficiency determination, multiplicity dependence of the track detection efficiency, trigger efficiency and luminosity measurement are thus assumed to be fully correlated.The uncertainties from the branching fractions cancel partially and the remaining uncertainty is 2.2% on the cross-section ratio.The uncertainties due to calibration sample sizes are considered uncorrelated.The statistical uncertainty on the cross-section ratio is dominant.The ratio is predicted to be 0.94 ± 0.30 for SPS [10], and 0.282 ± 0.027 for DPS, according to Eq. 1 and the measured single-J/ψ and single-ψ(2S) cross-sections [50][51][52] assuming that the σ eff is the same for these two processes.The measured value is much smaller than the SPS prediction and consistent with the DPS prediction, indicating a prominent DPS contribution to charmonium pair production.
The cross-section ratio is also measured as a function of ∆y, ∆ϕ, p di-ψ T and y di-ψ , as shown in Figure 4.The symbol di-ψ indicates J/ψ-J/ψ or J/ψ-ψ(2S).The results are also presented in Tables 7-10 in Appendix A. No clear dependence of the cross-section ratio on the kinematic variables is seen given the large uncertainties so far.

Conclusion
The cross-section of associated J/ψ-ψ(2S) production in pp collisions at √ s = 13 TeV is measured to be 4.5±0.7 (stat)±0.3(syst) nb for J/ψ and ψ(2S) mesons with p T < 14 GeV/c and 2.0 < y < 4.5, using a data sample corresponding to an integrated luminosity of 4.2 fb −1 collected by the LHCb experiment.The differential cross-sections of J/ψ-ψ(2S) production are measured as functions of ∆y, ∆ϕ, p J/ψ-ψ(2S) T , y J/ψ-ψ(2S) and m J/ψ-ψ(2S) , and compared to the NLO* CS predictions for the SPS process.The obtained results are consistent with the theory given the presence of large theoretical uncertainties, although the DPS contribution is not subtracted from the measurement.The PRA+NRQCD predictions for SPS are larger than the data at small m J/ψ-ψ(2S) and consistent with them at large m J/ψ-ψ(2S) .The cross-section ratio between J/ψ-ψ(2S) and J/ψ-J/ψ production is also measured, and found to be significantly smaller than the SPS prediction and consistent with the DPS expectation.Previously [17][18][19][20], a prominent DPS contribution to the charmonium pair production was established mainly from the study of the kinematic relation of the two charmonia states.This is confirmed in this paper in a novel way.

Table 1 :
Summary of relative systematic uncertainties on the measurement of J/ψ-ψ(2S) production cross-section.The total systematic uncertainty is the quadratic sum of the individual contributions.