New structures in the J/$\psi$J/$\psi$ mass spectrum in proton-proton collisions at $\sqrt{s}$ = 13 TeV

A search is reported for near-threshold structures in the J/$\psi$J/$\psi$ invariant mass spectrum produced in proton-proton collisions at $\sqrt{s}$ = 13 TeV from data collected by the CMS experiment, corresponding to an integrated luminosity of 135 fb$^{-1}$. Three structures are found, and a model with quantum interference among these structures provides a good description of the data. A new structure is observed with a significance above 5 standard deviations at a mass of 6638 $^{+43}_{-38}$ (stat) $^{+16}_{-31}$ (syst) MeV. Another structure with even higher significance is found at a mass of 6847 $^{+44}_{-28}$ (stat) $^{+48}_{-20}$ (syst) MeV, which is consistent with the X(6900) resonance reported by the LHCb experiment and confirmed by the ATLAS experiment. Evidence for another new structure, with a local significance of 4.7 standard deviations, is found at a mass of 7134 $^{+48}_{-25}$ (stat) $^{+41}_{-15}$ (syst) MeV. Results are also reported for a model without interference, which does not fit the data as well and shows mass shifts up to 150 MeV relative to the model with interference.


1
The prospect of exotic hadrons-states other than qq or qqq-dates back to Gell-Mann's seminal 1964 paper [1].However, interest in exotics, both theoretical and experimental, wavered for decades.The discovery of the charmoniumlike X(3872) state by the Belle Collaboration in 2003 [2] propelled exotics from speculative chimeras to the frontier of hadronic physics [3].While many tetraquark candidates containing heavy quarks (c, b) are now known, questions still abound, including which of these are truly exotic hadrons [3] and, if they are bound states, what is their internal structure (e.g., molecules, bound states of diquarks, etc. [4][5][6]).
Beginning with the J/ψ discovery in 1974 [7,8], heavy quarkonia brought clarity to the quark model.Similarly, exotic states composed entirely of heavy quarks may offer analogous insights.Exotic states decaying into the J/ψJ/ψ and ΥΥ channels are particularly promising experimentally because of the efficient triggering and reconstruction of muonic channels.The possibility of structure appearing in these channels is enhanced by an empirical pattern first observed with light mesons: the mass spectra of two vector states with the same isospin (e.g., ϕϕ and ωϕ) have near-threshold enhancements, whereas systems of two vector mesons of different isospin (e.g., ϕρ and ωρ) have no structure [9].Such a behavior has been corroborated by the enhancements observed in J/ψω [10] and J/ψϕ [11] spectra.
Here we report on the J/ψJ/ψ invariant mass spectrum from proton-proton (pp) collisions as recorded with the CMS detector, where the J/ψ mesons are reconstructed from µ + µ − pairs.The data sample corresponds to an integrated luminosity of 135 fb −1 [39][40][41] at a center-ofmass energy of 13 TeV.
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 end cap sections.Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.A more detailed description of the CMS detector can be found in Ref. [42].
Several Monte Carlo (MC) samples are used to model the detector mass resolution and efficiency of signals.As simulations of tetraquarks states are not available and the analysis is independent of signal kinematic details, we use the production of other particles to simulate the signal, with the masses defined appropriately.We use J P = 0 + mesons for a variety of masses decaying to J/ψJ/ψ as the signal and simulate their production via gluon-gluon interactions.We model signal particles by χ b0 (1P) production using the PYTHIA event generator (version 8.230) [43].As an alternate production model, signal is also simulated by generating a sample with "Higgs bosons" decaying to ZZ ( * ) , where the Z ( * ) is redefined as a J/ψ, using the JHUGen event generator (version 7.40) [44,45].To understand the detector mass resolution, the samples are generated with a negligible natural width.Major sources of background, i.e., selected events that are not from a resonant decay to two J/ψ candidates, are expected to originate from two genuine J/ψ decays arising either from a single parton-parton collision, i.e., nonresonant single-parton scattering (NRSPS) [46][47][48][49][50][51], or a pair of parton-parton interactions in a single pp collision, i.e., double-parton scattering (DPS) [52][53][54][55][56][57][58].The NRSPS contribution is expected to dominate the DPS contribution near the J/ψJ/ψ threshold, with the DPS contribution dominating at masses above 11 GeV.To simulate these backgrounds we use the PYTHIA generator.The next-to-leading order CASCADE generator [59] and the next-to-next-to-leading order HELAC-Onia generator [47,48,60,61] are also used to simulate alternative NRSPS shapes for the estimation of the systematic uncertainty in the background.Direct production, as well as feed-down processes where J/ψ mesons arise as decay products, are included.Feed-down processes from double-charmonium resonances [38] are generated using the PYTHIA generator separately.Generated events are processed through the full CMS detector simulation, which utilizes the GEANT4 toolkit [62].
Events of interest are selected using a two-tiered trigger system [63].The first level L1, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events.The second level, known as the high-level trigger (HLT), consists of a farm of processors running a faster version of the full event reconstruction software, which reduces the event rate to around 1 kHz.The data collected in 2016 correspond to an integrated luminosity of 36.3 fb −1 , and were obtained with an L1 trigger that required at least three muon candidates to balance the requirements of high efficiency for the signal and a low trigger rate.The HLT required that the pseudorapidity of each muon satisfies |η| < 2.5 and that a muon pair candidate is composed of oppositely charged muon candidates with an invariant mass between 2.95 and 3.25 GeV, a distance of closest approach between the two muons less than 0.5 cm, and a fit to a common vertex with a χ 2 probability greater than 0.5%.The data collected in 2017-2018 correspond to an integrated luminosity of 98.6 fb −1 and utilized an L1 trigger requiring at least three muons with at least two muons having p T > 3 GeV and one of the two muons having p T > 5 GeV, and at least one oppositely charged pair of muons with invariant mass below 9 GeV.The HLT criteria replicated those of 2016 with the additional requirement that the two muons with invariant mass between 2.95 and 3.25 GeV each have p T > 3.5 GeV.
The off-line selection criteria were finalized before unblinding the signal region (defined as invariant mass below 7.8 GeV) and are based in part on an examination of data collected in 2011 and 2012.It requires at least four muons, each of them having p T > 2 GeV, |η| < 2.4, and satisfying the "soft muon identification" requirements [64].Pairs of oppositely charged muons are used to form J/ψ candidates, and these are required to have p T > 3.5 GeV.To confirm the HLT requirement, each muon pair is fit to a vertex, requiring a fit probability greater than 0.5%.Combinations of four muons from two J/ψ candidates are fit to a common vertex, which must have a fit probability greater than 0.5%.A mass constrained vertex fit is applied to the oppositely charged muon pairs, in which the muon pair's mass is constrained to the known J/ψ mass m J/ψ [65], and the probability is required to be greater than 0.1%.The two constrained J/ψ candidates are then fit to a common vertex and its fit probability must exceed 0.1%.An event may have multiple J/ψ-pair candidates.If both charge permutations of the same four muons have masses within 3 standard deviations (3σ) of the known J/ψ mass, we select the pairing that minimizes (∆m 1 /σ m 1 ) 2 + (∆m 2 /σ m 2 ) 2 , where ∆m is the difference between the reconstructed dimuon mass and the J/ψ mass, and σ m is the calculated uncertainty of the reconstructed dimuon mass (affects 0.2% of events).All J/ψJ/ψ candidates are kept if they arise from more than four muons (affects 0.2% of events).
The J/ψJ/ψ invariant mass spectrum of the 8651 selected candidate pairs below 9 GeV is shown in Fig. 1.To better constrain the smooth background continuum, the full analysis range extends up to an invariant mass of 15 GeV (the corresponding mass spectrum is shown in the Supplemental Material [66]).[GeV] 1.The J/ψJ/ψ invariant mass spectrum in the range up to 9 GeV, with fits consisting of three signal functions (BW 1 , BW 2 , and BW 3 ) and a background model (see text).Left: the fit without interference.Right: the fit that includes interference, where "Interfering BWs" refers to the total contribution of all the interfering amplitudes and their cross terms.For clarity, only the sum of the three background components (NRSPS + DPS + BW 0 ) is shown on the plots.The lower portion of the plots shows the pulls, i.e., the number of standard deviations (statistical uncertainties only) that the binned data differ from the fit.
Using an unbinned extended likelihood method, we fit the data from the threshold 2m J/ψ to 15 GeV.Signals are represented by relativistic S-wave Breit-Wigner (BW) functions [65,67,68] convolved with resolution functions, which are sums of two Gaussian functions with the same mean value.The detector mass resolution varies from about 10 MeV at a mass of 6500 MeV to 18 MeV at 7300 MeV.The BW functions are not modified by acceptance or efficiency corrections because these vary only slowly, by less than 1.3% over the natural width of the X(6900).
The NRSPS and DPS backgrounds are parametrizations of distributions of the corresponding simulated events.The NRSPS component includes an exponential function where one of its parameters is a free parameter in the fit.The effect of feed-down from heavier mass states such as the ψ(2S) is accounted for as a systematic uncertainty.Combinatorial backgrounds arise when one or more muon candidates are not from a J/ψ meson decay and from hadrons that are misidentified as a muon or misreconstructed.The shapes of combinatorial backgrounds are well modeled by the NRSPS + DPS parametrizations, and the residual effects are accounted for in the systematic uncertainties.
Aside from the NRSPS and DPS background components, a component to model a threshold enhancement is included, as was done in Ref. [14].This excess may be due to a resonance, but other processes could be responsible: coupled-channel interactions [69], triangle singularities [70], Pomeron exchange [36], or simply an inadequate NRSPS model.This region could also include feed-down from higher mass tetraquarks that are only partially reconstructed.Because of these uncertainties, we regard this enhancement as an additional background component and model it with a BW function (BW 0 ) with free mass and width, which provides a good ad hoc description of this feature.
We begin with a background-only fit (NRSPS + DPS + BW 0 ) and then add signal contributions one at a time to the fit, with the yields of each signal and background component being free parameters of the fit.In this first series of fits, the additional components are added without interference.A signal BW is kept as long as its local significance exceeds 3σ, with the local significance calculated from the difference in log-likelihood between including and not including the component and only accounting for statistical uncertainties.Three resonance structures, la-beled BW 1 , BW 2 , and BW 3 in order of increasing mass, are found to be statistically significant.The local significance of these peaks, calculated from the log-likelihood differences between the full fit and the fit with the BW function of interest removed, are 6.5σ, 9.4σ, and 4.1σ, respectively.Figure 1 (left) shows this fit, with numerical results in Table I.The BW function widths are all much larger than the detector mass resolution.Over the full mass range of 6-15 GeV, the proportions of the background components to the total are about 58%, 25%, and 9% for NRSPS, DPS, and BW 0 , respectively.We quantify the goodness of fit by the probability of the fit's χ 2 (over the full mass range with the bin size of Fig. 1) and the number of degrees of freedom.Because the signal region is only a small fraction of the fit range, the impact of deviations from the model in this region is diluted.
As an alternative, we also compute a probability using only the bins below 7.8 GeV, which we refer to as the signal-region χ 2 probability.
The full χ 2 fit probability is 98%.However, the dips in the data around 6750 and 7150 MeV are poorly described: the signal-region χ 2 probability is 9%.Agreement may be improved by introducing interference, as was done by LHCb [14].
Our construction of interference models is based on the three structures found in the nointerference fit.We consider interference between pairs of ("two-way"), or all three ("threeway") components.The three-way interference is implemented with a term proportional to |r 1 exp(iϕ 1 )BW 1 + BW 2 + r 3 exp(iϕ 3 )BW 3 | 2 , where r 1,3 and ϕ 1,3 are the relative magnitudes and phases of BW 1,3 with respect to BW 2 [65].We report the results from the three-way interference model, which shows a signal-region χ 2 probability of 65%, whereas all two-way models have probabilities below 30%.The three-way interference fit is shown in Fig. 1 (right), and the fit parameters are listed in Table I.
In this fit, the local significances of all three structures increase with respect to the model without interference.For the structure with the lowest significance (BW 3 ), a value of 4.7σ is obtained.The global significance of BW 3 is estimated by generating MC pseudoexperiments and determining the probability that a statistical fluctuation yields a "signal" whose local significance equals or exceeds that of BW 3 in the windows 7.05-7.8GeV for mass and 50-260 MeV for width.The BW 3 phase is unconstrained.The resulting global significance for BW 3 is 3.4σ.
The interference between different resonances is motivated by the idea that the states could have the same quantum numbers and be coherently produced.One can also consider interference between the resonances and the NRSPS background.We consider this scenario less probable as the background is likely a mixture of J PC states.We have investigated two models with this sort of interference.The first model is identical to Model II from Ref. [14] with a signal BW function for the X(6900), an auxiliary BW function around 6700 MeV, the NRSPS and DPS background, and interference between the NRSPS background and auxiliary BW function.This did not provide a good fit to our data as shown in the Supplemental Material [66].The second model starts from the three-way interference model and adds interference between the three signal BW functions and the NRSPS background.This model did not significantly improve the fit quality.
While interference between resonances is one possible mechanism, other models may also be able to reproduce the dips.For example, inspired by predictions of dense spectroscopy of tetracharm states [15][16][17][18][19][20][21], representing individual peaks by multiple overlapping narrow resonances could provide good fits by enabling the BW functions to fall more deeply into the dips.
Systematic uncertainties for masses and widths are determined by varying aspects of, and inputs to, the fits.For a given source, the largest deviation from a parameter's nominal value is taken as its systematic uncertainty.The sources of systematic uncertainty that have been considered are as follows: different BW function shapes (P and D wave, and alternative values for the "interaction radius" in the Blatt-Weisskopf barrier factor [67,68]); alternative DPS parametrizations, such as those obtained by mixing data events (artificial J/ψJ/ψ mass spectrum formed by two J/ψs from different events); NRSPS parametrizations from different simulations (CASCADE and HELAC-Onia) or floating otherwise fixed parameters individually; momentum scale (based on the shift in the unconstrained fitted J/ψ meson mass from the world average value [65]); detector mass resolution (slightly different resolution models based on the PYTHIA or JHUGen generators); combinatorial background shape (the result of altering the parameterization of the combinatorial background); efficiency corrections (the difference between not applying and applying efficiency corrections based on the PYTHIA or JHUGen generator); and the effects of including various feed-down components from hypothetical heavier tetraquarks [38].The effects of feed-down components in the interference model can produce asymmetric uncertainties, so an asymmetric uncertainty is assigned for this source.For other sources, symmetric uncertainties are assigned.The principal systematic effects are summarized in Table II.The total uncertainties are their sum in quadrature and also appear in Table I.
Including feed-down in the fit affects signal structures, especially for the BW 1 , because of their overlap in the J/ψJ/ψ mass distribution.Compared to the no-interference model, the greater complexity and increased parameter correlations of the interference model result in larger systematic uncertainties.
The impact of systematic uncertainties on the local significances of the structures was checked in the interference fit model by using discrete sets of individual alternative hypotheses and recomputing the significances.The systematic uncertainties introduce no appreciable change.
The BW 2 parameters are within 2σ, when comparing the analogous no-interference or interference models, of the X(6900) values reported by LHCb.For their model I, consisting of a BW function for the X(6900) signal plus two auxiliary threshold BW functions, they reported an X(6900) mass of 6905 ± 11 ± 7 MeV with a natural width of 80 ± 19 ± 33 MeV, and 6886 ± 11 ± 11 MeV and 168 ± 33 ± 69 MeV for model II [14].However, because the two experiments use different fit models, quantitative comparisons between them are not straightforward, and thus we also fit our data with the two LHCb models.Neither LHCb model provides a good description of the CMS data; nevertheless, our X(6900) parameters are again compatible.More information about these comparisons, including where some small discrepancies are present, is summarized in the Supplemental Material [66].
Our measured masses appear compatible with recent calculations of the cccc spectrum [21,71], which would indicate that these three structures may be a family of radial excitations of the same J PC .This is the case for both no-interference and interference masses, albeit for different theoretical models.
In summary, the J/ψJ/ψ invariant mass spectrum has been presented.The data were collected with the CMS detector from pp collisions at √ s = 13 TeV and correspond to an integrated luminosity of 135 fb −1 .Three structures are found in the J/ψJ/ψ invariant mass spectrum.The spectrum is better described by a model with interference among three resonances.Two new structures, tentatively named X(6600) and X(7100), are found with masses of 6638 +43 −38 (stat) +16 −31 (syst) and 7134 +48 −25 (stat) +41 −15 (syst) MeV; and the X(6900) structure observed by LHCb is confirmed with a mass of 6847 +44 −28 (stat) +48 −20 (syst) MeV.The local statistical significances of these peaks are, for increasing mass, 7.9, 9.8, and 4.7 standard deviations.Numeric results are provided in the HEPData record [72].
Note added.-Recently, the ATLAS Collaboration released a Letter [73] confirming the X(6900) structure in the J/ψJ/ψ spectrum and also reported near-threshold excesses, including a possible feature in the ψ(2S)J/ψ channel.

