Study of the doubly charmed tetraquark $T_{cc}^+$

An exotic narrow state in the $D^0D^0\pi^+$ mass spectrum just below the $D^{*+}D^0$ mass threshold is studied using a data set corresponding to an integrated luminosity of 9 fb$^{-1}$ acquired with the LHCb detector in proton-proton collisions at centre-of-mass energies of 7, 8 and 13 TeV. The state is consistent with the ground isoscalar $T^+_{cc}$ tetraquark with a quark content of $cc\bar{u}\bar{d}$ and spin-parity quantum numbers $\mathrm{J}^{\mathrm{P}}=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 confirmed by the $D^0\pi^+$ mass distribution. The mass of the resonance and its coupling to the $D^{*}D$ system are analysed. Resonance parameters including the pole position, scattering length, effective range and compositeness are measured to reveal important information about the nature of the $T^+_{cc}$ 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 (q 1q2 ) and baryons (q 1 q 2 q 3 ) have been actively discussed since the birth of the quark model [1][2][3][4][5][6][7]. Since the discovery of the χ c1 (3872) state [8] many tetraquark and pentaquark candidates, listed in Table 1, have been observed [9][10][11][12][13][14][15][16][17][18]. For all but the X 0 (2900) and X 1 (2900) states the minimal quark content implies the presence of either a cc or bb quark-antiquark pair. The masses of many tetra-and pentaquark states are close to mass thresholds, e.g. D ( * ) D ( * ) or B ( * ) B ( * ) , where D ( * ) or B ( * ) represents a hadron containing a charm or beauty quark, respectively. Therefore, these states are likely to be hadronic molecules [15,[19][20][21] where colour-singlet hadrons are bound by residual nuclear forces, such as the exchange of a pion or ρ meson [22], 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 cc quarks is the Ξ ++ cc (ccu) baryon, a long-lived, weakly-decaying compact object [23,24]. The recently observed X(6900) structure in the J/ψJ/ψ mass spectrum [25] belongs to both categories simultaneously. Its proximity to the χ c0 χ c1 threshold could indicate a molecular structure [26,27]. 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 [28][29][30][31]. The existence and properties of Q 1 Q 2q1q2 states with two heavy quarks and two light antiquarks have been widely discussed for a long time [60][61][62][63][64][65]. 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 Q 1 Q 2 form a point-like color-antitriplet object, analogous to an antiquark, and Table 1: Tetra-and pentaquark candidates and their plausible valence quark content.

States
Quark content X 0 (2900), X 1 (2900) [32,33]  as a result the Q 1 Q 2q1q2 system has similar degrees of freedom for its light quarks as an antibaryon with a single heavy quark, e.g. the Λ − c or Λ 0 b antibaryons. The beauty quark is considered heavy enough to sustain the existence of a bbud state that is stable with respect to the strong and electromagnetic interactions with a mass of about 200 MeV below the BB * mass threshold. In the case of the bcud and ccud 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 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 [9], with cd and cu quark content, respectively. 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 with ccud quark content is reported in Ref. [98].
In the work presented here, the properties of the T + cc state are studied by constructing a dedicated amplitude model that accounts for the D * + D 0 and D * 0 D + decay channels. In addition, the mass spectra of other DD ( * ) and opposite-sign DD ( * ) combinations are explored. Furthermore, production-related observables, such as the event multiplicity and transverse momentum (p T ) 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 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 [99,100] 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
T + cc signal in the D 0 D 0 π + mass spectrum The D 0 D 0 π + final state is reconstructed using the D 0 → K − π + decay channel with two D 0 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. [101][102][103][104] and described in detail in Methods. The background not originating from true D 0 mesons is subtracted using an extended unbinned maximum-likelihood fit to the two-dimensional distribution of the masses of the two D 0 candidates from selected D 0 D 0 π + combinations, see Methods and Supplementary Fig. 1(a). The obtained D 0 D 0 π + mass distribution for selected D 0 D 0 π + combinations is shown in Fig. 1. 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 corresponds to the T + cc → D 0 D 0 π + 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 F BW with a Blatt-Weisskopf form factor [105,106] is used in Ref. [98] as the natural resonance profile. That function, while sufficient to  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 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.
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 [107][108][109][110][111][112][113][114][115][116][117]. A unitarised Breit-Wigner profile F U , described in Methods Eq. (48), is used in this analysis. The function F U , is built under two main assumptions.
Assumption 1. The newly observed state has quantum numbers J P = 1 + and isospin I = 0 in accordance with the theoretical expectation for the T + cc ground state.
Assumption 2. The T + cc state is strongly coupled to the D * D channel, which is motivated by the proximity of the T + cc mass to the D * D mass threshold. The derivation of the F U profile relies on the assumed isospin symmetry for the T + cc → D * D decays and the coupled-channel interaction of the D * + D 0 and D * 0 D 0 system as required by unitarity and causality following Ref. [118]. The resulting energy-dependent width of the T + cc state accounts explicitly for the T + cc → D 0 D 0 π + , T + cc → D 0 D + π 0 and T + cc → D 0 D + γ decays. The modification of the D * meson lineshape [119] due to contributions from triangle diagrams [120] to the final-state interactions is neglected. Similarly to the F BW profile, the F U function has two parameters: the peak location m U , 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 T + cc → D * D decay. The detector mass resolution, R, 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, that accounts for a small residual difference between simulation and data [54, 121,122]. The root mean square of the resolution function is around 400 keV/c 2 .
A study of the D 0 π + mass distribution for selected 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 D 0 D 0 π + combinations contain a true D * + meson. Therefore, the background component is parameterised with a product of the two-body phase space function Φ D * + D 0 [123] and a positive polynomial function P n , convolved with the detector resolution function R where n denotes the order of the polynomial function, n = 2 is used in the default fit. The D 0 D 0 π + mass spectrum with non-D 0 background subtracted is shown in Fig. 1 with the result of the fit using a model based on the F U signal profile overlaid. The fit gives a signal yield of 186 ± 24 and a mass parameter relative to the D * + D 0 mass threshold, δm U of −359 ± 40 keV/c 2 . The statistical significances of the observed T + cc → D 0 D 0 π + signal and for the δm U < 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 |g|. With increasing |g|, 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 |g|, the F U signal profile exhibits a scaling property similar to the Flatté function [121,124,125]. The parameter |g| effectively decouples from the fit model, and the model resembles the scattering-length approximation [108]. The likelihood profile for the parameter |g| is shown in Fig. 2, where one can see a plateau at large values. At small values of the |g| parameter, |g| < 1 GeV, the likelihood function is independent of |g| because the resonance is too narrow for the details of the F U signal profile to be resolved by the detector. The lower limits on the |g| parameter of |g| > 7.7 (6.2) GeV at 90 (95) % confidence level (CL) are obtained as the values where the difference in the negative log-likelihood −∆ log L is equal to 1.35 and 1.92, respectively. Smaller values for |g| are further used for systematic uncertainty evaluation.
The mode relative to the D * + D 0 mass threshold, δm, and the full width at half maximum (FWHM), w, for the F BW profile are found to be δm = −361 ± 40 keV/c 2 and w = 47.8 ± 1.9 keV/c 2 , to be compared with those quantities determined for the F U signal profile of δm = −279 ± 59 keV/c 2 and w = 409 ± 163 keV/c 2 . 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 F U 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 F BW function. The obtained mean and root mean square (RMS) values for the parameters δm BW and Γ BW over the ensemble are shown in Table 2. The mass parameter δm BW agrees well with the value determined from data [98]. 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 F BW profile is analysed with a model based on the F U function. The obtained mean and RMS values for the δm U 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 explains 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.

Systematic uncertainties
Systematic uncertainties for the δm U 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 to 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: • 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 [126], a modified Gaussian function with symmetric power-law tails on both sides of the distribution [127,128], a generalised symmetric Student's t-distribution [129,130], a symmetric Johnson's S U distribution [131,132], and a modified Novosibirsk function [133].
• 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 [54, 121, 122, 134-136] 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 * + → 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 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 D 0 D 0 π + combinations without an intermediate D * + meson, a more general family of background models is tested where Φ D 0 D 0 π + denotes the three-body phase-space function [137]. The functions B 0 , B 1 , B 3 and B nm 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 F U signal profile. These coupling constants are calculated from the known branching fractions of the D * → Dπ and D * → Dγ decays [9], the measured natural width of the D * + meson [9,138] and the derived value for the natural width of the D * 0 meson [93,108,139]. To assess the associated systematic uncertainty, a set of alternative models built around the F U profiles, obtained with coupling constants varying within their calculated uncertainties, is studied.
• Unknown value of the |g| parameter. In the baseline fit the value of the |g| parameter is fixed to a large value. To assess the effect of this constraint the fit is repeated using the value of |g| = 8.08 GeV, that corresponds to −2∆ log L = 1 for the most conservative likelihood profile for |g| that accounts for the systematic uncertainty. The change of 7 keV/c 2 of the δm U 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 [140]. The accuracy of the procedure has been checked using 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 [141]. This uncertainty is propagated for the parameters of interest using simulated samples, with momentum scale corrections of (1 ± δα) 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 [142,143]. The amount of the material traversed in the tracking system by a charged particle is known to 10% accuracy [144]. 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 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 [9]. This procedure produces negligible uncertainties for the δm U parameter due to imprecise knowledge of the D 0 mass. However, the small uncertainty of 2 keV/c 2 for the known D * + − D 0 mass difference [9,138,145] directly affects the values of these parameters and is assigned as corresponding systematic uncertainty.
For the lower limit on the parameter |g|, 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: |g| > 5.1 (4.3) GeV at 90 (95) % CL.

Results
Studying the D 0 π + mass distribution for T + cc → D 0 D 0 π + decays allows testing the hypothesis that the T + cc → D 0 D 0 π + decay proceeds through an intermediate off-shell D * + meson. The background-subtracted D 0 π + mass distribution for selected D 0 D 0 π + candidates with the D 0 D 0 π + mass with respect to the D * + D 0 mass threshold, δm D 0 D 0 π + , below zero is shown in Fig. 3. Both D 0 π + combinations are included in this plot. The two-dimensional distribution of the mass of one D 0 π + combination versus the mass of another D 0 π + combination is presented in Supplementary Fig. 10.
A fit is performed to this distribution with a model containing signal and background components. The signal component is derived from the A U amplitude, see Methods Eq. (49), and is convolved with a detector resolution for the D 0 π + mass. This detector resolution function is modelled with a modified Gaussian function with power-law tails on both sides of the distribution [127,128] and parameters taken from simulation. Similarly to the correction used for the D 0 D 0 π + mass resolution function R, the width of the Gaussian function is corrected by a factor of 1.06 which is determined by studying large samples of D * + → D 0 π + decays. The root mean square of the resolution function is around 220 keV/c 2 . The shape of the background component is derived from data for δm D 0 D 0 π + > 0.6 MeV/c 2 . The fit results are overlaid in Fig. 3. The background component vanishes in the fit, and the D 0 π + spectrum is consistent with the hypothesis that the T + cc → D 0 D 0 π + 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 T + cc signal to the D * + D 0 mass threshold, and the small energy release in the D * + → D 0 π + decay, the D 0 D 0 mass distribution from the T + cc → D 0 D 0 π + decay forms a narrow peak just above the D 0 D 0 mass threshold. In a similar way, a peaking structure in the D + D 0 mass spectrum just above the D + D 0 mass threshold is expected from T + cc → D + D 0 π 0 and T + cc → D + D 0 γ decays, both proceeding via off-shell intermediate D * + D 0 and D * 0 D + states. The D 0 D 0 and D + D 0 final states are reconstructed and selected similarly to the D 0 D 0 π + final state, where the D + → K − π + π − decay channel is used. The background-subtracted D 0 D 0 and D + D 0 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 F DD described in Methods Eqs. (71) and (72) and obtained via integration of the matrix elements for the T + cc → DDπ/γ decays with the F U profile, and a background component, parameterised as a product of the two-body phase space function Φ DD and  a positive linear function P 1 . The fit results are overlaid in Fig. 4(top). The signal yields in the D 0 D 0 and D + D 0 spectra are found to be 263 ± 23 and 171 ± 26, respectively. The statistical significance of the observed T + cc → D 0 D 0 X and T + cc → D + D 0 X signals, where X stands for non-reconstructed pions or photons, is estimated using Wilks' theorem [146] and is found to be in excess of 20 and 10 standard deviations, respectively. The relative yields for the signals observed in the D 0 D 0 π + , D 0 D 0 and D 0 D + mass spectra agree with the expectations of the model described in Methods where the decay of an isoscalar T + cc state via the D * D channel with an intermediate off-shell D * meson is assumed. The observation of the near-threshold signals in the D 0 D 0 and D + D 0 mass spectra, along with the signal shapes and yields, all agree with the isoscalar T + cc hypothesis for the narrow signal observed in the D 0 D 0 π + mass spectrum. However, an alternative interpretation could be that this state is the I 3 = 0 component of aT cc isotriplet (T 0 cc ,T + cc , T ++ cc ) with ccuu, ccud and ccdd quark content, respectively. Assuming that the observed  peak corresponds to theT + cc component and using the estimates for theT cc mass splitting from Methods Eqs. (86) and (87), the masses of theT 0 cc andT ++ cc states are estimated to be slightly below the D 0 D * 0 and slightly above the D + D * + mass thresholds, respectively: With these mass assignments, assuming equal production of all threeT cc components, theT 0 cc state would be an extra narrow state that decays into the D 0 D 0 π 0 and D 0 D 0 γ final states via an off-shell D * 0 meson. These decays would contribute to the nar-row near-threshold enhancement in the D 0 D 0 spectrum, and increase the signal in the D 0 D 0 mass spectrum by almost a factor of three. TheT ++ cc state would decay via an on-shell D * + mesonT ++ cc → D + D * + , therefore it could be a relatively wide state, with width up to a few MeV [147]. Therefore, it would manifest itself as a peak with a moderate width in the D + D 0 π + mass spectrum with a yield comparable to that of theT + cc → D 0 D 0 π + decays. In addition, it would contribute to the D + D 0 mass spectrum, tripling the contribution from theT + cc decays. However, due to the larger mass of theT ++ cc state and its larger width, this contribution should be wider, making it more difficult to disentangle from the background. Finally, theT ++ cc state would make a contribution to the D + D + spectrum with a yield similar to the contribution from T + cc → D 0 D + π 0 /γ decays to the D 0 D + spectrum, but wider. The mass spectra for D + D + and D + D 0 π + 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 + D 0 π + mass spectra is therefore a strong argument in favour of the isoscalar nature of the observed peak in the D 0 D 0 π + mass spectrum. The interference between two virtual channels for the T + cc → D 0 D 0 π + 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 T + cc → D + D 0 π 0 decays gives the contribution of 45% for the interference between the (D * + → D + π 0 ) D 0 and (D * 0 → D 0 π 0 ) D + channels. For T + cc → D + D 0 γ decays the role of the interference between the (D * + → D + γ) D 0 and (D * 0 → D 0 γ) D + channels is estimated by equating to zero the F + F * 0 and F * + F 0 terms in Methods Eqs. (46) and (47). The interference contribution is found to be 33%.
Using the model described earlier and results of the fit to the D 0 D 0 π + mass spectrum, the position of the amplitude poleŝ in the complex plane, responsible for the appearance of the narrow structure in the D 0 D 0 π + mass spectrum is determined. The pole parameters, mass m pole and width Γ pole , are defined through the pole locationŝ as The pole locationŝ is a solution of the equation where A II U (s) denotes the amplitude on the second Riemann sheet defined in Methods Eq. (65). For large coupling |g| the position of the resonance pole is uniquely determined by the parameter δm U , i.e. the binding energy and the width of the D * + meson. Figure 5 shows the complex plane of the δ √ s variable, defined as All possible positions of the pole for |g| m D 0 + m D * + 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 * + D 0 mass threshold), the width gets narrower; and when the parameter δm U approaches zero, the pole touches the D 0 D * + cut and moves to the other complex sheet, i.e. the state becomes virtual. For smaller values of |g|, the pole is located between the limiting curve and the s = 0 line. The pole parameters are found to be where the first uncertainty is due to the δm U parameter and the second is due to the unknown value of the |g| parameter. The peak is well separated from the D * + D 0 threshold in the D 0 D 0 π + mass spectrum. Hence, as for an isolated narrow resonance, the parameters of the pole are similar to the visible peak parameters, namely the mode δm and FWHM w.
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 valie. For example, as shown earlier the data are consistent with a wide range of FHWM w values for the signal profile. Therefore the pole width Γ pole is based mainly on the T + cc amplitude model and the value of the m U parameter determined from the fit to the D 0 D 0 π + mass spectrum.
A study of the behaviour of the A U (s) amplitude in the vicinity of the D * + D 0 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 [148], i.e.
where k is the wave number. For δ √ s −Γ D * + the inverse amplitude from Eq. (48) matches Eq. (13) up to a scale parameter obtained numerically, see Methods Eq. (75). 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 [149]; however, in this case the non-zero imaginary part is related to the lower threshold, T + cc → D 0 D 0 π + , 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 [15], For the A U amplitude the effective range r is non-positive and proportional to |g| −2 , see Methods Eq. (77). Its value is consistent with zero for the baseline fit. An upper limit on the −r value is set as 0 ≤ −r < 11.9 (16.9) fm at 90 (95)% CL .
corresponding to the openings of the D 0 D + γ, D 0 D + π 0 and D 0 D 0 π + decay channels, are outside of the displayed region.
The Weinberg compositeness criterion [150,151] 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.
Another estimate of the characteristic size is obtained from the value of the binding energy ∆E. Within the interpretation of the T + cc state as a bound D * + D 0 molecular-like state, the binding energy is ∆E = −δm U . The characteristic momentum scale γ [15] is estimated to be where µ is the reduced mass of the D * + D 0 system. This value of the momentum scale in turn corresponds to a characteristic size R ∆E of the molecular-like state, which is consistent with the R a estimate from the scattering length. For high-energy hadroproduction of a state with such a large size, R a or R ∆E , one expects a strong dependency of the production rate on event multiplicity, similar to that observed for the χ c1 (3872) state [152]. The background-subtracted distribution of the number of tracks reconstructed in the vertex detector, N tracks , is shown in Fig. 6(left) together with the distributions for low-mass D 0 D 0 pairs with m D 0 D 0 < 3.87 GeV/c 2 and low-mass D 0 D 0 pairs with mass 3.75 < m D 0 D 0 < 3.87 GeV/c 2 . The former is dominated by pp → ccX production, while the latter is presumably dominated by the double parton scattering process [101,153]. The chosen interval for D 0 D 0 pairs includes the region populated by the χ c1 (3872) → D 0 D 0 π 0 /γ 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 [152]. It is noteworthy that the track multiplicity distribution for the T + cc state differs from that of the low-mass D 0 D 0 pairs, in particular no suppression at large multiplicity is observed. A p-value for the consistency of the track multiplicity distributions for T + cc production and low-mass D 0 D 0 pairs is found to be 0.1%. It is interesting to note that the multiplicity distribution for T + cc production and the one for D 0 D 0 -pairs with 3.75 < m D 0 D 0 < 3.87 are consistent with a corresponding p-value of 12%. The similarity between T + cc production, which is inherently a single parton scattering process, and the distribution for process dominated by a double parton scattering is surprising.
The transverse momentum spectrum for the T + cc state is compared with those for the low-mass D 0 D 0 and D 0 D 0 pairs in Fig. 6(right). The p-values for the consistency of the p T spectra for the T + cc state and low-mass D 0 D 0 pairs are 1.4%, and 0.02% for low-mass D 0 D 0 pairs. More data are needed for further conclusions.
The background-subtracted D 0 D 0 mass distribution in a wider mass range is shown in Fig. 7 together with a similar distribution for D 0 D 0 pairs. In the D 0 D 0 mass spectrum the near-threshold enhancement is due to χ c1 (3872) → D 0 D 0 π 0 and χ c1 (3872) → D 0 D 0 γ decays via intermediate D * 0 mesons [104]. This structure is significantly wider than the structure in the D 0 D 0 mass spectrum from T + cc → D 0 D 0 π + decays primarily due to the larger natural width and smaller binding energy for the χ c1 (3872) state [121,122]. With more data, and with a better understanding of the dynamics of χ c1 (3872) → D 0 D 0 π 0 /γ decays, and therefore of the corresponding shape in the D 0 D 0 mass spectrum, it will be possible to estimate the relative production rates for the T + cc and χ c1 (3872) states. Background-subtracted D 0 D 0 π + and D 0 D + mass distributions together with those for D 0 D 0 π + and D 0 D − are shown in Supplementary Figs. 3 and 4.

