Study of the doubly charmed tetraquark [Formula: see text].

Quantum chromodynamics, the theory of the strong force, describes interactions of coloured quarks and gluons and the formation of hadronic matter. Conventional hadronic matter consists of baryons and mesons made of three quarks and quark-antiquark pairs, respectively. Particles with an alternative quark content are known as exotic states. Here a study is reported of an exotic narrow state in the D0D0π+ mass spectrum just below the D*+D0 mass threshold produced in proton-proton collisions collected with the LHCb detector at the Large Hadron Collider. The state is consistent with the ground isoscalar [Formula: see text] tetraquark with a quark content of [Formula: see text] and spin-parity quantum numbers JP = 1+. Study of the DD mass spectra disfavours interpretation of the resonance as the isovector state. The decay structure via intermediate off-shell D*+ mesons is consistent with the observed D0π+ mass distribution. To analyse the mass of the resonance and its coupling to the D*D system, a dedicated model is developed under the assumption of an isoscalar axial-vector [Formula: see text] state decaying to the D*D channel. Using this model, resonance parameters including the pole position, scattering length, effective range and compositeness are determined to reveal important information about the nature of the [Formula: see text] state. In addition, an unexpected dependence of the production rate on track multiplicity is observed.

Introduction

Hadrons with quark content other than that seen in mesons () and baryons (q1q2q3) have been actively discussed since the birth of the quark model[1-8]. Since the discovery of the χc1(3872) state[9] many tetraquark and pentaquark candidates, listed in Table 1, have been observed[10-19]. For all but the X0(2900) and X1(2900) states the minimal quark content implies the presence of either a or quark-antiquark pair. The masses of many tetra- and pentaquark states are close to mass thresholds, e.g., or , where D(*) or B(*) represents a hadron containing a charm or beauty quark, respectively. Therefore, these states are likely to be hadronic molecules[16,20-22] where colour-singlet hadrons are bound by residual nuclear forces, such as the exchange of a pion or ρ meson[23], similar to electromagnetic van der Waals forces attracting electrically neutral atoms and molecules. These states are expected to have a spatial extension significantly larger than a typical compact hadron. Conversely, the only hadron currently observed that contains a pair of c quarks is the (ccu) baryon, a long-lived, weakly-decaying compact object[24,25]. The recently observed X(6900) structure in the J/ψJ/ψ mass spectrum[26] belongs to both categories simultaneously. Its proximity to the χc0χc1 threshold could indicate a molecular structure[27,28]. Alternatively, it could be a compact object, where all four quarks are within one confinement volume and each quark interacts directly with the other three quarks via the strong force[29-32].

Table 1

Tetra- and pentaquark candidates and their plausible valence quark content.

States	Quark content
X0(2900), X1(2900)[148,149]	c¯dus¯
χc1(3872)[9]	cc¯qq¯
Zc(3900)[150–154], Zc(4020)[155,156], Zc(4050)[157], X(4100)[158], Zc(4200)[159], Zc(4430)[160–163], Rc0(4240)[162]	cc¯ud¯
Zcs(3985)[164], Zcs(4000), Zcs(4220)[165]	cc¯us¯
χc1(4140)[166–169], χc1(4274), χc0(4500), χc0(4700)[169], X(4630), X(4685)[165], X(4740)[96]	cc¯ss¯
X(6900)[26]	cc¯cc¯
Zb(10610), Zb(10650)[170]	bb¯ud¯
Pc(4312)[171], Pc(4380)[172], Pc(4440), Pc(4457)[171], Pc(4357)[173]	cc¯uud
Pcs(4459)[174]	cc¯uds

The existence and properties of states with two heavy quarks and two light antiquarks have been widely discussed for a long time[33-38]. In the limit of large masses of the heavy quarks the corresponding ground state should be deeply bound. In this limit, the two heavy quarks Q1Q2 form a point-like colour-antitriplet object, analogous to an antiquark, and as a result, the system has similar degrees of freedom for its light quarks as an antibaryon with a single heavy quark, e.g., the or antibaryons. The beauty quark is considered heavy enough to sustain the existence of a state that is stable with respect to the strong and electromagnetic interactions with a mass of about 200 MeV/c2 below the BB* mass threshold. In the case of the and systems, there is currently no consensus in the literature whether such states exist and if their natural widths are narrow enough to allow for experimental observation. The theoretical predictions for the mass of the ground state with spin-parity quantum numbers JP = 1+ and isospin I = 0, denoted hereafter as , relative to the D*+D0 mass thresholdlie in the range −300 < δm < 300 MeV/c2 [39-70], where and denote the known masses of the D*+ and D0 mesons[10], with and quark content, respectively. The observation of a narrow state in the D0D0π+ mass spectrum near the D*+D0 mass threshold, compatible with being a tetraquark state with quark content is reported in Ref. [71].

In the work presented here, the properties of the state are studied by constructing a dedicated amplitude model that accounts for the D*+D0 and D*0D+ decay channels. In addition, the mass spectra of other DD(*) and opposite-sign combinations are explored. Furthermore, production-related observables, such as the event multiplicity and transverse momentum (pT) spectra that are sensitive to the internal structure of the state, are discussed. This analysis is based on proton–proton (pp) collision data, corresponding to an integrated luminosity of 9 fb−1, collected with the LHCb detector at centre-of-mass energies of 7, 8 and 13 TeV. The LHCb detector[72,73] 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.

Results

signal in the D0D0π+ mass spectrum

