Observation of an Excess of Dicharmonium Events in the Four-Muon Final State with the ATLAS Detector

A search is made for potential cc ¯ c ¯ c tetraquarks decaying into a pair of charmonium states in the four muon final state using proton-proton collision data at ﬃﬃﬃ s p ¼ 13 TeV, corresponding to an integrated luminosity of 140 fb − 1 recorded by the ATLAS experiment at LHC. Two decay channels, J= ψ þ J= ψ → 4 μ and J= ψ þ ψ ð 2 S Þ → 4 μ , are studied. Backgrounds are estimated based on a hybrid approach involving Monte Carlo simulations and data-driven methods. Statistically significant excesses with respect to backgrounds dominated by the single parton scattering are seen in the di-J= ψ channel consistent with a narrow resonance at 6.9 GeVand a broader structure at lower mass. A statistically significant excess is also seen in the J= ψ þ ψ ð 2 S Þ channel. The fitted masses and decay widths of the structures are reported.

Beyond the conventional mesons ( q) and baryons ( or q q q), exotic hadrons composed of four ( q q) or five quarks ( q) are also allowed under color confinement.The  (3872) particle discovered by Belle in 2003 was the first tetraquark (TQ) candidate [1], and was followed by a series of further candidates designated as X, Y, and Z states [2].In 2020, LHCb observed a narrow  (6900) structure in the di-/ channel [3].The structure could be interpreted as a tetraquark with four charm quarks,   c c [4][5][6][7][8][9][10][11].An additional enhancement closer to the di-/ mass threshold was also observed in the LHCb data.Since the 6.9 GeV LHCb resonance is above the /+(2S) mass threshold, a structure in the /+(2S) channel is also possible.Both channels are investigated by ATLAS in a quite different phase space region from LHCb, and the new channel of /+(2S) provides more information for di-charmonium excesses.For example in some predictions, the two channels are coupled via Pomeron exchange between the two charmonia, and  (6900) is dynamically produced [12].
A search in the 4 final state produced through the di-/ and /+(2S) channels is carried out, using 140 fb −1 of LHC proton-proton ( ) data collected by the ATLAS experiment at a center-of-mass energy of √  = 13 TeV between the years 2015 and 2018.Only the data where all detector systems are functional and recording high-quality data are used.The ATLAS detector [13] covers nearly the entire solid angle around the collision point 1 with layered tracking detectors, calorimeters and muon chambers.The muon and tracking systems are of particular importance in the reconstruction of charmonia.The inner tracking detector (ID) consists of a silicon pixel detector, a silicon microstrip detector and a transition radiation tracker.The muon spectrometer (MS) surrounds the calorimeters and consists of three large superconducting air-core toroids with eight coils each, a system of tracking chambers, and detectors for triggering.Muons are reconstructed using information from the ID and MS systems.
The data sample was collected with triggers requiring either two muons with invariant mass compatible with / or (2S) mesons (mass in the range of [2.5, 4.3] GeV), or three muons containing at least one such di-muon pair [33,34].Combinations of triggers with different prescales [35] depending on the run period are used to give the largest acceptance.The trigger efficiency for the  (6900) relative to the offline selection is about 72%.Dedicated ATLAS offline software is used to reconstruct the charmonium and 4 candidates in each event recorded by the triggers.In each event containing at least four muons with two opposite-charge pairs, the ID tracks are fit to a common vertex.Afterwards, each vertex of the two pairs is refit with a / or (2S) mass constraint [36].The resolution of the TQ mass with these mass constraints ( 4 ) is about 0.33% for  (6900).
The loose identification selection criteria [37] are required for all muon candidates.Depending on the muon trigger thresholds and muon identification requirements, different muon momenta on the four muons are required.Several requirements are imposed on the following variables to further suppress the background: The SPS and DPS backgrounds contain two prompt charmonia and are modeled by MC simulations.Because the event generator does not reproduce the data distributions well, kinematic corrections are derived from two dedicated control regions.Since SPS (DPS) events are characterized by two charmonia which are nearby (distant) in  −  space, the control region is defined with a 4 mass sideband within [7.5, 12] GeV ([14, 24.5] GeV) without the Δ requirement.The corrections are implemented by assigning event weights to MC simulations such that distributions of kinematic variables such as di-/  T , Δ and Δ between charmonia, and the lowest muon  T match the data in the control regions.The corrections are then applied to all mass regions and the SPS modelling is validated in the Δ > 0.25 control region.
The non-prompt background also contains two charmonia, albeit originating from -hadron decays.These typically contain a decay vertex that is displaced from the primary   interaction.This background is also modeled using MC simulation, but normalized and validated by dedicated control regions obtained by reversing the vertex quality requirements as shown in Table 1.Events from prompt single charmonium production and non-resonant di-muon production are collectively called "Others", and have at least one charmonium candidate containing random combinations of mostly fake muons.Fake muons are tracks, typically charged hadrons, that are misidentified as muons.A data-driven method is used because MC simulations do not accurately estimate this kind of background.A fake muon control region from data is used to model the Others background, which is defined by requiring that one charmonium candidate contains a track that is not reconstructed as a muon candidate, with all the other requirements kept unchanged.Events in the charmonium mass sidebands are used for both the normalization and shape corrections for Others.
Unbinned maximum likelihood fits are performed to extract the signal information from data in the 4 mass spectra.The likelihood used for the fit is where L  (L   ) is the likelihood in the signal (control) region, ì  are the parameters of interest,   are nuisance parameters (NP) which account for systematic uncertainties shared between the two regions.Each NP has a Gaussian distribution constraint with a subsidiary measurement  ′  , a mean   and a width set to   = 1 by construction.Only the background yields in the control regions are used in simultaneous fits with the signal regions.Background yields in the two regions are related by a transfer factor obtained from MC predictions and data-driven estimations, with systematic variations for both components.
In the di-/ channel, the feed-down normalized by Eq. 1 is included as an additional background, and two fit models are considered.In model A, the signal probability density function in L  consists of three interfering S-wave Breit-Wigner (BW) resonances multiplied with a phase space factor and convolved with a mass resolution function, which gives where   (Γ  ()) are the masses (widths) of resonances,   are complex numbers representing the relative magnitudes and phases ( 1 is fixed to unity with zero phase for this purpose), Γ  () = Γ        , where  (  ) is the momentum of one charmonium in the rest frame of the di-charmonium system at the invariant mass equal to  (  ) [38], and  is the mass resolution function.The   terms are ordered by the subscripts.In model B, two resonances are considered.The first one interferes with the SPS background, while the second is standalone.The signal+SPS probability distribution function gives where () and  are the SPS background amplitude and phase relative to the resonance at  0 (| ()| 2 reproduces the non-interfering SPS background from the MC prediction).In this model, the control region becomes irrelevant and is excluded from the likelihood given in Eq. 2.
Models A and B are analogous to models I and II of the LHCb study [3], respectively.However, interferences between the signal resonances are introduced in model A, which is not done in the analysis by LHCb.The number of resonances in model A starts from one and increases to three with the fit quality gradually improving.A 4th resonance is added only for systematics, as the fit quality does not improve appreciably.For comparison, a two-resonance model with interference, and a three-resonance model without interferences, are also tried.It is found that when compared with model A, these models are excluded with a confidence level of more than 95% based on toy MC studies.
In the /+(2S) channel, two fit models are also considered.Model  assumes that the same interfering resonances observed in the di-/ channel also decay into /+(2S), in addition to a standalone fourth resonance in this channel.The signal probability distribution function gives where the parameters of the first three resonances, whose contribution appears as a structure just above the  / +   (2) mass threshold, are fixed to the values from the fit to the di-/ channel.In contrast, model  assumes a single resonance in this channel (i.e., without the  0,1,2 terms in Eq. 5).
The systematic uncertainties are classified into those affecting exclusively normalizations, and those affecting the mass spectrum shape as well.Only the latter are relevant, since the signal and background normalizations are freely floating parameters.The systematic uncertainties in  4 , with and without the muon momentum calibration corrections, are treated as resolution uncertainties.Because the resolution of the  4 is mass dependent and a constant mass is used in the nominal fit, resolutions in different mass ranges are treated as systematic uncertainties.A shape uncertainty is assigned to account for bin-to-bin fluctuations from the limited MC sample size for backgrounds.In the SPS background, a Pythia model parameter uncertainty from pT0timesMPI [28], which controls the suppression of the soft double charmonia production, is assigned, and its nominal value is tuned to data in the SPS control region.A shape uncertainty in the background due to residual di-charmonium  T mismodeling is applied.Based on toy MC studies, biases from the fit in the resonance parameters are also considered as systematic uncertainties.
The P and D-wave BW functions are substituted for the S-wave for resonances away from the threshold to estimate systematic uncertainties due to different orbital angular momentum assumptions 4 .Systematic shape variations in the  (6900) in the di-/ channel, and in the second resonance for the /+(2S) channel due to the Δ and muon  T requirements are considered as well.In the di-/ channel, a 4th resonance around 7.2 GeV (hinted by the LHCb analysis) is added to the fit, and the feed-down background normalizations are varied according to the uncertainties in /+(2S).The transfer factor uncertainty is dominated by the SPS model parameter, so it is not treated as a separate NP.In the /+(2S) channel, the uncertainty in a transfer factor between the signal and control regions, and a shape uncertainty derived from the non-prompt region due to Others (shape inconsistency), are included.Interference between the 4th resonance and the other ones are included in systematic uncertainties.In model , systematic uncertainties on the lower resonance shape from the di-/ channel model A fit are also included.Other systematic uncertainties such as the parton PDF and Pythia parameters affect signal and background normalizations only, and are not incorporated in the fits.
The 4 mass spectra fit to data in the two channels are shown in Figure 1.The fitted masses and widths of resonances are given in Table 2.Both the significance of all resonances, and the one for  (6900) alone, far exceed 5 5 .The mass of the third resonance,  2 , is consistent with the LHCb mass.Although both the models A and B describe the data well, the broad structure at the lower mass could result from other physical effects, such as the feed-down from higher di-charmonium resonances, e.g.,   c c →     ′ → // where the soft photons are not reconstructed.In the /+(2S) channel, the signal significance with signal shape parameters of model  () fixed to their best-fit values is 4.7 (4.3).In the fit with model , the significance of the second resonance alone is found to be 3.0.±20% ± 12% 5 The asymptotic formula based on the profile likelihood ratio,  = √︂ 2 ln , is used to calculate the overall significance, where  is the signal yield and  are NPs [39].Similarly for  (6900) alone,  = √︂ 2 ln is used.In the calculations, the signal shape parameters are all fixed to their best-fit values.In conclusion, the results of a search for potential  c c tetraquarks decaying into a pair of / charmonium states, or into a / and (2S), in the 4 final state are presented based on   collisions data collected by the ATLAS experiment at √  = 13 TeV corresponding to an integrated luminosity of 140 fb −1 .A significant excess of events (far exceeding 5) in data above the expected background is observed in the di-/ channel.Analogous to LHCb observations, a broad structure at lower mass and a resonance around 6.9 GeV are observed.A three-resonance model with interferences, or a model with the lower broad structure interfering with the SPS background, describes the excess better than models with fewer interfering resonances or with no interferences.In the /+(2S) channel, a 4.7 excess of events is observed when considering a model involving two resonances, one of which is near the 6.9 GeV threshold.In both channels, details of the lower-mass structure cannot be discerned directly from the data, and other interpretations (e.g.multiple non-interfering resonances, reflection effects and threshold enhancements) cannot be excluded.More data are required to better characterize the excesses observed in both channels.
We thank CERN for the very successful operation of the LHC, as well as the support staff from our institutions without whom ATLAS could not be operated efficiently.

