Observation of structure in the $J/\psi$-pair mass spectrum

Using proton-proton collision data at centre-of-mass energies of $\sqrt{s} = 7$, $8$ and $13\mathrm{\,TeV}$ recorded by the LHCb experiment at the Large Hadron Collider, corresponding to an integrated luminosity of $9\mathrm{\,fb}^{-1}$, the invariant mass spectrum of $J/\psi$ pairs is studied. A narrow structure around $6.9\mathrm{\,GeV/}c^2$ matching the lineshape of a resonance and a broad structure just above twice the $J/\psi$ mass are observed. The deviation of the data from nonresonant $J/\psi$-pair production is above five standard deviations in the mass region between $6.2$ and $7.4\mathrm{\,GeV/}c^2$, covering predicted masses of states composed of four charm quarks. The mass and natural width of the narrow $X(6900)$ structure are measured assuming a Breit--Wigner lineshape.


Introduction
The strong interaction is one of the fundamental forces of nature and it governs the dynamics of quarks and gluons. According to quantum chromodynamics (QCD), the theory describing the strong interaction, quarks are confined into hadrons, in agreement with experimental observations. The quark model [1,2] classifies hadrons into conventional mesons (qq) and baryons (qqq or qqq), and also allows for the existence of exotic hadrons such as tetraquarks (qqqq) and pentaquarks (qqqqq). Exotic states provide a unique environment to study the strong interaction and the confinement mechanism [3]. The first experimental evidence for an exotic hadron candidate was the χ c1 (3872) state observed in 2003 by the Belle collaboration [4]. Since then a series of novel states consistent with containing four quarks have been discovered. Recently, the LHCb collaboration observed resonances interpreted to be pentaquark states [5][6][7][8]. All hadrons observed to date, including those of exotic nature, contain at most two heavy charm (c) or bottom (b) quarks, whereas many QCD-motivated phenomenological models also predict the existence of states consisting of four heavy quarks, i.e. T Q 1 Q 2 Q 3 Q 4 , where Q i is a c or a b quark . Theoretically, the interpretation of the internal structure of such states usually assumes the formation of a diquark (Q 1 Q 2 ) and an antidiquark (Q 3 Q 4 ) attracting each other. Application of this diquark model successfully predicts the mass of the doubly charmed baryon Ξ ++ cc [34, 35] and helps to explain the relative rates of bottom baryon decays [36]. Tetraquark states comprising only bottom quarks, T bbbb , have been searched for by the LHCb and CMS collaborations in the Υ µ + µ − decay [37,38], with the Υ state consisting of a bb pair. However, the four-charm states, T cccc , have not yet been studied in detail experimentally. A T cccc state could disintegrate into a pair of charmonium states such as J/ψ mesons, with each consisting of a cc pair. Decays to a J/ψ meson plus a heavier charmonium state, or two heavier charmonium states, with the heavier states decaying subsequently into a J/ψ meson and accompanying particles, are also possible. Predictions for the masses of T cccc states vary from 5.8 to 7.4 GeV/c 2 [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24], which are above the masses of known charmonia and charmonium-like exotic states and below those of bottomonium hadrons. This mass range guarantees a clean experimental environment to identify possible T cccc states in the J/ψ-pair (also referred to as di-J/ψ) invariant mass (M di-J/ψ ) spectrum.
In proton-proton (pp) collisions, a pair of J/ψ mesons can be produced in two separate interactions of gluons or quarks, named double-parton scattering (DPS) [39][40][41], or in a single interaction, named single-parton scattering (SPS) [42][43][44][45][46][47][48][49]. The SPS process includes both resonant production via intermediate states, which could be T cccc tetraquarks, and nonresonant production. Within the DPS process, the two J/ψ mesons are usually assumed to be produced independently, thus the distribution of any di-J/ψ observable can be constructed using the kinematics from single J/ψ production. Evidence of DPS in pp collisions has been found in studies at the Large Hadron Collider (LHC) experiments [50-54]. The LHCb experiment has measured the di-J/ψ production in pp collisions at centre-ofmass energies of √ s = 7 [55] and 13 TeV [56]. The DPS contribution is found to dominate the high M di-J/ψ region, in agreement with expectation.
In this paper, fully charmed tetraquark states T cccc are searched for in the di-J/ψ invariant mass spectrum, using pp collision data collected by LHCb at √ s = 7, 8 and 13 TeV, corresponding to an integrated luminosity of 9 fb −1 . The two J/ψ candidates in a pair are reconstructed through the J/ψ → µ + µ − decay, and are labelled randomly as either J/ψ 1 or J/ψ 2 .