The D0D0π+ final state is reconstructed using the D0 → K−π+ decay channel with two D0 mesons and a pion all produced promptly in the same pp collision. The inclusion of charge-conjugated processes is implied throughout the paper. The selection criteria are similar to those used in Refs. [74-77] and described in detail in Methods. The background not originating from true D0 mesons is subtracted using an extended unbinned maximum-likelihood fit to the two-dimensional distribution of the masses of the two D0 candidates from selected D0D0π+ combinations, see Methods and Supplementary Fig. 1a. The obtained D0D0π+ mass distribution for selected D0D0π+ combinations is shown in Fig. 1.

Fig. 1

Distribution of D0D0π+ mass.

Distribution of D0D0π+ mass where the contribution of the non-D0 background has been statistically subtracted. The result of the fit described in the text is overlaid. 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.

An extended unbinned maximum-likelihood fit to the D0D0π+ mass distribution is performed using a model consisting of signal and background components. The signal component corresponds to the decay and is described as the convolution of the natural resonance profile with the detector mass resolution function. A relativistic P-wave two-body Breit–Wigner function with a Blatt–Weisskopf form factor[78,79] is used in Ref. [71] as the natural resonance profile. That function, while sufficient to reveal the existence of the state, does not account for the resonance being in close vicinity of the D*D threshold. To assess the fundamental properties of resonances that are close to thresholds, advanced parametrisations ought to be used[80-90]. A unitarised Breit–Wigner profile , described in Methods Eq. (47), is used in this analysis. The function is built under two main assumptions.

Assumption 1

The newly observed state has quantum numbers JP = 1+ and isospin I = 0 in accordance with the theoretical expectation for the ground state.

Assumption 2

The state is strongly coupled to the D*D channel, which is motivated by the proximity of the mass to the D*D mass threshold.

The derivation of the profile relies on the assumed isospin symmetry for the decays and the coupled-channel interaction of the D*+D0 and D*0D0 system as required by unitarity and causality following Ref. [91]. The resulting energy-dependent width of the state accounts explicitly for the , and decays. The modification of the D* meson lineshape[92] due to contributions from triangle diagrams[93] to the final-state interactions is neglected. Similarly to the profile, the function has two parameters: the peak location mU, defined as the mass value where the real part of the complex amplitude vanishes, and the absolute value of the coupling constant g for the decay.

The detector mass resolution, , is modelled with the sum of two Gaussian functions with a common mean, and parameters taken from simulation, see Methods. The widths of the Gaussian functions are corrected by a factor of 1.05, which accounts for a small residual difference between simulation and data[94-96]. The root mean square (RMS) of the resolution function is around 400 keV/c2.

A study of the D0π+ mass distribution for selected D0D0π+ combinations in the region above the D*0D+ mass threshold and below 3.9 GeV/c2 shows that approximately 90% of all D0D0π+ combinations contain a true D*+ meson. Therefore, the background component is parameterised with a product of the two-body phase-space function [97] and a positive polynomial function P, convolved with the detector resolution function where n denotes the order of the polynomial function, n = 2 is used in the default fit.

The D0D0π+ mass spectrum with non-D0 background subtracted is shown in Fig. 1 with the result of the fit using a model based on the signal profile overlaid. The fit gives a signal yield of 186 ± 24 and a mass parameter relative to the D*+D0 mass threshold, δmU of −359 ± 40 keV/c2. The statistical significances of the observed signal and for the δmU < 0 hypothesis are determined using Wilks’ theorem to be 22 and 9 standard deviations, respectively.

The width of the resonance is determined by the coupling constant g for small values of . With increasing , the width increases to an asymptotic value determined by the width of the D*+ meson, see Methods and Supplementary Fig. 7. In this regime of large , the signal profile exhibits a scaling property similar to the Flatté function[94,98,99]. The parameter effectively decouples from the fit model, and the model resembles the scattering-length approximation[81]. The likelihood profile for the parameter is shown in Fig. 2, where one can see a plateau at large values. At small values of the parameter, , the likelihood function is independent of because the resonance is too narrow for the details of the signal profile to be resolved by the detector. The lower limits on the parameter of at 90% (95%) confidence level (CL) are obtained as the values where the difference in the negative log-likelihood is equal to 1.35 and 1.92, respectively. Smaller values for are further used for systematic uncertainty evaluation.

Fig. 2

Likelihood profile for the parameter.

Likelihood profile for the absolute value of the coupling constant g from the fit to the background-subtracted D0D0π+ mass spectrum with a model based on the signal profile.

The mode relative to the D*+D0 mass threshold, , and the full width at half maximum (FWHM), , for the profile are found to be and , to be compared with those quantities determined for the signal profile of and . They appear to be rather different. Nonetheless, both functions properly describe the data given the limited sample size, and accounting for the detector resolution, and residual background. To quantify the impact of these experimental effects, two ensembles of pseudoexperiments are performed. Firstly, pseudodata samples are generated with a model based on the profile. The parameters used here are obtained from the default fit, and the size of the sample corresponds to the size of data sample. Each pseudodata sample is then analysed with a model based on the function. The obtained mean and RMS values for the parameters δmBW and ΓBW over the ensemble are shown in Table 2. The mass parameter δmBW agrees well with the value determined from data[71]. The difference for the parameter ΓBW does not exceed one standard deviation. Secondly, an ensemble of pseudodata samples generated with a model based on the profile is analysed with a model based on the function. The obtained mean and RMS values for the δmU parameter over an ensemble are also reported in Table 2. These values agree well with the result of the default fit to data. The results of these pseudoexperiments explain the seeming inconsistency between the models and illustrate the importance of an accurate description of the detector resolution and residual background given the limited sample size.

Table 2

Mean and root mean square (RMS) values for the δmBW, ΓBW and δmU parameters obtained from pseudoexperiments produced as a consistency check.

		Pseudoexperiments
Parameter		Mean	RMS	Data
δmBW	keV/c^2	−301	50	−273 ± 61[71]
ΓBW	keV	222	121	410 ± 165[71]
δmU	keV/c^2	−378	46	−359 ± 40

Systematic uncertainties

Systematic uncertainties for the δmU parameter are summarised in Table 3 and described in greater detail below. The systematic uncertainty related to the fit model is studied using pseudoexperiments with a set of alternative parameterisations. For each alternative model, an ensemble of pseudoexperiments is performed with parameters obtained from a fit to data. A fit with the baseline model is performed on each pseudoexperiment, and the mean values of the parameters of interest are evaluated over the ensemble. The absolute values of the differences between these mean values and the corresponding parameter values obtained from the fit to data are used to assess the systematic uncertainty due to the choice of the fit model. The maximal value of such differences over the considered set of alternative models is taken as the corresponding systematic uncertainty. The following sources of systematic uncertainty related to the fit model are considered:

Table 3

Systematic uncertainties for the δmU parameter.

Source	σ_δm_UkeV/c^2
Fit model
Resolution model	2
Resolution correction factor	2
Background model	2
Coupling constants	1
Unknown value of g	−0+7
Momentum scaling	3
Energy loss	1
D*+ − D0 mass difference	2
Total	−6+9

The total uncertainty is calculated as the sum in quadrature of all components.

Imperfect knowledge of the detector resolution model. To estimate the associated systematic uncertainty a set of alternative resolution functions is tested: a symmetric variant of an Apollonios function[100], a modified Gaussian function with symmetric power-law tails on both sides of the distribution[101,102], a generalised symmetric Student’s t-distribution[103,104], a symmetric Johnson’s SU distribution[105,106], and a modified Novosibirsk function[107].

A small difference in the detector resolution between data and simulation. A correction factor of 1.05 is applied to account for known discrepancies in modelling the detector resolution in simulation. This factor was studied for different decays[94-96,108-110] and found to lie between 1.0 and 1.1. For decays with relatively low momentum tracks, this factor is close to 1.05, which is the nominal value used in this analysis. This factor is also cross-checked using large samples of D*+ → D0π+ 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 alternatives.

Parameterisation of the background component. To assess the associated systematic uncertainty, the order of the positive polynomial function of Eq. (2) is varied. In addition, to estimate a possible effect from a small contribution from three-body D0D0π+ combinations without an intermediate D*+ meson, a more general family of background models is testedwhere denotes the three-body phase-space function[111]. The functions B0, B1, B3 and with n ≤ 2, m ≤ 1 are used as alternative models for the estimation of the systematic uncertainty.

Values of the coupling constants for the D* → Dπ and D* → Dγ decays affecting the shape of the signal profile. These coupling constants are calculated from the known branching fractions of the D* → Dπ and D* → Dγ decays[10], the measured natural width of the D*+ meson[10,112] and the derived value for the natural width of the D*0 meson[66,81,113]. To assess the associated systematic uncertainty, a set of alternative models built around the profiles, obtained with coupling constants varying within their calculated uncertainties, is studied.

Unknown value of the parameter. In the baseline fit the value of the parameter is fixed to a large value. To assess the effect of this constraint the fit is repeated using the value of , that corresponds to for the most conservative likelihood profile for that accounts for the systematic uncertainty. The change of 7 keV/c2 of the δmU parameter is assigned as the systematic uncertainty.

The calibration of the momentum scale of the tracking system is based upon large samples of B+ → J/ψK+ and J/ψ → μ+μ− decays[114]. The accuracy of the procedure has been checked using fully reconstructed B decays together with two-body ϒ(nS) and decays and the largest deviation of the bias in the momentum scale of δα = 3 × 10−4 is taken as the uncertainty[115]. This uncertainty is propagated for the parameters of interest using simulated samples, with momentum scale corrections of applied. Half of the difference between the obtained peak locations is taken as an estimate of the systematic uncertainty.

In the reconstruction the momenta of the charged tracks are corrected for energy loss in the detector material using the Bethe–Bloch formula[116,117]. The amount of the material traversed in the tracking system by a charged particle is known to have 10% accuracy[118]. To assess the corresponding uncertainty the magnitude of the calculated corrections is varied by ±10%. Half of the difference between the obtained peak locations is taken as an estimate of the systematic uncertainty due to energy loss corrections.

The mass of D0D0π+ combinations is calculated with the mass of each D0 meson constrained to the known value of the D0 mass[10]. This procedure produces negligible uncertainties for the δmU parameter due to imprecise knowledge of the D0 mass. However, the small uncertainty of 2 keV/c2 for the known D*+ − D0 mass difference[10,112,119] directly affects the values of these parameters and is assigned as a corresponding systematic uncertainty.

For the lower limit on the parameter , only systematic uncertainties related to the fit model are considered. For each alternative model the likelihood profile curves are built and corresponding 90 and 95% CL lower limits are calculated using the procedure described above. The smallest of the resulting values is taken as the lower limit that accounts for the systematic uncertainty: at 90 (95%) CL.

Results

Studying the D0π+ mass distribution for decays allows testing the hypothesis that the decay proceeds through an intermediate off-shell D*+ meson. The background-subtracted D0π+ mass distribution for selected D0D0π+ candidates with the D0D0π+ mass with respect to the D*+D0 mass threshold, , below zero is shown in Fig. 3. Both D0π+ combinations are included in this plot. The two-dimensional distribution of the mass of one D0π+ combination versus the mass of another D0π+ combination is presented in Supplementary Fig. 10.

Fig. 3

Mass distribution for D0π+ pairs.

Mass distribution for D0π+ pairs from selected D0D0π+ candidates with a mass below the D*+D0 mass threshold with non-D0 background subtracted. The overlaid fit result is described in the text. The background component vanishes in the fit. 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.

A fit is performed to this distribution with a model containing signal and background components. The signal component is derived from the amplitude, see Methods Eq. (48), and is convolved with a detector resolution for the D0π+ mass. This detector resolution function is modelled with a modified Gaussian function with power-law tails on both sides of the distribution[101,102] and parameters taken from simulation. Similarly to the correction used for the D0D0π+ mass resolution function , the width of the Gaussian function is corrected by a factor of 1.06 which is determined by studying large samples of D*+ → D0π+ decays. The RMS of the resolution function is around 220 keV/c2. The shape of the background component is derived from data for . The fit results are overlaid in Fig. 3. The background component vanishes in the fit, and the D0π+ spectrum is consistent with the hypothesis that the decay proceeds through an intermediate off-shell D*+ meson. This in turn favours the 1+ assignment for the spin-parity of the state.

Due to the proximity of the observed signal to the D*+D0 mass threshold, and the small energy release in the D*+ → D0π+ decay, the D0D0 mass distribution from the decay forms a narrow peak just above the D0D0 mass threshold. In a similar way, a peaking structure in the D+D0 mass spectrum just above the D+D0 mass threshold is expected from and decays, both proceeding via off-shell intermediate D*+D0 and D*0D+ states. The D0D0 and D+D0 final states are reconstructed and selected similarly to the D0D0π+ final state, where the D+ → K−π+π− decay channel is used. The background-subtracted D0D0 and D+D0 mass distributions are shown in Fig. 4 (top), where narrow structures are clearly visible just above the DD thresholds. Fits to these distributions are performed using models consisting of two components: a signal component FDD described in Methods Eqs. (70) and (71) and obtained via integration of the matrix elements for the decays with the profile, and a background component, parameterised as a product of the two-body phase-space function ΦDD and a positive linear function P1. The fit results are overlaid in Fig. 4 (top). The signal yields in the D0D0 and D+D0 spectra are found to be 263 ± 23 and 171 ± 26, respectively. The statistical significance of the observed and signals, where X stands for non-reconstructed pions or photons, is estimated using Wilks’ theorem[120] and is found to be in excess of 20 and 10 standard deviations, respectively. The relative yields for the signals observed in the D0D0π+, D0D0 and D0D+ mass spectra agree with the expectations of the model described in Methods where the decay of an isoscalar state via the D*D channel with an intermediate off-shell D* meson is assumed.

Fig. 4

Mass distributions for selected D0D0, D+D0, D+D+ and D+D0π+ combinations.

(Top) DD and DDπ+ mass distributions for selected (left) D0D0 and (right) D+D0 candidates with the non-D background subtracted. The overlaid fit results are described in the text. For visibility the is stacked on top of the component. (Bottom) Mass distributions for selected (left) D+D+ and (right) D+D0π+ candidates with the non-D background subtracted. The vertical coloured band indicates the expected mass for the hypothetical state. The overlaid fit results with background-only functions are described in the text. 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.

The observation of the near-threshold signals in the D0D0 and D+D0 mass spectra, along with the signal shapes and yields, all agree with the isoscalar hypothesis for the narrow signal observed in the D0D0π+ mass spectrum. However, an alternative interpretation could be that this state is the I3 = 0 component of a isotriplet (,,) with , and quark content, respectively. Assuming that the observed peak corresponds to the component and using the estimates for the mass splitting from Methods Eqs. (85) and (86), the masses of the and states are estimated to be slightly below the D0D*0 and slightly above the D+D*+ mass thresholds, respectively:

With these mass assignments, assuming equal production of all three components, the state would be an extra narrow state that decays into the D0D0π0 and D0D0γ final states via an off-shell D*0 meson. These decays would contribute to the narrow near-threshold enhancement in the D0D0 spectrum, and increase the signal in the D0D0 mass spectrum by almost a factor of three. The state would decay via an on-shell D*+ meson ; therefore, it could be a relatively wide state, with a width up to a few Me[121]. Therefore, it would manifest itself as a peak with a moderate width in the D+D0π+ mass spectrum with a yield comparable to that of the decays. In addition, it would contribute to the D+D0 mass spectrum, tripling the contribution from the decays. However, due to the larger mass of the state and its larger width, this contribution should be wider, making it more difficult to disentangle from the background. Finally, the state would make a contribution to the D+D+ spectrum with a yield similar to the contribution from decays to the D0D+ spectrum, but wider. The mass spectra for D+D+ and D+D0π+ combinations are shown in Fig. 4 (bottom).

Neither distribution exhibits any narrow signal-like structure. Fits to these spectra are performed using the following background-only functions:

The results of these fits are overlaid in Fig. 4 (bottom). The absence of any signals in the D+D+ and D+D0π+ mass spectra is therefore a strong argument in favour of the isoscalar nature of the observed peak in the D0D0π+ mass spectrum.

The interference between two virtual channels for the decay, corresponding to two amplitude terms, see Methods Eq. (35), is studied by setting the term proportional to C in Methods Eq. (39) to be equal to zero. This causes a 43% reduction in the decay rate, pointing to a large interference. The same procedure applied to the decays gives the contribution of 45% for the interference between the and channels. For decays the role of the interference between the and channels is estimated by equating to zero the and terms in Methods Eqs. (45) and (46). The interference contribution is found to be 33%.

Using the model described earlier and results of the fit to the D0D0π+ mass spectrum, the position of the amplitude pole in the complex plane, responsible for the appearance of the narrow structure in the D0D0π+ mass spectrum is determined. The pole parameters, mass mpole and width Γpole, are defined through the pole location as

The pole location is a solution to the equationwhere denotes the amplitude on the second Riemann sheet defined in Methods Eq. (64). For large coupling the position of the resonance pole is uniquely determined by the parameter δmU, i.e., the binding energy and the width of the D*+ meson. Figure 5 shows the complex plane of the variable, defined as

Fig. 5

Complex plane of the variable.

