Associated production of prompt $J/\psi$ and $\mathit{\Upsilon}$ mesons in $pp$ collisions at $\sqrt{s}=13\,\mathrm{TeV}$

The associated production of prompt $J/\psi$ and $\mathit{\mathit{\Upsilon}}$ mesons in $pp$ collisions at a centre-of-mass energy of $\sqrt{s}=13\,\mathrm{TeV}$ is studied using LHCb data, corresponding to an integrated luminosity of $4\,\mathrm{fb}^{-1}$. The measurement is performed for $J/\psi$ ($\mathit{\Upsilon}$) mesons with a transverse momentum $p_{\mathrm{T}}<10\,(30)\,\mathrm{GeV}/c$ in the rapidity range $2.0<y<4.5$. In this kinematic range, the cross-section of the associated production of prompt $J/\psi$ and $\mathit{\Upsilon}(1S)$ mesons is measured to be $133 \pm 22 \pm 7 \pm 3 \, \mathrm{pb}$, with a significance of $7.9\,\sigma$, and that of prompt $J/\psi$ and $\mathit{\Upsilon}(2S)$ mesons to be $76\pm 21 \pm 4 \pm 7 \, \mathrm{pb}$, with a significance of $4.9\,\sigma$. The first uncertainty is statistical, the second systematic, and the third due to uncertainties on the used branching fractions. This is the first observation of the associated production of $J/\psi$ and $\mathit{\Upsilon}(1S)$ in proton-proton collisions. Differential cross-sections are measured as functions of variables that are sensitive to kinematic correlations between the $J/\psi$ and $\mathit{\Upsilon}(1S)$ mesons. The effective cross-sections of the associated production of prompt $J/\psi$ and $\mathit{\Upsilon}$ mesons are obtained and found to be compatible with measurements using other particle productions.


