Amplitude analysis of the $B^+\to D^+D^-K^+$ decay

Results are reported from an amplitude analysis of the $B^+\to D^+D^-K^+$ decay. The analysis is carried out using LHCb proton-proton collision data taken at $\sqrt{s}=7,8,$ and $13$ TeV, corresponding to a total integrated luminosity of 9 fb$^{-1}$. In order to obtain a good description of the data, it is found to be necessary to include new spin-0 and spin-1 resonances in the $D^-K^+$ channel with masses around 2.9 GeV$/c^2$, and a new spin-0 charmonium resonance in proximity to the spin-2 $\chi_{c2}(3930)$ state. The masses and widths of these resonances are determined, as are the relative contributions of all components in the amplitude model, which additionally include the vector charmonia $\psi(3770)$, $\psi(4040)$, $\psi(4160)$ and $\psi(4415)$ states and a nonresonant component.


Introduction
Decays of B mesons to multibody final states involving two open-charm mesons and a strange meson, henceforth labelled B → DDK decays, proceed at quark level through b → ccs transitions and comprise a relatively large fraction of the total width of the B mesons. Their branching fractions have been measured previously [1][2][3][4], but few studies of their resonant structure exist. Such analyses are valuable as a means to study resonant structure in both DD and charm-strange systems. Conventional cc charmonium states can produce resonant structures in a neutral DD system, but it is now known that exotic charmonium-like states, which can decay to both neutral and charged DD combinations, also exist [5][6][7]. Conventional resonances can also be observed in charged DK systems, containing charm and antistrange (cs) quarks. 1 There is no previous experimental evidence of exotic hadrons containing a charm and a strange quark (cs), and the possible existence of such states has not been widely discussed in the theoretical literature, although some predictions do exist [8][9][10].
In the B + → D + D − K + decay, resonances in the D − K + channel must have minimal quark contentcdsu and hence would be exotic, as would doubly charged D + K + states. Since conventional resonances can only contribute in the D + D − channel, this B decay stands to provide a clean environment to study charmonium states and to address open questions concerning cc resonant structure, in particular to identify and determine the properties of spin-0 and spin-2 states [11][12][13]. Properties of the vector charmonium states are better known from studies of their production in e + e − collisions, but improved knowledge of their rates of production in B + decays will aid characterisation of the cc contribution in B + → K + µ + µ − decays [14,15]. A more detailed discussion of the current knowledge of charmomium spectroscopy, as relevant to the B + → D + D − K + decay, is given in Sec. 7.1.
No prior study of B + → D + D − K + resonant structure has been published, but a few previous amplitude analyses of other B → DDK decays exist. The Belle collaboration analysed the resonant structure of the B + → D 0 D 0 K + decay [2], while Dalitz-plot analyses of both the B + → D 0 D 0 K + and B 0 → D 0 D − K + final states have been performed by the BaBar collaboration [16]. The signal yields in these previous measurements ranged from about 400 to just under 2000, with relatively high background levels giving a maximum signal purity of 40%. Contributions from the vector ψ(3770) and ψ(4160) charmonium states, and the D * s2 (2573) + and D * s1 (2700) + charm-strange resonances, were determined. A large nonresonant contribution to the B 0 → D 0 D − K + decay was also found.
In this paper the first amplitude analysis of the B + → D + D − K + decay is described. The analysis is based on LHCb proton-proton (pp) collision data taken at √ s = 7, 8, and 13 TeV, corresponding to a total integrated luminosity of 9 fb −1 . In Secs. 2 and 3, the dataset and candidate selection are described. The procedure to determine the signal and background yields, using a fit to the B-candidate invariant-mass spectrum, is presented in Sec. 4. The amplitude modelling formalism used is detailed in Sec. 5, and a description of the selection efficiency and residual background modelling is given in Sec. 6. The development of the model itself follows in Sec. 7, with results given in Sec. 8. Sources of systematic uncertainties that affect the measurements are described in Sec. 9. Studies of the significance of various features in the model are presented in Sec. 10, and a summary that observed in events containing selected B + → D + D − K + candidates. Good agreement is seen between the simulated samples and data for the variables used in the analysis.
The momentum scale is calibrated using samples of J/ψ → µ + µ − and B + → J/ψK + decays collected concurrently with the data sample used for this analysis [36,37]. The relative accuracy of this procedure is estimated to be 3 × 10 −4 using samples of other fully reconstructed b hadrons, Υ and K 0 S mesons.