Complex plane of the variable. The dashed red line shows the allowed region for large values. The filled red circle indicates the best estimate for the pole location and the filled regions show 1σ and 2σ confidence regions. Open blue and green circles show the branch points corresponding to the D*+D0 and D*0D+ channels, respectively, and the corresponding blue and green lines indicate branch cuts. Three other branch points at of , and , corresponding to the openings of the D0D+γ, D0D+π0 and D0D0π+ decay channels, are outside of the displayed region.

All possible positions of the pole for are located on a red dashed curve in Fig. 5. The behaviour of the curve can be understood as follows: with an increase of the binding energy (distance to the D*+D0 mass threshold), the width gets narrower; and when the parameter δmU approaches zero, the pole touches the D0D*+ cut and moves to the other complex sheet, i.e., the state becomes virtual. For smaller values of , the pole is located between the limiting curve and the ℜ s = 0 line. The pole parameters are found to bewhere the first uncertainty is due to the δmU parameter and the second is due to the unknown value of the parameter. The peak is well separated from the D*+D0 threshold in the D0D0π+ mass spectrum. Hence, as for an isolated narrow resonance, the parameters of the pole are similar to the visible peak parameters, namely the mode and FWHM .

The systematic uncertainties quoted here do not account for the possibility that any of the underlying assumptions on which the model is built are not valid. For example, as shown earlier the data are consistent with a wide range of FHWM values for the signal profile. Therefore the pole width Γpole is based mainly on the amplitude model and the value of the mU parameter determined from the fit to the D0D0π+ mass spectrum.

A study of the behaviour of the amplitude in the vicinity of the D*+D0 mass threshold leads to the determination of the low-energy scattering parameters, namely the scattering length, a, and the effective range, r. These parameters are defined via the coefficients of the first two terms of the Taylor expansion of the inverse non-relativistic amplitude[122], i.e.,where k is the wave number. For the inverse amplitude from Eq. (47) matches Eq. (13) up to a scale parameter obtained numerically, see Methods Eq. (74). The value of the scattering length is found to be

Typically, a non-vanishing imaginary part of the scattering length indicates the presence of inelastic channels[123]; however, in this case the non-zero imaginary part is related to the lower threshold, , and is determined by the width of the D*+ meson. The real part of the scattering length a is negative indicating attraction. This can be interpreted as the characteristic size of the state[16],

For the amplitude the effective range r is non-positive and proportional to , see Methods Eq. (76). Its value is consistent with zero for the baseline fit. An upper limit on the −r value is set as

The Weinberg compositeness criterion[124,125] makes use of the relation between the scattering length and the effective range to construct the compositeness variable Z,for which Z = 1 corresponds to a compact state that does not interact with the continuum, while Z = 0 indicates a composite state formed by compound interaction. Using the relation between r and ∣g∣ from Methods Eq. (76), one finds for large values of . The default fit corresponds to large values of , and thus, Z approaching to zero. A non-zero value of Z would require a smaller value of , i.e., smaller resonance width, see Supplementary Fig. 7. The following upper limit of the compositeness parameter Z is set:

Another estimate of the characteristic size is obtained from the value of the binding energy ΔE. Within the interpretation of the state as a bound D*+D0 molecular-like state, the binding energy is ΔE = − δmU. The characteristic momentum scale γ[16] is estimated to bewhere μ is the reduced mass of the D*+D0 system. This value of the momentum scale in turn corresponds to a characteristic size RΔ of the molecular-like state,which is consistent with the R estimate from the scattering length.

For high-energy hadroproduction of a state with such a large size, R or RΔ, one expects a strong dependency of the production rate on event multiplicity, similar to that observed for the χc1(3872) state[126]. The background-subtracted distribution of the number of tracks reconstructed in the vertex detector, Ntracks, is shown in Fig. 6 (left) together with the distributions for low-mass pairs with and low-mass D0D0 pairs with mass . The former is dominated by production, while the latter is presumably dominated by the double-parton scattering process[74,127]. The chosen interval for pairs includes the region populated by the decays; however, this contribution is small, see Fig. 7. The χc1(3872) production cross-section is suppressed with respect to the conventional charmonium state ψ(2S) at large track multiplicities[126]. It is noteworthy that the track multiplicity distribution for the state differs from that of the low-mass pairs, in particular, no suppression at large multiplicity is observed. A p value for the consistency of the track multiplicity distributions for production and low-mass pairs is found to be 0.1%. It is interesting to note that the multiplicity distribution for production and the one for D0D0-pairs with are consistent with a corresponding p value of 12%. The similarity between production, which is inherently a single parton scattering process, and the distribution for process dominated by a double-parton scattering is surprising.

Fig. 6

Track multiplicity and transverse momentum distributions.

(Left) Background-subtracted distributions for the multiplicity of tracks reconstructed in the vertex detector for (red circles) signal, low-mass (blue open squares) and (green filled diamonds) D0D0 pairs. The binning scheme is chosen to have an approximately uniform distribution for pairs. The distributions for the and D0D0 pairs are normalised to the same yields as the signal. (right) Background-subtracted transverse momentum spectra for (red circles) signal, (blue open squares) low-mass and (green filled diamonds) D0D0 pairs. The binning scheme is chosen to have an approximately uniform distribution for pairs. The distributions for the and D0D0 pairs are normalised to the same yields as signal. For better visualisation, the points are slightly displaced from the bin centres. For better visualisation, the points are slightly displaced from the bin centres. Uncertainties on the data points are statistical only and represent one standard deviation.

Fig. 7

Mass distributions for D0D0 and candidates.

Background-subtracted D0D0 and mass distributions. The near-threshold enhancement in the D0D0 channel corresponds to partially reconstructed decays, while in the channel the threshold enhancement corresponds to partially reconstructed decays. The mass distribution is zero-suppressed for better visualisation. 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.