Introduction
Hadroproduction of heavy quarkonia has been extensively studied to probe quantum chromodynamics (QCD) [1]. Various theoretical treatments have been proposed that model the production of quark-antiquark pairs QQ and their subsequent hadronisation into quarkonia. Promising approaches include the colour-singlet model (CSM) [2][3][4] and the non-relativistic QCD (NRQCD) framework [5]. The CSM considers only the coloursinglet intermediate QQ state with the same J P C quantum numbers as the final-state quarkonium, while the NRQCD framework considers both colour-singlet and colouroctet QQ states. The NRQCD framework assumes process-independent long-distance matrix elements (LDMEs) to describe the transition probabilities from the intermediate QQ states to the final-state quarkonium. Quarkonium hadroproduction cross-sections as a function of the transverse momentum (p T ) are successfully described by NRQCD calculations [6][7][8][9][10][11]. However, there are still challenges in providing a coherent description of all heavy quarkonium hadroproduction data, e.g. there is no unified description of the cross-sections and polarisation of the vector quarkonium states; and the inconsistency between the LDMEs of the η c and J/ψ states-which are related to the heavy-quark spin symmetry-remains to be understood [12][13][14][15][16][17][18][19].
In hadron collisions, a pair of heavy quarkonia can be produced through either singleparton scattering (SPS), or double-parton scattering (DPS), depending on the number of parton-parton scatterings in the process. A representative Feynman diagram of each process is presented in Fig. 1. In the SPS process, two quarkonia are produced in the same parton-parton interaction, leading to kinematic correlations between them. In contrast, two parton-parton interactions occur simultaneously and independently in the DPS process, producing two quarkonia with weak kinematic correlations. While the SPS process provides new insights into the heavy quarkonium production mechanism according to the left Feynman diagram in Fig. 1 [20][21][22][23], the DPS process probes the parton distribution function inside the colliding hadrons and provides valuable information on the hadronic wave functions describing correlations among partons [1]. Taking the J/ψ and Υ associated production (referred to as J/ψ-Υ production throughout) as an example 1 , the DPS cross-section σ DPS (J/ψ-Υ ) can be estimated from the J/ψ and Υ cross-sections, σ(J/ψ) and σ(Υ ), as where σ eff is the effective cross-section parameter encoding the transverse overlapping region between partons in the protons [24][25][26]. The effective cross-section is expected to be universal based on the assumptions that the two sets of colliding partons are uncorrelated and the longitudinal and transverse components of the parton distribution function factorise. Measurements in different processes for the effective cross-section are mostly within the range of 5 to 35 mb [27,28]. Further study is needed to test its universality [29]. A range of associated production processes of heavy quarkonia has been studied experimentally [20,28,[30][31][32][33][34][35][36][37][38][39][40][41][42]. Heavy-quarkonium-pair production was first observed by the NA3 collaboration in the study of J/ψ-pair production in π − -induced interactions [ In this paper, J/ψ-Υ production in pp collisions at a centre-of-mass energy of √ s = 13 TeV is studied for J/ψ (Υ ) mesons with a transverse momentum p T < 10 (30) GeV/c in the rapidity range 2.0 < y < 4.5 using data collected by LHCb in the years 2016, 2017, and 2018. The data correspond to an integrated luminosity of L = 4.18 ± 0.08 fb −1 . The J/ψ and Υ mesons are reconstructed in the dimuon final state. The polarisation of the J/ψ and Υ mesons is assumed to be zero. The polarisation of J/ψ-Υ pairs is not yet measured, but all the LHC measurements so far indicate a small polarisation for quarkonia [52][53][54][55][56][57].
The online event selection is performed by a trigger [60], which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. For this analysis, data were collected with varying trigger conditions particularly for the hardware stage, which requires two muons with p T (µ 1 ) × p T (µ 2 ) > (1.3) 2 , (1.5) 2 or (1.8) 2 GeV 2 /c 2 . In the software trigger stage, at least two muons with p T > 300 MeV/c, p > 6 GeV/c, and dimuon invariant mass m(µ + µ − ) > 2.7 GeV/c 2 are required. Candidate J/ψ and Υ mesons are built from two oppositely charged tracks that are identified as muons and form a good vertex. These muons from the J/ψ (Υ ) decays are required to be in the kinematic region 2 < η < 5, p T > 650 (3000) MeV/c, and p > 6 GeV/c. The decay vertices of the J/ψ and Υ mesons have to be consistent with the same PV. The invariant mass of the J/ψ candidates is required to be within [3.0, 3.2] GeV/c 2 and that of the Υ candidates within [9,11] GeV/c 2 .
Two independent simulated samples of J/ψ and Υ decays are used to model effects of the detector acceptance and the imposed selection requirements for J/ψ and Υ decays respectively. Primary pp collisions are generated using Pythia 8 [61] with a specific LHCb configuration [62]. Decays of unstable particles are simulated using EvtGen [63], in which the final-state radiation is generated with PHOTOS [64]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [65] as described in Ref. [66]. All simulated events are reconstructed and selected in the same way as the data.

Cross-section determination
The J/ψ-Υ production cross-section is measured according to where N cor is the efficiency-corrected signal yield of the J/ψ-Υ associated production, B J/ψ→µ + µ − and B Υ →µ + µ − are the branching fractions of the J/ψ → µ + µ − and Υ → µ + µ − decays [67], respectively, and L = 4.18 ± 0.08 fb −1 is the integrated luminosity determined with the van-der-Meer-scan method [68,69]. The total detection efficiency ϵ tot for the offline selected J/ψ-Υ candidates is factorised as where ϵ acc is the geometrical acceptance, ϵ rec&sel is the reconstruction and selection (not including PID requirements) efficiency of all the candidates within the geometrical acceptance, ϵ PID is the PID selection efficiency of candidates reconstructed and passing other selection requirements, and ϵ trig is the trigger selection efficiency of candidates passing all selection requirements. The first three terms on the right-hand side of Eq. 3 are computed by multiplying the individual efficiencies for J/ψ and Υ candidates, i.e., where i denotes acc, rec&sel, or PID. For the trigger selection, it is required that either the J/ψ or Υ candidate meets the trigger requirements specified in Sec. 2; therefore the trigger efficiency ϵ trig can be factorised as  Each efficiency component is determined separately for J/ψ and Υ mesons as a function of their p T and rapidity y. The acceptances, the reconstruction and selection efficiencies, and the trigger efficiencies are estimated using simulation, while the PID efficiencies are obtained using a data-driven method [70].
To extract the signal yield, an unbinned extended maximum-likelihood fit is performed to the two-dimensional J/ψ-Υ invariant-mass distribution. Defining the signal distribution to be S and combinatorial-background distribution B, the components in the model used to fit the data fall into three categories, each corresponding to the following cases: • a true J/ψ (Υ ) candidate associated with a combinatorial Υ (J/ψ) background candidate (S J/ψ × B Υ and B J/ψ × S Υ ), • and a purely combinatorial background (B J/ψ × B Υ ).
The two-dimensional probability density function for each component factorises into the product of two one-dimensional distributions for the J/ψ and Υ candidates. The signal invariant-mass distribution of the J/ψ meson and that of each Υ meson are described by a combination of a double-sided Crystal Ball (DSCB) function [71] and a Gaussian function with a shared mean. The DSCB power-law tail parameters, the ratio of the width of the Gaussian function to that of the Gaussian core in the DSCB function, and the relative fraction of the DSCB to the Gaussian function are fixed to values estimated from simulation. The ratios between the means of the Υ Gaussian functions are fixed to the ratios of the known Υ masses. As the mass resolution of signals is linearly dependent on invariant mass, the width ratios of the Gaussian functions are also fixed to the known mass ratios [67]. The mean value and the width of the Gaussian function for the Υ (1S) are free parameters in the fit.
Each combinatorial background is modelled by an exponential function. The invariantmass distributions of the J/ψ and Υ candidates, together with the fit projections, are shown in Fig. 2. Associated J/ψ-Υ (1S) production is observed with a significance of 7.9 σ, and evidence for J/ψ-Υ (2S) production with a significance of 4.9 σ is seen. The total significance for the associated production of J/ψ with a combination of Υ (1S), Υ (2S), and Υ (3S) states is determined to be 8.5 σ. The significance is estimated by comparing the likelihood of the fit with and without the considered component in the fit, with relevant systematic uncertainties taken into consideration. The raw yields of signals from the fit and their significances are summarised in Table 1.
Based on the fit result, a signal weight w i is calculated with the sPlot method [72] for each J/ψ-Υ candidate, using the J/ψ and Υ invariant masses as discriminating variables. These weights are used to determine the efficiency corrected yield N cor = j w j /ϵ tot, j , where j runs over all the J/ψ-Υ candidates in data. The total efficiency for the j th candidate, ϵ tot, j , is defined in Eq. 3 and determined as a function of the p T and y of the J/ψ and Υ candidate. The corrected yield, N cor , is shown in Table 1.
The signal extracted from the fit also contains J/ψ mesons produced from b-hadron decays (denoted J/ψ-from-b). The fraction of J/ψ-from-b decays is determined to be 2.1% using the prompt J/ψ and J/ψ-from-b production cross-sections measured by LHCb in pp collisions at √ s = 13 TeV [73]. The difference between the average reconstruction and selection efficiencies of the two components is evaluated from simulation; the selection criteria constraining J/ψ and Υ decay vertices to the PV favour the prompt J/ψ candidates and suppress the J/ψ-from-b candidates. With current measurements of σ eff , the J/ψ-Υ associated production is predicted to be dominated by the DPS process. It is therefore assumed that in events with J/ψ and Υ candidates the relative cross-section of the J/ψ-from-b production with respect to the prompt J/ψ production is the same as the relative cross-section of J/ψ-from-b in inclusive J/ψ events. An uncertainty is assigned for systematic effects due to this assumption. The J/ψ-from-b component is subtracted when calculating the total cross-section of the associated J/ψ-Υ production.

Systematic uncertainties
Systematic uncertainties from various sources are studied and summarised in Table 2. The total relative systematic uncertainty is 5.5% for J/ψ-Υ (1S) and 10.3% for J/ψ-Υ (2S) production, corresponding to the quadratic sum of all the contributions discussed below.
The uncertainty introduced by the imperfect modelling of the signal shape is estimated by fitting the data with an alternative signal model. The alternative shape is the shape of the simulated signal distribution, convolved with a Gaussian function to account for the different invariant-mass resolution observed between data and simulation. The 2.0% relative change of the signal yield is assigned as the systematic uncertainty for both  The track reconstruction efficiency as a function of muon p and η is estimated with J/ψ → µ + µ − decays in a calibration data sample using a tag-and-probe technique [74]. The difference between data and simulation is used to calibrate the reconstruction efficiency estimated using simulated signal decays. The uncertainties on these calibrations are propagated to the efficiency-corrected yield using a pseudoexperiment Monte Carlo method. The uncertainty is 0.9% for both the J/ψ-Υ (1S) and J/ψ-Υ (2S) signals. An additional uncertainty of 0.8% per track is assigned to account for the difference in detector occupancy between data and simulation [74]. A total relative uncertainty of 3.3% is attributed to the tracking efficiency calibration.
The muon PID efficiency is calculated with J/ψ → µ + µ − decays using a data-driven tag-and-probe method [70]. The efficiency is estimated in bins of p T and η of the probe track, for which the binning scheme introduces an uncertainty. Alternative binning schemes are used to estimate variations of the PID efficiency. The relative change with respect to the default value, 0.5%, is quoted as the systematic uncertainty.
The trigger efficiency estimated from simulation is compared with that from data with a data-driven method [46,75]. A relative difference of 1.0% is quoted as the systematic uncertainty.
The uncertainties introduced by the mismodelling of the radiative decay tails and the requirement on the vertex fit quality is estimated to be 2% and 0.6%, respectively, following Ref. [73,76] by comparing data and simulation.
The fraction of J/ψ-from-b is calculated from the p T -and y-dependent J/ψ-from-b fraction in single J/ψ production and the relative detection efficiency between J/ψ-from-b and prompt production, convolved with the J/ψ kinematic distributions in simulation. The uncertainty of J/ψ kinematic distributions used to estimate the fraction and the possible difference of J/ψ-from-b fraction between inclusive J/ψ production and associated J/ψ-Υ production are sources of systematic uncertainties. The first uncertainty is evaluated using the maximum and minimum values of the inclusive J/ψ-from-b fraction in the studied kinematic range of this analysis. The second source is estimated conservatively by varying the effective cross-section σ eff of the Υ -B hadron associated production by 50%. In total, this systematic uncertainty is estimated to be 1.1%. The effect of events where J/ψ and Υ mesons originate from different PVs is estimated to be negligible.
The ratio of the cross-sections of J/ψ-Υ (2S) and J/ψ-Υ (1S) multiplied by the respective branching ratio is calculated to be where the first uncertainty is statistical and the second systematic. The result is consistent with the measurement of Υ (1S, 2S) production [76] and of the production of Υ (1S, 2S) mesons associated with an open charm hadron [38], as is expected by the DPS assumption. Using the DPS cross-section of the J/ψ-Υ associated production and the promptproduction cross-sections of single J/ψ and Υ mesons, the effective cross-section is calculated according to Eq. 1. The DPS cross-sections are estimated by subtracting the SPS cross-sections, 20 +52 −15 pb for J/ψ-Υ (1S) and 8 +22 −6 pb for J/ψ-Υ (2S), from the associatedproduction cross-sections. The SPS cross-sections are adapted from the calculations with the NRQCD framework in Ref. [20]. The prompt-production cross-sections of J/ψ and Υ are taken from Ref. [73] and Ref. [76], respectively. The effective cross-section is σ eff (J/ψ-Υ (1S)) = 26 ± 5 ± 2 +22 −3 mb, σ eff (J/ψ-Υ (2S)) = 14 ± 5 ± 1 +7 −1 mb, where the first uncertainties are statistical, the second uncertainties are systematic, and the last terms are uncertainties due to theory calculations. The results of σ eff are broadly compatible with measurements using hadroproduction of other particles, as shown in Fig. 3. The effective cross-sections for J/ψ-Υ (1S) and J/ψ-Υ (2S) production are consistent.
The distributions of some observables related to both J/ψ and Υ (1S) mesons are studied to probe the kinematic correlation between the two mesons. These variables include the J/ψ-Υ (1S) invariant mass, the transverse momentum of the J/ψ-Υ (1S) system, the relative azimuthal angle, and the rapidity difference between J/ψ and Υ (1S) mesons. The distributions of these variables and the transverse momentum distributions of the J/ψ and Υ (1S) mesons are presented in Fig. 4, after correcting for the total detection efficiency and subtracting the background contribution. The large uncertainties in the correlation distributions are due to the small size of the signal data sample, which dominate over systematic uncertainties. These distributions are compared with pseudoexperiments corresponding to different contributions of the DPS and SPS processes. Since the J/ψ and Υ mesons are produced independently from two different parton-parton scatterings in the DPS processes, the kinematic distribution of each meson in the associated production  [73,76]. Black points with error bars represent the data distributions with statistical uncertainties. Red solid lines represent distributions of the DPS pseudodata. Green dashed lines and blue dotted lines represent the DPS distribution together with the lower limit and the upper limit of the SPS distributions, respectively. is expected to be similar to that in the single meson production. As a result, using the measured p T and y distributions in single-J/ψ and -Υ (1S) production [73,76], the DPS distributions of J/ψ and Υ (1S) kinematic variables in the associated production can be modelled. The distributions for the SPS production are taken from calculations in Ref. [20]. Three samples of pseudodata are generated, differing in the relative SPS to DPS fractions. For one sample, only DPS production is considered. For the other two samples, the SPS fraction is calculated using either the upper or lower limit of the predicted cross-section in Ref. [20], and the DPS cross-section is the LHCb measurement of the total cross-section with the SPS cross-section subtracted. As shown in Fig. 4, the J/ψ and Υ (1S) p T distributions and the kinematic correlations between the J/ψ and Υ (1S) mesons in the three pseudodata samples are all found to be consistent with measurements.

Summary
Using pp collision data collected by the LHCb experiment at a centre-of-mass energy of 13 TeV corresponding to an integrated luminosity of 4.18 ± 0.08 fb −1 , a signal for J/ψ-Υ (1S) production is observed with a significance of 7.9 σ, and evidence is found for J/ψ-Υ (2S) production. Within the fiducial region 2.0 < y < 4.5 and p T < 10 (30) GeV/c for J/ψ (Υ ) mesons, the cross-section for J/ψ-Υ (1S) production is measured to be 133 ± 22 (stat) ± 7 (syst) ± 3 (B) pb and that for J/ψ-Υ (2S) production 76 ± 21 (stat) ± 4 (syst) ± 7 (B) pb. The effective cross-sections for J/ψ-Υ (1S) and J/ψ-Υ (2S) production are calculated and found to be compatible with previous measurements in various processes. Kinematic correlations between the J/ψ and Υ (1S) mesons are studied, showing results consistent with both the hypotheses that the associated production is dominated by the DPS process or has DPS and SPS contributions. In the future, further data from LHCb combined with improved theory calculations will help constrain the associated quarkonium production mechanism.