A Fit projections in full range
The J/ψJ/ψ invariant mass spectrum for the full range of the fit is shown in Fig. A.1 for the no-interference and the interference models.

FIG.
A.1.The J/ψJ/ψ invariant mass spectrum with the no-interference fit (upper) and the interference fit (lower) in the full fit range (see text for model details).The "Interference BWs" curve is the total contribution of all the interference amplitudes and their cross terms.The lower portion of the plots shows the pulls, i.e., the number of standard deviations (statistical uncertainties only) that the binned data differ from the fit.

B Fitting with the LHCb models
In order to make a direct comparison with the results from the LHCb Collaboration [14] we also used their two principal models (I and II) to fit the CMS data.The details are provided here.
Model I consists of two auxiliary relativistic Breit-Wigner (BW) functions below the X(6900) BW.The parameter values for these two BWs are not reported in Ref. [14], but the second auxiliary BW is in the same mass region as the CMS "BW 1 " and so we extend this labeling to the LHCb models.
The Model I fit is shown in Fig. B.1 (upper), with numerical results listed in Table B.1.The fit quality is poor with a signal-region χ 2 probability of 0.9%.The no-interference fit using the CMS model with three structures-and roughly the same number of J/ψJ/ψ candidates after the final selection as LHCb-has a signal-region χ 2 probability of 9%.The J/ψJ/ψ invariant mass spectrum with fits using LHCb models [14]: Model I (nointerference) (upper) and Model II (interference) (lower).Model details are given in the text.The "Interfering components" is the total contribution of BW 1 , NRSPS, and their interference.The lower portion of the plots shows the pulls between the binned data and the fit.The 6750 MeV dip is now well described: the pulls in Fig. B.1 (lower) are well behaved in the vicinity of the X(6900) peak.However, systematic deviations are evident elsewhere in the pulls, and the signal-region χ 2 probability of 0.8% for this fit is actually lower than for Model I.
FIG.B.1.The J/ψJ/ψ invariant mass spectrum with fits using LHCb models[14]: Model I (nointerference) (upper) and Model II (interference) (lower).Model details are given in the text.The "Interfering components" is the total contribution of BW 1 , NRSPS, and their interference.The lower portion of the plots shows the pulls between the binned data and the fit.

TABLE . I
. Summary of the fit results for the J/ψJ/ψ invariant mass distribution.The mass m and natural width Γ for both the no-interference model and the interference model, and the signal yields N for the no-interference model, are given for the three signal structures.The dual uncertainties are the statistical followed by the systematic components, and single uncertainties are only statistical.

TABLE .
II. Dominant contributions to the systematic uncertainties in masses and widths, in MeV, for the two fits.The "Total uncertainty" is the quadratic sum of all individual components, including the unlisted nondominant contributions.M BW 2 M BW 3 Γ BW 1 Γ BW 2 Γ BW 3

TABLE . B
[14]ummary of results for the two LHCb fit models from LHCb and CMS: the mass m and natural width Γ, in MeV, are given for the principal structures (the Breit-Wigner function at the threshold is not included in Model I), and LHCb values not reported in Ref.[14]are marked as "-".Single uncertainties are statistical only.Model II 6741 ± 6 288 ± 16 6886 ± 11 ± 11 168 ± 33 ± 69 CMS Model II 6736 ± 38 439 ± 65 6918 ± 10 187 ± 40 Model II has, aside from the signal BW, a single auxiliary BW (here labeled as 'BW 1 '), which interferes with the NRSPS background.This interference was introduced to better model the dip around 6750 MeV.Our fit results for Model II are shown in Fig. B.1 (lower) and Table B.1.