Discussion
The exotic narrow tetraquark state T + cc observed in Ref. [98] 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 D 0 π + mass distribution indicates that the T + cc → D 0 D 0 π + decay proceeds via an intermediate off-shell D * + meson. Together with the proximity of the state to the D * D 0 mass threshold this favours the spin-parity quantum numbers J P to be 1 + . Narrow near-threshold structures are observed in the D 0 D 0 and D 0 D + mass spectra with high significance. These are found to be consistent with originating from off-shell T + cc → D * D decays followed by the D * → Dπ and D * → Dγ decays. No signal is observed in the D + D 0 π + 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 J P = 1 + ccud-tetraquark ground state. A dedicated unitarised three-body Breit-Wigner amplitude is built on the assumption of strong isocalar  coupling of the axial-vector T + cc 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 T + cc state, relative to the D * + D 0 mass threshold, is determined to be where the first uncertainty is statistic and the second systematic. The lower limit on the absolute value of the coupling constant of the T + cc state to the D * D system is |g| > 5.1 (4.3) GeV at 90 (95) % CL .
Using 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 be a = − (7.16 ± 0.51) + i (1.85 ± 0.28) fm , −r < 11.9 (16.9) fm at 90 (95)% CL , The characteristic size calculated from the binding energy is R ∆E = 7.49 ± 0.42 fm. This value is consistent with the estimation from the scattering length, R a = 7.16±0.51 fm. Both R ∆E and R a 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 D 0 D 0 π + threshold, atŝ = m pole − i 2 Γ pole , where where the first uncertainty accounts for statistical and systematic uncertainties for the δm U parameters, and the second is due to the unknown value of the |g| parameter. The pole position, scattering length, effective range and compositeness form a complete set of observables related to the T + cc → D 0 D 0 π + reaction amplitude, which are crucial for inferring the nature of the T + cc tetraquark. Unlike in the prompt production of the χ c1 (3872) state, no suppression of the T + cc production at high track multiplicities is observed relative to the low-mass D 0 D 0 pairs. The observed similarity with the multiplicity distribution for the low-mass D 0 D 0 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. D 0 → K − π + π + π − , detailed studies of the properties of this new state and its production mechanisms could be possible.
In conclusion, the T + cc tetraquark observed in D 0 D 0 π + decays is studied in detail, using a unitarised model that accounts for the relevant thresholds by taking into account the D 0 D 0 π + and D 0 D + π 0 (γ) decay channels with intermediate D * resonances. This model is found to give an excellent description of the D 0 π + mass distribution in the T + cc → D 0 D 0 π + decay and of the threshold enhancements observed in the D 0 D 0 and D 0 D + spectra. Together with the absence of a signal in the D 0 D + and D + D 0 π + mass distributions this provides a strong argument for interpreting the observed state as the isoscalar T + cc tetraquark with spin-parity J P = 1 + . The precise T + cc 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 T + cc tetraquark. In addition, the counter-intuitive dependence of the production rate on track multiplicity will pose a challenge for theoretical explanations.