Detector and data set
The LHCb detector is designed to study particles containing b or c quarks at the LHC. It is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, described in detail in Refs. [57,58]. The online event selection is performed by a trigger, 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. At the hardware stage, events are required to have at least one muon with high momentum transverse to the beamline, p T . At the software stage, two oppositely charged muon candidates are required to have high p T and to form a common vertex. Events are retained if there is at least one J/ψ candidate selected by both the hardware and software trigger stages. Imperfections in the description of the magnetic field and misalignment of subdetectors lead to a bias in the momentum measurement of charged particles, which is calibrated using reconstructed J/ψ and B + mesons [59], with well-known masses.
Simulated samples are used to model the signal properties, including the invariant mass resolution and the reconstruction efficiency. In the simulation, pp collisions are generated using Pythia [60] with a specific LHCb configuration [61]. Decays of unstable particles are described by EvtGen [62], in which final-state radiation is generated using Photos [63]. The interaction of the generated particles with the detector and its response are implemented using the Geant4 toolkit [64], as described in Ref. [65].

Candidate selection
In the offline selection, two pairs of oppositely charged muon candidate tracks are reconstructed, with each pair forming a vertex of a J/ψ candidate. Each muon track must have p T > 0.65 GeV/c and momentum p > 6 GeV/c. The J/ψ candidates are required to have a dimuon invariant mass in the range 3.0 < M µµ < 3.2 GeV/c 2 . A kinematic fit is performed for each J/ψ candidate constraining its vertex to coincide with a primary pp collision vertex (PV) [66]. The requirement of a good kinematic fit quality strongly suppresses the contamination of di-J/ψ candidates stemming from feed-down of b-hadrons, which decay at displaced vertices. The four muon tracks in a J/ψ-pair candidate are required to originate from the same PV, reducing to a negligible level the number of pile-up candidates with the two J/ψ candidates produced in separated pp collisions. Fake di-J/ψ candidates, comprising two muon-track candidates reconstructed from the same real particle, are rejected by requiring muons of the same charge to have trajectories separated by an angle inconsistent with zero. For events with more than one reconstructed di-J/ψ candidate, accounting for about 0.8% of the total sample, only one pair is randomly chosen.  The di-J/ψ signal yield is extracted by performing an extended unbinned maximumlikelihood fit to the two-dimensional distribution of J/ψ 1 and J/ψ 2 invariant masses, µµ ), as displayed in Fig. 1, where projections of the data and the fit result are shown. For both J/ψ candidates, the signal mass shape is modelled by a Gaussian kernel with power-law tails [67]. Each component of combinatorial background, consisting of random combinations of muon tracks, is described by an exponential function. The total di-J/ψ signal yield is measured to be (33.57 ± 0.23) × 10 3 , where the uncertainty is statistical.
The di-J/ψ transverse momentum (p di-J/ψ T ) in SPS production is expected to be, on average, higher than that in DPS [48]. The high-p di-J/ψ T region is thus exploited to select a data sample with enhanced SPS production, which could include contributions from T cccc states. Two different approaches are applied. In the first approach (denoted as p di-J/ψ T -threshold), J/ψ-pair candidates are selected with the requirement p di-J/ψ T > 5.2 GeV/c, which maximises the statistical significance of the SPS signal yield. In the second approach (denoted as p di-J/ψ T -binned), di-J/ψ candidates are categorised into six p di-J/ψ T intervals with boundaries {0, 5, 6, 8, 9.5, 12, 50} GeV/c, defined to obtain equally populated bins of SPS signal events in the M di-J/ψ range between 6.2 and 7.4 GeV/c 2 . This mass region covers the predicted masses of T cccc states decaying into a J/ψ pair. For both scenarios, the DPS yield in the T cccc signal region is extrapolated from the high-M di-J/ψ region using the wide-range distribution constructed from available double-differential J/ψ cross-sections [68][69][70]. The high-M di-J/ψ region is chosen such that the SPS yield is negligible compared to DPS. The SPS yield is obtained by subtracting the DPS contribution from the total number of J/ψ-pair signals. Fig. 2. The di-J/ψ mass is calculated by constraining the reconstructed mass of each J/ψ candidate to its known value [71]. The spectrum shows a broad structure just above twice the J/ψ mass threshold ranging from 6.2 to 6.8 GeV/c 2 (dubbed threshold enhancement in the following) and a narrower structure at about 6.9 GeV/c 2 , referred to hereafter as X(6900). There is also a hint of another structure around 7.2 GeV/c 2 , whereas there are no evident structures at higher invariant mass. Several cross-checks are performed to investigate the origin of these structures and to exclude that they are experimental artifacts. The threshold enhancement and the X(6900) structure become more pronounced in higher p di-J/ψ T intervals, and they are present in subsamples split according to different beam or detector configurations for data collection. The structures are not caused by the experimental efficiency, since the efficiency variation across the whole M di-J/ψ range is found to be marginal. Residual background, in which a muon track is reused or at least one J/ψ candidate is produced from a b-hadron decay, is observed to have no structure. The possible contribution of J/ψ pairs from Υ decays is estimated to be negligible and distributed uniformly in the M di-J/ψ distribution. In Fig. 2, the M di-J/ψ distribution for background pairs with M (1),(2) µµ in the range 3.00 − 3.05 GeV/c 2 or 3.15 − 3.20 GeV/c 2 is also shown, with the yield normalised by interpolating the background into the J/ψ signal region, which accounts for around 15% of the total candidates. There is no evidence of structures in the M di-J/ψ distribution of background candidates.