Selection
Data samples collected in pp collisions during the Run 1 (2011 and 2012) and Run 2 (2015-2018) data-taking periods of the Large Hadron Collider are used, corresponding to integrated luminosities of 3 fb −1 and 6 fb −1 , respectively. Signal B + candidates are built from sets of well-reconstructed pions and kaons, where intermediate charm mesons are reconstructed via the D + → K − π + π + decay. The final-state particles are ensured to be well displaced from the interaction point by requiring that their χ 2 IP with respect to any PV be greater than 4, where χ 2 IP is defined as the difference in the vertex-fit χ 2 of a given PV reconstructed with and without the particle under consideration. The PV that fits best to the flight direction of the B candidate is taken as the associated PV. All charged final-state particles are required to have momentum greater than 1 GeV/c and transverse momentum above 0.1 GeV/c. At least one of them must have momentum greater than 10 GeV/c and transverse momentum exceeding 1.7 GeV/c, whilst also having an impact parameter with respect to the B-candidate's associated PV of at least 0.1 mm. The D-candidates' invariant masses are required to lie within 20 MeV/c 2 of the known D ± mass [38] and their decay vertices must be well reconstructed. The reconstructed momentum and the vector between production and decay vertices are required to be well aligned for both B and D candidates. The flight time (distance significance) from the associated PV for the B-(D)-meson candidates is required to exceed 0.2 ps (6). Finally, PID information is employed to aid identification of final-state K and π mesons.
A boosted decision tree (BDT) [39,40] algorithm, implemented in the TMVA toolkit [41], is employed to separate signal from background. The boosting algorithm assigns weights during training both to correct for classification error and to prioritise uniformity in the Dalitz plot variables. The signal sample for the training consists of correctly reconstructed simulated B + → D + D − K + candidates and the background sample is composed of candidates from the data samples where the B-candidate mass exceeds 5.6 GeV/c 2 . No evidence of overtraining is observed. Candidates are retained if the BDT response exceeds a threshold, chosen to maximise the product of signal significance and sample purity, S 2 / (S + B) 3/2 , where S and B are the expected signal and background yields in the range 5.265 GeV/c 2 < m(D + D − K + ) < 5.295 GeV/c 2 . The invariant mass is calculated from a kinematic fit in which the masses of the charm-meson candidates are fixed to the known D ± mass value and the B meson is constrained to originate from its associated PV. Given the variations in hardware and software trigger criteria, separate BDT classifiers are developed for Run 1 and Run 2 data. The variables entering the BDT are: the χ 2 of the reconstructed B-meson decay vertex; the angle between the B-meson flight direction from the associated PV and its reconstructed momentum; the χ 2 IP of the B-and D-meson candidates and of the final state pions and kaons; the ratio of the flight distance, parallel to the beampipe, of each of the D ± candidates to its uncertainty; and the PID variables of the final-state K and π mesons. Decays of B + mesons to the same set of final-state pions and kaons, having only one or no intermediate D ± mesons or where final-state particles are associated to the wrong D meson, are a potentially important source of background since they produce a peak in the reconstructed D + D − K + invariant-mass distribution. To suppress these backgrounds, vetoes are imposed on narrow invariant-mass structures formed between specific pairs of final-state pions and kaons where the two particles originate from different D ± mesons or the pair involves the kaon produced directly in the B-meson decay. In addition, the two D ± mesons are required to be displaced significantly, parallel to the beampipe, from their production vertex. These requirements are efficient for the B + → D + D − K + signal, and examination of the sidebands of the reconstructed D ± invariant-mass distributions, illustrated in Fig. 1, confirms that there is negligible residual background contamination from this source.
The fraction of events containing more than one reconstructed candidate is measured to be below 1%. All such candidates are retained.

B-candidate invariant-mass fit
An extended maximum-likelihood fit is applied to the m(D + D − K + ) distribution shown in Fig. 2, for candidates in the range between 5.22 and 5.60 GeV/c 2 . The selected candidates in this region are predominantly from signal with a small amount of combinatorial background. There is no significant contribution from partially reconstructed B decays, which appear at lower m(D + D − K + ) values.
The probability density function (PDF) used to model the B + → D + D − K + signal component consists of a double-sided Crystal Ball (DSCB) function [42], having tails on opposite sides of the peak in order to describe the asymmetric power-law tails of the distribution due to detector resolution and final-state radiation. An exponential function accounts for combinatorial background. In the simultaneous fit to each year of the Run 1 and Run 2 datasets, the mean and width of the signal component's Gaussian core are allowed to vary separately for the two periods, and the parameters of the DSCB tails are fixed to their values obtained in fits to simulated samples. The sample purities are very high, so if the background yield falls below 0.01 candidates for one subset of the data, the background component is removed for that subset and the fit re-run to ensure stability. The fit projection is shown in Fig. 2, the yields of the included components are given in Table 1, and the values of the varying parameters are recorded in Table 2.
Of the 1374 candidates to which the invariant-mass fit is applied, 1260 have a value of m(D + D − K + ) within 20 MeV/c 2 of the known B + mass, which is the window applied for the amplitude analysis. Within this signal window, the purity is greater than 99.5%.