Experimental setup
The LHCb detector [99,100] 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 [154]. 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 [155] with a specific LHCb configuration [156]. Decays of unstable particles are described by EvtGen [157], in which final-state radiation is generated using Photos [158]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [159] as described in Ref. [160].

Event selection
The D 0 D 0 , D 0 D + and D 0 D 0 π + final states are reconstructed using the D 0 → K − π + and D + → K − π + π + decay channels. The selection criteria are similar to those used in Refs. [101][102][103][104]. Kaons and pions are selected from well-reconstructed tracks within the acceptance of the spectrometer that are identified using information from the ringimaging 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 D 0 and D + candidates, referred to as D hereafter. The resulting D candidates are required to have good vertex quality, mass within ±65 MeV/c 2 and ±50 MeV/c 2 of the known D 0 and D + masses [9], 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 D 0 and D + candidates consistent with originating from a common primary vertex are combined to form D 0 D 0 and D 0 D + candidates. The resulting D 0 D 0 candidates are combined with a pion 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 * + mass by more than 155 MeV/c 2 . For each D 0 D 0 , D 0 D + and D 0 D 0 π + candidate a kinematic fit [161] 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 D 0 D 0 π + , 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 [101]. 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 + D 0 π + , D 0 D 0 , D 0 D − and D 0 D 0 π + 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 D 0 D 0 π + , D 0 D 0 and D 0 D + 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 [162], which is based on an extended unbinned maximum-likelihood fit to these two-dimensional distributions with the function described in Ref. [101]. This function consists of four components: • a component corresponding to genuine D 1 D 2 pairs and described as a product of two signal functions, each parameterised with a modified Novosibirsk function [133]; • 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 D 0 candidates as discriminating variables. The D 0 D 0 π + mass distributions for each of the subtracted background components are presented in Supplementary Fig. 2, where fit results with background-only functions B 10 , defined in Eq. (3) are overlaid.
Resolution model for the D 0 D 0 π + mass In the vicinity of the D * + D 0 mass threshold the resolution function R for the D 0 D 0 π + 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/c 2 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 [54, 121,122]. The root mean square of the resolution function is around 400 keV/c 2 .
Matrix elements for T + cc → DDπ/γ decays Assuming isospin symmetry, the isoscalar vector state T + cc that decays into the D * D final state can be expressed as Therefore, the S-wave amplitudes for the T + cc → D * + D 0 and T + cc → D * 0 D + decays have different signs where g is a coupling constant, T + cc is the polarisation vector of the T + cc particle and D * 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 T + cc masses significantly above the D * D threshold, higher-order waves also need to be considered. The amplitudes for the D * → Dπ decays are written as where f denotes a coupling constant, and p D stands for the momentum of the D meson. The amplitude for the D * → Dγ decays is where h denotes a coupling constant, µ stands for the magnetic moment for D * → Dγ transitions, p D * and p γ are the the D * -meson and photon momenta, respectively, and γ is the polarisation vector of the photon. The three amplitudes for T + cc → πDD and T + cc → γDD decays are where s ij = p 2 ij = (p i + p j ) 2 and the F 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 T + cc 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 T + cc state is assumed to be produced unpolarized, therefore the squared absolute value of the decay amplitudes with pions in the final state, averaged over the initial spin-state are E = (pp 12 )(pp 13 )(p 2 p 12 )(p 3 p 13 ) ss 12 s 13 and λ (x, y, z) stands for the Källén function [123]. The additional factor of 2! in the denominator of Eq. (39) is due to the presence of two identical particles (D 0 ) 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 is The 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 [9], the measured natural width of the D * + meson [9,138] and the derived value for the natural width for the D * 0 meson [93,108,139]. The magnetic moment µ + is taken to be 1 and the ratio of magnetic moments µ 0 /µ + is calculated according to Refs. [163][164][165].

