Observation of triple J/$\psi$ meson production in proton-proton collisions at $\sqrt{s}$ = 13 TeV

Protons consist of three valence quarks, two up-quarks and one down-quark, held together by gluons and a sea of quark-antiquark pairs. Collectively, quarks and gluons are referred to as partons. In a proton-proton collision, typically only one parton of each proton undergoes a hard scattering - referred to as single-parton scattering - leaving the remainder of each proton only slightly disturbed. Here, we report the study of double- and triple-parton scatterings through the simultaneous production of three J/$\psi$ mesons, which consist of a charm quark-antiquark pair, in proton-proton collisions recorded with the CMS experiment at the Large Hadron Collider. We observed this process - reconstructed through the decays of J/$\psi$ mesons into pairs of oppositely charged muons - with a statistical significance above five standard deviations. We measured the inclusive fiducial cross section to be 272 $^{+141}_{-104}$ (stat) $\pm$ 17 (syst) fb, and compared it to theoretical expectations for triple-J/$\psi$ meson production in single-, double- and triple-parton scattering scenarios. Assuming factorization of multiple hard-scattering probabilities in terms of single-parton scattering cross sections, double- and triple-parton scattering are the dominant contributions for the measured process.


1
High-energy particle accelerators are unique tools to study the structure of matter at the shortest distances. The most powerful accelerator today is the CERN Large Hadron Collider (LHC) that has so far collided beams of protons with energies of up to 6.5 TeV, resulting in center-ofmass energies of up to 13 TeV. Protons are used in energy-frontier colliders because they are relatively easy to accelerate and keep in a circular orbit to enable high collision rates. However, they are not fundamental particles and, in fact, have a complicated quantum mechanical structure. Protons consist of three quarks, two up-type and one down-type, and gluons, which hold together the three valence quarks, as well as of a "sea" of virtual quark-antiquark pairs, which are fundamental elements of the quantum vacuum. All these components are referred to as "partons". In the standard picture of proton-proton (pp) collisions, typically only a few partons undergo a hard scattering with one another, with the remainders of each proton only slightly disturbed in the process. As the collision energy increases, the densities of gluons and sea quarks probed inside each proton grow rapidly. Thus, at high enough energies, more than one pair of partons can undergo a hard scattering in a single pp collision, leading to the simultaneous and independent production of two or more particles with transverse momentum (p T ) and/or mass (m) above a few GeV. Double-parton scatterings (DPS) were first observed at the CERN Intersecting Storage Rings some 35 years ago [1,2] and have been a subject of intense theoretical and experimental investigations ever since [3]. Numerous DPS processes have been studied in measurements with many combinations of pairs of heavy and/or high-p T produced particles [3]. Triple-parton scatterings (TPS), where three hard parton interactions take place simultaneously, have only been proposed for study recently [4] and have not been experimentally explored yet.
Studies of n-parton scattering (NPS) processes are important to elucidate the complicated inner structure of the proton and its evolution with energy [5,6]. Many of these features, e.g., the parton density profile in the plane transverse to the colliding beams, as well as the various correlations (in position, momentum, flavor, color, spin, etc.) between individual partons, are very difficult to calculate theoretically, and can only be mapped out through experimental studies of NPS in different systems and for different numbers n of simultaneous scatterings [3]. Measurements of such processes not only allow for a deeper understanding of the proton structure, but are also of relevance at the LHC to predict backgrounds in rare standard model processes, and in searches for new physics, in final states where multiple heavy particles are produced [7,8]. This work presents the experimental study of the simultaneous production of three massive particles originating in NPS, and the observation of triple-J/ψ meson production. The study of TPS via triple-J/ψ production provides input for further theoretical and experimental progress in understanding the NPS dynamics.
In the simplest approach, ignoring any correlations between the individual partons, the probability to produce n high-p T particles in a given pp collision must be proportional to the n thproduct of probabilities to independently produce each of them. Thus, the probability to produce two or three high-p T particles in DPS and TPS scales with the square and cube, respectively, of the corresponding single-parton scattering (SPS) probabilities [9]. The occurrence of DPS and TPS processes is therefore more likely for final states with large SPS cross sections, such as quarkonia states (e.g., J/ψ and Υ mesons), than for rarer heavier particles such as, e.g., electroweak (EW) bosons [10]. In the DPS case, the cross section to produce, e.g., two charmonium mesons ψ 1 Here, m is a combinatorial factor to avoid double counting, m = 1 (2) if ψ 1 = ψ 2 (ψ 1 = ψ 2 ), and σ eff,DPS is an effective cross section that, in a purely geometric approach, can be determined from the pp transverse overlap [9]. A smaller value of σ eff,DPS , which is proportional to the average (squared) transverse separation of the partons participating in the two hard scatterings, implies larger DPS yields. for processes involving pairs of high-p T jets and/or EW bosons, which are found to lie in the range σ eff,DPS ≈ 10-20 mb [13][14][15][16][17][18][19]. This discrepancy has been mostly explained as evidence of parton correlations in the collision not accounted for in the purely geometrical approaches [20]. In addition, significantly lower σ eff,DPS ≈ 3-10 mb values have been extracted from measurements of quarkonium pair production (J/ψJ/ψ [21-25], J/ψΥ [26], and ΥΥ [8, 27]) that have been interpreted as due to the different dominant species (mostly gluons for quarkonia, and quarks for EW bosons) in the parton distribution functions (PDFs) probed in the different scatterings [3], but can be also attributed in some cases to poorly controlled subtractions of SPS contributions [10].
The study of TPS via triple-J/ψ production can help solve all the issues mentioned above. The equivalent of Eq. (1) for the production of three charmonium mesons in a TPS process reads where m = 1, 3, or 6 (depending on whether all three, two, or none of the ψ i states are identical).
In the absence of parton correlations, the effective cross section σ eff,TPS is closely related to its DPS counterpart via σ eff,TPS = κ σ eff,DPS , with κ of order unity. A value of κ = 0.82 ± 0.11 has been derived in [4] for a variety of proton transverse parton profiles. A theoretical study of the production of three prompt-J/ψ mesons [28], based on the nonrelativistic quantum chromodynamics (NRQCD) approach at leading-order (LO) accuracy as implemented in the HELAC-ONIA code [29,30], has demonstrated that the pure SPS contributions are negligible compared to those from DPS and TPS. Namely, the upper left diagram of Fig. 1 is irrelevant compared to the two other diagrams in the left column of the figure. The experimental measurement of pp → J/ψ J/ψ J/ψ X is thus a golden channel for the study of TPS and, in addition, provides an alternative extraction of σ eff,DPS , thereby shedding light on the underlying dynamics of hard NPS. The production of J/ψ states can also proceed nonpromptly through the decay of a beauty-quark (b) hadron. Notwithstanding a small branching fraction, B b→J/ψ X ≈ 1% [31], the cross section to produce bb pairs is large at the LHC, σ(pp → bb X) ≈ 0.5 mb [4]. The contributions of such processes to inclusive triple-J/ψ production are schematically shown in Fig. 1 (diagrams to the right of the vertical dashed line).
This Letter reports the observation of the simultaneous production of three J/ψ mesons in pp collisions. The analysis is based on a data sample collected at √ s = 13 TeV by the CMS experiment, corresponding to an integrated luminosity of 133 fb −1 . The J/ψ mesons are reconstructed in their dimuon decay mode over a fiducial phase space in transverse momenta and (pseudo)rapidities (p µ,J/ψ T , |η µ |, and |y J/ψ |) defined to maximize the signal purity and the detector acceptance and efficiency. The analysis of the 6µ final state offers a very clean experimental signature for inclusive triple-J/ψ production, comprising prompt and nonprompt components. The Methods section provides more details on the experimental setup and event reconstruction. The signal yield is extracted with a three-dimensional unbinned extended maximum likelihood fit of the m µ + µ − distributions of all J/ψ candidates in the event over the 2.9 < m µ + µ − < 3.3 GeV range. The expected J/ψ mass peaks are modeled with a Gaussian function with mean fixed to their nominal value (m J/ψ = 3.097 GeV) [31] and the root-mean-square (RMS) width fixed to the resolution derived from the MC simulation (σ m ≈ 30 MeV). Given the very low number of events passing the selection, the mass mean and RMS width of the J/ψ mesons cannot be left as free parameters in the fit. The dimuon background is described with an exponential function [23, 32-35]. The fit has eight free parameters for the yields given by the combination of each of the three J/ψ candidates as being either signal or background. The extracted signal yield (red shaded areas in the m µ + µ − distributions of Fig. 2) corresponds to N 3J/ψ sig triple-J/ψ events, with 1.0 +1.4 −0.8 background events. The statistical significance of the signal is evaluated using various methods. From the likelihood ratio of two fits (background-only imposing N 3J/ψ sig = 0, and the default signal-plus-background), with the standard asymptotic formula [36] assuming that the conditions to apply Wilks' theorem [37] are satisfied, a significance of 6.7 standard deviations (std. dev.) is obtained. The significance derived assuming a Poisson counting experiment yields 5.8 std. dev., and it is found to be above 5.5 std. dev. by using MC pseudoexperiments.
To cross-check the size of the combinatorial background derived from the fit, two tests are carried out. First, the fit is repeated over the extended dimuon mass range [2.5-3.3] GeV for the two subleading J/ψ mesons (dotted curves in Fig. 2). This mass range corresponds to an asym- metric window of about [−20σ m , +7σ m ] around the J/ψ nominal mass that, as aforementioned, covers lower dimuon masses where the background, if any, should be larger. The obtained signal yield is fully consistent with the default result. A second test is performed whereby the opposite-sign (OS) requirement is removed to allow also for same-sign dimuon combinations (µ ± µ ± ) for the two subleading pairs. After applying the rest of the selection criteria of the default analysis, no triplet events containing same-sign muon pairs are observed.
In order to estimate the average prompt and nonprompt contributions in the triple-J/ψ events, the proper decay-length of each J/ψ is calculated as T | is the transverse distance between the J/ψ decay vertex and the PV ( r is the vector from the PV to the J/ψ vertex). The experimental resolution on L J/ψ is ≈20 µm. Prompt J/ψ mesons are defined as those having L J/ψ < 60 µm, a choice that reduces the nonprompt contribution by more than one order in magnitude [33]. A signal-only subset of events is defined based on all J/ψ candidates found in a narrower invariant mass region, within ±3 std. dev. around m J/ψ , than that used in the default signal extraction. In this range, five triple-J/ψ events are found that can be classified as two events being consistent with the "2 nonprompt + 1 prompt J/ψ" hypothesis, and one event each in the "1 nonprompt + 2 prompt J/ψ", "3 nonprompt J/ψ", and "3 prompt J/ψ" categories. Such a classification is confirmed with a second method where prompt and nonprompt weights are extracted with the SPLOT technique [38], exploiting an unbinned maximum likelihood fit of the L J/ψ distributions.
The cross section for inclusive triple-J/ψ production measured in the fiducial region defined in Table 1 is the number of signal events, L int the total integrated luminosity, B J/ψ→µ + µ − the dimuon branching fraction, and = trig id reco the total efficiency composed of trigger, identification, and reconstruction components. The J/ψ muon identification and reconstruction efficiencies are extracted with the tag-and-probe method using correction factors from the large inclusive J/ψ → µ + µ − data samples used in Ref.
[39]. Since they depend on the p T and η of the muons, they propagate into the final cross section via two-dimensional maps, yielding id reco = 0.78. The trigger efficiency is found to be trig = 0.84 from a study of the MC samples.
The impact on the extracted cross section of the choice of functions used to reproduce the shapes of the signal and background dimuon invariant masses is studied. For the signal, the Gaussian distribution is changed to a Crystal-Ball function [40] as well as to a Gaussian function with the RMS width left to vary in the fit. The background shape is changed from For all muons For all J/ψ mesons p T > 6 GeV and |y| < 2.4 2.9 < m µ + µ − < 3.3 GeV the default exponential to first-and zeroth-order polynomials. The relative differences in the cross sections obtained from the alternative modeling for signal and background are 0.8 and 3.4%, respectively, and are assigned as corresponding systematic uncertainties. Uncertainties arising from the muon reconstruction and identification efficiencies are derived by allowing the tag-and-probe correction factors for each (p T , η) bin to vary within their precision, and checking the effect on the extracted cross section. The maximal variation observed is ±1.0%. Varying the relative composition of double-and single-J/ψ meson production in the MC event sample used for the determination of the trigger efficiency leads to a 3.4% propagated uncertainty. Uncertainties of 1.6% and 3.0% are added from the integrated luminosity measurement [41-43], and from the simulated signal sample size, respectively. The uncer- propagates into a 1.7% uncertainty in the cross section. The total systematic uncertainty of the measured cross section is 6.2%, obtained by adding all individual sources in quadrature (Extended Data Table 5). The measured cross section for triple-J/ψ production, within the fiducial region defined in Table 1, is The total inclusive triple-J/ψ cross section is expected to correspond to the sum of the contributions from the SPS, DPS, and TPS processes schematically shown in Fig. 1, each of which contains various combinations of prompt (p) and nonprompt (np) J/ψ contributions, Under the simplest assumption of factorization of multiple hard-scattering probabilities in terms of SPS cross sections, the DPS and TPS contributions to triple-J/ψ production (last row of Eq. (3)) can be written through Eqs. (1) and (2) as a combination of products of single-and double-J/ψ SPS cross sections as and with combinatorial prefactors m 1 = 2/2 = 1, m 2 = 3/3! = 1/2, and m 3 = 1/3! = 1/6. Therefore, from the eight individual SPS cross sections for single-, double-, and triple-J/ψ production, one can determine the total triple-J/ψ production cross section via Eqs. (3)(4)(5). The values of the relevant SPS cross sections, each within the fiducial phase space defined in Table 1, are computed as described next and listed in Table 2. Table 2: Theoretical J/ψ, J/ψ + J/ψ, and J/ψ + J/ψ + J/ψ cross sections. Predictions for single-, double-, and triple-J/ψ production cross sections in SPS processes, which pass the fiducial criteria listed in Table 1, derived from the HELAC-ONIA (HO) and MADGRAPH5 aMC@NLO (MG5NLO) matrix element calculators, complemented with the PYTHIA 8 (PY8) parton shower, as described in the text.
SPS single-J/ψ production SPS double-J/ψ production SPS triple-J/ψ production The values of the SPS single, double, and triple prompt-J/ψ cross sections (σ  Table 2  , except for the single prompt-J/ψ predictions that have a better precision because they are determined with an explicit fit of the NRQCD predictions to the LHC data [46] and have an associated 10% uncertainty of experimental origin. All these sources are treated as uncorrelated and the corresponding uncertainties are added in quadrature. Using Eqs. (4,5) with the SPS cross sections listed in Table 2, and assuming that the effective DPS and TPS cross sections are related by σ eff,TPS = (0.82 ± 0.11) σ eff,DPS [4] in a baseline approach that ignores parton correlations, one can extract the value of the effective DPS cross section that yields the experimentally measured σ 3J/ψ tot value. Following such a procedure, the value σ eff,DPS = 2.7 +1.4 −1.0 (exp) +1.5 −1.0 (theo) mb is derived, where the first uncertainty is due to the experimental σ 3J/ψ tot precision and the second is due to the propagation of all sources of theoretical uncertainties in the ingredients of Eqs. (3)(4)(5).
The inclusive triple-J/ψ theoretical cross sections and yields for each individual process contributing to the total production are listed in Table 3. The expected contributions from SPS, DPS, and TPS processes to the total triple-J/ψ cross section amount to about 6, 74, and 20%, respectively. This confirms the conclusion of Ref.
[28] that triple-J/ψ production is a golden channel to study DPS and TPS, with minimal SPS contamination. The largest contributors to the triple-J/ψ cross section are σ amounting to about 7% each, and σ 3np SPS representing about 4% of the total production. In terms of prompt and nonprompt contributions, the theoretical expectation for the production of three promptly produced J/ψ mesons is ≈5% of the total yield, whereas the percentage expected for Table 3: Expected contributions to triple-J/ψ production. Central values of the predictions for triple-J/ψ production cross sections (in fb) and yields from SPS, DPS (for σ eff,DPS = 2.7 mb), and TPS (for σ eff,TPS = 0.82 σ eff,DPS = 2.2 mb) processes, and their total sum. The values are given in columns for combinations of n prompt and (3 − n) nonprompt J/ψ mesons, the last column giving their corresponding sums. The DPS and TPS results are derived via Eqs. (3-5) from the SPS cross sections listed in Table 2   three nonprompt J/ψ mesons is ≈45%. The remaining half of the triple-J/ψ events are expected to be due to the combination of J/ψ mesons produced promptly and from beauty hadron decays. This result is consistent, within the large statistical uncertainties of the present data set, with the combination of prompt and nonprompt J/ψ mesons derived from the decay length of each dimuon candidate.
In Fig. 3 The effective cross sections obtained from quarkonium measurements favor a smaller value of σ eff,DPS ≈ 3-10 mb compared to the σ eff,DPS ≈ 10-20 mb derived from harder or heavier final states. Such an apparent process-dependent σ eff,DPS value is suggestive of different parton transverse profiles, and/or correlations present, probed inside the proton at varying fractional momenta, given by x = At midrapidity (|η| < 2.5), quarkonia are produced mostly in gluon-gluon scatterings carrying a fraction x ≈ 5 × 10 −4 of the proton momentum, whereas mostly quarks with x ≈ 10 −2 participate in the production of EW bosons. The fact that LHCb measurements of double-quarkonia and quarkonia-pluscharm [25, 61] at forward rapidities (η ≈ 2-4.5), processes that originate in parton scatterings with asymmetric fractional momenta x 1 ≈ 10 −4 and x 2 ≈ 10 −2 , lead to values of σ eff,DPS ≈ 15 mb, larger than those measured at midrapidity for similar final states, seems to confirm the dependence of the effective DPS cross section on the relevant parton species and x fractions probed.    [26] D0 Collaboration, "Evidence for simultaneous production of J/ψ and Υ mesons", Phys.

Experimental setup
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Forward calorimeters extend the η coverage provided by the barrel and endcap detectors. Muons are measured over the range |η| < 2.4 in gas-ionization detectors, embedded in the steel flux-return yoke outside the solenoid, made using three technologies: drift tubes, cathode strip chambers, and resistive plate chambers. A more detailed description of the CMS detector can be found in Ref. [64].
Events of interest are selected using a two-tiered trigger system. The first trigger level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz with a fixed latency of about 4 µs [65]. The second level (or high-level trigger, HLT) consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage [66]. The present analysis employs an HLT that requires three muons, each having p T > 3.5 GeV for |η| < 1.2 (barrel) or p T > 2. . Triple-J/ψ events are simulated by combining events from two different MC samples with varying relative fractions. The first MC sample contains three singly produced J/ψ mesons, whereas the second sample combines single-and double-J/ψ meson production, following Ref. [28]. Since the J/ψ mesons originating in both processes have slightly different p T spectra, the variation of their relative fractions in the simulation allows to effectively scan the triple-muon trigger efficiency in the more uncertain region close to the trigger p T thresholds. The generated events are then processed through a detailed GEANT4 simulation [69] of the CMS detector response.

Event reconstruction and selection
Muons are reconstructed by combining information from the silicon tracker and the muon system [39]. The matching between tracks reconstructed in each of the subsystems proceeds either outside-in, starting from a track in the muon system, or inside-out, starting from a track provided by the silicon tracker. Matching muons to tracks measured in the silicon tracker leads to a relative p T resolution of 1% in the barrel and 3% in the endcaps for muons with p T up to 100 GeV [39]. The candidate vertex with the largest value of summed physics-object p 2 T in the event is taken to be the primary pp interaction vertex (PV) [70]. Only charged particles that are either directly participating in the PV determination, or are closest to the PV along the z direction and not used in any other PV identification (namely, they are consistent with secondary decays from the closest PV), are used in the analysis. Simulation studies show that this procedure efficiently suppresses any potential contamination from charged particles produced in any other pp collision taking place in the same bunch crossing (pileup).
The offline data analysis starts by selecting events with six or more reconstructed muons, each passing the p T and η criteria implemented at the HLT. The muons are combined into OS pairs, and are considered for further study if their invariant mass is 2.9 < m µ + µ − < 3.3 GeV (corresponding to about ±7 times the J/ψ mass resolution discussed below) and originate from a common vertex with a probability greater than 0.5%, thereby reducing the random pairing of muons originating from other background sources in the same event, as discussed below. All selected muon pairs are further required to share the same PV (including the possibility that they originate from secondary vertices associated to a common PV) to eliminate accidental combinations of muons from different pp collisions in the same or neighboring bunch crossings. The analysis thereby includes prompt-J/ψ mesons coming directly from the pp interaction (or from feed-down decays of promptly produced and decayed resonances), and nonprompt ones coming from the decays of beauty hadrons. In order to ensure high purity and reconstruction efficiency, the J/ψ candidates are required to have p T > 6 GeV and |y| < 2.4.
After applying the selection criteria discussed above, six triple-J/ψ events are observed. No alternative pairings among the OS muons in these six events are found to satisfy the analysis requirements. In all events, the J/ψ meson with highest p T ("leading" J/ψ) is found to correspond to the muon pair that passed the online trigger selection. In inclusive J/ψ measurements [23, 32-35], a continuum background is present in the m µ + µ − distribution due to random combinations of OS muons originating from semileptonic beauty-and charm-quark hadron decays, b → µ + X and c → µ + X, and Drell-Yan events, which pass the trigger and data analysis criteria. To monitor the sideband background population, the mass windows corresponding to the two subleading J/ψ mesons are extended well below the J/ψ mass region. The corresponding dimuon invariant mass distributions are shown in Fig. 2 ordered (left to right) by decreasing p T of the µ + µ − system. Only one additional background event is found in the extended mass region indicated by the dotted curves, confirming that the combinatorial continuum is suppressed by the requirement of having three reconstructed J/ψ candidates in the same event.

Event kinematic properties
Detailed information of the kinematic properties of the J/ψ candidates passing the triple-J/ψ selection criteria are shown in Extended Data Table 4. The kinematic distributions of all six J/ψ triplets do not show any local peak, which could be indicative of, e.g., a resonance decaying into three J/ψ mesons, but are distributed featureless (with large point-to-point statistical uncertainties) over triple-J/ψ mass, p T , and η.
All sources of systematic uncertainty in the triple-J/ψ cross section measurement are listed in Extended Data Table 5. Table 4: Dimuon invariant mass, proper decay-length, transverse momentum, rapidity, and azimuthal angle of each of the three J/ψ candidates measured in the six triple-J/ψ events passing our selection criteria. (The last J/ψ candidate in the last row, with invariant mass m J/ψ,3 = 2.94 GeV, is likely a background event).

Data availability
Tabulated results are provided in the HEPData record for this analysis [71]. Release and preservation of data used by the CMS Collaboration as the basis for publications is guided by the CMS policy as stated in CMS data preservation, re-use and open access policy.

Code availability
The CMS core software is publically available at https://github.com/cms-sw/cmssw.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid and other centers for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC, the CMS detector, and the supporting computing infrastructure provided by the following funding agencies: