Model-Independent Study of Structure in B + → D + D − K + Decays

The only anticipated resonant contributions to B þ → D þ D − K þ decays are charmonium states in the D þ D − channel. A model-independent analysis, using LHCb proton-proton collision data taken at center-of-mass energies of ﬃﬃﬃ s p ¼ 7 , 8, and 13 TeV, corresponding to a total integrated luminosity of 9 fb − 1 , is carried out to test this hypothesis. The description of the data assuming that resonances only manifest in decays to the D þ D − pair is shown to be incomplete. This constitutes evidence for a new contribution to the decay, potentially one or more new charm-strange resonances in the D − K þ channel with masses around 2 . 9 GeV =c 2 .

The B þ → D ðÃÞþ D ðÃÞ− K þ family of decays offers unique opportunities to study charmonium states. The constrained environment of B-meson decays allows the masses, widths, and quantum numbers of such states to be determined using amplitude analysis techniques, with low backgrounds from other processes. In particular, resonances in the D ðÃÞ− K þ or D ðÃÞþ K þ channels would be manifestly exotic, having minimal quark contentcdus or cdus, respectively. While many exotic hadrons containing cc or bb quarks have recently been observed [1][2][3], there is to date no significant evidence of the existence of exotic hadrons with open flavor, i.e., with nonzero strangeness, charm, or beauty quantum numbers. Studies of B þ → D ðÃÞþ D ðÃÞ− K þ decays are therefore expected to help resolve open questions regarding charmonium spectroscopy [4,5]. In addition, measurements of these processes have been proposed as a method to aid characterization of the cc contribution in B þ → K þ μ þ μ − decays [6,7].
The branching fractions of B þ → D ðÃÞþ D ðÃÞ− K þ decays have been measured [8,9], but no prior analyses of their resonant structure exist [10]. Recent studies have shown that extremely pure samples of these decays can be obtained using LHCb data [9] with yields much larger than those available at previous experiments.
A model-dependent study of the resonant structure in B þ → D þ D − K þ decays [11], carried out in parallel to this work, has revealed structure in the D − K þ invariantmass spectrum that cannot be described by reflections of charmonium resonances. This highly surprising observation, along with the limited current knowledge of the charmonium spectrum in this mass range, particularly among spin-0 and spin-2 states, motivates the study of this decay using a model-independent approach as presented in this Letter. This method is particularly useful when applied to three-body decays where resonances are only expected to form between one pair of the final-state particles, such that the decay kinematics are described through one mass and one angular variable. Unexpected exotic, contributions to the decay process manifest as high-order moments in the distribution of the angular variable, as has been demonstrated by the use of the method to identify exotic resonances contributing to B 0 → ψð2SÞK þ π − [12], Λ 0 b → J=ψpK − [13], and B 0 → J=ψK þ π − [14] decays.
The model-independent analysis of the B þ → D þ D − K þ decay involves consideration of the distribution of the variable hðD þ D − Þ defined as the cosine of the D þ D − helicity angle, i.e., the angle between the momenta of the K þ and D − particles in the D þ D − rest frame. A description of the B þ → D þ D − K þ Dalitz plot is obtained by decomposing the hðD þ D − Þ distribution in terms of Legendre polynomials. The decomposition is done within slices of the D þ D − invariant mass, mðD þ D − Þ, thereby accounting for the two degrees of freedom in the B þ → D þ D − K þ decay kinematics. The description can be projected onto the other invariant-mass distributions in order to identify regions where exotic contributions are needed, and the significance of such deviations can be quantified. If only D þ D − resonances contribute, the projections will be well described using only low-order moments up to twice the maximum spin of the charmonium resonances present. If peaking contributions from other channels enter, higherorder moments will be required. The narrower the structure, the higher the order that will be needed. Consequently, a description employing only low-order moments will be incomplete.
This method is applied to a sample of B þ → D þ D − K þ candidates selected from LHCb proton-proton (pp) collision datasets, corresponding to integrated luminosities of 3 fb −1 recorded during 2011 and 2012 (Run 1) and 6 fb −1 from 2015 to 2018 (Run 2). The data sample, selection criteria, background, and efficiency modeling are identical to those in the amplitude analysis of the same process, described in detail in Ref. [11] and briefly summarized here. The LHCb detector [15,16] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5 that is designed for the study of particles containing b or c quarks. Simulation, produced with the software packages described in Refs. [17][18][19][20], is used to model the effects of the detector acceptance and the imposed selection requirements. The online event selection is performed by a trigger [21], 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 and which identifies a two-, three-, or fourtrack secondary vertex by means of a multivariate algorithm. The charm mesons are reconstructed via the D þ → K − π þ π þ decay. Reconstructed B þ → D þ D − K þ candidates that pass the trigger criteria are subjected to further selection requirements, including the use of a boosted decision tree (BDT) algorithm [22,23] to reduce the combinatorial background. Variables characterizing the particular topology of the decay (flight distance of the D mesons and displacement of the reconstructed intermediate-and final-state particles from the B production point) and particle identification information are used as inputs to the BDT algorithm. Specific requirements are imposed to suppress contributions from B decays involving one or no D mesons but having the same set of final-state pions and kaons as the signal decays; inspection of the sidebands of the D candidates' invariant-mass distributions confirms that any residual background from this source is at a negligible level.
An extended maximum-likelihood fit is applied to the invariant-mass distribution, mðD þ D − K þ Þ, of the selected candidates shown in Fig. 1(a). There are 1260 candidates inside the signal window of mðD þ D − K þ Þ within 20 MeV=c 2 of the known B þ mass [24] in which the sample purity is greater than 99.5% and the residual background is combinatorial in nature. The distribution of these candidates, which are retained for further analysis, in the Dalitz plot is shown in Fig. 1(b). The Dalitz-plot coordinates, m 2 ðD − K þ Þ and m 2 ðD þ D − Þ, are determined after refitting the candidate decays, imposing the constraints that the reconstructed B þ and D AE masses should match their known values and that the reconstructed B þ meson should originate at its associated primary pp interaction vertex. Charmonium resonances are clearly visible as horizontal bands in the Dalitz plot, but additional structure also appears to be present. A signal efficiency map is determined as a function of position in the Dalitz plot with simulation, where the particle identification response is calibrated using data control samples [25,26]. The efficiency is found to vary with mðD þ D − Þ at the AE10% level and to depend only weakly on hðD þ D − Þ.
The mðD þ D − Þ distribution is divided into slices of width 20 MeV=c 2 , which is large compared to the resolution but narrower than any expected structure. Within each slice, the distribution of the cosine of the helicity angle is decomposed according to the basis of Legendre polynomials. Including a factor to ensure normalization over the domain −1 to 1, these are given by In bin j of the mðD þ D − Þ distribution, the coefficient of the expansion at order k is referred to as the kth unnormalized moment, where the sum is over the N Data j candidates in that bin, w l is a weight assigned to each candidate to achieve a background subtraction and efficiency correction, and h l ðD þ D − Þ is the value of hðD þ D − Þ for candidate l. To probe whether charmonium resonances with spins up to and including J max account for the structures observed in the Dalitz plot, the expansion can be truncated at a given order, k max ¼ 2J max .
A simulated sample, generated uniformly in the Dalitz plot and weighted using the truncated expansion, is used in order to visualize the description of the mðD − K þ Þ and mðD þ K þ Þ distributions and to compare them to data. The weights applied to the simulated sample are where i indexes the generated candidates and N Sim j is the number of candidates in the simulation in bin j of the mðD þ D − Þ spectrum, centered on m j ðD þ D − Þ.
The significance of any deviation between the truncated Legendre polynomial description and the data can be assessed using pseudoexperiments. They are generated according to a probability density function (PDF) constructed as a function of mðD þ D − Þ and hðD þ D − Þ given a hypothesis H regarding k max , The binned PDF P½m j ðD þ D − Þ is given by where N is a normalization factor. The PDF P½hðD þ D − ÞjH; m j ðD þ D − Þ is a function of the moments and Legendre polynomial functions, reproducing the helicity angle dependence in bin j of mðD þ D − Þ, Since reflections of exotic contributions to the D − K þ or D þ K þ channels would produce complicated structure in the ½mðD þ D − Þ; hðD þ D − Þ plane, the most sensitive modelindependent test statistic is based on the PDF for mðD − K þ Þ or mðD þ K þ Þ. The PDF P½hðD þ D − ÞjH; m j ðD þ D − Þ is projected onto mðD − K þ Þ or mðD þ K þ Þ by generating candidates uniformly in the ½mðD þ D − Þ; hðD þ D − Þ plane and assigning a weight to each according to Eq. (4). A representation of P½mðD − K þ ÞjH or P½mðD þ K þ ÞjH is then obtained by filling a histogram of mðD − K þ Þ or mðD þ K þ Þ with these weighted candidates, respectively. A test statistic is constructed to discriminate between the hypothesis, H 0 , that only D þ D − resonances contribute up to order k max and the hypothesis that allows for contributions from higher-order moments to describe higher-spin or exotic contributions, H 1 . The test statistic, formulated in terms of determining the significance of deviations in the D − K þ channel, has the form [27] t where P½m l ðD − K þ ÞjH is the value of the PDF in the bin of mðD − K þ Þ where candidate l is found, s l is the signal weight effecting a background subtraction [28], and I H is a normalization factor computed by Monte Carlo integration, where ϵ l is the efficiency appropriate for candidate l. The distributions in the D þ D − invariant mass mðD þ D − Þ of the first nine unnormalized moments defined in Eq. (2) are computed for the selected candidates and are shown in Fig. 2. Significant structure is visible at mðD þ D − Þ ≈ 3.8 GeV=c 2 up to and including the second moment, and not at higher orders, as expected for a contribution from the spin-1 resonance ψð3770Þ. In the vicinity of the χ c2 ð2PÞ resonance near mðD þ D − Þ ¼ 3.9 GeV=c 2 , significant structure appears at order 2 and persists, albeit weakly, at order 4. This is found, in the model-dependent analysis [11], to be due to the presence of both spin-0 and spin-2 charmonia in this region. Structure at low mðD þ D − Þ in the first moment indicates interference between S and P waves and, similarly, that around mðD þ D − Þ ¼ 3.9 GeV=c 2 in the third moment could indicate interference between P and D waves. Structure apparent at all orders for mðD þ D − Þ > 4.1 GeV=c 2 -though having large uncertainties at orders above 5-could indicate reflection from a structure in another two-body combination.
In order to test how well the B þ → D þ D − K þ Dalitz plot can be described using a truncated sum over mðD þ D − Þ moments, a sample of 10 7 B þ → D þ D − K þ decays is generated uniformly in the ½mðD þ D − Þ; hðD þ D − Þ plane. Weights are applied according to Eq. (3), and the resulting distribution of the weighted sample is compared to that for the candidates selected from the LHCb data. In the first instance, k max is set to a high value of 29 in the construction of weights to allow all but the smallest of fluctuations in data to be captured. The comparison between the generated decays and the data sample is shown in Fig. 3. The excellent agreement, limited only by statistical fluctuations that can generate structure to arbitrarily high moments, in the mðD − K þ Þ and mðD þ K þ Þ invariant-mass distributions is also to be expected, given the high value of k max .
The effect of truncating the sum over moments at a lower value is explored. A value of k max ¼ 4 is chosen under the assumption that only resonances with spin up to 2 appear in the D þ D − channel, since production of high-spin resonances in B-meson decays is suppressed and no evidence for a contribution with spin-3 or higher is seen in either Fig. 2 or the model-dependent analysis [11]. Figure 4 shows the comparison between the weighted generated sample and the data. A prominent discrepancy is apparent around mðD − K þ Þ ¼ 2.9 GeV=c 2 . No narrow regions of disagreement are evident in the D þ K þ spectrum.
The significance of the discrepancy in the mðD − K þ Þ distribution between the data and the weighted generated sample in Fig. 4(a) is evaluated using the test statistic defined in Eq. (7). An ensemble of pseudoexperiments, in which each dataset has the same size as the real dataset, is prepared according to the PDF defined in Eq. (6), where k max is taken to be 4. The tiny background contribution is ignored, which introduces negligible uncertainty due to the high purity of the selected B þ → D þ D − K þ sample. For each pseudoexperiment, a new efficiency map is generated to incorporate the systematic uncertainty arising from the limited size of the simulated sample. This ensemble of nearly 260 000 pseudoexperiments allows determination of the distribution of the test statistic under the hypothesis, H 0 , that only D þ D − resonances up to spin-2 are present, as shown in Fig. 5. The value of the test statistic obtained from data, t Data , allows the H 0 hypothesis to be rejected at the 99.994% level, corresponding to a significance of 3.9 Gaussian standard deviations (σ). The impact of allowing moments up to order 6 is investigated with a smaller ensemble of nearly 35 000 pseudoexperiments; the significance of the discrepancy remains above 3.7σ.
In summary, a model-independent technique has been employed to confirm whether or not the observed mðD − K þ Þ distribution in B þ → D þ D − K þ decays reconstructed in the LHCb data sample can be explained in terms of reflections from charmonium resonances alone. It is found that the intermediate structure of the decay cannot be described using only D þ D − resonances of spin up to 2. The significance of the disagreement in the mðD − K þ Þ distribution is 3.9σ and is most apparent in the region mðD − K þ Þ ¼ 2.9 GeV=c 2 . This discrepancy could be explained by a new, manifestly exotic, charm-strange resonance decaying to the D − K þ final state. The outcome of this model-independent study therefore supports the results of the amplitude analysis of the same data [11], where both new spin-0 and spin-1 components are included  FIG. 4. Comparison between data (points with error bars) and a weighted generated sample (filled histogram) as a function of (a) mðD − K þ Þ and (b) mðD þ K þ Þ, where the weights account for the Legendre polynomial moments of orders up to and including 4. The uncertainty on the weighted shape (dark band) is also shown.
We express our gratitude to our colleagues in the CERN accelerator departments for the excellent performance of the LHC. We thank the technical and administrative staff at the LHCb institutes. We acknowledge support from CERN and from the national agencies: CAPES, CNPq, FAPERJ, and