The transverse momentum spectrum for the state is compared with those for the low-mass and D0D0 pairs in Fig. 6 (right). The p values for the consistency of the pT spectra for the state and low-mass pairs are 1.4%, and 0.02% for low-mass D0D0 pairs. More data are needed for further conclusions.

The background-subtracted D0D0 mass distribution in a wider mass range is shown in Fig. 7 together with a similar distribution for pairs. In the D0 mass spectrum the near-threshold enhancement is due to and decays via intermediate D*0 mesons[77]. This structure is significantly wider than the structure in the D0D0 mass spectrum from decays primarily due to the larger natural width and smaller binding energy for the χc1(3872) state[94,95]. With more data, and with a better understanding of the dynamics of decays, and therefore of the corresponding shape in the D0 mass spectrum, it will be possible to estimate the relative production rates for the and χc1(3872) states. Background-subtracted D0D0π+ and D0D+ mass distributions together with those for and D0D− are shown in Supplementary Figs. 3 and 4.

Discussion

The exotic narrow tetraquark state observed in Ref. [71] is studied using a dataset corresponding to an integrated luminosity of 9 fb−1, collected by the LHCb experiment in pp collisions at centre-of-mass energies of 7, 8 and 13 TeV. The observed D0π+ mass distribution indicates that the decay proceeds via an intermediate off-shell D*+ meson. Together with the proximity of the state to the D*D0 mass threshold, this favours the spin-parity quantum numbers JP to be 1+. Narrow near-threshold structures are observed in the D0D0 and D0D+ mass spectra with high significance. These are found to be consistent with originating from off-shell decays followed by the D* → Dπ and D* → Dγ decays. No signal is observed in the D+D0π+ mass spectrum, and no structure is observed in the D+D+ mass spectrum. These non-observations provide a strong argument in favour of the isoscalar nature for the observed state, supporting its interpretation as the isoscalar JP = 1+ -tetraquark ground state. A dedicated unitarised three-body Breit–Wigner amplitude is built on the assumption of strong isoscalar coupling of the axial-vector state to the D*D channel. This assumption is supported by the data; however, alternative models are not excluded by the distributions studied in this analysis. Probing alternative models and the validity of the underlying assumptions of this analysis will be a subject for future studies.

Using the developed amplitude model, the mass of the state, relative to the D*+D0 mass threshold, is determined to bewhere the first uncertainty is statistic and the second systematic. The lower limit on the absolute value of the coupling constant of the state to the D*D system isUsing the same model, the estimates for the scattering length a, effective range r, and the compositeness, Z are obtained from the low-energy limit of the amplitude to beThe characteristic size calculated from the binding energy is RΔ = 7.49 ± 0.42 fm. This value is consistent with the estimation from the scattering length, R = 7.16 ± 0.51 fm. Both RΔ and R correspond to a spatial extension significantly exceeding the typical scale for heavy-flavour hadrons. Within this model the resonance pole is found to be located on the second Riemann sheet with respect to the D0D0π+ threshold, at , wherewhere the first uncertainty accounts for statistical and systematic uncertainties for the δmU parameters, and the second is due to the unknown value of the parameter. The pole position, scattering length, effective range and compositeness form a complete set of observables related to the reaction amplitude, which are crucial for inferring the nature of the tetraquark.

Unlike in the prompt production of the χc1(3872) state, no suppression of the production at high track multiplicities is observed relative to the low-mass pairs. The observed similarity with the multiplicity distribution for the low-mass D0D0 production process, that is presumably double-parton-scattering dominated, is unexpected. In the future with a larger dataset and including other decay modes, e.g., D0 → K−π+π+π−, detailed studies of the properties of this new state and its production mechanisms could be possible.

In conclusion, the tetraquark observed in D0D0π+ decays is studied in detail, using a unitarised model that accounts for the relevant thresholds by taking into account the D0D0π+ and D0D+π0(γ) decay channels with intermediate D* resonances. This model is found to give an excellent description of the D0π+ mass distribution in the decay and of the threshold enhancements observed in the D0D0 and D0D+ spectra. Together with the absence of a signal in the D0D+ and D+D0π+ mass distributions this provides a strong argument for interpreting the observed state as the isoscalar tetraquark with spin-parity JP = 1+. The precise mass measurement will rule out or improve on a considerable range of theoretical models on heavy quark systems. The determined pole position and physical quantities derived from low-energy scattering parameters reveal important information about the nature of the tetraquark. In addition, the counter-intuitive dependence of the production rate on track multiplicity will pose a challenge for theoretical explanations.

Methods

Experimental setup

The LHCb detector[72,73] 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, the impact parameter, is measured with a resolution of (15 + 29/pT) μm, where pT 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[128]. 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.

Simulation

Simulation is required to model the effects of the detector acceptance, resolution, and the efficiency of the imposed selection requirements. In the simulation, pp collisions are generated using PYTHIA[129] with a specific LHCb configuration[130]. Decays of unstable particles are described by EVTGEN[131], in which final-state radiation is generated using PHOTOS[132]. The interaction of the generated particles with the detector, and its response, are implemented using the GEANT4 toolkit[133,134] as described in Ref. [135].

Event selection