Investigation of the J/ψ-pair invariant mass spectrum
To remove background pairs that have at least one background J/ψ candidate, the sPlot weighting method [72] is applied, where the weights are calculated from the fit to the two-dimensional (M µµ ) distribution. The background-subtracted di-J/ψ spectra in intervals, which are investigated by weighted unbinned maximum-likelihood fits [73]. The M di-J/ψ distribution of signal events is expected to be dominated by the sum of the nonresonant SPS (NRSPS) and DPS production, which have smooth shapes (referred to as continuum in the following). The DPS continuum is described by an empirical function and its yield determined by extrapolation from the M di-J/ψ > 12 GeV/c 2 region, which is dominated and well described by the DPS distribution. The continuum NRSPS is modelled by a two-body phase-space distribution multiplied by an exponential function determined from the data. The combination of continuum NRSPS and DPS does not provide a good description of the data, as is illustrated in Fig. 3(a). The M di-J/ψ spectrum in the data is tested against the hypothesis that only NRSPS and DPS components are present in the range 6.2 < M di-J/ψ < 7.4 GeV/c 2 (null hypothesis) using a χ 2 test statistic. Pseudoexperiments are generated and fitted according to the null hypothesis, and the fraction of these fits with a χ 2 value exceeding that in the data is converted into a significance. Considering the sample in the p di-J/ψ T > 5.2 GeV/c region, the null hypothesis is inconsistent with the data at 3.4 standard deviations (σ). A test performed simultaneously in the aforementioned six p di-J/ψ T regions yields a discrepancy of 6.0 σ with the null hypothesis. A higher value is obtained in the latter case as more detailed information on the p The structures in the M di-J/ψ distribution can have various interpretations. There may be one or more resonant states T cccc decaying directly into a pair of J/ψ mesons, or T cccc states decaying into a pair of J/ψ mesons through feed-down of heavier quarkonia, for example T cccc → χ c (→ J/ψγ)J/ψ where the photon escapes detection. In the latter case, such a state would be expected to peak at a lower M di-J/ψ position, close to the di-J/ψ mass threshold, and its structure would be broader compared to that from a direct decay. This feed-down is unlikely an explanation for the narrow X(6900) structure. Rescattering of two charmonium states produced by SPS close to their mass threshold may also generate a narrow structure [74][75][76][77]. The two thresholds close to the X(6900) structure could be formed by χ c0 χ c0 pairs at 6829.4 MeV/c 2 and χ c1 χ c0 pairs at 6925.4 MeV/c 2 , respectively. Whereas a resonance is often described by a relativistic Breit-Wigner (BW) function [71], the lineshape of a structure with rescattering effects taken into account is more complex. In principle, resonant production can interfere with NRSPS of the same spin-parity quantum numbers (J PC ), resulting in a coherent sum of the two components and thus a modification of the total M di-J/ψ distribution.
Two different models of the structure lineshape providing a reasonable description of the data are investigated. The X(6900) lineshape parameters and yields are derived from fits to the p di-J/ψ T -threshold sample. Simultaneous p di-J/ψ T -binned fits are also performed as a cross-check and the variation of lineshape parameters is considered as a source of systematic uncertainties. Due to its low significance, the structure around 7.2 GeV/c 2 has been neglected.
In model I, the X(6900) structure is considered as a resonance, whereas the threshold enhancement is described through a superposition of two resonances. The lineshapes of these resonances are described by S-wave relativistic BW functions multiplied by a two-body phase-space distribution. The experimental resolution on M di-J/ψ is below 5 MeV/c 2 over the full mass range and negligible compared to the widths of the structures. The projections of the p di-J/ψ T -threshold fit using this model are shown in Fig. 3(b). The mass, natural width and yield are determined to be m[X(6900)] = 6905 ± 11 MeV/c 2 , Γ[X(6900)] = 80 ± 19 MeV and N sig = 252 ± 63, where biases on the statistical uncertainties have been corrected using a bootstrap method [78]. The goodness of fit is studied using a χ 2 test statistic and found to be χ 2 /ndof = 112.7/89, corresponding to a probability of 4.6%. The fit is also performed assuming the threshold enhancement as due to a single wide resonance (see Supplementary Material); the fit quality is found significantly poorer and thus this model is not further investigated.
A comparison between the best fit result of model I and the data reveals a tension around 6.75 GeV/c 2 , where the data shows a dip. In an attempt to describe the dip, model II allows for interference between the NRSPS component and a resonance for the threshold enhancement. The coherent sum of the two components is defined as where A and φ are the magnitude and phase of the nonresonant component, relative to the BW lineshape for the resonance, assumed to be independent of M di-J/ψ , and f nr (M di-J/ψ ) is an exponential function. The interference term in Eq. (1) is then added incoherently to the BW function describing the X(6900) structure and the DPS description. The fit to the p di-J/ψ T -threshold sample with this model has a probability of 15.5% (χ 2 /ndf = 104.7/91), and its projections are illustrated in Fig. 3(c). In this case, the mass, natural width and yield are determined to be m[X(6900)] = 6886 ± 11 MeV/c 2 , Γ[X(6900)] = 168 ± 33 MeV and N sig = 784 ± 148. A larger X(6900) width and yield are preferred in comparison to model I. Here it is assumed that the whole NRSPS production is involved in the interference with the lower-mass resonance despite that there may be several components with different quantum numbers in the NRSPS and more than one resonance in the threshold enhancement. Material. An additional model describing the dip and the X(6900) structure simultaneously by using the interference between the NRSPS and a BW resonance around 6.9 GeV/c 2 is also considered, however the fit quality is significantly poorer, as illustrated in the Supplementary Material. Alternative lineshapes, other than the BW, may also be possible to describe these structures and will be the subject of future studies.
The increase of the likelihood between the fits with or without considering the X(6900)

