The decays of b hadrons into final states involving two open-charm hadrons form a large family of topologically similar processes that include many intermediate states such as charmonia, highly excited D (s) states, and possible exotic hadrons. The Dalitz plot distributions of B 0 ! D 0 D K + , B + ! D 0 D 0 K + and B + ! D + D K + decaysfoot_0 have already been explored by the Belle [1], BaBar [2] and LHCb collaborations [3,4]. In these studies, the discovery of the charm-strange meson D s1 (2700) + , the charmonium-like state c0 (3930), and the open-charm tetraquark state X 0,1 (2900), were reported, prompting many theoretical investigations into the internal structure of these states [5].
The decays B + ! D D + s ⇡ + and B 0 ! D 0 D + s ⇡ are yet to be explored. They are ideal to study excited D mesons (D ⇤⇤ ) with natural spin-parity, to test isospin symmetry in the charged and neutral D⇡ resonances, and to test quantum chromodynamics (QCD) predictions [6]. The D ⇤ (2007) 0 , D ⇤ (2010) + , D ⇤ 0 (2300), and D ⇤ 2 (2460) mesons are already well-established. The D ⇤ 1 (2600) 0 and D ⇤ J (3000) 0 mesons were recently discovered in the inclusive proton-proton (pp) collisions and in B decays [7], while their charged isospin partners have not been observed, although some measurements suggest their existence [8]. These states could also be explored in B ! DD + s ⇡ decays. Figure 1 shows the Feynman diagrams of the dominant tree-level amplitudes contributing to the two decays.
Studies of B ! DD + s ⇡ decays also provide an excellent opportunity to search for exotic hadrons decaying into the D + s ⇡ and DD + s final states. The discoveries of the D ⇤ s0 (2317) + [9] and D s1 (2460) + [10] states prompted speculation that they may have a tetraquark component [6,7]. No evidence for isospin partners has been found in explicit searches [11,12], but if they exist they should contribute to the B ! DD + s ⇡ decays. The D0 collaboration claimed evidence for an X(5568) state [13,14], which however was not confirmed by other experiments [15][16][17][18]. An open-charm tetraquark state with four di↵erent quark flavors, analogous to the X(5568) state, is predicted by the diquark-antidiquark model [19,20], and can be investigated in D + s ⇡ ± final states. In particular, the D + s ⇡ + channel presents an attractive potential according to lattice QCD calculation [21], which can form a tetraquark resonance. Some theoretical studies on the X 0,1 (2900) state suggest searching for a potential doubly charged charm-strange tetraquark candidate, together with its neutral isospin partner, in the D + s ⇡ + [csu d] and D + s ⇡ [csūd] final states [22][23][24][25][26][27]. In addition, searches for possible DD + s resonances are well-motivated by the recent observations of the open-strange hidden-charm tetraquark state Z cs (3985) decaying into D ⇤ D + s + DD ⇤+ s at BESIII [28,29], as well as the Z cs (4000) and Z cs (4220) states decaying into KJ/ at LHCb [30].
In this paper, an amplitude analysis of B + ! D D + s ⇡ + and B 0 ! D 0 D + s ⇡ decays is presented for the first time, revealing the contributions of D⇡ resonances in the two decays, and allowing searches for possible exotic states. As the two channels are closely related by isospin, three di↵erent fit scenarios are performed: a fit performed independently on the two decay channels is called the separate fit; a simultaneous fit of the B 0 ! D 0 D + s ⇡ and B + ! D D + s ⇡ + decays, by assuming that all the D⇡ resonances in the two decays are isospin-related, is denoted the simultaneous D⇡ fit; a simultaneous fit, in which all the parameters of D⇡ states and potential D + s ⇡ or DD + s resonances are shared between the two decays, is called the full simultaneous fit. The separate fit and simultaneous D⇡ fit are discussed in this paper, while the full simultaneous fit, which is considered as the default result of this analysis, is described in Ref. [31]. The analysis is based on pp collision data collected using the LHCb detector, corresponding to a total integrated luminosity of 3 fb 1 at center-of-mass energies p s = 7, 8 TeV, referred to as Run 1, and 6 fb 1 at p s = 13 TeV, referred to as Run 2. The paper is organized as follows. A brief introduction of LHCb detector and reconstruction and simulation procedures is provided in Sec. 2. The event selection criteria are shown in Sec. 3, and signal and background yields are determined in Sec. 4. The formalism of amplitude analysis is summarized in Sec. 5. The signal e ciency and background models in amplitude analysis are studied in Sec. 6. The separate fit result is presented in Sec. 7, while the result from simultaneous D⇡ fit is provided in Sec. 8. The systematic uncertainties are evaluated in Sec. 9. All the results are summarized in Sec. 10.
The LHCb detector [32,33] 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 siliconstrip 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. 2 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. Di↵erent types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors. 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.
At the hardware trigger stage, events are required to have a muon with high p T or a hadron, photon or electron with high transverse energy in the calorimeters. For hadrons, the transverse energy threshold is 3.5 GeV. The software trigger requires a two-, three-or four-track secondary vertex with a significant displacement from any primary pp interaction vertex. At least one charged particle must have a transverse momentum p T > 1.6 GeV and be inconsistent with originating from a PV. A multivariate algorithm [34,35] is used for the identification of secondary vertices consistent with a decay of a b hadron.
Simulation is used to model the e↵ects of the detector acceptance and the imposed selection requirements. In the simulation, pp collisions are generated using Pythia 8 [36,37] with a specific LHCb configuration [38]. Decays of unstable particles are described by EvtGen [39], in which final-state radiation is generated using Photos [40]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [41] as described in Ref. [42]. The underlying pp interaction is reused multiple times, with an independently generated signal decay for each [43].
The particle identification (PID) response for charged tracks in the simulated samples is corrected based on special samples of D ⇤+ ! D 0 ⇡ + , D 0 ! K ⇡ + decays. For each PID response of a track, the unbinned four-dimensional probability density functions (PDF) for the data, p data (x|p T , ⌘, N tr ), and for the simulated samples p sim (x|p T , ⌘, N tr ) are extracted based on a kernel density estimation [44], where x is the PID response, p T and ⌘ are the transverse momentum and pseudorapidity of the track, and N tr is the number of tracks in the event. The cumulative distribution functions for the data P data (x|p T , ⌘, N tr ) and for the simulated samples P sim (x|p T , ⌘, N tr ) are determined, and the corrected PID response in the simulated samples is evaluated by transforming the x sim into x corr with
During the transformation, the N tr distribution in the simulated samples is scaled by a factor to match the same distribution in the corresponding datasets. The PID response in the simulated samples shows good agreement with that in the data samples after the correction.
The momentum scale is calibrated using control samples of J/ ! µ + µ and B + ! J/ K + decays collected concurrently with the data samples used for this analysis [45,46]. The relative uncertainty on the momentum scale is 3 ⇥ 10 4 .
In LHCb, trigger decisions are associated with reconstructed particles. Selection requirements can therefore be made on the trigger selection itself and on whether the decision was due to the signal candidate, other particles produced in the pp collision, or a combination of both.
In the analysis, the B 0 ! D 0 D + s ⇡ and B + ! D D + s ⇡ + candidates are formed using charged kaon and pion candidates, in which the D 0 candidates is reconstructed through the D 0 ! K + ⇡ and D 0 ! K + ⇡ ⇡ ⇡ + decays, D + through the D + ! K ⇡ + ⇡ + decays, and D + s through the D + s ! K + K ⇡ + decays. The invariant mass of D 0 ⇡ is required to be larger than 2.05 GeV in order to veto the contribution from the B 0 ! D ⇤ D + s decay, which is not the focus of this analysis. To consider potential variations of signal e ciencies and background distributions, selections are designed and optimized separately for the six datasets, namely Run 1 and Run 2 datasets with three reconstruction channels,
, and B + ! D D + s ⇡ + decay. Loose requirements on the p, p T , PID response, and minimum 2 IP of the charged tracks are first applied to improve track quality and remove tracks originating directly from the pp collision. Here 2 IP is defined as the di↵erence in the vertex-fit 2 of a given PV reconstructed with and without the considered track. To separate B and D (s) candidates from random combinations of tracks directly produced in pp collisions, loose requirements on the invariant mass and on the 2 of the flight distance of D candidates with respect to the associated PV are imposed, where the associated PV is defined as the PV yielding the smallest 2 IP for the considered B candidate. The quantity cos ✓ dir for a B candidate, where ✓ dir is the angle between the momentum direction and the vector from the associated PV to the vertex of the candidate, is required to be close to unity.
A boosted decision tree (BDT) classifier [47,48] implemented in the TMVA toolkit [49] is used to further suppress the combinatorial background. The BDT classifier depends on the 2 IP of the B, D, D s candidates and the final tracks, the 2 of the B decay vertex and its cos ✓ dir , the PID response of all the final tracks, and the signed significance of the separation of D (s) and B vertices parallel to the beam pipe (s
z FD ). For the BDT classifier training, the signal samples are the simulated signal candidates, and the backgrounds are B sideband candidates in data with B invariant mass within [5500, 6950] MeV. The requirement on the BDT response is determined by maximising S 2 / (S + B) 3 2 , where the S and B are the expected signal and background yields in the signal mass window |M (B) m(B)| < 20 MeV, which is 2.5 to 3 times wider than the mass resolution for the di↵erent channels, where the M (B) and m(B) are the reconstructed and known masses of the corresponding B meson [7], respectively.
Two kinds of misidentified (misID) backgrounds are vetoed through additional requirements. A background to the B + ! D D + s ⇡ + signal occurs at the mass threshold of the D + s ⇡ + spectrum due to e + ⇡ + misidentification. It is removed by a stringent requirement on the PID response of the companion pion. For
misID candidates are vetoed by tightening PID requirements of the companion ⇡ for candidates with invariant mass within ±30 MeV of the B 0 s or ⇤ 0 b mass [7], after replacing the mass hypothesis of the companion ⇡ to K or p.
A track-swapped background is found to peak in the B signal region where the
is swapped with a ⇡ in the D 0 decays. As the momentum of the ⇡ in D ⇤ decays is almost zero in the D ⇤ rest frame, this type of background is removed by requiring M (K + ⇡ ⇡ ⇡ + ) M (K + ⇡ ⇡ + ) > 160 MeV in the D 0 decays, to veto those ⇡ with negligible momentum.
Backgrounds with the same final-state tracks but with one or even zero intermediate charmed mesons, called non-double-charm (NDC) backgrounds, are suppressed by the requirements on s D (s) z FD . After the selection, the invariant mass distributions of all possible two-body and three-body combinations of the final charged tracks in D (s) sideband samples are further checked, and all the visible narrow structures, namely K ⇤ (892) 0 , , D 0 , D + and D + s particles, are vetoed, to further suppress the NDC backgrounds. After applying the full o✏ine selection, the D and D s masses are required to lie within ±15 MeV around their known mass [7], corresponding to two to three times the detector resolution. Roughly 1% to 2% of events contain more than one B candidate; one is retained at random in these cases.
To improve the resolution of B candidate invariant mass distributions, a fit based on the Kalman filter method [50], which contains topological and kinematic information of the decay chain, is applied. By updating the four-momentum of all the final-state tracks, the invariant mass of D and D s candidates is constrained to the their known masses [7], and the updated invariant mass distributions of B candidates are used to determine the signal and background yields. The four-momenta of charged tracks from another fit, which additionally constrains B candidate mass to the known B mass, are used to calculate kinematic variables used in the amplitude fit.
Figure 2 shows the invariant mass distributions of the B 0 and B + candidates in each dataset after the application of the selection requirements. An extended unbinned maximumlikelihood fit to the data in the mass range [5230,5630] MeV is used to extract the signal and background yields, which are used later in the amplitude fit.
The total PDF comprises a signal PDF and an exponential function to describe the distribution of combinatorial background. The signal PDF is a double-sided Crystal Ball (DSCB) function [51] which consists of a Gaussian kernel and independent tail parameters on both the left and right sides to model e↵ects such as the detector resolution and final-state radiation. In the fit, the signal and background yields, the parameters of the exponential function, and those of the Gaussian kernel in the DSCB function are allowed to vary independently in each dataset, while the tail parameters are fixed to the values obtained from a fit to the corresponding simulated sample.
In the D 0 D + s ⇡ mode, an additional DSCB function is included in the fit model to describe the singly Cabibbo-suppressed B 0 s ! D 0 D + s ⇡ contribution, whose yield is determined from data. The width of the Gaussian kernel is shared with the B 0 signal and the mean value µ B 0 s is defined as
, where the µ B 0 is the mean value of the Gaussian kernel of the B 0 signal PDF, and m(B 0 ) and m(B 0 s ) are the known masses of the B 0 and B 0 s mesons [7]. The tail parameters are fixed to the same values as the B 0 signal.
The results of the fit to the data samples are shown in Fig. 2. The values of the fitted parameters are listed in Table 1, and Table 2 summarises the signal yields, the number of candidates, and the purity inside the signal mass window used for the amplitude analysis.
The amplitude formalism and fit method of three-body B ! abc decays, where abc denotes any sequence of D, D + s and ⇡ states, is established in this Section.
The amplitude of the three-body B ! DD s ⇡ decays is constructed following the isobar formalism [52][53][54], which is a coherent sum of quasi two-body amplitudes, either resonant or nonresonant,
where c i is a complex parameter for the i-th contribution that is determined from data, x denotes variables calculated from the four-momenta of the final-state particles and ⇥ i is a set of parameters used to describe the i-th lineshape. The amplitude of the i-th quasi two-body decay to a and b (a, b represent any pair of D, D s , ⇡ mesons) is where T (✓ ab ) describes the angular distribution which depends on the spin J of the intermediate resonant state R(ab). The helicity angle, ✓ ab , is defined as the angle between the R(ab) momentum direction in the B rest frame, and the momentum direction of a as determined in the R(ab) rest frame. The definitions of T (✓ ab ) up to J = 4 are
The function f (M ab ) is the lineshape of the R(ab) resonance where M ab is the invariant mass of the pair. The complex relativistic Breit-Wigner (RBW) function is used as the default lineshape,
where p is the momentum of particle a in the rest frame of the resonance R(ab), and q denotes the momentum of R(ab) in the B rest frame. The mass-dependent running width is
where m 0 and 0 are the mass and width of the resonance, respectively. The quantities p 0 and q 0 are these momenta evaluated when M = m 0 . The orbital angular momentum between R(ab) and c is denoted by L 1 , while L 2 refers to the orbital angular momentum between particles a and b. Conservation of angular momentum implies that
The Blatt-Weisskopf form-factor [55] F (M, L) is parameterized as
where z(M ) = pd, z 0 = p 0 d, and d stands for the radial parameter, which is taken to be 3.0 GeV 1 by default for all resonances. Nonresonant (NR) contributions are parameterized using an exponential function,
where ↵ is the slope parameter that is allowed to vary in the fit, and m 2 min (D⇡) = 4 GeV 2 is an approximation of the M 2 (D⇡) lower threshold.
The RBW functions do not provide an adequate description of overlapping resonant states. Furthermore, the latest experimental [56] and theoretical [57] studies show that a simple BW lineshape is not su cient for the D ⇤ 0 (2300) resonance. A quasi-modelindependent (qMI) parameterization [56] is used for the D⇡ S-wave, where the M (D⇡) range is divided into k slices, and the line shape is replaced by a set of complex coe cients assigned to each slice, each free to vary in the fit. The real and imaginary parts are independently interpolated using cubic splines. Details are given in Sec. 7.
The normalized PDF for the signal is expressed as
The normalization factor, I sig (⇥), which is obtained by integrating over the phase space using simulated samples after full selection and thus including ✏(x) implicitly, can be expressed as
The signal e ciency, ✏(x), is obtained as described in Sec. 6. The w j , also described in Sec. 6, are applied to the simulated events to correct for discrepancies between simulated samples and data, where the subscript j runs over all the events of the sample. The resolutions of the Dalitz plot coordinates are much smaller than the widths of the narrowest structures and thus related e↵ects can be neglected. The total PDF is given by
where f sig and f bkg are the fractions of signal and background contributions, determined from the fit to the m(DD s ⇡) invariant mass distributions, and P norm bkg (x) is the normalized background PDF described in Sec. 6.
An amplitude fit is performed, minimising the unbinned negative log-likelihood
The signal e ciency map and background map are obtained separately for di↵erent samples. For the B 0 ! D 0 D + s ⇡ decay, the di↵erent LHC Run (1 or 2) and reconstruction channels (B 0 ! D 0
) are fitted simultaneously. For the B + ! D D + s ⇡ + decay, the Run 1 and Run 2 datasets are simultaneously fitted. Where isospin symmetry is imposed, the fit is performed simultaneously on all datasets.
The fit fraction F i for a given contribution i is calculated from the fitted parameters, ⇥ 0 , and is defined as
The fit fractions do not necessarily add up to 1 due to interference e↵ects between components. The interference term between any pair of components is defined as
and thus we have
6 Signal e ciency and background models
The amplitude analysis is only sensitive to signal e ciency variations across the Dalitz plot, not the absolute e ciency. These are extracted for each dataset as a function of position in the square Dalitz plot (SDP), whose coordinates are defined by
Here m(D⇡) is the invariant mass of the D⇡ combination, m min D⇡ and m max D⇡ are the kinematic limits of m(D⇡) in B ! DD s ⇡ decays, and ✓(D⇡) is the D⇡ helicity angle.
The e ciency maps across the SDP are evaluated from the simulated samples by kernel density estimation [44] and are shown in Fig. 3. The tracking e ciency and the e ciency of the trigger requirements have been corrected using control samples in data.
The background distributions over the phase space are estimated using candidates in the B mass sidebands between [5450, 6000] MeV for the B 0 ! D 0 D + s ⇡ decays and [5400, 6000] MeV for the B + ! D D + s ⇡ + decays. The requirement on the BDT response is relaxed to increase the number of events in these regions as no significant change in the shapes of the distributions of sideband samples is observed. A Gaussian process extrapolation method [58] is applied to extrapolate the background Dalitz plot distribution into the B-meson signal region, in order to account for the correlations between the Dalitz variables and the invariant mass of the B candidates. The SDP distributions of the background shape for each dataset are shown in Fig. 4. The broad structures in the SDP distributions are related to the accumulation of combinatorial background with low pion momentum.
Conventional resonances are only expected to decay to D 0 ⇡ and D ⇡ + final states in the B 0 ! D 0 D + s ⇡ and the B + ! D D + s ⇡ + decays, respectively. The straightforward way to perform the amplitude analysis is by including all the D⇡ resonances with natural spin-parity, as listed in Table 3, in the fit model. It is called a model-dependent (MD) description. The result with the MD description is shown in Sec. 7.7. The description of the D⇡ S-wave distributions is improved by introducing a 0 + D⇡ qMI description, which accounts for both the broad spin-0 D⇡ states and nonresonant components. The result with the qMI description is considered as the default, which is further discussed in this section.
The basic fit model is defined by considering all known D ⇤⇤ states with natural spinparity [7], as listed in Table 3, except for the broad D ⇤ 0 (2300) state. Their masses and widths are fixed to their default values. As described in Sec. 5.1, a qMI description of the D⇡ S-wave is used [56], with 11 spline points. 3 The first and last points are outside of the invariant mass range, and their amplitudes are fixed to zero. The other points are each assigned a complex coe cient that varies freely in the fit. Moreover, as the D ⇤ (2007) 0 mass is lower than the D ⇡ + mass threshold, the q 0 value in Eq. 6 would be imaginary. The q 0 value in this case is taken as the value calculated from D ⇤ (2007) 0 ! D 0 ⇡ 0 rather than D ⇤ (2007) 0 ! D ⇡ + in the default model. The D ⇤ J (3000) state was first observed in the D ⇤ ⇡ + decay mode, and its spin has not been determined yet [8]. A similar structure has been seen in B ! D + ⇡ ⇡ decays [56], with J P = 2 + . In this analysis di↵erent J P hypotheses for the D ⇤ J (3000) state are tested, either 1 , 2 + , 3 , or 4 + . In each case its mass and width are fixed to the corresponding default values [7]. The test results favor J P = 4 + , which is used as the default.
The fit results, where only D⇡ resonances are included and the two decays are considered independently, are given in Figs. 5 and6 for
s ⇡ decays are combined when plotting here and subsequently. A peaking structure at about 2.9 GeV is visible in the M (D s ⇡) distribution of each decay, and is not well described by the included contributions. Furthermore, the addition of further D ⇤ states up to J P = 4 + does not resolve the discrepancy. The normalized residuals of the fits are shown in Fig. 7. In the plot, the pull value in bin i is defined as (N i sig N i fit )/ q N i sig , where the N i sig and N i fit are the number of signal candidates and number of expected candidates from the fit result. The 2 /ndf , where ndf is the number of degrees of freedom, is 78.2/35 for B 0 ! D 0 D + s ⇡
To improve the description of the M (D + s ⇡) distributions for the two decays, an additional D + s ⇡ state is added to each decay, whose mass and width are free parameters, and di↵erent J P assignments are tested. No relationship is assumed for the two D + s ⇡ states. Both states with J P = 0 + give the best description of the data, while the D + s ⇡ states with the other spin-parity are disfavored compared to the 0 + hypothesis (see Sec. s ⇡ resonances, which are named as T a cs0 (2900) 0 and T a cs0 (2900) ++ following the convention in Ref. [59], are evident.
In the M (D + s ⇡ ) and M (D + s ⇡ + ) distributions, both the peaks near 2.9 GeV and the dips near 3.0 GeV are better described by the presence of the new states and their interference with the existing D ⇤ states. The masses and widths of the T a cs0 (2900) 0 and T a cs0 (2900) ++ states are listed in Table 4. Fit fractions are given in Tables 5 and6 for B 0 ! D 0 D + s ⇡ and B + ! D D + s ⇡ + decays, respectively. These results include the systematic uncertainties and corrections of fit bias, which are described in Sec. 9. The amplitudes and phases of the complex coe cients of the resonant contributions, relative to those of D ⇤ 2 (2460), are also displayed in Tables 5 and6. The two-dimensional pull plots are given in Fig. 9. The 2 /ndf is 43.2/31 and 63.0/31 for B 0 ! D 0 D + s ⇡ and B + ! D D + s ⇡ + decays, respectively. The distributions of Legendre polynomial weighted moments, together with the fit results with and without T a cs0 (2900) 0 and T a cs0 (2900) ++ states, are shown in Appendix A; these also suggest the existence of the new exotic states. The above model with a new 0 + T a cs0 (2900) is set as the default fit model.
Mass ( GeV) Width ( GeV) T a cs0 (2900) 0 2.879 ± 0.017 ± 0.018 0.153 ± 0.028 ± 0.020 T a cs0 (2900) ++ 2.935 ± 0.021 ± 0.013 0.143 ± 0.038 ± 0.025
Numerous additional resonances are tested to verify the stability of the default fit result and to better understand the system. First, the exotic D + s ⇡ states with other spin-parity hypotheses are tested. The LL values, with definition LL = NLL Default NLL Other , are summarized in Table 7. It is clear that at least one D + s ⇡ exotic state is needed to improve the value of NLL of the fit to each decay channel. The default model has the best fit quality, while the model with a spin-1 D + s ⇡ state also provides reasonably good description. The data are tested against the spin 0 and spin 1 hypotheses for the T a cs0 (2900) 0 and T a cs0 (2900) ++ states. The results are shown in Sec. 7.6 and 0 + is favored.
Secondly, the existence of extra D⇡, D + s ⇡ and DD + s states with natural spin-parity up to 3 is explored when including the 0 + T a cs0 (2900) 0 and T a cs0 (2900) ++ states. The Table 7: Tested fit models and the corresponding LL value. masses and widths of the additional resonances are allowed to vary freely in the fit. When considering an additional spin-0 or spin-1 D + s ⇡ state for each decay, the mass of the new state converges near the D + s ⇡ mass threshold. The resulting changes in the value of NLL are insignificant, corresponding to statistical significance of less than 2 , and these states are not included in the default model.
The D⇡ states with natural spin-parities have been well investigated by the amplitude analyses of B 0 ! D 0 ⇡ + ⇡ [60], B + ! D ⇡ + ⇡ + [56] and other topologically similar B-meson decays in LHCb [61][62][63] with large yields. No extra D⇡ state is expected to be observed in this analysis, which is consistent with the D⇡ results in Table 7. Additional DD + s exotic states with natural spin-parities are also found to be disfavored, which is consistent with the previous results [28][29][30], where only 1 + Z cs states are observed.
With the T a cs0 (2900) 0 and T a cs0 (2900) ++ states in place, the excited D states are investigated. Some tension in the measured mass of the D ⇤ 2 (2460) state between inclusive results [8] and those obtained in amplitude fits [56,[60][61][62] 2 (2460) 0 state. These results are in better agreement with earlier measurements in amplitude analyses. However they are also consistent with those obtained from inclusive results within 3 when only considering statistical uncertainties [7].
The charged isospin partners of the D ⇤ 1 (2600) 0 and D ⇤ J (3000) 0 states have not yet been observed, their significances are estimated in the B 0 ! D 0 D + s ⇡ decays based on the default fit model by fixing their masses and widths to the known values [7] and assuming the D ⇤ J (3000) 0 spin-parity to be 4 + . The statistical significance of the D ⇤ 1 (2600) , and D ⇤ J (3000) resonances is estimated to be 4.8 and 2.2 , respectively. When the mass and width of the D ⇤ 1 (2600) resonance are allowed to vary in the fit, they are determined to be m 0 = (2640 ± 51) MeV, 0 = (122 ± 35) MeV which are consistent with the default values. The masses and widths of the D ⇤ J (3000) state cannot be determined due to the limited sample size. While the significance of these states is small they are still included in the default fit for a conservative evaluation of the exotic contributions. The nature of the D ⇤ s0 (2317) + state is still in debate. Some theoretical models interpret the D ⇤ s0 (2317) + state as an isoscalar [cqsq] tetraquark state, and suggest searching for the isotriplet partners in the D + s ⇡ + and D + s ⇡ final states [64][65][66][67][68][69]. The B 0 ! D 0 D + s ⇡ and B + ! D D + s ⇡ + decays are ideal for such studies. However, in Sec. 7.3, additional D + s ⇡ exotic states with freely varying mass and width values, under di↵erent spin-parity hypotheses, are found to be insignificant. By assuming that the masses of the neutral and doubly charged partners are the same as that of the D ⇤ s0 (2317) + state, the upper limit on fit fractions with three di↵erent scenarios is evaluated. The first scenario is that the natural width of the new D ⇤ s0 (2317) state is ignored. A Gaussian function is used to describe the lineshape of the new state, of which the width represents the detector resolution. The second is that a Breit-Wigner function is added to model the J P = 0 + , D + s ⇡ state, of which the width is set to 3.8 MeV, the current upper limit on the D ⇤ s0 (2317) + width [7]. No resolution e↵ect is considered. The third is that an additional 0 + D + s ⇡ Breit-Wigner function, the width of which is set to be the same as the T a cs0 (2900) 0 and T a cs0 (2900) ++ states, is included in the fit model. The upper limits on the fit fractions of neutral and doubly charged D ⇤ s0 (2317) with di↵erent hypotheses at 90% confidence level (C.L.), which are all less than 1%, are summarized in Table 8.
Pseudoexperiments are carried out to determine the significance of the T a cs0 (2900) 0 and T a cs0 (2900) ++ states, accounting for the look-elsewhere e↵ect [30,70]. The pseudoexperiments use data generated according to the fit results without the new D + s ⇡ resonances, and the yield of each generated sample follows a Poisson distribution whose mean is the yield in the corresponding dataset. The pseudodatasets are fitted with the model with and without the D + s ⇡ contributions. The di↵erence in the value of 2 LL between these two fits is obtained, and fitted with a 2 PDF. The distributions and fit results are shown in Fig. 10 for the two decays. The numbers of degrees of freedom after considering the look-elsewhere e↵ect are found to be 7.39 ± 0.17 and 6.93 ± 0.17 for the T a cs0 (2900) 0 and T a cs0 (2900) ++ states, respectively, with the corresponding significances estimated to be 7.3 and 5.3 . After accounting for the systematic uncertainties discussed in Sec. 9, these are reduced to 6.6 and 4.8 .
The spin-parity values for the T a cs0 (2900) 0 and T a cs0 (2900) ++ states are also determined using pseudoexperiments. For each decay, 500 pseudoexperiments are generated based on the fit results with the 0 + D + s ⇡ state included, while another 500 pseudoexperiments are generated according to the fit results using a model assuming the 1 spin hypothesis. Each pseudodataset is fitted in the same way as for data. The 2 LL between the two fits is
The blue solid and black dashed histograms are the distributions of the 2 LL for the pseudoexperiments generated based on the fit results with 0 + or 1 D + s ⇡ exotic state, respectively. The purple vertical line shows the 2 LL value for the data fitted with the new D + s ⇡ exotic state under the J P = 0 + and J P = 1 hypotheses. The red curve is the result of a fit to the black dashed histogram with a Gaussian function. calculated, and the distributions of 2 LL for the two sets of pseudoexperiments are shown in Fig. 11. The 2 LL distribution of the pseudoexperiments with the 1 D + s ⇡ hypothesis is fitted with a Gaussian function, and the significance of the data disfavoring the 1 hypothesis is evaluated to be 4.3 for the T a cs0 (2900) 0 state and 4.2 for the T a cs0 (2900) ++ state when no isospin relationship is imposed on the D⇡ or D s ⇡ components.
The Argand diagrams [7] of the T a cs0 (2900) 0 and T a cs0 (2900) ++ states are shown in Fig. 12. The Breit-Wigner function follows a counter-clockwise circular path on the complex plane while contributions which are not genuine resonances are expected to have a di↵erent shape. Seven spline points on M D + s ⇡ near the measured mass of the T a cs0 (2900) 0 and T a cs0 (2900) ++ states (m ± 1.5 ) are used to model these regions instead of Breit-Wigner functions. The complex parameters of all points are allowed to vary in the fit. The fitted parameters of the spline points, together with the D + s ⇡ lineshape are shown in Fig. 12. The spline lineshape shows similar behavior as the Breit-Wigner distribution, which confirms the resonant character of the two new exotic states.
Instead of modeling the D⇡ S-wave with the qMI description, the amplitude fit with a fully model-dependent (MD) description is carried out. In the MD description, the 0 + qMI model is replaced by the RBW of the D ⇤ 0 (2300) component together with a 0 + nonresonant component, while the parameters of all the other resonances are set to be the same. The obtained parameters of the T a cs0 (2900) 0 and T a cs0 (2900) ++ states are summarized in Table 9 and are consistent with those determined using the qMI model, within statistical uncertainties. Figures 13 and14 show the fit results of the MD description of the B 0 ! D 0 D + s ⇡ and B + ! D D + s ⇡ + decays, respectively. The 2 /ndf is 79.1/44 for the B 0 ! D 0 D + s ⇡ decay and 114.7/44 for the B + ! D D + s ⇡ + decay. As a comparison, the Argand diagrams of the MD and qMI 0 + D⇡ components of the two Table 9: Masses, widths and fit fractions of the T a cs0 (2900) 0 and T a cs0 (2900) ++ states obtained from the MD fit.
LL Mass (GeV) Width (GeV) Fraction (%) T a cs0 (2900) 0 70.1 2.871 ± 0.012 0.135 ± 0.025 3.0 ± 0.5 T a cs0 (2900) ++ 33.2 2.922 ± 0.014 0.161 ± 0.033 2.0 ± 0.
In the default fit model, all known D ⇤⇤ states with natural spin-parity [7] are included. However, their fit fractions, except those of the D ⇤ 2 (2460) and D ⇤ 1 (2600) states, are consistent with 0, as shown in Tables 5 and6. Moreover, the parameters of the qMI 0 + D⇡ spline points in the higher D⇡ mass region have large uncertainties due to the smaller sample size. To improve the precision and stability of the fit results, a simultaneous fit of the B 0 ! D 0 D + s ⇡ and B + ! D D + s ⇡ + decays is performed, as the two decays are related by isospin symmetry. In the simultaneous fit, all complex parameters of the D ⇤⇤ states are shared, except for the D ⇤ (2007) 0 and D ⇤ (2010) states allowing for small isospin symmetry breaking e↵ects near the D⇡ mass threshold. The simultaneous fit results with only D⇡ states are shown in Ref. [31]. The description of the mass spectra is consistent with the separate fit as shown in Sec. 7.1, which supports the feasibility to
Mass ( GeV) Width ( GeV) T a cs0 (2900) 0 2.892 ± 0.014 ± 0.015 0.119 ± 0.026 ± 0.013 T a cs0 (2900) ++ 2.921 ± 0.017 ± 0.020 0.137 ± 0.032 ± 0.017 perform simultaneous fit for the two decays.
The fit results are shown in Fig. 16. As separate fit, neutral and doubly charged D + s ⇡ states are needed to describe the data well, where their parameters are set to be di↵erent. The masses, widths and fit fractions of the T a cs0 (2900) 0 and T a cs0 (2900) ++ states after considering the systematic uncertainties and possible fit bias, which are summarized in Tables 10 and11, show good agreement with the separate fit result, with significant improvement on the relative statistical uncertainties. The total 2 /ndf is evaluated to be 140.3/89.
The ndf of the T a cs0 (2900) 0 and T a cs0 (2900) ++ states in the simultaneous D⇡ fit model are evaluated to be 7.29 ± 0.18 and 8.57 ± 0.17 in the same way as described in Sec. 7.6, with the corresponding significances estimated to be 9.0 and 7.4 , respectively. After accounting for the systematic uncertainties discussed in Sec. 9, these are reduced to 8.0 and 6.5 . The constraints on the D⇡ contributions using isospin symmetry lead to higher significance, as expected.
To estimate the isospin-breaking e↵ects between the T a cs0 (2900) 0 and T a cs0 (2900) ++ states, the mass di↵erence, M = M (T a cs0 (2900) ++ ) M (T a cs0 (2900) 0 ), and width di↵erence, = (T a cs0 (2900) ++ ) (T a cs0 (2900) 0 ), are evaluated to be 28 ± 20 ± 12 MeV and 15 ± 39 ± 16 MeV, respectively, where the first and second uncertainties are statistical and systematic. The statistical uncertainties in M and are evaluated using pseudoexperiments to account for the correlations, and some of the systematic uncertainties are canceled. The masses and widths of the two exotic states are consistent with each other within 1 , which confirms that they are related by isospin symmetry.
The full simultaneous fit, where the parameters of the T a cs0 (2900) 0 and T a cs0 (2900) ++ states are shared in the fit, is also performed, and described in a separate Letter [31].
The sources of systematic uncertainty fall into two categories: experimental and those related to the amplitude model. In the first category there are e↵ects related to the fixed signal yields of B candidates, the models of the background distributions, and the signal e ciency computation. Those arising from the amplitude model are mainly due to the fixed parameters of the model. The total systematic uncertainty is found by summing these in quadrature. The signal yields in the amplitude analysis are taken from the results of the fits to the invariant mass distributions of B candidates. To determine the systematic uncertainty, the signal yield of each dataset is varied according to a Gaussian distribution whose width corresponds to a signal-yield uncertainty that includes uncertainties due to the modelling of the invariant mass distribution. The amplitude fit is repeated with the new signal yields and the RMS value of each fit parameter is taken as the systematic uncertainty.
Backgrounds are modelled using a Gaussian process extrapolation method [58] according to sideband distributions. To evaluate the associated systematic uncertainty, the background model is replaced by the result of a kernel density estimation [44] applied to the Dalitz-plot distributions of the sideband samples. The deviations of the fit parameters from the default result are taken as the associated systematic uncertainties.
Knowledge of the signal e ciency variation over the Dalitz plot is limited by four e↵ects: uncertainty in the PID response, trigger e ciency calibration uncertainties, signal e ciency determination, and simulation sample size. The systematic uncertainty due to the PID calibration is estimated by regenerating the PID responses with a perturbed kernel density estimation [44], extracting the e ciencies, and repeating all the fitting procedures. For the trigger e ciency calibration e↵ect, a conservative systematic uncertainty is determined by repeating all the fitting procedures with the signal e ciency maps without any trigger e ciency correction. The systematic uncertainty related to the signal e ciency determination is estimated by performing the amplitude analysis with the e ciency maps obtained using an alternative kernel density estimation. The deviations of the fit parameters in each fit are taken to be the systematic uncertainties.
The systematic uncertainty due to the limited size of the simulated samples is determined by generating 200 samples following the bootstrap method [71], where the simulated candidates after all the selection criteria are allowed to be picked multiple times. The e ciencies are then extracted from each bootstrapped sample and applied in the amplitude analysis. The RMS of each fit parameter in these fits is taken as the systematic uncertainty.
The fixed parameters in the amplitude analysis include the radius d in the Blatt-Weisskopf form factor, and the parameters of the D⇡ lineshapes. The systematic uncertainty associated to d is evaluated by setting d to 1.5 GeV 1 and 4.5 GeV 1 . The largest di↵erence compared to the default results is assigned as the systematic uncertainty. The fixed D⇡ parameters consist of the masses and widths of all the D ⇤ states in Table 3, the choice of the qMI spline points, the constant q 0 in the D ⇤ (2007) 0 state lineshape, and the spin hypothesis of D J (3000). The amplitude fits are performed several times, with one or more parameters changed. The fixed masses and widths of all the D ⇤ states are allowed to vary one at a time in the fit but are constrained by a Gaussian function within uncertainties in their default values [7]. The positions of the qMI spline points are chosen empirically in the default fit, and shifted by ±10 MeV in the amplitude analysis. The systematic uncertainty from the number of the qMI spline points is also explored, by adding a new point at M (D⇡) = 2.35 GeV to try to improve the model description near the D ⇤ 0 (2300) mass region, or removing the point at M (D⇡) = 2.6 GeV as there is no 0 + D⇡ state observed in this region. The q 0 of the D ⇤ (2007) 0 state is taken as the q 0 of the D ⇤ (2007) 0 ! D 0 ⇡ 0 decay in the default fit, and replaced by a value calculated from the D ⇤ (2007) 0 e↵ective mass [56]. The spin of the D J (3000) state is found to be 4 + in the default fit, and altered to 2 + , the result measured in Ref. [56]. The deviations of the results between each fit and the default fit, are summed in quadrature, and taken as the systematic uncertainty.
Possible fit biases are investigated using pseudoexperiments, and used to correct the results. For each dataset, 500 samples are generated according to the default fit results, where the yield of each is sampled from a Poisson distribution whose mean value is the number of B candidates in the corresponding dataset. The fit parameters of the pseudodata are then extracted using the default fit model. The residual distributions (µ pseudo µ default ) and pull distributions (µ pseudo µ default )/ pseudo are extracted from the pseudoexperiments and default fit result, and used to correct the fit results. Here the µ pseudo and pseudo are the mean values and uncertainties in the pseudoexperiments. Both the residual distributions and pull distributions of each parameter are fitted with a Gaussian function. The mean value of the residual distribution from the Gaussian fit is used to correct the mean value of the parameter, while the width of the pull distribution is used to scale the statistical uncertainty. The results are summarized in Tables 4,5, and 6. The systematic uncertainties of the simultaneous D⇡ fit are also evaluated in the same way, and summarized in Tables 10 and11.
Amplitude analyses of B 0 ! D 0 D + s ⇡ and B + ! D D + s ⇡ + decays are performed for the first time, using LHCb pp collision data taken at center-of-mass energies of 7, 8 and 13 TeV, corresponding to a total integrated luminosity of 9 fb 1 . In total, signal yields of 4009 ± 70 and 3750 ± 64 candidates are obtained from the B 0 ! D 0 D + s ⇡ and B + ! D D + s ⇡ + decays, respectively. When all known D⇡ resonances with spin-parities of 1 , 2 + , 3 and 4 + [7] are included, along with a qMI spline model to describe the 0 + D⇡ distributions, the results show that the D + s ⇡ invariant-mass distributions are not well described. To improve the model description, a 0 + D + s ⇡ resonance is added to each decay mode. The masses and widths of the two resonances are determined to be T a cs0 (2900) 0 : M = (2.879 ± 0.017 ± 0.018) GeV, = (0.153 ± 0.028 ± 0.020) GeV, T a cs0 (2900) ++ : M = (2.935 ± 0.021 ± 0.013) GeV, = (0.143 ± 0.038 ± 0.025) GeV, where the first uncertainty is statistical and the second systematic. The significances, accounting for the look-elsewhere e↵ect and systematic uncertainties, of the exotic T a cs0 (2900) 0 and T a cs0 (2900) ++ states are 6.6 and 4.8 , respectively. A simultaneous D⇡ amplitude fit assuming isospin symmetry in the B 0 ! D 0 D + s ⇡ and B + ! D D + s ⇡ + decays is also performed to provide better control on the contributions from D⇡ resonances, especially the 0 + D⇡ spline model, and to improve the precision of the measured parameters of exotic states. The masses and widths of the two resonances in the simultaneous D⇡ fit are measured to be T a cs0 (2900) 0 : M = (2.892 ± 0.014 ± 0.015) GeV, = (0.119 ± 0.026 ± 0.013) GeV, T a cs0 (2900) ++ : M = (2.921 ± 0.017 ± 0.020) GeV, = (0.137 ± 0.032 ± 0.017) GeV, with the significances evaluated to be 8.0 and 6.5 for the T a cs0 (2900) 0 and T a cs0 (2900) ++ states, including systematic uncertainties. The mass and width di↵erences between T a cs0 (2900) ++ and T a cs0 (2900) 0 are evaluated to be M = (28 ± 20 ± 12) MeV, = (15 ± 39 ± 16) MeV, based on simultaneous D⇡ amplitude fit, and consistent with zero. A simultaneous fit with the parameters of the D + s ⇡ exotic states shared is also performed, and described in a separate Letter [31]. All the results of the di↵erent fit scenarios show good agreement. This is the first observation of an isospin triplet of manifestly exotic mesons with four di↵erent quark flavors. The masses and widths of the two states are consistent with the X 0 (2900) and X 1 (2900) states [3,4], but have an opposite strangeness number. No hint of a DD + s structure is observed in the analysis. With the significantly larger data samples that will be collected by the upgraded LHCb detector in the coming years, the nature of the isospin triplet of exotic mesons, and the existence of the possible D + s ⇡ exotic states with J P = 1 in the same region, will be further explored.
Moments analysis is useful to suggest possible resonant structures in the decay. The formulation and the results are provided in this section.
The Legendre polynomial of a certain order k, as expressed in Eq. 17, is used to weight the data,
For example, when focusing on the resonant structures of the D⇡ channel, the variable x in this case would be the helicity variable in that decay chain, namely cos ✓ D⇡ D . Therefore, the total amplitude modified by the order-k Legendre polynomial can be expressed as
where w i is the original weight for the data point i.
By analytical calculation, the relationship between the Legendre-weighted total amplitude hY k i and combination of di↵erent orders of partial waves [56] can be bridged. Considering the existence of partial waves with the first J orders of orbital momentum, only hY k i of k up to 2J are nonzero. The weighted distributions can be visualized on the M 2 (D⇡) axis, which can be helpful to distinguish the structures from di↵erent ordered partial waves. The moments on the two axes, M 2 (D⇡), M 2 (D s ⇡), are shown up to the eighth order. The results, which are shown in Figs. 17 o Università di Padova, Padova, Italy p Università di Perugia, Perugia, Italy q Scuola Normale Superiore, Pisa, Italy r Università di Pisa, Pisa, Italy s Università della Basilicata, Potenza, Italy t Università di Roma Tor Vergata, Roma, Italy u Università di Siena, Siena, Italy v Università di Urbino, Urbino, Italy w Universidad de Alcalá, Alcalá de Henares , Spain † Deceased
Charge conjugation is implied throughout this paper.
Natural units with ~= c = 1 are used throughout this paper.
[1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.9, 3.4] GeV