Observation of an exotic narrow doubly charmed tetraquark

Conventional hadronic matter consists of baryons and mesons made of three quarks and quark-antiquark pairs, respectively. The observation of a new type of hadronic state, a doubly charmed tetraquark containing two charm quarks, an anti-$u$ and an anti-$d$ quark, is reported using data collected by the LHCb experiment at the Large Hadron Collider. This exotic state with a mass of about 3875 MeV$/c^2$ manifests itself as a narrow peak in the mass spectrum of $D^0D^0\pi^+$ mesons just below the $D^{*+}D^0$ mass threshold. The near threshold mass together with a strikingly narrow width reveals the resonance nature of the state.

Quantum chromodynamics (QCD), the theory of the strong force, describes interactions of coloured quarks and gluons and the formation of hadronic matter, the so-called mesons and baryons. While QCD makes precise predictions at high energies, the theory has difficulties describing the interactions of quarks in hadrons from first principles due to the highly non-perturbative regime at the corresponding energy scale. Hence, the field of hadron spectroscopy is driven by experimental discoveries that sometimes are unexpected, which could lead to changes in the research landscape. Along with conventional mesons and baryons, made of a quark-antiquark pair (q 1 q 2 ) and three quarks (q 1 q 2 q 3 ), respectively, particles with an alternative quark content, known as exotic states, have been actively discussed since the birth of the constituent quark model [1][2][3][4][5][6][7][8]. The discussion has been revived by recent observations of numerous tetraquark q 1 q 2 q 3 q 4 and pentaquark q 1 q 2 q 3 q 4 q 5 candidates . Due to the closeness of their masses to known particle-pair thresholds [37,38], many of these states are likely to be hadronic molecules [39][40][41][42] where colour-singlet hadrons are bound by residual nuclear forces similar to the electromagnetic van der Waals forces attracting electrically neutral atoms and molecules. An ordinary example of a hadronic molecule is the deuteron formed by a proton and a neutron. On the other hand, an interpretation of exotic states as compact multiquark structures is also possible [43].
All exotic hadrons observed so far predominantly decay via the strong interaction and their decay widths vary from a few to a few hundred MeV. A discovery of a long-lived exotic state, stable with respect to the strong interaction, would be intriguing. A hadron with two heavy quarks Q and two light antiquarks q, i.e. Q 1 Q 2 q 1 q 2 , is a prime candidate to form such a state [44][45][46][47][48][49]. In the limit of a large heavy-quark mass the two heavy quarks Q 1 Q 2 form a point-like heavy colour-antitriplet object, that behaves similarly to an antiquark, and the corresponding state should be bound. It is expected that the b quark is heavy enough to sustain the existence of a stable bbud state with its binding energy of about 200 MeV with respect to the sum of masses of the pseudoscalar, B − or B 0 , and vector, B * − or B * 0 , beauty mesons that defines the minimal mass for the strong decay to be allowed. In the case of the bcud and ccud systems, there is currently no consensus as to whether such states exist and are narrow enough to be detected experimentally. The similarity of the ccud tetraquark state and the Ξ ++ cc baryon containing two c quarks and a u quark, leads to a relationship between the properties of the two states. In particular, the measured mass of the Ξ ++ cc baryon with quark content ccu [50-52] implies that the mass of the ccud tetraquark is close to the sum of masses of D 0 and D * + mesons with quark content of cu and cd, respectively, as suggested in Ref. [53]. Theoretical predictions for the mass of the ccud ground state with spin-parity quantum numbers J P = 1 + and isospin I = 0, denoted hereafter as T + cc , relative to the D * + D 0 mass threshold lie in the range −300 < δm < 300 MeV/c 2 , where m D * + and m D 0 denote the known masses of the D * + and D 0 mesons [38]. Lattice QCD calculations also do not provide a definite conclusion on the existence of the T + cc state and its binding energy [73,[85][86][87]. The observation of the Ξ ++ cc baryon [50, 51] and of a new exotic resonance decaying to a pair of J/ψ mesons [29] by the LHCb experiment, motivates the search for the T + cc state. In this Letter, the observation of a narrow state in the D 0 D 0 π + mass spectrum near the D * + D 0 mass threshold compatible with being a T + cc tetraquark state is reported. Throughout this Letter, charge conjugate decays are implied. The study is based on proton-proton (pp) collision data collected with the LHCb detector at the Large Hadron Collider (LHC) at CERN at centre-of-mass energies of 7, 8 and 13 TeV, corresponding to integrated luminosity of 9 fb −1 . The LHCb detector [88,89] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks and is further described in Methods. Pseudorapidity η is defined as − log tan θ 2 , where θ is a polar angle of the track relative to the proton beam line.
The D 0 D 0 π + final state is reconstructed by selecting events with two D 0 mesons and a positively charged pion all produced at the same pp interaction point. Both D 0 mesons are reconstructed in the D 0 → K − π + decay channel. The selection criteria are similar to those used in Ref. [90]. To subtract background not originating from two D 0 candidates an extended unbinned maximum-likelihood fit to the two-dimensional distribution of the masses of the two D 0 candidates is performed. The corresponding procedure, together with the selection criteria, is described in detail in Methods. To improve the δm mass resolution and to make the determination insensitive to the precision of the D 0 meson mass, the mass of the D 0 D 0 π + combinations is calculated with the mass of each D 0 meson constrained to the known value [38]. The resulting D 0 D 0 π + mass distribution for selected D 0 D 0 π + combinations is shown in Fig. 1. A narrow peak near the D * + D 0 mass threshold is clearly visible.
An extended unbinned maximum-likelihood fit to the D 0 D 0 π + mass distribution is performed using a model consisting of signal and background components. The signal component is described by the convolution of the detector resolution with a resonant shape, which is modelled by a relativistic P-wave two-body Breit-Wigner function modified by a Blatt-Weisskopf form factor with a meson-radius parameter of 3.5 GeV −1 . The use of a P-wave resonance is motivated by the expected J P = 1 + quantum numbers for the T + cc state. A two-body decay structure T + cc → AB is assumed with m A = 2m D 0 and m B = m π + , where m π + stands for the known mass of the π + meson. Several alternative prescriptions are used for evaluation of systematic uncertainties. Despite its simplicity, the model serves well to quantify the existence of the T + cc state and to measure its properties, such as the position and the width of the resonance. A follow-up study [91] investigates the underlying nature of the T + cc state, expanding on the modelling of the signal shape and determining its physical properties. The detector resolution is modelled by the sum of two Gaussian functions with a common mean, where the additional parameters are taken from simulation (see Methods) with corrections applied [32,92,93]. The root mean square of the resolution function is around 400 keV/c 2 . A study of the D 0 π + mass distribution for D 0 D 0 π + combinations in the region above the D * 0 D + mass threshold and below 3.9 GeV/c 2 shows that approximately 90% of all random D 0 D 0 π + combinations contain a genuine D * + meson. Based on this observation, the background component is parametrised by the product of a two-body phase space function and a positive second-order polynomial. The resulting function is convolved with the detector resolution.
The fit results are shown in Fig. 1, and the parameters of interest, namely the signal yield, N , the mass parameter of the Breit-Wigner function relative to the D * + D 0 mass threshold, δm BW ≡ m BW − (m D * + + m D 0 ), and the width parameter, Γ BW , are listed in Table 1. The statistical significance of the observed T + cc → D 0 D 0 π + signal is estimated using Wilks' theorem to be 22 standard deviations. The fit suggests that the mass parameter of the Breit-Wigner shape is slightly below the D * + D 0 mass threshold. The statistical significance of the hypothesis δm BW < 0 is estimated to be 4.3 standard deviations.
Yield/(200 keV/c 2 ) Figure 1: Distribution of D 0 D 0 π + mass. Distribution of D 0 D 0 π + mass where the contribution of the non-D 0 background has been statistically subtracted. The result of the fit with the two-component function described in the text is overlaid. The D * + D 0 and D * 0 D + thresholds are indicated with the vertical dashed lines. The horizontal bin width is indicated on the vertical axis legend. Inset shows a zoomed signal region with fine binning scheme, Uncertainties on the data points are statistical only and represent one standard deviation, calculated as a sum in quadrature of the assigned weights from the background-subtraction procedure.
To validate the presence of the signal component, several additional cross-checks are performed. The data are categorised according to data-taking periods including the polarity Table 1: Parameters obtained obtained from the fit to the D 0 D 0 π + mass spectrum. Signal yield, N , Breit-Wigner mass relative to D * + D 0 mass threshold, δm BW , and width, Γ BW , are listed. The uncertainties are statistical only.

Parameter
Value of the LHCb dipole magnet and the charge of the T + cc candidates. Instead of statistically subtracting the non-D 0 background, the mass of each D 0 → K − π + candidate is required to be within a narrow region around the known mass of the D 0 meson [38]. The results are found to be consistent among all samples and analysis techniques. Furthermore, dedicated studies are performed to ensure that the observed signal is not caused by kaon or pion misidentification, doubly Cabibbo-suppressed D 0 → K + π − decays and D 0 D 0 oscillations, decays of charm hadrons originating from beauty hadrons, or artefacts due to the track reconstruction creating duplicate tracks.
Systematic uncertainties for the δm BW and Γ BW parameters are summarised in Table 2 and described below. The largest systematic uncertainty is related to the fit model and is studied using pseudoexperiments with alternative parametrisations of the D 0 D 0 π + mass shape. Several variations in the fit model are considered: changes in the signal model due to the imperfect knowledge of the detector resolution, an uncertainty in the correction factor for the resolution taken from control channels, parametrisation of the background component and the additional model parameters of the Breit-Wigner function. The model uncertainty related to the assumption of J P = 1 + quantum numbers of the state is estimated and listed separately. The results are affected by the overall detector momentum scale, which is known to a relative precision of δα = 3 × 10 −4 [94]. The corresponding uncertainty is estimated using simulated samples where the momentum-scale is modified by factors of (1 ± δα). In the reconstruction, the momenta of charged tracks are corrected for energy loss in the detector material, the amount of which is known with a relative uncertainty of 10. The resulting uncertainty is assessed by varying the energy loss correction by ±10%. As the mass of the D 0 D 0 π + combinations is calculated with the mass of each D 0 meson constrained to the known value of the D 0 mass, the δm BW parameter is insensitive to the precision of the D 0 mass. However, the small uncertainty of 2 keV/c 2 for the D * + − D 0 mass difference [38] directly affects the δm BW value. The corresponding systematic uncertainty is added.
In summary, using the full dataset corresponding to an integrated luminosity of 9 fb −1 , collected by the LHCb experiment at centre-of-mass energies of 7, 8 and 13 TeV, a narrow peak is observed in the mass spectrum of D 0 D 0 π + candidates produced promptly in pp collisions. The statistical significance of the peak is overwhelming. Using the Breit-Wigner parametrisation, the location of the peak relative to the D * + D 0 mass threshold, δm BW , and the width, Γ BW , are determined to be where the first uncertainty is statistical, the second systematic and the third is related to the J P quantum numbers assignment. The measured δm BW value corresponds to a mass of approximately 3875 MeV. This is the narrowest exotic state observed to date [37,38]. The minimal quark content for this state is ccud. Two heavy quarks of the same flavour make it manifestly exotic, i.e. beyond the conventional pattern of hadron formation found in mesons and baryons. Moreover, a combination of the near-threshold mass, narrow decay width and its appearance in prompt hadroproduction demonstrates its genuine resonance nature. The measured mass and width are consistent with the expected values for a T + cc isoscalar tetraquark ground state with quantum numbers J P = 1 + . The precision of the mass measurement with respect to the corresponding threshold is superior to those of all other exotic states, which will give better understanding of the nature of exotic states. A dedicated study of the reaction amplitudes for the T + cc → D 0 D 0 π + and T + cc → D 0 D + π 0 (γ) decays that uses the isospin symmetry for the T + cc → D * D transition [91] yields insights on the fundamental resonance properties, like the pole position, the scattering length and the effective range. The observation of this ccud tetraquark candidate close to the D * + D 0 threshold provides strong support for the theory approaches and models that predict the existence of a bbud tetraquark stable with respect to the strong and electromagnetic interactions.

Experimental setup
The LHCb detector [88,89] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary pp collision vertex (PV), the impact parameter (IP), is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors [95]. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers. 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. The trigger selection algorithms are primarily based on identifying key characteristics of beauty and charm hadrons and their decay products, such as high p T final state particles, and a decay vertex that is significantly displaced from any of the pp interaction vertices in the event.

Simulated samples
Simulation is required to model the effects of the detector acceptance, resolution and the imposed selection requirements. In the simulation, pp collisions are generated using Pythia with a specific LHCb configuration [96]. Decays of unstable particles are described by EvtGen [97]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit as described in Ref. [98].

Selection
The selection of D 0 → K − π + candidates and D 0 D 0 π + combinations is similar to those used in Ref. [90]. Kaon and pion candidates are selected from well-reconstructed tracks within the acceptance of the spectrometer. Particle identification is provided using information from the ring-imaging Cherenkov detectors. Kaons and pions that have transverse momenta larger than 250 MeV/c and are inconsistent with being produced in a pp interaction vertex are combined together to form D 0 candidates. The resulting D 0 candidates are required to have good vertex quality, mass within ±65 MeV/c 2 of the known D 0 mass [38] (mass resolution for the D 0 → K − π + signal is 7 MeV/c 2 ), transverse momentum larger than 1 GeV/c, decay time larger than 100 µm/c and a momentum direction that is consistent with the vector from the primary to the secondary vertex. Selected pairs of D 0 candidates consistent with originating from a common primary vertex are then combined with pion candidates of the same charge as the pions from the D 0 → K − π + decay candidates to form D 0 D 0 π + candidates. At least one of the two D 0 π + combinations is required to have good vertex quality and mass not exceeding the known D 0 mass by more than 155 MeV/c 2 . For each D 0 D 0 π + candidate a kinematic fit [99] is performed. This fit requires both D 0 candidates and a pion to originate from the same primary vertex. A requirement on the quality of this fit is applied to further suppress combinatorial background and reduce the background from D 0 candidates produced in two independent pp interactions or in decays of beauty hadrons [90]. To suppress background from kaon and pion candidates reconstructed from one common track, all track pairs of the same charge are required to have opening angle inconsistent with being zero and mass of the combination to be inconsistent with the sum of masses of the two constituents.

Non-D 0 background subtraction
The two-dimensional distribution of the mass of one D 0 candidate versus the mass of the other D 0 candidate from selected D 0 D 0 π + combinations is used to subtract the background from fake D 0 candidates. The procedure employs the sPlot technique [100], where an extended unbinned maximum-likelihood fit to the two-dimensional distribution is performed. The signal is described using a modified Novosibirsk function and the background is modelled by a product of an exponential function and a positive polynomial function [90]. Each candidate is assigned a positive weight for being signal-like or a negative weight for being background-like, with the masses of the two D 0 candidates as the discriminating variables. All candidates are then retained and the weights are used in the further analysis for statistical subtraction of non-D 0 background.

Contributions from the D 0 D 0 oscillations
A hypothetical narrow charmonium-like state decaying into the D 0 D 0 π + final state, followed by the D 0 → D 0 transition or doubly Cabibbo-suppressed decay D 0 → K − π + , would produce a narrow signal in the reconstructed D 0 D 0 π + mass spectrum. If the observed narrow near-threshold peak in the reconstructed D 0 D 0 π + system is caused by the D 0 D 0 oscillations or doubly Cabibbo-suppressed decays, a much larger signal should be visible in the reconstructed D 0 D 0 π + mass spectrum at the same mass. No such structure is observed, see Fig. 9 in Ref. [91].

Systematic uncertainties
Several sources of systematic uncertainty on the mass δm BW and width Γ BW of the T + cc state have been evaluated. The largest systematic uncertainty is related to the fit model and is studied using a set of alternative parametrisations and pseudoexperiments. For each alternative model, an ensemble of pseudoexperiments is produced; each is generated using the model under consideration with parameters obtained from a fit to the data. A subsequent fit with the default model to each pseudoexperiment is performed and the mean values of the parameters of interest over the ensemble are evaluated. The absolute value of the difference between the ensemble mean and the value of the parameter obtained from the fit to the data sample is used to characterise the difference between the alternative model and the default model. The maximal value of such a difference over the considered set of alternative models is taken as the corresponding systematic uncertainty for the mass δm BW and width Γ BW of the T + cc state. The following sources of systematic uncertainties related to the fit model are considered: • Imperfect knowledge of detector resolution model: to estimate the associated systematic uncertainty alternative resolution functions are studied, namely a symmetric variant of an Apollonios function a modified Gaussian function with symmetric power-law tails on both sides of the distribution a generalised symmetric Student's t-distribution; a symmetric Johnson's S U distribution and a modified Novosibirsk function.
• Difference in detector resolution due to imperfect modelling: a correction factor of 1.05 for the resolution is applied for the default fit to account for such a difference. This factor was studied for several other decays measured with the LHCb detector and found to lie between 1.0 and 1.1 [92,93] For decays with relatively low-momentum tracks, this factor is close to 1.05. The factor is also cross-checked using large samples of D * + → D 0 π + decays, where a value of 1.06 is obtained. To assess the systematic uncertainty related to this factor, detector resolution models with correction factors of 1.0 and 1.1 are studied as alternative models.
• Parametrisation of the background component: to assess the associated systematic uncertainty, the order of the positive polynomial function used for the baseline fit is varied. In addition, to estimate the possible effect of a small contribution from D 0 D 0 π + combinations without an intermediate D * + meson, a three-body background component is added to the fit. This component is described by a product of the three-body phase space function and a positive linear or second-order polynomial function. The contribution from non-resonant D 0 D 0 π + background is negligible in the low-mass region due to the O(Q 2 ) scaling of the three-body phase-space factor near threshold.
• Model parameters of the Breit-Wigner function: alternative parametrisations include different choices for the decay structure, m A = m D 0 and m B = m D 0 +m π + ; the meson radius, 1.5 GeV −1 and 5 GeV −1 , and the orbital angular momentum between A and B particles, corresponding to S-and D-waves. The effect of the different decay structure and choice of meson radius is smaller than 1 keV/c 2 and 1 keV for the δm BW and Γ BW parameters, respectively. The parameters of interest are more sensitive to the choice of orbital angular momentum, in which the S-wave function gives larger δm BW and smaller Γ BW , and the D-wave function corresponds to smaller δm BW and larger Γ BW . As the S-wave and D-wave imply that the quantum numbers of the T + cc state differ from J P = 1 + , the corresponding systematic uncertainty is considered separately and is not included it in the total systematic uncertainty.
The calibration of the momentum scale of the tracking system is based upon large calibration samples of B + → J/ψK + and J/ψ → µ + µ − decays. The accuracy of the procedure has been checked using other fully reconstructed B decays together with two-body Υ(nS) and K 0 S decays and the largest deviation of the bias in the momentum scale of δα = 3 × 10 −4 is taken as the uncertainty [94]. This is then propagated to uncertainties for the parameters of interest using simulated samples, where momentum scale corrections of (1 ± δα) are applied. Half of the difference between the peak locations obtained with 1 + δα and 1 − δα corrections applied to the same simulated sample is taken as an estimate of the systematic uncertainty due to the momentum scale. The main contribution to this uncertainty is due to the bachelor pion track, since the D 0 mass constraint reduces the contributions from the kaon and pion tracks originating from D 0 meson decays.
In the reconstruction the momenta of the charged tracks are corrected for the energy loss in the detector material. The energy loss corrections are calculated using the Bethe-Bloch formula. The amount of material traversed in the tracking system by a charged particle is known to 10% accuracy. To assess the corresponding uncertainty, the magnitude of the calculated corrections is varied by ±10%. Half of the difference between the peak locations obtained with +10% and −10% corrections applied to the same simulated sample is taken as an estimate of the systematic uncertainty due to the energy loss corrections.
The mass of the D 0 D 0 π + combinations is calculated with both D 0 candidate masses constrained to the known D 0 meson mass [38]. This procedure removes the uncertainty on the δm BW parameter related to imprecise knowledge of the D 0 mass. In contrast, the small uncertainty of 2 keV/c 2 for the known D * + −D 0 mass difference [38] directly affects the δm BW value and therefore is assigned as the corresponding systematic uncertainty.
in writing software, calibrating sub-systems, operating the detector and acquiring data and finally analysing the processed data.

Competing Interests Statement
The authors declare no competing interests.

Correspondence and requests for materials
Correspondence and requests for materials should be addressed to I. Belyaev Ivan.Belyaev@itep.ru.

Data Availability Statement
LHCb data used in this analysis will be released according to the LHCb external data access policy, that can be downloaded from http://opendata.cern.ch/record/410/files/LHCb-Data-Policy.pdf.
The raw data in all of the figures of this manuscript can be downloaded from https://cds.cern.ch/record/2780001, where no access codes are required. In addition, the unbinned background-subtracted data, shown in Fig. 1 have been added to the HEPData record at https://www.hepdata.net/record/ins1915457.

Code Availability Statement
LHCb software used to process the data analysed in this manuscript is available at GitLab repository https://gitlab.cern.ch/lhcb. The specific software used in data analysis is available at Zenodo repository DOI:10.5281/zenodo.5595937.