Structure
Significance p di-J/ψ T -threshold p di-J/ψ T -binned Any structure beyond NRSPS plus DPS 3.4 σ 6.0 σ Threshold enhancement plus X(6900) 6.4 σ 6.9 σ Threshold enhancement 6.0 σ 6.5 σ X(6900) 5.1 σ 5.4 σ NRSPS and DPS continuum. The significance is determined from both p di-J/ψ T -threshold and p di-J/ψ T -binned fits, and summarised in Table 1. The results are above 5 σ for the two structures, with slightly higher significance for the p di-J/ψ T -binned case. Systematic uncertainties on the measurements of the mass and natural width of the X(6900) structure are reported in Table 2. They include variations of the results obtained by: including an explicit component in the M di-J/ψ fits for the J/ψ combinatorial background rather than subtracting it using the weighting method (sPlot weights in Table 2); convolving the M di-J/ψ fit functions with a Gaussian function of 5 MeV/c 2 width to account for the invariant mass resolution (Experimental resolution); modelling the threshold structure using an alternative Gaussian function with asymmetric power-law tails, or fitting in a reduced M di-J/ψ range excluding the threshold structure (Threshold structure shape); using alternative functions to describe the NRSPS component and varying the DPS yield (NRSPS plus DPS modelling); allowing the relative phase in the NRSPS component to vary linearly with M di-J/ψ (NRSPS phase) using an alternative P -wave BW function for the X(6900) structure and varying the hadron radius in the BW function from 2 to 5 GeV −1 [X(6900) shape]; obtaining results from a simultaneous fit to the M di-J/ψ distributions in the six p di-J/ψ T bins (Cut on p di-J/ψ T ); including an explicit contribution for J/ψ mesons from b-hadron feed-down (b-hadron feed-down) or adding a BW component for the 7.2 GeV/c 2 structure (Structure at 7.2 GeV/c 2 ). The total uncertainties are determined to be 7 MeV/c 2 and 33 MeV for the mass and natural width, respectively, without considering any interference, and 11 MeV/c 2 and 69 MeV when the interference between NRSPS and the threshold structure is introduced.
For the scenario without interference, the production cross-section of the X(6900) structure relative to that of all J/ψ pairs (inclusive), times the branching fraction B(X(6900) → J/ψJ/ψ), R, is determined in the pp collision data at √ s = 13 TeV. The measurement is obtained for both J/ψ mesons in the fiducial region of transverse momentum below 10 GeV/c and rapidity between 2.0 and 4.5. An event-by-event efficiency correction is performed to obtain the signal yield at production. The residual contamination from b-hadron feed-down is subtracted from inclusive J/ψ-pair production following Ref. [70]. The systematic uncertainties on the X(6900) yield are estimated in a similar way to that for the mass and natural width, while other systematic uncertainties mostly cancel in the ratio. The production ratio is measured to be R = [1.1 ± 0.4 (stat) ± 0.3 (syst)]% without any p di-J/ψ T requirement and R = [2.6 ± 0.6 (stat) ± 0.8 (syst)]% for p di-J/ψ T > 5.2 GeV/c.