Amplitude analysis formalism
The distribution of B + → D + D − K + decays across the Dalitz plot is fitted using the Laura ++ software package [43]. Generic details of the formalism and its implementation in the analysis of LHCb data can be found in the literature [44][45][46]; only aspects specifically relevant to the current analysis are described here. The PDF used to fit the Dalitz plot structure of the selected candidates is composed of signal and background contributions and is a function of position in the B-decay phase space, x. It includes dependence upon model parameters such as mass, width, or spin of individual components in the signal model. The fit procedure maximises the likelihood, where N sig and N bg are the signal and background yields obtained from the invariant-mass fit, respectively, N c is the total number of candidates in the data sample, and P sig,bg ( x j ) are the PDFs for candidate j, which differ for Run 1 and Run 2 data since different efficiency and background models are employed. Gaussian constraints with parameters µ c and σ c are applied to the values of model parameters, p c , such as the masses or widths of intermediate resonances given in Sec. 7. The background PDF, P bg ( x), is an empirical shape used to represent the residual combinatorial background that enters the selected sample of B + → D + D − K + candidates, and is described further in Sec. 6. The signal PDF is given by where N is a normalisation factor that ensures the integral of P sig ( x) over the Dalitz plot ( x) is unity, and total ( x) is the total efficiency for the selected candidates, described further in Sec. 6. The signal amplitude, A sig ( x), is constructed according to the isobar formalism [47][48][49] and contains a coherent sum of resonant and nonresonant amplitudes, where the sum runs over the components in the model, indexed by j. The c j factors are complex coefficients that multiply the complex amplitudes F j ( x), which contain information about the dynamics of each component in the amplitude model. For a D + D − resonance, for example, where R and T describe the invariant-mass and angular dependence of the amplitude, and the X functions are Blatt-Weisskopf barrier factors. The invariant-mass dependence, R (m(D + D − )), is given by a relativistic Breit-Wigner function for all resonant contributions and the angular terms, T ( p, q), are constructed using the non-relativistic Zemach tensor formalism [50,51]. Nonresonant contributions are described with a lineshape that includes an exponential form factor, with alternative models also considered during the model-building and determination of systematic uncertainties. The momenta p and q are those of the third particle (not involved in the resonance) and one of the particles produced in the resonance decay, respectively, both evaluated in the rest frame of the resonance.
The choice of which of the particles produced in the resonance decay is taken to define q corresponds to a convention for the definition of the helicity angle of the resonance. The helicity angle is defined to be, in the rest frame of the resonance, the angle between one of the two particles produced in the resonance decay and the third particle. In this study, the choice is: • θ(D + D − ) is the angle between the K + and D − particles, in the D + D − rest frame, • θ(D + K + ) is the angle between the D − and K + particles, in the D + K + rest frame, and • θ(D − K + ) is the angle between the particles D + and K + , in the D − K + rest frame.
The square Dalitz plot (SDP) provides a useful representation of the phase space. The large B + mass means that resonant structure is often found close to the edge of the regular Dalitz plot, and the SDP provides greater granularity in exactly these regions. Moreover, the SDP aligns a rectangular grid with the edges of the phase space, avoiding edge effects associated with rectangular binning of the regular Dalitz plot.
The two degrees of freedom used to define the SDP are the variables m (D + D − ) and θ (D + D − ), which are defined as where m min,max D + D − are the minimum and maximum kinematically allowed values of m(D + D − ) (equal to 2m D + and m B + − m K + , respectively). With these definitions both m and θ are bounded in the range 0-1.
The complex coefficients, c j in Eq. (3), depend on choices of phase convention and normalisation. In order to be able to compare results between different analyses, it is therefore helpful to report the convention-independent fit fractions, which are defined as the integral of the absolute value of the amplitude squared for each component, j, divided by that of the coherent matrix-element squared for all intermediate contributions, Interference between amplitudes in the coherent sum within A sig ( x) can cause the sum of the fit fractions to depart from unity. This deviation can be quantified by means of interference fit fractions, Interference effects between different partial waves in the same two-body combination cancel when integrated over the helicity angle, due to the angular terms having the form of Legendre polynomials, which form an orthogonal basis.