The D0D0, D0D+ and D0D0π+ final states are reconstructed using the D0 → K−π+ and D+ → K−π+π+ decay channels. The selection criteria are similar to those used in Refs. [74-77]. Kaons and pions are selected from well-reconstructed tracks within the acceptance of the spectrometer that are identified using information from the ring-imaging Cherenkov detectors. The kaon and pion candidates that have transverse momenta larger than 250 MeV/c and are inconsistent with being produced at a pp interaction vertex are combined together to form D0 and D+ candidates, referred to as D hereafter. The resulting D candidates are required to have good vertex quality, mass within ±65 and ±50 MeV/c2 of the known D0 and D+ masses[10], respectively, 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 secondary vertex. Selected D0 and D+ candidates consistent with originating from a common primary vertex are combined to form D0D0 and D0D+ candidates. The resulting D0D0 candidates are combined with a pion to form D0D0π+ candidates. At least one of the two D0π+ combinations is required to have good vertex quality and mass not exceeding the known D*+ mass by more than 155 MeV/c2. For each D0D0, D0D+ and D0D0π+ candidate a kinematic fit[136] is performed. This fit constrains the mass of the D candidates to their known values and requires both D mesons, and a pion in the case of D0D0π+, to originate from the same primary vertex. A requirement is applied to the quality of this fit to further suppress combinatorial background and reduce background from D candidates produced in two independent pp interactions or in the decays of beauty hadrons[74]. To suppress background from kaon and pion candidates reconstructed from a common track, all track pairs of the same charge are required to have an opening angle inconsistent with zero and the mass of the combination must be inconsistent with the sum of the masses of the two constituents. For cross-checks additional final states D+D+, D+D0π+, , D0D− and are reconstructed, selected and treated in the same way.

Non-D background subtraction

Two-dimensional distributions of the mass of one D candidate versus the mass of the other D candidate from selected D0D0π+, D0D0 and D0D+ combinations are shown in Supplementary Fig. 1. These distributions illustrate the relatively small combinatorial background levels due to fake D candidates. This background is subtracted using the sPlot technique[137], which is based on an extended unbinned maximum-likelihood fit to these two-dimensional distributions with the function described in Ref. [74].

This function consists of four components:

a component corresponding to genuine D1D2 pairs and described as a product of two signal functions, each parameterised with a modified Novosibirsk function[107];

two components corresponding to combinations of one of the D mesons with combinatorial background, described as a product of the signal function and a background function, which is parameterised with a product of an exponential function and a positive first-order polynomial;

a component corresponding to pure background pairs and described by a product of exponential functions and a positive two-dimensional non-factorisable second-order polynomial function.

Based on the results of the fit, 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 D0 candidates as discriminating variables. The D0D0π+ mass distributions for each of the subtracted background components are presented in Supplementary Fig. 2. where fit results with background-only functions , defined in Eq. (3) are overlaid.

Resolution model for the D0D0π+ mass

In the vicinity of the D*+D0 mass threshold the resolution function for the D0D0π+ mass is parametrised with the sum of two Gaussian functions with a common mean. The widths of the Gaussian functions are σ1 = 1.05 × 263 keV/c2 and σ2 = 2.413 × σ1 for the narrow and wide components, respectively, and the fraction of the narrow Gaussian is α = 0.778. The parameters α and σ1,2 are taken from simulation, and σ1,2 are corrected with a factor of 1.05 that accounts for a small difference between simulation and data for the mass resolution[94-96]. The RMS of the resolution function is around 400 keV/c2.

Matrix elements for decays

Assuming isospin symmetry, the isoscalar vector state that decays into the D*D final state can be expressed as

Therefore, the S-wave amplitudes for the and decays have different signswhere g is a coupling constant, is the polarisation vector of the particle and is the polarisation vector of D* meson, and the upper and lower Greek indices imply the summation in the Einstein notation. The S-wave (corresponding to orbital angular momentum equal to zero) approximation is valid for a near-threshold peak. For masses significantly above the D*D threshold, higher-order waves also need to be considered. The amplitudes for the D* → Dπ decays are written aswhere f denotes a coupling constant, and pD stands for the momentum of the D meson. The amplitude for the D* → Dγ decays iswhere h denotes a coupling constant, μ stands for the magnetic moment for D* → Dγ transitions, and p are the D*-meson and photon momenta, respectively, and ϵ is the polarisation vector of the photon. The three amplitudes for and decays arewhere and the functions that denote the Breit–Wigner amplitude for the D* mesons are

A small possible distortion of the Breit–Wigner shape of the D* meson due to three-body final-state interactions is neglected in the model. The impact of the energy-dependence of the D* meson self-energy is found to be insignificant. The decays of the state into the D+D+π− final state via off-shell D*0 → D+π− decays are highly suppressed and are not considered here. The last terms in Eqs. (35) and (36) imply the same amplitudes with swapped momenta.

The state is assumed to be produced unpolarised; therefore, the squared absolute value of the decay amplitudes with pions in the final state, averaged over the initial spin state arewhereand stands for the Källén function[97]. The additional factor of 2! in the denominator of Eq. (39) is due to the presence of two identical particles (D0) in the final state. The squared absolute values of the decay amplitude with a photon in the final state, averaged over the initial spin state isThe coupling constants f and h for the D* → Dπ and D* → Dγ decays are calculated using Eqs. (31)–(34), from the known branching fractions of the D* → Dπ and D* → Dγ decays[10], the measured natural width of the D*+ meson[10,112] and the derived value for the natural width for the D*0 meson[66,81,113]. The magnetic moment μ+ is taken to be 1 and the ratio of magnetic moments μ0/μ+ is calculated according to Refs. [138-140].

Unitarised Breit–Wigner shape