Unitarised Breit-Wigner shape
A unitarised three-body Breit-Wigner function is defined as where f ∈ {D 0 D 0 π + , D 0 D + π 0 , D 0 D + γ} denotes the final state. The decay matrix element for each channel integrated over the three-body phase space is denoted by where |M f | 2 is defined by Eqs.
where Φ D * D (s) denotes the two-body phase-space function, the constants c 1 , c 2 and c 3 are chosen to ensure the continuity of the functions f (s), and a value of √ s * = 3.9 GeV/c 2 is used. The functions f (s) are shown in Supplementary Fig. 5. The complex-valued widtĥ Γ(s) is defined via the self-energy function Σ(s) [166] im UΓ (s) ≡ |g| 2 Σ(s) , where |g| 2 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 [167]: The real part of the self-energy function is computed using Kramers-Kronig dispersion relations with a single subtraction [168,169], where the Cauchy principal value (p.v.) integral over tot (s) is understood as and s f denotes the threshold value for the channel f . The subtraction is needed since the integral tot (s)/s ds diverges. The term ξ(m 2 U ) in Eq. (57) corresponds to the choice of subtraction constant such that A U (m 2 U ) = 0. The function ξ(s) is shown in Supplementary Fig. 6.
Alternatively, the isoscalar amplitude A U is constructed using the K-matrix approach [170] with two coupled channels, D * + D 0 and D * 0 D + . 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 [118] and expressed in a symbolical way in Supplementary Eq. (1). where suppressed D * 0 → D + π − transition is neglected. The D * + D 0 ↔ D * 0 D + 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. (55) and (57), are related as Similar to the Flatté function [124] for large values of the |g| parameter, the F U signal profile exhibits a scaling property [121,125]. For large values of the |g| parameter the width approaches asymptotic behaviour, see Supplementary Fig. 7. The unitarised three-body Breit-Wigner function F U for T + cc → D 0 D 0 π + decays with parameters m U and |g| 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 pole where m and w are the mode and full width at half maximum, respectively.