Efficiency and background models
The absolute efficiency is not needed for the amplitude analysis but the variation of the efficiency across the Dalitz plot must be accounted for. Efficiency variations as a function of position in the Dalitz plot are evaluated using simulated samples. Four contributing factors are considered: The geometrical efficiency, geom , quantifies the probability for all final state particles to be within the LHCb detector acceptance. This efficiency is found not to vary significantly across the phase space. The efficiencies of the trigger requirements, trig , and that of the reconstruction, reco , all with respect to the preceding step, do however have significant dependence on Dalitz-plot position. The BDT, which dominates the offline selection criteria and is designed to minimise induced efficiency variations across the Dalitz plot, behaves as expected with offline being approximately independent of position in phase space. The total efficiency, total , is shown in Fig. 3 as a function of position in both standard Dalitz plot and SDP. Smooth functions are obtained by kernel estimation [35] LHCb (a) and the model obtained using the SDP is used in the analysis to avoid edge effects. Given the differences between Run 1 and Run 2 data for every element of Eq. (9), separate efficiency maps are used for the two data-taking periods. The residual combinatorial background contribution, though small, is accounted for in the fit. A model is derived from candidates in the high B-candidate mass sideband, between 5.35 and 5.69 GeV/c 2 . In order to increase the sample size available for this modelling, the BDT requirement is relaxed by an amount that is seen not to influence the distribution of the background candidates in the Dalitz plot significantly. A kernel estimation procedure is applied to the selected background candidates to reduce the impact of statistical fluctuations. Due to the different selections applied to Run 1 and Run 2 data, both online and offline, separate background models are obtained for each. The background candidates in the regular Dalitz plot are shown in Fig. 4, along with the derived background model as a function of SDP position, obtained using a kernel density estimation [35].   Table 3.
Contributions to the S wave can be expected, but there are few previous experimental results on scalar cc resonances. The Belle collaboration [53] has reported the observation of a χ c0 (3860) state, 2 seen as a DD resonance in the process e + e − → J/ψDD, where the J P C = 0 ++ hypothesis is favoured over the 2 ++ hypothesis at the level of 2.5 σ. This resonance is yet to be confirmed, and there could be other states or nonresonant S-wave DD contributions. The PDG listing [38] includes a X(3915) state, with J P C = 0 ++ or 2 ++ , seen produced in γγ collisions by the Belle [54] and BaBar [55] collaborations (and also possibly in B → X(3915)K decays [56,57]) and decaying to the J/ψω final state -it has not been seen in the DD final state. It appears that this structure may be caused by the χ c2 (3930) state [58], which has also been seen by BaBar to be produced in γγ collisions [59] and has been studied more recently and precisely by LHCb in pp collisions [52]. However, the existence of both spin-0 and spin-2 states near 3930 MeV/c 2 [13,60] is not excluded. At higher mass, the χ c0 (4500) and χ c0 (4700) states have been seen as J/ψφ resonances in Table 3: Components which may appear in the D + D − spectrum of B + → D + D − K + decays, and their properties as given by the Particle Data Group (PDG) [38]. For the ψ(3770) mass and the mass/width of both the χ c2 (3930) and X(3842), the values in Ref. [52] are used.
an LHCb amplitude analysis of B + → J/ψφK + decays [61,62], with masses and widths M = 4506 +16 −19 MeV/c 2 , Γ = 92 ± 29 MeV and M = 4704 +17 −26 MeV/c 2 , Γ = 120 ± 50 MeV respectively. Given that their quantum numbers have been measured as J P C = 0 ++ , these could in principle be seen in B + → D + D − K + decays, but since their composition is unclear it is difficult to make any prediction as to whether this is likely or not.
A larger number of vector cc states have been observed, since these can be produced directly in e + e − collisions. The ψ(3770), ψ(4040), ψ(4160) and ψ(4415) states are all well established and known to decay to DD, therefore all might be expected to appear in B + → D + D − K + decays. The ψ(3770) and ψ(4160) resonances were included in the previous BaBar [16] and Belle [2] amplitude analyses of the B + → D 0 D 0 K + decay, while ψ(4040) and ψ(4415) components were additionally included in an LHCb amplitude analysis of B + → K + µ + µ − decays [63] but found not to contribute significantly. The ψ(4260) state, originally called Y (4260), was observed by the BaBar collaboration through radiative return in e + e − production to the J/ψπ + π − final state [64]. Subsequently confirmed by CLEO, Belle and BESIII collaborations [65][66][67], including through direct e + e − production, it has not been observed in the DD final state, nor is there convincing evidence for its production in B decays. The only ψ(4260) decays to be observed to date contain a J/ψ meson in the final state, although a ψ(4230) state with similar mass and width (M = 4218 + 5 − 4 MeV/c 2 , Γ = 59 + 12 − 10 MeV) has been seen by BESIII to be produced in e + e − collisions in the χ c0 ω, h c π + π − and ψ(2S)π + π − final states [68][69][70]. It is sufficient to consider one of the two as a candidate contribution to the B + → D + D − K + Dalitz plot; the ψ(4260) is used as it is considered to be better established in the PDG 2019 listings. 3 Two further vector states, the ψ(4360) and ψ(4660), have been seen in radiative return from e + e − collisions to the ψ(2S)π + π − final state by the BaBar and Belle collaborations [71,72]. Moreover, a BESIII scan of the energy dependence of the e + e − → J/ψπ + π − cross-section [67] suggests that the structure around 4260 MeV/c 2 is composed of two states: one with M = 4222.0 In the PDG 2019 edition, the results for the first are included in the averages of the properties of the ψ(4260), while those for the second are included in the ψ(4360) averages. Both the ψ(4360) and ψ(4660) are considered unlikely to be present in B + → D + D − K + decays since they have never previously been observed to either be produced in B decays or to decay to DD final states. They are therefore not included in Table 3.
In the D wave, the χ c2 (3930) state has recently been studied by LHCb in pp collisions [52], leading to significant improvement in the knowledge of its properties. However, its quantum numbers are assumed, and while previous analyses have indicated a preference for a spin-2 particle in this mass range [59,73] it is not experimentally excluded that the measured structure is spin-0 or, at least, has a spin-0 contribution. Therefore, it is important to determine the spin of the χ c2 (3930) resonance in this analysis.
Finally, a candidate for the spin-3 ψ 3 (1 3 D 3 ) charmonium state, the X(3842), has recently been observed by LHCb decaying to DD [52]. Its quantum numbers have not been measured, but its properties fit the expectation for that state. Production of spin-3 states in B-meson decays is suppressed, especially when there is little phase space available, and therefore this state is not expected to contribute at a significant level in , for Run 1 and Run 2 respectively. The Dalitz-plot coordinates are determined after refitting the candidate decays, imposing the constraints that the reconstructed B + and D ± masses should match their known values and that the reconstructed B + meson should originate at its associated primary vertex. This improves the resolution of the Dalitz-plot coordinates; for example, the m(D + D − ) resolution is reduced from 10-13 MeV/c 2 to 1.5-3.5 MeV/c 2 , depending upon position in the Dalitz plot. As the resolution is much smaller than the width of the narrowest resonance considered in the analysis, it is neglected in the amplitude fit. A simultaneous fit of the Run 1 and Run 2 datasets is carried out with separate efficiency maps, background models, and fixed signal yields for the two samples. All other model parameters are shared.