A unitarised three-body Breit–Wigner function is defined aswhere denotes the final state. The decay matrix element for each channel integrated over the three-body phase-space is denoted bywhere is defined by Eqs. (39)–(46) and the unknown coupling constant g is taken out of the expression for . For large values of s, in excess of s*, such as , the functions ϱ(s) are defined aswhere denotes the two-body phase-space function, the constants c1,c2 and c3 are chosen to ensure the continuity of the functions ϱ(s), and a value of is used. The functions ϱ(s) are shown in Supplementary Fig. 5. The complex-valued width is defined via the self-energy function Σ(s)[141]where is again factored out for convenience. The imaginary part of Σ(s) for real physical values of s is computed through the optical theorem as half of the sum of the decay probability to all available channels[142]:The real part of the self-energy function is computed using Kramers–Kronig dispersion relations with a single subtraction[143,144],where the Cauchy principal value (p.v.) integral over ϱtot(s) is understood asand s denotes the threshold value for the channel f. The subtraction is needed since the integral ∫ϱtot(s)/s ds diverges. The term in Eq. (56) corresponds to the choice of subtraction constant such that . The function ξ(s) is shown in Supplementary Fig. 6.

Alternatively, the isoscalar amplitude is constructed using the K-matrix approach[145] with two coupled channels, D*+D0 and D*0D+. The relation reads:where a production vector P and an isoscalar potential K are defined as

The propagation matrix G describes the D*D → D*D rescattering via the virtual loops including the one-particle exchange process[91] and expressed in a symbolical way in Supplementary Eq. (1). where suppressed D*0 → D+π− transition is neglected. The D*+D0 ↔ D*0D+ rescattering occurs due to non-diagonal element of the K-matrix (contact interaction) and non-diagonal elements of the G matrix (long-range interaction). The matrix G and the self-energy function Σ(s) from Eqs. (54) and (56), are related as

Similar to the Flatté function[98] for large values of the parameter, the signal profile exhibits a scaling property[94,99]. For large values of the parameter the width approaches asymptotic behaviour, see Supplementary Fig. 7. The unitarised three-body Breit–Wigner function for decays with parameters mU and obtained from the fit to data is shown in Supplementary Fig. 8. The inset illustrates the similarity of the profile with the single-pole profile in the vicinity of the polewhere and are the mode and FWHM, respectively.

Analytic continuation

Equation (56) defines Σ(s) and the amplitude for real values of s. Analytic continuation to the whole first Riemann sheet is calculated aswhere the integral is understood as in Eq. (58). The search for the resonance pole requires knowledge of the amplitude on the second Riemann sheet denoted by . According to the optical theorem[142], the discontinuity of the inverse amplitude across the unitarity cut is given by :

For the complex-values s, the analytic continuation of ϱtot(s) is needed: the phase-space integral in Eq. (49) is performed over a two-dimensional complex manifold (see discussion on the continuation in Ref. [146]):where the limits of the second integral represent the Dalitz plot borders[111],The integration is performed along straight lines connecting the end points in the complex plane.

DD spectra from decays

The shapes of the D0D0 and D+D− mass spectra from decays are obtained via integration of the expressions from Eqs. (39)–(46) over the s and s12 variables with the amplitude squared, , from Eq. (47):where the lower and upper integration limits for s12 at fixed s and s23 are given in Eq. (66). The function f(s) is introduced to perform a smooth cutoff of the long tail of the profile. Cutoffs are chosen to suppress the profile for regions , where and are the mode and FWHM for the distribution. Two cutoff functions f(s) are studied:

A Gaussian cutoff is defined as

A power-law cutoff defined as

Fits to the background-subtracted D0D0π+ mass spectrum using a signal profile of the form show that the parameter δmU is insensitive to the choice of cutoff function when and σ≥1 MeV/c2. The power-law cutoff function with parameters and σ = 1 MeV/c2 is chosen. The shapes for the D0D0 and D+D0 mass distributions are defined as

Low-energy scattering amplitude

The unitarized Breit–Wigner amplitude is formally similar to the low-energy expansion given by Eq. (13) once the factor is divided out

The function iϱtot(s) matches ik up to a slowly varying energy factor that can be approximated by a constant in the threshold region. The proportionality factor w has the dimension of an inverse mass and is found by matching the decay probability to the two-body phase-space expression:where c1 is a coefficient computed in Eq. (50). The comparison of and that validates the matching is shown in Supplementary Fig. 9.

The inverse scattering length is defined as the value of the amplitude in Eq. (72) at the D*+D0 threshold:

The imaginary part is fully determined by the available decay channels, while the real part depends on the constant adjusted in the fit. The quadratic term, k2 in Eq. (72), corresponds to the linear correction ins since k2 = (s − sth)/4 for the non-relativistic case. Hence, the slope of the linear term in the amplitude is related to the effective range as follows:

Mass splitting for the isotriplet

While the degrees of freedom of the light diquark for the isoscalar state are similar to those for the state, for the isotriplet (,,) the light diquark degrees of freedom would be similar to those for the (anti)triplet. Assuming that the difference in the light quark masses, the Coulomb interaction of light quarks in the diquark, and the Coulomb interaction of the light diquark with the c-quark are responsible for the observed mass splitting in the Σc isotriplet, the masses for the Σc states can be written aswhere mΣ is a common mass parameter; the second and third terms describe the contribution from the light quark masses, mu and md, into the mass splitting; terms proportional to a describe Coulomb interactions of light quarks in the diquark; terms proportional to b describe the Coulomb interactions of the diquark with the c-quark; and qq denotes the charge of the q-quark. Similar expressions can be written for the isotriplet:where is the common mass parameter, and qcc = 2qc is the charge of a cc diquark. Using the known masses of the light quarks and Σc states[10] and taking and , the mass splitting for the isotriplet is estimated to beThe validity of this approach is tested by comparing the calculated mass splitting between Σbp and Σbm states of −6.7 ± 0.7 MeV/c2 with the measured value of −5.1 ± 0.2 MeV/c2[10]. Based on the small observed difference, an additional uncertainty of 0.8 MeV/c2 is added in quadrature to the results from Eqs. (83) and (84), and finally one gets

These results agree with the assigned uncertainty with results based on a more advanced model from Ref. [147].

Supplementary information

Editorial Assessment Report