Appendix
Signal events are simulated with the event generator JHU [41] and CTEQ6L1 PDF [42], or with Pythia and the NNPDF23LO PDF.Feed-down backgrounds from the /+(2S) channel to di-/ are included.A natural width of 100 MeV is assumed for all the resonances with no interference between them.An extensive software suite [43] is used in data simulation [44], in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.The MC simulated events are weighted to reproduce the same number of   interactions per bunch crossing (pileup) and trigger conditions as occur in data.
The / and 4 mass distributions of data and predictions in the signal regions of the two channels before the fits are shown in Figure 2(a,b,c).Similar structures were also observed by CMS [45].Similar 4 mass distributions in the control regions are shown in Figure 3. Distributions in the SPS and DPS control regions without the Δ requirement between the charmonia are shown in Figure 4-5.The systematic uncertainties in the fitted masses and widths of the highest resonances in models A and  of the two channels are summarized in Table 3.

(c)
Figure 2: The / mass spectrum with 6.7 GeV <  4 <7.1 GeV (a) and the 4 mass spectrum (b) in the signal region in the di-/ channel, and the similar mass spectrum in the /+(2S) channel (c).The signal from the  (6900) is scaled to match data around 6.9 GeV.The bars and shaded areas represent uncertainties of data and predictions in each bin, respectively.The arrows in the lower panel indicate that the ratio of data to prediction is out of range in that bin.