Model development
Models which reproduce the Dalitz plot distribution of the data are developed by considering resonances listed in Table 3 and additional resonant and nonresonant components. The ψ(3770) → D + D − and χ c2 (3930) → D + D − resonances, which are both clearly seen in the data, are taken as a starting point. Further components are included in the model if they cause a significant reduction in the negative log-likelihood obtained from the fit to data, while not causing instabilities in the fit or producing excessively large inference effects and hence a sum of fit fractions far from 100%. The complex coefficients associated with all resonant or nonresonant components are allowed to vary freely, with the exception of that for the ψ(3770) → D + D − component, which is fixed to unit length along the real axis to serve as a reference amplitude. The masses and widths of contributing resonances are all allowed to vary, though Gaussian constraints, with parameters corresponding to the central values and uncertainties in Table 3, are applied to those of the ψ(3770), ψ(4040), ψ(4160), and ψ(4415) states. It is observed that significantly better agreement between the model and the data is obtained when including a spin-0 DD component that overlaps with the χ c2 (3930) state, labelled χ c0 (3930). The presence of a spin-0 component in this χ cJ (3930) region may mean that previous measurements of the mass and width of the χ c2 (3930) state, based on an assumption of a single resonance, are not reliable. Therefore, the masses and widths of both the spin-0 and spin-2 components are allowed to vary freely.
It is found that the inclusion of at least one nonresonant component is essential to obtain a good fit to data. A number of parameterisations are considered, including the case of completely uniform Dalitz plot density and modulation of the nonresonant amplitude by either polynomial or exponential form factors, and the possibility of a spin-1, instead of spin-0, angular term. A quasi-model-independent partial wave description of the S wave, as used for example in Refs. [45,46,74], is also attempted, but is not viable with the current sample sizes. In all cases, parameters associated with the nonresonant model are allowed to vary freely in the fit to data.
For each configuration, the minimisation is repeated 100 times, randomising the starting parameters at each iteration. The minimisation that is consistently found to yield the best likelihood value is selected. In order to assess the fit quality, a χ 2 computation is performed, with an adaptive binning scheme ensuring a minimum of 20 candidates in each bin. The associated number of degrees of freedom is determined using an ensemble of pseudoexperiments generated at the fit minimum. The goodness of fit is assessed using this figure of merit as well as the change in negative log-likelihood value between different It is concluded that a satisfactory description of the data cannot be obtained without including one or more components that model structure in m(D − K + ) explicitly. The same conclusion is reached with a model-independent analysis, as described in Ref. [17].

Baseline model including D − K + resonances
The simplest way to account for the m(D − K + ) structure is by adding resonances to the model. Analysis of the current data sample cannot, however, exclude the possibility that hadronic effects such as rescattering may be important, in particular given the observation that the structure appears near to the D * K * threshold. More detailed investigations of plausible explanations for the observed structure will require new theoretical models to be developed and larger data samples to be analysed.
The baseline model includes the same components as in Sec. 8.1, but adds both spin-1 and spin-0 D − K + resonances. An exponential S-wave lineshape in the D − K + channel remains the best description of the nonresonant contribution. The projections of the Dalitz plot, with fit results superimposed, are shown in Fig. 10. In Appendix A, the results are compared to the helicity-angle distributions in eight bins of the invariant-mass distribution of each pair of particles. A comparison to the distributions of the angular moments (defined in Ref. [17]) of each pair of particles is made in Appendix B. The results for the fit parameters and the fit fractions for each component are shown in Tables 4 and 5, where X 1 (2900) and X 0 (2900) are used to label the new spin-1 and spin-0 D − K + states, respectively. These results include systematic uncertainties, the evaluation of which is  described in Sec. 9. The coefficient of the nonresonant exponential lineshape is found to be (0.08 ± 0.05) (GeV 2 /c 4 ) −1 , where the uncertainty is statistical only. The interference fit fractions are given in Table 6, with their statistical and systematic uncertainties.  As described in Sec. 7.1, DD resonant structure has previously been observed in Table 5: Lineshape parameters for the χ c0,2 (3930) and X 0,1 (2900) resonances determined from the fit. The first uncertainty is statistical and the second is the sum in quadrature of all systematic uncertainties.

Resonance
Mass (GeV/c 2 ) Width (MeV) χ c0 (3930) 3.9238 ± 0.0015 ± 0.0004 17.4 ± 5.1 ± 0.8 χ c2 (3930) 3.9268 ± 0.0024 ± 0.0008 34.2 ± 6.6 ± 1.1 X 0 (2900) 2.866 ± 0.007 ± 0.002 57 ± 12 ± 4 X 1 (2900) 2.904 ± 0.005 ± 0.001 110 ± 11 ± 4 the χ cJ (3930) region, however it has usually been assumed to arise from the χ c2 (3930) resonance. The mass and helicity-angle distributions of candidates in this region, shown in Fig. 11, clearly demonstrate that both spin-0 and spin-2 contributions are necessary. The masses and widths of these two components are completely free to vary in the fit; they are found to have consistent masses while the fit prefers a narrower width for the spin-0 state. If both spin-0 and spin-2 states are present at the same mass, one would generically expect the spin-0 state to be broader since its decay to a D + D − pair is in S wave, as compared to D wave for the spin-2 state, and therefore is not suppressed by any angular momentum barrier. This expected pattern is seen in some explicit calculations of the properties of the χ cJ (2P ) states [11], however the observed pattern is consistent with other theoretical predictions [13]. Moreover, the fitted χ c0 (3930) parameters are consistent with those of the X(3915) state.
The χ c0 (3930) state is the only component in the D + D − S wave in the baseline model. The broad χ c0 (3860) state, reported by the Belle collaboration [53], has been included in alternative fit models but is disfavoured. Fits in which additional S-wave structure is introduced e.g. through a nonresonant component, have been attempted but tend to destabilise the fit, which is understood as a consequence of there being too much freedom in the S wave. In fact the nonresonant component in the D − K + projection covers most of the m(D + D − ) range, as can be seen in Fig. 10 top row, but only allows a small contribution at low m(D + D − ) values.    4415) contributions, together with reflections from the D − K + structures. Inclusion of the ψ(4260) resonance was also considered during the model-building process, but its inclusion together with the ψ(4160) state leads to fit instabilities, due to the similarity of their masses and widths. Between the two, a slight preference was visible in negative log-likelihood value for the ψ(4160) component.
The impact of the X 1 (2900) and X 0 (2900) states on the agreement between the data and the model is highlighted in Fig. 12(a) by restricting the phase space to exclude low mass charmonium resonances in the same way as in Fig. 9. The need for both spin-1 and spin-0 components is seen in the helicity-angle distribution shown in Fig. 12(b).