Analytic continuation
Equations (57) define Σ(s) and the amplitude A U (s) for real values of s. Analytic continuation to the whole first Riemann sheet is calculated as where the integral is understood as in Eq. (59). The search for the resonance pole requires knowledge of the amplitude on the second Riemann sheet denoted by A II U . According to the optical theorem [167], the discontinuity of the inverse amplitude across the unitarity cut is given by i |g| 2 tot (s): For the complex-values s, the analytic continuation of tot (s) is needed: the phase space integral in Eq. (50) is performed over a two-dimensional complex manifold D (see discussion on the continuation in Ref. [171]): where the limits of the second integral represent the Dalitz plot borders [137], The integration is performed along straight lines connecting the end points in the complex plane.
DD spectra from T + cc → DDπ/γ decays where the lower and upper integration limits for s 12 at fixed s and s 23 are given in Eq. (67). The function f C (s) is introduced to perform a smooth cutoff of the long tail of the T + cc profile. Cutoffs are chosen to suppress the profile for regions | √ s − m| w, where m and w are the mode and FWHM for the F U (s) distribution. Two cut-off functions f C (s) are studied: 2. A power-law cut off f P C (s) defined as Fits to the background-subtracted D 0 D 0 π + mass spectrum using a signal profile of the form F U (s) × f C (s) show that the parameter δm U is insensitive to the choice of cut-off function when x c ≥ m D * 0 + m D + and σ c ≥ 1 MeV/c 2 . The power-law cut-off function f P C (s) with parameters x c = m D * 0 + m D + and σ c = 1 MeV/c 2 is chosen. The shapes for the D 0 D 0 and D + D 0 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 1 2 |g| 2 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 c 1 is a coefficient computed in Eq. (51). The comparison of A −1 NR and A −1 U × 2w/ |g| 2 that validates the matching is shown in Supplementary Fig. 9.
The inverse scattering length is defined as the value of the amplitude in Eq. (73) at the D * + D 0 threshold: The imaginary part is fully determined by the available decay channels, while the real part depends on the constant ξ(m 2 U ) adjusted in the fit. The quadratic term, k 2 in Eq. (73), corresponds to the linear correction in s since k 2 = (s − s th )/4 for the non-relativistic case. Hence, the slope of the linear term in the A −1 U amplitude is related to the effective range as follows: Mass splitting for theT cc isotriplet While the degrees of freedom of the light diquark for the isoscalar T + cc state are similar to those for the Λ − c state, for theT cc isotriplet (T 0 cc ,T + cc ,T ++ cc ) the light diquark degrees of freedom would be similar to those for the Σ c (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 as where m Σ is a common mass parameter; the second and third terms describe the contribution from the light quark masses, m u and m d , 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 q q denotes the charge of the q-quark. Similar expressions can be written for theT cc isotriplet: where mT cc is the common mass parameter, q q = −q q and q cc = 2q c is the charge of a cc diquark. Using the known masses of the light quarks and Σ c states [9] and taking a = a and b = b, the mass splitting for theT cc isotriplet is estimated to be The validity of this approach is tested by comparing the calculated mass splitting between Σ + b and Σ − b states of −6.7 ± 0.7 MeV/c 2 with the measured value of −5.1 ± 0.2 MeV/c 2 [9]. Based on the small observed difference, an addition uncertainty of 0.8 MeV/c 2 is added in quadrature to the results from Eqs. (84) and (85), and finally one gets These results agree within the assigned uncertainty with results based on a more advanced model from Ref. [172].

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 Figs. 1, 3 and 4 have been added to the HEPData record at https://www.hepdata.net/record/ins1915358.

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.

Author Contribution Statement
All contributing authors, as listed at the end of this manuscript, have contributed to the publication, being variously involved in the design and the construction of the detector, in writing software, calibrating sub-systems, operating the detector and acquiring data and finally analysing the processed data.  Candidates/(5    Fig. 3: Mass distributions for D 0 D 0 π + and D 0 D 0 π − candidates. Background-subtracted D 0 D 0 π + and D 0 D 0 π − mass distributions. 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.     [MeV] Supplementary Fig. 9: Comparison of the A U and A NR amplitudes. The real and imaginary parts of the inverse A U and A NR amplitudes. The yellow band correspond to the pole position and vertical dashed lines show the D * + D 0 and D * 0 D + mass thresholds.  Fig. 10: Two-dimensional D 0 π + mass distribution. Background-subtracted two-dimensional D 0 π + mass distribution for the D 0 D 0 π + events with δm D 0 D 0 π + ≤ 0. Dashed vertical and horizontal lines indicate the known D * + mass. Red and dashed blue lines show the boundary corresponding to δm D 0 D 0 π + = 0 and δm D 0 D 0 π + = −359 keV/c 2 , respectively.

Propagation matrix G
The propagation matrix G describes the D * D → D * D rescattering via the virtual loops including the one-particle exchange process and expressed in a symbolical way as where suppressed D * 0 → D + π − transition is neglected.      Fig. 10: Two-dimensional D 0 π + mass distribution. Background-subtracted two-dimensional D 0 π + mass distribution for the D 0 D 0 π + events with δm D 0 D 0 π + ≤ 0. Dashed vertical and horizontal lines indicate the known D * + mass. Red and dashed blue lines show the boundary corresponding to δm D 0 D 0 π + = 0 and δm D 0 D 0 π + = −359 keV/c 2 , respectively.

Propagation matrix G
The propagation matrix G describes the D * D → D * D rescattering via the virtual loops including the one-particle exchange process and expressed in a symbolical way as where suppressed D * 0 → D + π − transition is neglected.