(c)
Figure 4: The 4 mass spectrum within [7.5, 24.5] GeV and without the Δ requirement (a),  T of the di-charmonium in the SPS control region with 7.5 GeV <  4 < 12.0 GeV (b), and Δ between the charmonia in the DPS control region with 14.0 GeV <  4 < 24.5 GeV (c), in the di-/ channel.

Figure 1 :
Figure 1: The fit to the mass spectra in the signal regions in the di-/ (a,b) and /+(2S) (c,d) channels.Fit results for models A (a), B (b),  (c) and  (d) are shown.The purple dash-dotted lines represent the components of individual resonances, and the green short dashed ones represent the interferences among them.

Figure 3 :
Figure 3: The 4 mass spectra in the control regions with Δ ≥ 0.25 in the di-/ (a) and /+(2S) (b) channels.The bars and shaded areas represent uncertainties of data and predictions in each bin, respectively.The arrows in the lower panel indicate that the ratio of data to prediction is out of range in that bin.

Table 2 :
The fitted masses and natural widths (in GeV), and relative uncertainties of signal yields (Δ/) in the di-/ and /+(2S) channels.The results of both the models are given in each channel.The first uncertainties are statistical while the second ones are systematic.

Table 3 :
Different sources of systematic uncertainty in the mass and natural width (in MeV) of the third (second) resonance in model A () of the di-/ (/+(2S)) channel.