Other models
Numerous variations in the composition of the decay amplitude are considered in the process of establishing the baseline model. These include consideration of one or two states with different spins in the χ cJ (3930) region, and zero, one or two states in the X(2900) region, as well as the inclusion of a contribution from the X(3842) state (assumed to be spin 3). The impact of these different model choices on the negative log-likelihood resulting from the fit is summarised in Table 7. Models with two components with the same spin in the same two-body combination, and with freely varying masses and widths, tend to make the fit unstable and are therefore not included. Similarly, variations in the description of the nonresonant component that destabilise the fit are not included as the obtained negative log-likelihood values are not reliable.
Among the models with variations to the description of the χ cJ (3930) region, those including a spin-1 state (denoted ψ(3930)) are considered unlikely since any vector state in this region would have been seen by previous experiments, as discussed in Sec. 7.1. Moreover, including such a state in the model, either by itself or together with a χ c2 (3930) state, has a large impact on other components of the model. The X 1 (2900) component moves to higher mass and much broader width, with the nonresonant lineshape also changing significantly. These models are therefore excluded from Table 7. The model with  Among the variations in the D − K + channel, the need for two states is clear from the improvement in the NLL and χ 2 values. Noting the proximity to the D * K * threshold, a model with spin-0 and spin-2 states is theoretically well motivated. However, when the masses and widths of the states are allowed to vary freely in the fit, the spin-2 component takes an extremely large (> 500 MeV) width, effectively becoming a nonresonant spin-2 component. While this may be due to residual imperfections in the model (discussed below), this configuration cannot be considered further in the current analysis and is therefore excluded from Table 7. Studies of larger data samples may help to shed light on whether it is possible to describe the structure in m(D − K + ) with spin-0 and spin-2 components. A model with spin-1 and spin-2 D − K + resonances gives comparable, but less favourable, goodness-of-fit indicators to the baseline model.
The model with the inclusion of the X(3842) state, assumed to be spin-3, demonstrates that there is no significant contribution from that component. This supports the assumption, made in Ref. [17], that only states of spin up to 2 are present in B + → D + D − K + decays. Fits with this model are, for simplicity, made neglecting resolution effects since this is done for all other fits. If the narrow X(3842) state were present in the data it would be necessary to account for resolution effects properly, but the fit neglecting them is sufficient to confirm qualitatively the absence of this contribution at any significant level.

Residual imperfections in the baseline model
The goodness of fit is visualised using the binned normalised residual distribution in  Fig. 14, which shows a clear asymmetry most likely originating from interference between the ψ(3770) P-wave state and S-wave D + D − structure. Since the baseline model has only very limited S wave in this region, the asymmetry observed in data cannot be reproduced in the model. This disagreement can also be seen in some other projections, for example at high m(D − K + ) in the projection of the whole Dalitz plot (Fig. 10).
The second of the aforementioned regions of data-model disagreement corresponds to low values of m 2 (D + K + ). No particular disagreement is seen in other projections of this region, and therefore it is not considered a source of concern. There does seem to be some disagreement at high m(D + D − ) values (Fig. 10), but this does not make a large contribution to the χ 2 value. While the region around the ψ(4415) resonance does not appear to be perfectly modelled in the projection, it is probable that at least some of this is statistical, since a very sharp structure at m(D + D − ) ∼ 4.47 GeV/c 2 seems unlikely to be physical.
In summary, while the baseline model does not perfectly reproduce the observed Dalitz-plot distribution, it gives the best description of the currently available data, with a stable fit, among a large range of considered models. Analysis of a larger sample in future will be of great interest to resolve issues associated with the imperfections of the baseline model, as will improved knowledge of D + D − and D − K + structures that may be obtained by analysis of other systems.

