Evidence for X(3872) in PbPb collisions and studies of its prompt production at $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV

The first evidence for X(3872) production in relativistic heavy ion collisions is reported. The X(3872) production is studied in lead-lead (PbPb) collisions at a center-of-mass energy of $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV per nucleon pair, using the decay chain X(3872) $\to$ J$/\psi\, \pi^+\pi^- \to$ $\mu^+\mu^-\pi^+\pi^-$. The data were recorded with the CMS detector in 2018 and correspond to an integrated luminosity of 1.7 nb$^{-1}$. The measurement is performed in the rapidity and transverse momentum ranges $|y|$ $\lt$ 1.6 and 15 $\lt$ $p_\mathrm{T}$ $\lt$ 50 GeV$/c$. The significance of the inclusive X(3872) signal is 4.2 standard deviations. The prompt X(3872) to $\psi$(2S) yield ratio is found to be $\rho^\mathrm{PbPb} = $ 1.08 $\pm$ 0.49 (stat) $\pm$ 0.52 (syst), to be compared with typical values of 0.1 for pp collisions. This result provides a unique experimental input to theoretical models of the X(3872) production mechanism, and of the nature of this exotic state.


1
The X(3872), also known as χ c1 (3872), is an exotic particle that was first observed by the Belle Collaboration [1], and then confirmed and studied by other experiments at electronpositron [2, 3] and hadron colliders [4][5][6][7][8]. The quantum numbers of the X(3872) have been narrowed down by the CDF [9], and later determined to be J PC = 1 ++ by the LHCb [10] Collaborations. However, the nature of this particle is still not fully understood and interpretations in terms of conventional charmonium (a bound state of charm-anticharm quarks), D * (2010) 0 D 0 molecules [11], tetraquark states [12], and their admixture [13] have been proposed. The production and survival of the X(3872) in a quark-gluon plasma (QGP), a deconfined state of quarks and gluons [14,15], or after the QGP, in a hadronic phase, is expected to depend upon the X(3872)'s internal structure [16,17]. Thus, the recent large data set of lead-lead (PbPb) collisions at a center-of-mass energy of √ s NN = 5.02 TeV per nucleon pair, delivered by the Large Hadron Collider (LHC) at CERN at the end of 2018, opened new opportunities to probe the nature of this exotic state [18][19][20].
It is expected that in relativistic heavy ion collisions, the formation of the QGP could enhance or suppress the production of the X(3872) particle. Coalescence mechanisms could enhance the X(3872) production yield [16,19]. These mechanisms can be modeled via the overlap of the density matrix of the constituents in an emission source with the Wigner function of the produced particle [21]. Therefore, the enhancement of the X(3872) production in the QGP would depend on the spatial configuration of the exotic state. Moreover, a longer distance between the quarks and antiquarks that constitute the state could also lead to a higher X(3872) dissociation rate, similar to that from the mechanism of quarkonium suppression in heavy ion collisions [22]. Therefore, the study of the X(3872) state in the QGP may be used as a tool to distinguish a compact tetraquark configuration with a radius ∼0.3 fm from a molecular state with a radius greater than 1.5 fm [23]. Such a measurement would be complementary to the recent evidence for the radiative decay X(3872) → ψ(2S)γ in proton-proton (pp) collisions reported by LHCb Collaboration [24], which does not support a pure D * (2010) 0 D molecular interpretation. In addition, measurements of prompt X(3872) production could provide an interesting test of the statistical hadronization model, which assumes that the produced matter is in thermodynamic equilibrium at the phase transition to hadrons [25,26].
In this Letter, the first evidence for X(3872) production in PbPb collisions at √ s NN = 5.02 TeV is reported. The PbPb sample corresponds to an integrated luminosity of 1.7 nb −1 . The X(3872) candidates are reconstructed through the decay chain X(3872) → J/ψ π + π − → µ + µ − π + π − , and are measured in the 15 < p T < 50 GeV/c and |y| < 1.6 kinematic region. At the LHC energies, the inclusive X(3872) yields in pp and PbPb collisions contain a significant nonprompt contribution coming from b hadron decays [8]. The nonprompt X(3872) component is related to the medium-modified beauty hadron production in heavy ion collisions, which is out of the scope of this paper. Here, we focus on the prompt component, from charm quark fragmentation, for which the ratio ρ i (i is pp or PbPb) between the corrected yields of X(3872) and ψ(2S) mesons (where the ψ(2S) is reconstructed with the same final-state particles in order to reduce systematic uncertainties) is presented: The ratios in pp and PbPb collisions are connected to the nuclear modification factors R X(3872) AA and R ψ(2S) AA (the meson yield ratio in nucleus-nucleus and pp interactions normalized by the number of inelastic nucleon-nucleon collisions) via the following relation: ( Events of interest were selected in real-time using the CMS two-tiered trigger system: the first level (L1), composed of custom hardware processors [28], and the high-level trigger (HLT), consisting of a farm of processors running a version of the full event reconstruction software optimized for fast processing [29]. The selection required the presence of two muon candidates, with at least one muon reconstructed in the outer muon spectrometer, and one muon reconstructed using information from both the outer muon spectrometer and the inner tracker. The dimuon candidate invariant mass is required to be 1 < m µµ < 5 GeV/c 2 . For the offline analysis, events have to pass a set of selection criteria designed to reject events from background processes (beam-gas collisions, beam scraping events and electromagnetic interactions) as described in Ref. [36]. Events are required to have at least one reconstructed primary interaction vertex formed by two or more tracks, with a distance from the center of the nominal interaction point of less than 15 cm along the beam axis. The shapes of the clusters in the pixel detector have to be compatible with those expected from particles produced by a PbPb collision [37]. In order to select hadronic collisions, the PbPb events are also required to have at least two towers (i.e., a geometrically defined group of calorimeter cells) in each of the HF detectors with total energy deposits of more than 4 GeV per tower. This analysis is restricted to events within centrality 0-90%, for which the hadronic event selection is fully efficient. Multiple-collision events (pileup) have a negligible effect on the measurement, since the average number of additional collisions per bunch crossing is approximately 0.002.
Dedicated PbPb X(3872) and ψ(2S) Monte Carlo (MC) simulated samples were generated in order to estimate the acceptance and selection efficiencies, to study the background components, and to evaluate systematic uncertainties. The PYTHIA 8 v212 [38] Tune CP5 [39] was used to generate the X(3872) and ψ(2S) signals at √ s NN = 5.02 TeV. It was assumed that X(3872) and ψ(2S) are unpolarized. Since the X(3872) cannot be generated by PYTHIA, the χ c1 (1P) particle is used, with a modified mass of 3.8716 GeV/c 2 [40]. The χ c1 (1P) has the quantum numbers J PC = 1 ++ , identical to those of the X(3872). The X(3872) particle was forced to decay into J/ψ π + π − (assuming the ρ resonance dominates the pion pair spectrum [7, 8]), followed by the J/ψ meson decaying into two muons. Final-state radiation was generated using PHOTOS 2.0 [41]. The χ c1 (1P) → J/ψ π + π − decay is generated with EVTGEN. The samples with prompt (fragmenting in charm quarks) and nonprompt (originating from b hadron decays) ψ(2S)/X(3872) production are generated separately. Each PYTHIA event is embedded in a PbPb collision event generated with HYDJET 1.8 [42], which is tuned to reproduce global event properties such as the charged-hadron p T spectrum and particle multiplicity. The X(3872) signal is extracted in the following steps. Each muon candidate must be matched to a triggered muon and have p µ T > 3.5 GeV/c in the interval |η µ | < 1.2, p µ T > (5.47-1.89 |η µ |) GeV/c in the interval 1.2 < |η µ | < 2.1, or p µ T > 1.5 GeV/c in the forward region 2.1 < |η µ | < 2.4. Two muons of opposite sign, with an invariant mass within ±150 MeV/c 2 of the world-average J/ψ meson mass [40] are selected to reconstruct a J/ψ candidate. The opposite-sign muon pairs are fitted with a common vertex constraint and are kept if the χ 2 probability of the fit is greater than 1%, thus reducing the background from charm and beauty hadron semileptonic decays. The X(3872) and ψ(2S) candidates are built by combining the J/ψ candidates with two additional tracks, which have p T > 0.9 GeV/c, |η| < 2.4 and pass a high purity selection, assumed to be produced by two pions. Only candidates that have 15 < p T < 50 GeV/c and |y| < 1.6 are considered. Then, a kinematic fit to the J/ψ π + π − system is performed, requiring that the four tracks originate from a common vertex and forcing the mass of the dimuon pair to be equal to the nominal J/ψ mass [40]. The selection is further optimized, using a boosted decision tree (BDT) algorithm [43]. The X(3872) decay vertex probability, the radial distance between the pion and the J/ψ candidate momentum vectors ( where φ is the azimuthal angle), and the p T of each pion are used in the BDT algorithm to distinguish signal and combinatorial background formed by random combinations of tracks. For the multivariate training, the X(3872) signal sample is taken from simulation, and the background sample consists of data from the sidebands of the X(3872) meson peak (i.e., 0.07 < |m − m X(3872) | < 0.128 GeV/c 2 ). The nominal BDT selection is chosen in order to maximize the statistical significance of X(3872), defined as S/ (S + B) where S and B are the numbers of X(3872) and background in the signal region, respectively. Because there is no reliable theoretical calculation available for the X(3872) production, an estimation of the X(3872) cross-section is obtained from data, by extracting the signal yield in data when applying a tight BDT cut. This is used in conjunction with the MC-calculated efficiency to obtain the yield for each BDT cutoff value. The BDT selection determined in this way is applied to entries in the whole invariant mass range from 3.62 to 4 GeV/c 2 in data.
The raw inclusive yields of X(3872) and ψ(2S) are extracted by an extended unbinned maximumlikelihood fit. A double-Gaussian function with a common mean but independent widths is used to model the signal component for each of the X(3872) and ψ(2S) peaks. This was preferred to a single-Gaussian or a Breit-Wigner function since it described better (i.e., superior χ 2 of the fit) the signal shape in MC simulations. For describing the combinatorial background, mostly produced by the random combination of a J/ψ candidate with tracks that are not coming from X(3872) or ψ(2S) decay, a 4th-order polynomial is used, which gives the best fit in terms of χ 2 per degrees of freedom and stability during all studies. For the signal, only the magnitude of the two peaks is left free in the fit, the rest (the mean and widths of the two Gaussian functions, as well as their relative contribution to the signal yield in either X(3872) or ψ(2S) peaks) are set to the values derived from simulation. The five parameters of the combinatorial background are all allowed to float. The invariant mass range considered for the fit is 3.62 to 4 GeV/c 2 . The invariant mass fits for both the inclusive and nonprompt samples, with BDT selection optimized for X(3872), are shown in Fig. 1. The significance of the inclusive X(3872) signal against background-only hypothesis is 4.2 standard deviations. The systematic uncertainty (described below) contributing to this significance is the one related to the X(3872) invariant mass fit. After performing a likelihood scan for each alternative signal and background shape considered, the significance was calculated as the square root of the logarithm of the profile likelihood ratio where the signal is zero, with the smallest value obtained being chosen among all scans.
The contribution from b hadron decays is subtracted from the inclusive result using the The results of the unbinned maximum-likelihood fits for the signal+background and background alone are also shown by the solid and dashed lines, respectively. The pull distribution is represented by the shaded bars. The X(3872) peak mass resolution, σ X(3872) , calculated at the half-maximum of the signal-shape distribution, is also listed for reference.
"pseudo-proper" decay length l xy , defined as the distance in the transverse plane, L xy , between the vertex formed by the 4-tracks and the primary vertex, corrected by the transverse Lorentz boost of the candidate: l xy = L xy m J/ψππ c/| p T |. The prompt-component fraction ( f prompt ) is estimated using a cutoff-based method in the following way. Since the l xy of the prompt component was found in MC studies to be smaller than 0.1 mm, the prompt fraction, f prompt , can be derived from: i) the raw inclusive yield, N incl , obtained from the fit to the invariant mass distributions of all candidates, shown in Fig. 1, and ii) N b-enr , the "benriched" yield, obtained from a fit to the invariant mass distribution only containing candidates that passed the selection l xy > 0.1 mm, also shown in Fig. 1. In addition, N b-enr has to be corrected to account for nonprompt candidates, with l xy < 0.1 mm, that have been missed: f b-enr nonprompt = N nonprompt (l xy > 0.1 mm)/N nonprompt . The correction was obtained from simulation. The raw prompt fraction is then calculated as: separately for X(3872) and ψ(2S) states. The corrected yield of the prompt component can then be derived as: where i is X(3872) or ψ(2S), α is the acceptance, reco is the candidate reconstruction efficiency and sel is the candidate selection efficiency. Since the two states are reconstructed in the same decay channel and are relatively close in mass, their corresponding α reco values are similar.
The choice of the BDT optimization criteria results in sel being higher for the X(3872) than for the ψ(2S).
The measurement of ρ PbPb is affected by several sources of systematic uncertainty, arising from the candidate selection, invariant mass fit, and efficiency corrections. To estimate the systematic uncertainty associated with the BDT selection, the BDT cutoff values are varied within a range that allows a robust invariant mass fit procedure (i.e., signal statistical significance larger than 2), and for each variation all factors in Eq. (4) are recalculated, separately for X(3872) and ψ(2S). The maximum difference of the final ρ PbPb value from the nominal result (40%) is quoted as the systematic uncertainty. The relatively large ρ PbPb uncertainty associated with BDT cutoff value is the convolution of mainly two causes: the BDT variables distribution differences in data and MC samples for the X(3872) meson, and the statistical limitation of the signal in data. The largest differences (∼2 standard deviations) between data and MC samples are in the distributions of the p T of the pions, and the radial distance between the pion and the J/ψ candidate momentum vectors.
The uncertainty in the invariant mass fit (8.0%) is calculated by adding in quadrature the maximum deviations from the nominal result to that found using two alternative fitting functions for both signal and background. For the signal, one variation consists of using a triple-Gaussian function, while for the other the signal width of the nominal fit is allowed to float to account for the resolution difference between data and MC. Other choices for the signal shape (e.g., one-Gaussian function) were not considered because of their poor-quality fits. For the background, the fit function is changed once to a third-order polynomial (as an exponential function or lower-order polynomials could not describe the data), and the fit range is also changed from 3.62-4 GeV/c 2 to 3.62-3.9 GeV/c 2 to exclude the right-hand shoulder.
The efficiency corrections obtained from simulation are sensitive to how well the p T spectrum of the X(3872) and ψ(2S) candidates is modeled. The uncertainty related to the simulated p T shape is evaluated by comparing the reconstruction and selection efficiencies calculated using the default PYTHIA MC sample, with another MC sample in which the p T distributions of X(3872) and ψ(2S) are tuned to reproduce the extracted X(3872) and ψ(2S) p T and y spectra obtained in data, by performing mass fits in bins of X(3872) and ψ(2S) p T and y. The p T and y spectra of the alternative MC samples are allowed to vary within the statistical uncertainties in data. The mean of the differences between efficiencies from the alternative MC samples and the default PYTHIA MC due to the variation of p T and y spectra, which is 13%, is quoted as the systematic uncertainty.
The uncertainties in the trigger efficiency, in the muon reconstruction and identification are evaluated using single muons from J/ψ meson decays in both simulated and collision data, with the tag-and-probe method [44,45]. This combined uncertainty is found to be negligible, below 1%. Scale factors, calculated as the ratio of data to simulated efficiencies as a function of p µ T and η µ , are applied to each dimuon pair on a muon-by-muon basis. The uncertainties of the scaling factors from tag-and-probe studies are quoted as systematic uncertainties.
To estimate the uncertainty in the prompt fraction arising from potential differences between the resolution in data and simulation, a template fit of the l xy distribution in data is performed using prompt and nonprompt l xy templates from simulation. Data are binned in l xy , and an invariant mass fit is performed to extract the inclusive yield in each l xy bin. This backgroundsubtracted l xy distribution is then fitted using a two component fit, which includes the prompt and nonprompt l xy templates from simulation. The widths of the simulated DCA distributions are varied by a floating scale factor, and the best simulated smearing scale factor to match data is determined by minimizing the χ 2 of the two-component fit. The difference between the ratio of the prompt fractions of X(3872) to ψ(2S) using the template fit method and the nominal result (8.1%) is quoted as a systematic uncertainty.
When calculating the uncertainties in the ratio of the acceptance-corrected yields of prompt X(3872) production over prompt ψ(2S) production, the uncertainties of X(3872) and ψ(2S) yields are assumed to be independent except for the systematic uncertainties from muon re-construction, efficiencies, and prompt fractions.
The ratio ρ PbPb between the prompt X(3872) and ψ(2S) mesons is shown in Fig. 2, together with ρ pp measured as a function of p T . The pp data were measured at √ s = 7 and 8 TeV, in the |y| < 1.2 and |y| < 0.75 intervals, respectively [7, 8, 10]. The 7 TeV result was derived using the CMS Collaboration published ratio of the inclusive yields [7] and prompt fractions [7, 46]. From Fig. 2 it is clear that the prompt ρ pp does not depend significantly on collision energy or rapidity. In pp collisions at √ s = 8 TeV, in the kinematic range of 16 < p T < 22 GeV/c, the ρ pp measured by the ATLAS Collaboration is 0.106 ± 0.008 (stat) ± 0.004 (syst) [8]. This is to be compared to the prompt ρ PbPb measured in this Letter, ρ PbPb = 1.08 ± 0.49 (stat) ± 0.52 (syst).
In the interval 15 < p T < 20 GeV/c, the yield of the prompt ψ(2S) in PbPb collisions was reported to be significantly suppressed with respect to pp collisions, R ψ(2S) AA = 0.142 ± 0.061 (stat) ± 0.020 (syst) [47]. This leads, using Eq. (2), to an R X(3872) AA central value larger than 1 (i.e., enhancement of the prompt X(3872) yield in PbPb compared to pp collisions). However, the uncertainties are such that R X(3872) AA is compatible with 1 within ∼1 standard deviation, and with R ψ(2S) AA within ∼2 standard deviations. Thus, it is possible that in PbPb collisions, the prompt X(3872) yield has either no suppression with respect to pp collisions, or as much suppression as the ψ(2S) state. The much larger data sample expected in Run 3 at the LHC will answer whether R X(3872) AA is different from R ψ(2S) AA and significantly above unity. It may answer whether the ψ(2S) meson production (a bound state of a c and c quarks, with r ∼ 0.9 fm) [48], is affected differently by the medium produced in PbPb collisions than the X(3872) state (that could be made of c, c, u, and u quarks, with a radius of r ∼ 0.3 fm or r > 1.5 fm), the difference in both size and quark content playing a role into their production mechanisms. It will also answer whether the X(3872) prompt state production is different in PbPb collisions compared to pp collisions. The question whether X(3872) is a tetraquark or a molecule cannot yet be answered, because of the statistical limitation of the data, and the disagreement among models. For example, while the AMPT transport model [16] predicts R molecule AA >> R tetraquark AA with R molecule AA > 1, the TAMU transport model [17] predicts R molecule AA ∼ R tetraquark AA /2 (albeit, considering only the X(3872) from regeneration processes).
In summary, the first evidence for X(3872) production in heavy ion collisions is presented using lead-lead collisions at a center-of-mass energy of √ s NN = 5.02 TeV per nucleon pair, recorded with the CMS detector. The X(3872) state is reconstructed using the decay chain X(3872) → J/ψ π + π − → µ + µ − π + π − . The measurement is performed for transverse momentum values of the X(3872) of 15 < p T < 50 GeV/c and rapidity |y| < 1.6. The significance of the inclusive X(3872) signal is 4.2 standard deviations. The ratio ρ PbPb between the prompt X(3872) and ψ(2S) yields times their branching fractions into J/ψ π + π − is found to be 1.08 ± 0.49 (stat) ± 0.52 (syst), to be compared with typical values of 0.1 for pp collisions. This result provides a unique experimental input to theoretical models of the X(3872) production mechanism, and of the nature of this exotic state.

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     [9] CDF Collaboration, "Analysis of the quantum numbers J PC of the X(3872)", Phys. Rev.