Summary
In conclusion, using pp collision data at centre-of-mass energies of 7, 8 and 13 TeV collected with the LHCb detector, corresponding to an integrated luminosity of 9 fb −1 , the J/ψ-pair invariant mass spectrum is studied. The data in the mass range between 6.2 and 7.4 GeV/c 2 are found to be inconsistent with the hypothesis of NRSPS plus DPS continuum. A narrow structure, X(6900), matching the lineshape of a resonance and a broad structure next to the di-J/ψ mass threshold are found. The global significance of either the broad or the X(6900) structure is determined to be larger than five standard deviations. Describing the X(6900) structure with a Breit-Wigner lineshape, its mass and natural width are The X(6900) structure could originate from a hadron state consisting of four charm quarks, T cccc , predicted in various tetraquark models. The broad structure close to the di-J/ψ mass threshold could be due to a mixture of multiple four-charm states or have contributions from feed-down of four-charm states through heavier quarkonia. Other interpretations cannot presently be ruled out, for example the rescattering of two charmonium states produced close to their mass threshold. More data along with additional measurements, including determination of the spin-parity quantum numbers and p T dependence of the production cross-section, are needed to provide further information about the nature of the observed structure.

Supplementary Material
In the Supplementary Material, the J/ψ-pair mass distributions in bins of p di-J/ψ T are shown in Sec. A, the fits using several additional models to the J/ψ-pair mass spectrum are presented in Sec. B, and some supplemental information to the fit result of model II is given in Sec. C.
A J/ψ-pair mass distributions in bins of p di-J/ψ T 7000 8000 9000 B Additional fits to the J/ψ-pair mass spectrum Figure 6 shows the fits to the J/ψ-pair mass spectrum with (a) the threshold structure described by a single Breit-Wigner (BW) lineshape and (b) using a model that contains a single BW resonance interfering with the SPS continuum. The χ 2 /ndof of the two fits are 125.6/92 and 118.6/91, corresponding to a probability of 1.2% and 2.8%, respectively. Figure 7 shows the fit with an additional BW function introduced to describe the 7.2 GeV/c 2 structure, based on the model that contains two BW lineshapes for the threshold structure and a BW shape for the X(6900) structure on top of the NRSPS plus DPS continuum.  -threshold fit with an additional BW function introduced to describe the 7.2 GeV/c 2 structure, based on the model that contains two BW lineshapes for the threshold structure and a BW shape for the X(6900) structure on top of the NRSPS plus DPS continuum.

C Supplement to fit result of model II
In model II that contains a BW lineshape for the threshold structure interfering with the NRSPS, a BW shape for the X(6900) structure and the DPS continuum, the parameters of the lower-mass BW lineshape is determined to M = 6741 ± 6 (stat) MeV/c 2 and Γ = 288 ± 16 (stat) MeV. The systematic uncertainties on the mass and natural width are not studied. Due to the complex nature of the threshold structure, and the simple interference scenario considered, this study is not considered to claim a state with the parameters reported here.
Projections of the fit to the J/ψ-pair invariant mass spectra in bins of p di-J/ψ T assuming the interference between the threshold structure and the SPS continuum are shown in Fig. 8.   [31] Y. Chen and R. Vega-Morales, Golden probe of the di−Υ threshold, arXiv:1710.02738.