Systematic uncertainties
Systematic uncertainties arising from a variety of sources are investigated, and their impact on the model amplitudes, phases, and fit fractions is quantified. The effects on the masses and widths of resonances that are determined from the fit are also evaluated. Sources of systematic uncertainty are separated into those related to experimental effects and those related to model composition. The various systematic uncertainties on the complex coefficients and fit fractions are detailed in Table 8, while those on the masses and widths of resonances are given in Table 9.
The yield of the signal component in the amplitude fit is fixed according to the results of the invariant-mass fit. Repeats of the amplitude fit to data are performed where the signal yield is varied, each time being sampled from a Gaussian PDF centred at the value obtained from the invariant-mass fit, having a width equal to the statistical uncertainty on that yield. The RMS of the values of the fit parameters is taken as the systematic uncertainty. The magnitude of this uncertainty is negligible, and it is therefore omitted from Table 8.
The PDF used to model the signal component in the invariant-mass fit may be imperfect. A conservative estimate of the impact of mismodelling the signal shape is obtained by replacing the DSCB shape by a simple Gaussian function. The deviation of the fit parameters from their nominal values is taken as an estimate of the systematic uncertainty.
The size of the sideband sample limits knowledge of the residual background model in the amplitude fit. An ensemble of bootstrapped sideband data is prepared, from which an ensemble of background models is extracted. Repeated fits to data using the different models are performed, and the RMS of the fit parameters in the resulting ensemble of fit results is taken to represent the systematic uncertainty. This uncertainty is negligible, and is therefore omitted from Table 8.
The effect of the limited size of the simulated samples, used to determine the efficiency model, is quantified. A large ensemble of simulated samples is prepared by bootstrapping the original sample, such that variations within the ensemble are representative of statistical fluctuations expected for the size of that sample. For each variant the efficiency is obtained for Run 1 and Run 2 in the same way as for the nominal efficiency model. The fit to data is then repeated once per efficiency-model-variant, and the RMS of the values of the fit parameters is taken to represent the systematic uncertainty.
The PID response in the data is obtained from calibration samples. The systematic uncertainty incurred through this procedure principally arises from the kernel width used in the estimation of the PDFs. An alternative PID response is simulated using an alternative kernel estimation with changed width, and the efficiency models are regenerated. The fit to data is repeated with these alternative efficiency models in place and the absolute change in the fit parameters is taken as the systematic uncertainty. This uncertainty is omitted from Table 8 since it is negligible.
The hardware-level trigger decision is not expected to be perfectly modelled in the simulated samples. To estimate the impact of this mismodelling of this trigger, a correction obtained from data control samples is applied to the efficiency map. The fit is repeated with this alternative efficiency map and displacement in each parameter is computed. This procedure overestimates the effect, since the mismodelling only affects the efficiency for candidates triggered by hardware-level hadron requirement. Each displacement is therefore scaled according to the fraction of such candidates (64%) to evaluate the systematic uncertainty.
The default Blatt-Weisskopf barrier radii for the parent and intermediate resonances are set to 4.0 GeV −1 . To evaluate the systematic uncertainty arising from the fixed radii, the fit to data is repeated where the radius for each category -parent, charmonia or D − K + resonances -is sampled randomly from a Gaussian distribution centred at 4 GeV −1 and with a width of 1 GeV −1 , which is the approximate size of the uncertainty on the Blatt-Weisskopf barrier radii measured in comparable systems [44]. The RMS of the values of the fit parameters under these perturbations is taken to represent the systematic uncertainty, where the largest effect is seen when varying the Blatt-Weisskopf barrier radius of the charmonium resonances, which dominate the model. This is the largest systematic uncertainty for several of the parameters determined from the fit.
The baseline model includes contributions that are clearly established, but the true amplitude may include components that are not significant at the current level of precision and which are consequently omitted. In addition, the most appropriate way to model some of the components is not established, and mismodelling is a source of potential systematic uncertainty. While many possible model variations could be considered, including too many would lead to an artificial inflation of the uncertainty. Therefore this procedure is limited to specific variations in the partial waves where the modelling uncertainty is largest. With reference to the discussion in Sec. 7.1, these are • D + D − P wave: Inclusion of the ψ(4320) state, with fixed parameters [38].
These effects related to the composition of the amplitude model constitute the largest systematic uncertainty for many of the parameters determined from the fit. The statistical behaviour of the fit is investigated using pseudoexperiments, and the outcome of this study is used to correct the results of the fit to data as summarised in Table 8. The model obtained from the best fit to data is used to generate an ensemble of datasets. Each dataset includes the efficiency variation across the Dalitz plot and a background contribution, the yield of which is sampled for each pseudoexperiment from a Poisson distribution centred at the observed background yield in data. Separate datasets are generated for Run 1 and Run 2 data. The standard fit is then applied to each dataset, where the signal yield is fixed to the generated value. Both the residual, (P fit − P gen ), and normalised residual or "pull", (P fit − P gen ) /σ fit , are determined for the value P of each parameter, determined with uncertainty σ fit , in the fit to each dataset. The distribution of the residual for each fit parameter is fitted with a Gaussian function and the mean ("Bias") is used to correct the central value. The pull distribution for each fit parameter is also fitted with a Gaussian function, and the obtained width ("Pull width") is used to scale the reported statistical uncertainty for the parameter. For the fit fractions, which are calculated from the fitted complex coefficients, the width obtained from the fit of the distribution of the residuals with a Gaussian function is taken as the statistical uncertainty.

Significance of resonant structures
Pseudoexperiments are used to determine the significance of the D − K + structure. The pseudoexperiments are generated using an amplitude model where no D − K + resonances are included, with parameters obtained by fitting the data (see Sec. 8.1). For each dataset, the yields of the signal and background components are sampled from a Poisson distribution centred at the yields observed in the data, and the efficiency is applied to the signal component. Each dataset is fitted with both the model used for generation (H 0 ) and the baseline fit model (H 1 ) and the test statistic t = −2(log(L(H 1 ) − log(L(H 0 ))) is determined. The test statistic observed in data is compared to the distribution from the pseudoexperiments in Fig. 15(a), where the preference for the nominal hypothesis is overwhelming. These results confirm those of Sec. 8.3.
The significance of the X 1 (2900) and X 0 (2900) states in this amplitude analysis is much larger than the significance of exotic contributions obtained in the model-independent analysis of the same data sample [17]. This is expected since in the model-independent analysis the contributions from S, P and D waves in the D + D − system are independent in each m(D + D − ) bin, while in the amplitude analysis each partial wave is a continuous function of m(D + D − ) that is prescribed by the model. The amplitude analysis consequently has less freedom to absorb any structure in the m(D − K + ) distribution compared to the model-independent approach, unless explicit components are included to describe it, and correspondingly a higher significance is obtained.
A similar approach is taken to determine the significance of the presence of both spin-0 and spin-2 states in the χ cJ (3930) region. Three alternative configurations are considered, where these two components are replaced by a single resonance, having spin 0, 1, or 2. The results are shown in Fig. 15. The smallest, though still compelling, significance of the two state fit occurs when comparing to a single spin-1 resonance in the χ cJ (3930) region. Hence the need for two states in this region is clearly established. These results also    confirm those of Sec. 8.3, where issues with fits including a spin-1 state in the χ cJ (3930) region are discussed, leaving the configuration with spins 0 and 2 as the only candidate to describe the data.

Summary
The first amplitude analysis of the B + → D + D − K + decay has been carried out. The analysis is performed using LHCb pp collision data taken at √ s = 7, 8, and 13 TeV, corresponding to a total integrated luminosity of 9 fb −1 , from which a highly pure sample of 1260 signal candidates are selected.
It is not possible to describe the distribution across the Dalitz plot using only resonances in the D + D − system; this conclusion is supported by a model-independent analysis of the same data sample [17]. Reasonable agreement with the data is achieved by including new spin-0 and spin-1 resonances in the D − K + channel, described with Breit-Wigner lineshapes, the parameters of which are determined to be where the first uncertainties are statistical and the second systematic. While the significance of these contributions is overwhelming, and this model gives a good description of the data in this region, it cannot be ruled out that alternative models incorporating additional hadronic effects such as rescattering may also be able to accommodate these D − K + structures. Nonetheless, if the D − K + structures are interpreted as resonances, these results constitute the first clear observation of exotic hadrons with open flavour, and the first that do not contain a heavy quark-antiquark pair. More detailed investigations will require larger data samples and studies of additional decay modes. For example, it will be interesting to see if similar structures can be observed in B + → D − K + π + decays, where an analysis of a subset of the existing LHCb data sample [75] gave an indication of an excess -though not statistically significant -in the m(D − K + ) region where structure is now observed. The model also includes contributions from the ψ(3770), ψ(4040), ψ(4160) and ψ(4415) vector charmonia states. In addition, it is found necessary to include both spin-0 and spin-2 states in the χ cJ (3930) region, the parameters of which are determined from the fit to be Previous measurements of the properties of the χ c2 (3930) state have assumed a single state in this region and, in the light of these results, may be unreliable. There is no evidence for the χ c0 (3860) state reported by the Belle collaboration [53]. Further investigation and independent confirmation of these results concerning spin-0 and spin-2 charmonium states may be obtained in future by studies of B + → J/ψωK + decays. The size and purity of the sample demonstrates the potential impact of further studies of B → DDK decays in the LHCb dataset. In particular, the B + → D 0 D 0 K + mode is likely to shed further light on the production of charmonium states in B-meson decays, while analysis of B 0 → D 0 D − K + may provide crucial additional information on the D − K + structures. In both cases, however, contributions from D + s excitations decaying to D 0 K + will also need to be considered. The significantly larger sample anticipated to be collected by LHCb with an upgraded detector during Run 3 of the Large Hadron Collider also provides exciting prospects for further discoveries in this area.
to the communities behind the multiple open-source software packages on which we depend. Individual groups or members have received support from AvH Foundation (Germany); EPLANET, Marie Sk lodowska-Curie Actions and ERC (European Union); A*MIDEX, ANR, Labex P2IO and OCEVU, and Région Auvergne-Rhône-Alpes (France); Key Research Program of Frontier Sciences of CAS, CAS PIFI, and the Thousand Talents Program (China); RFBR, RSF and Yandex LLC (Russia); GVA, XuntaGal and GENCAT (Spain); the Royal Society and the Leverhulme Trust (United Kingdom).

Appendices A Helicity-angle distributions in slices of Dalitz-plot variables
To allow detailed inspection of the agreement between the result of the fit and the data, helicity angle distributions are shown in slices of the three invariant mass-squared combinations. Figure 16 defines the slices for these projections, with the helicity angle distributions themselves shown in Figs. 17-18.  Fig. 10.

B Angular moments
The angular moments of the data, in bins of m(D + D − ), are central to the modelindependent analysis presented in Ref. [17]. They also present a further way of checking the agreement between the result of the fit and the data. Moments 1-5, for each of m(D + D − ), m(D − K + ) and m(D + K + ) are presented in Fig. 19, with moments 6-9 in Fig. 20.