Dalitz plot analysis of B 0 s → D̄ 0 K − π þ decays

The resonant substructure of B0s→D¯0K−π+ decays is studied with the Dalitz plot analysis technique. The study is based on a data sample corresponding to an integrated luminosity of 3.0  fb−1 of pp collision data recorded by LHCb. A structure at m(D¯0K−)≈2.86  GeV/c2 is found to be an admixture of spin-1 and spin-3 resonances. The masses and widths of these states and of the D∗s2(2573)− meson are measured, as are the complex amplitudes and fit fractions for all the D¯0K− and K−π+ components included in the amplitude model. In addition, the D∗s2(2573)− resonance is confirmed to be spin 2.


Introduction
Several recent experimental discoveries have reinvigorated the field of heavy meson spectroscopy. Among the most interesting are the observations of the D * s0 (2317) − [1] and D s1 (2460) − [2] states below the DK and D * K thresholds, respectively, and hence narrow, in contrast to prior predictions. These are usually interpreted as two of the orbitally excited (1P) states, the other two being the long-established D s1 (2536) − and D * s2 (2573) − resonances, though the reason for the large mass splitting between the mesons below and above the D ( * ) K thresholds is not fully understood. Further interest in the field has been generated by the discovery of several D − sJ states with masses above that of the D * s2 (2573) − resonance through production in e + e − [3,4] or pp [5] collisions. A summary is given in Table 1.
An observation of a state with J P = 3 − would be a clear signature of that state being a member of the 1D family. Production of high-spin states is expected to be suppressed in B meson decay due to the angular momentum barrier [24], and indeed has never yet been observed. However, as the decays of high-spin resonances are suppressed for the same reason, they are expected to have relatively small widths, potentially enhancing their observability.
The Dalitz plot [25] analysis technique has proved to be a powerful tool for studies of charm meson spectroscopy. Analyses by the Belle [26,27] and BaBar [28] collaborations of B → Dππ decays have provided insight into the orbitally excited charm mesons. These results complement those on inclusive production of charm mesons [29][30][31] as the lower background allows broader states to be distinguished and the well-defined initial state allows the quantum numbers to be unambiguously determined. These advantages compensate to some extent for the smaller samples that are available from B meson decay compared to inclusive production.
Until now, few results on charm-strange meson spectroscopy have become available from Dalitz plot analyses, because the available samples of such mesons from B + and Table 1: Excited charm-strange states above the D * s2 (2573) − seen in D ( * ) K spectra by BaBar [4] in e + e − collisions and by LHCb [5] in pp collisions. Units of MeV/c 2 are implied. The first source of uncertainty is statistical and the second is systematic.

State
Mass B 0 decays are much smaller than those of non-strange charm mesons. An exception is a study of B + → D 0 D 0 K + decays by Belle [32], which produced the first observation of the D * s1 (2700) − meson and showed that it has J P = 1 − . Copious samples of charm-strange mesons are, however, available from decays of B 0 s mesons produced at high energy hadron colliders. These have been exploited to study the properties of the D s1 (2536) − [33] and D * s2 (2573) − [34] states produced in semileptonic B 0 s decays. Production of orbitally excited charm-strange mesons has also been seen in hadronic B 0 s decays [35]. In this paper, the first Dalitz plot analysis of the B 0 s → D 0 K − π + decay is presented. The D 0 meson is reconstructed through the K + π − decay mode, which is treated as flavourspecific i.e. the heavily suppressed B 0 s → D 0 K − π + , D 0 → K + π − contribution is neglected. The inclusion of charge conjugated processes is implied throughout the paper. Previously the resonant contribution from B 0 s → D 0 K * (892) 0 has been observed [36] and the inclusive three-body branching fraction has been measured [37]. In this work the contributions from excited charm-strange mesons and excited kaon states are separated from each other with the amplitude analysis technique. The results are important not only from the point-ofview of spectroscopy, but also as they will provide input to future studies of CP violation. In particular, the angle γ of the Cabibbo-Kobayashi-Maskawa Unitarity Triangle [38,39] can be determined from studies of CP violation in B 0 → D 0 K + π − decays [40][41][42]. In such analyses, B 0 s decays provide both an important control channel and a potential source of background (see, e.g., Ref. [43,44]).
The analysis is based on a data sample corresponding to an integrated luminosity of 3.0 fb −1 of pp collision data collected with the LHCb detector, approximately one third of which during 2011 when the collision centre-of-mass energy was √ s = 7 TeV and the rest during 2012 with √ s = 8 TeV. Amplitude analysis techniques have previously been used by LHCb to study B 0 and B 0 s meson decays to J/ψ K + K − [45,46] and J/ψ π + π − [47-50] final states, and to determine the quantum numbers of the X(3872) [51] and Z(4430) [52] resonances. This is, however, the first time that such an analysis has been performed by LHCb with a decay into a fully hadronic final state (i.e. without muons).
The paper is organised as follows. A brief description of the LHCb detector as well as reconstruction and simulation software is given in Sec. 2. The selection of signal candidates and the fit to the B 0 s candidate invariant mass distribution used to separate signal and background are described in Secs. 3 and 4, respectively. An overview of the Dalitz plot analysis formalism and a definition of the square Dalitz plot (SDP) are given in Sec. 5, and details of the implementation of the amplitude analysis are presented in Sec. 6. The evaluation of systematic uncertainties is described in Sec. 7. The results are given in Sec. 8, and a summary concludes the paper in Sec. 9. The highlights of the analysis are described in a shorter companion paper [53].

The LHCb detector
The LHCb detector [54] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector [55] 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 [56] placed downstream of the magnet. The combined tracking system provides a momentum measurement with a relative uncertainty that varies from 0.4% at low momentum, p, to 0.6% at 100 GeV/c, and an impact parameter (IP) measurement with a resolution of 20 µm for charged particles with large momentum transverse to the beamline, p T . Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors [57]. Photon, electron and hadron candidates are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers [58].
The trigger [59] 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. Signal candidates are accepted if one of the final state particles created a cluster in the hadronic calorimeter with sufficient transverse energy to fire the hardware trigger. Events that are triggered at the hardware level by another particle in the event are also retained. The software trigger requires a two-, three-or four-track secondary vertex with a large sum of the p T , of the tracks and a significant displacement from the primary pp interaction vertices (PVs). At least one track should have p T > 1.7 GeV/c and χ 2 IP with respect to any primary interaction greater than 16, where χ 2 IP is defined as the difference in χ 2 of a given PV reconstructed with and without the considered particle.
Simulated events are used to characterise the detector response to signal and certain types of background events. In the simulation, pp collisions are generated using Pythia [60] with a specific LHCb configuration [61]. Decays of hadronic particles are described by EvtGen [62], in which final state radiation is generated using Photos [63]. The interaction of the generated particles with the detector and its response are implemented using the Geant4 toolkit [64] as described in Ref. [65].

Selection requirements
The selection requirements are similar to those used in Refs. [37,66]. The B 0 → D 0 π + π − decay, which is topologically and kinematically similar to the signal mode, is used as a control channel to optimise the requirements and is not otherwise used in the analysis. A set of loose initial requirements is imposed to obtain a visible signal peak in the D 0 π + π − candidates. The tracks are required to be of good quality and to be above thresholds in p, p T and χ 2 IP , while the D 0 → K + π − candidate must satisfy criteria on its vertex quality (χ 2 vtx ) and flight distance from any PV and from the B candidate vertex. Only candidates with 1814 < m(K + π − ) < 1914 MeV/c 2 are retained. A requirement is also imposed on the output of a boosted decision tree that identifies D 0 mesons (with the appropriate final state) produced in b hadron decays (D 0 BDT) [67,68]. The B candidate must satisfy requirements on its invariant mass, χ 2 IP and on the cosine of the angle between the momentum vector and the line from the PV under consideration to the B vertex (cos θ dir ). A requirement is placed on the χ 2 of a kinematic fit [69] to the B decay hypothesis of the final state tracks. The four final state tracks are also required to satisfy particle identification (PID) requirements.
Further discrimination between signal and combinatorial background is achieved with a neural network [70]. The sPlot technique [71], with the B candidate mass as discriminating variable, is used to statistically separate signal from background among the remaining D 0 π + π − candidates. Signal and background weights obtained from this procedure are applied to the candidates, which are then used to train the network. A total of 16 variables are used in the network. They include the χ 2 IP of the four final state tracks and the following variables associated to the D 0 candidate: χ 2 IP ; χ 2 vtx ; the square of the flight distance divided by its uncertainty (χ 2 flight ); cos θ dir ; the output of the D 0 BDT. In addition, the following variables associated to the B candidate are included: p T ; χ 2 IP ; χ 2 vtx ; χ 2 flight ; cos θ dir . Information from the rest of the event is also included through variables that describe the p T asymmetry, A p T , and track multiplicity in a cone with half-angle of 1.5 units in the plane of pseudorapidity and azimuthal angle (measured in radians) [72] around the B candidate flight direction, with where the scalar sum is over the tracks contained in the cone excluding those associated with the signal B candidate. The input quantities to the neural network depend only weakly on position in the B decay Dalitz plot. A requirement imposed on the network output reduces the combinatorial background remaining after the initial selection by a factor of five while retaining more than 90 % of the signal. The B 0 s → D 0 K − π + candidates must satisfy all criteria applied to the D 0 π + π − sample with the exception of the PID requirement on the negatively charged "bachelor" track, i.e. that coming directly from the B 0 s decay, which is replaced with a requirement that preferentially selects kaons. The combined efficiency of the PID requirements on the four tracks in the final state is around 50 % and varies depending on the kinematics of the tracks, as described in detail in Sec. 6.2. The PID efficiency is determined using samples of D 0 → K − π + decays selected in data by exploiting the kinematics of the D * + → D 0 π + decay chain to obtain clean samples without using the PID information [57].
To improve the B 0 s candidate invariant mass resolution, a kinematic fit [69] is used to adjust the four-momenta of the tracks from the D 0 candidate so that their combined invariant mass matches the world average value for the D 0 meson [73]. To improve the resolution, track momenta are scaled [74,75] with calibration parameters determined by matching the measured peak of the J/ψ → µ + µ − decay to the known J/ψ mass [73]. An additional B 0 s mass constraint is applied in the calculation of the variables that are used in the Dalitz plot fit.
To remove potential background from D * ± decays, candidates are rejected if the difference between the invariant mass of the combination of the D 0 candidate and the π + bachelor and that of the D 0 candidate itself lies within ±2.5 MeV/c 2 of the nominal D * − -D 0 mass difference [73]. Candidates are also rejected if a similar mass difference calculated with the pion mass hypothesis applied to the kaon bachelor, satisfies the same criterion. Furthermore, it is required that the kaon from the D 0 candidate together with the bachelor kaon and the bachelor pion do not form an invariant mass in the range 1955-1980 MeV/c 2 to remove potential background from B 0 s → D − s π + decays. Potential background from B 0 s → D 0 D 0 decays [68] is removed by requiring that the pion and kaon originating directly from the B 0 s decay give an invariant mass outside the range 1835-1880 MeV/c 2 . At least one of the pion candidates is required to have no associated hits in the muon counters to remove potential background from B 0 → J/ψ K * 0 decays. Decays of B 0 s mesons to the same final state but without an intermediate charm meson are suppressed by the D 0 BDT criteria, and any surviving background from this source is removed by requiring that the D 0 candidate vertex is displaced by at least 1 mm from the B 0 s decay vertex. Figure 1 shows the D 0 candidate mass after the selection criteria are applied.
Signal candidates are retained for further analysis if they have an invariant mass in the range 5200-5900 MeV/c 2 . After all selection requirements are applied, fewer than 1 % of events with one candidate also contain a second candidate. Such multiple candidates are retained and treated in the same manner as other candidates; the associated systematic uncertainty is negligible.

Determination of signal and background yields
The signal and background yields are obtained from an extended unbinned maximum likelihood fit to the invariant mass distribution of B 0 s → D 0 K − π + candidates. In addition to signal decays and combinatorial background, the fit allows background contributions from other b hadron decays. The decay B 0 s → D * 0 K − π + , with D * 0 → D 0 π 0 or D 0 γ forms a partially reconstructed background that peaks at values below the B 0 s mass since the π 0 or γ is missed. Decays of B 0 mesons to the D 0 K − π + final state are Cabibbosuppressed, but may contribute a non-negligible background. Decays with similar topology Here the selection criteria have been modified to avoid biasing the distribution: the D 0 candidate invariant mass requirement has been removed, and the χ 2 of the kinematic fit is calculated without applying the D 0 mass constraint. and misidentified final state particles can also populate the mass region used in the fit. Studies using simulated background events show that contributions from B 0 → D ( * )0 π + π − and Λ 0 b → D ( * )0 pπ + [76] are expected, while background from B 0 (s) → D ( * )0 K + K − [77, 78] and Λ 0 b → D ( * )0 pK + is negligible. The signal and B 0 → D 0 K − π + shapes are each modelled with the sum of two Crystal Ball [79] functions which share a common mean and have tails on opposite sides. Studies using simulated events and the B 0 → D 0 π + π − control channel in data verify that this function gives an excellent description of the signal shape. The tail parameters are fixed to values determined from a fit to simulated signal decays. The mass difference between the peaks corresponding to B 0 and B 0 s decays is fixed to its known value [73]. The combinatorial background is modelled using a linear shape.

Parameter
Value Sec. 3, while those for (anti)protons are obtained from samples of Λ → pπ − decays. There are in total eleven free parameters determined by the fit: the peak position of the signal shape, the widths and the fraction of the shape contained within the narrower of the two Crystal Ball functions, the linear slope of the combinatorial background, and the yields of the six categories defined above. The results of the fit are shown in Fig. 2 and listed in Table 2. The fit gives a reduced χ 2 of 98.6/88 = 1.12. All yields are consistent with their expectations, based on measured or predicted production rates and branching fractions, and efficiencies or background rejection factors determined from simulations.
For the Dalitz plot analysis a signal region is defined as µ B 0 s ± 2.5σ 1 , where µ B 0 s and σ 1 are the peak position and core width of the signal shape, respectively, and are taken from

Component
Yield

Dalitz plot analysis formalism
The Dalitz plot [25] describes the phase-space of the three-body decay in terms of two of the three possible two-body invariant mass squared combinations. In B 0 s → D 0 K − π + decays, resonances are expected in the m 2 (D 0 K − ) and m 2 (K − π + ) combinations, and therefore this pair is a suitable choice to define the Dalitz plot axes. Given these two invariant mass squared combinations all other kinematic quantities can be uniquely determined.
The description of the complex amplitude is based on the isobar model [80][81][82], which describes the total amplitude as a coherent sum of amplitudes from resonant or nonresonant intermediate processes. As such the total amplitude is given by where c j are complex coefficients giving the relative contribution of each different decay channel. The resonance dynamics are contained within the F j m 2 (D 0 K − ), m 2 (K − π + ) terms, which are normalised such that the integral over the Dalitz plot of the squared magnitude of each term is unity and are composed of invariant mass and angular distributions.
where the functions R, X and T described below depend on parameters of the resonance such as its spin L, pole mass m 0 and width Γ 0 . In the case of a D 0 K − resonance, the π + is referred to as the "bachelor" particle and, since the B 0 s meson has zero spin, L is equivalently the orbital angular momentum between the resonance and the bachelor.
In Eq. (3), the function R m(D 0 K − ) is the resonance mass term (given e.g. by a Breit-Wigner shape), while p and q are the momenta of the bachelor particle and one of the resonance daughters, respectively, both evaluated in the rest frame of the resonance. The terms X(z), where z = | q | r BW or | p | r BW , are Blatt-Weisskopf barrier form factors [24], and are given by where z 0 represents the value of z when the invariant mass is equal to the pole mass of the resonance. The radius of the barrier, r BW , is taken to be 4.0 GeV −1 ≈ 0.8 fm [83]. The angular probability distribution terms, T ( p, q), are given in the Zemach tensor formalism [84,85] by which can be seen to have similar forms to the Legendre polynomials, P L (x), where x is the cosine of the angle between p and q (referred to as the "helicity angle"). The majority of the resonant contributions in the decay can have their mass terms described by the relativistic Breit-Wigner (RBW) function where the dependence of the decay width of the resonance on m is given by where the symbol q 0 denotes the value of q = | q | when m = m 0 . This shape can also describe so-called virtual contributions, from resonances with pole masses outside the kinematically accessible region of the Dalitz plot, with one modification: in the calculation of the parameter q 0 the pole mass is set to a value within the kinematically allowed range. This is accomplished with the ad-hoc formula where m max and m min are the upper and lower limits, respectively, of the kinematically allowed mass range. For virtual contributions only the tail of the RBW function enters the Dalitz plot. Due to the large phase-space available in three-body B meson decays, it is possible to have nonresonant amplitudes (i.e. contributions that are not associated with any known resonance, including virtual states) that are not, however, constant across the Dalitz plot.
A common approach to model nonresonant terms is to use an exponential factor [86], where α is a shape parameter that must be determined from the data. The RBW function is a very good approximation for narrow resonances well separated from any other resonant or nonresonant contribution in the same partial wave. This approximation is known to be invalid in the Kπ S-wave, since the K * 0 (1430) resonance interferes strongly with a slowly varying nonresonant term (see, for example, Ref. [87]). The so-called LASS lineshape [88] has been developed to combine these amplitudes, where cot δ B = 1 aq and where m 0 and Γ 0 are now the pole mass and width of the K * 0 (1430), and a and r are parameters that describe the shape. Most implementations of the LASS shape in amplitude analyses of B meson decays (e.g. Refs. [83,89]) have applied a cut-off to the slowly varying part close to the charm mass. The value of the cut-off used in this analysis is 1.7 GeV/c 2 .
In the absence of any reconstruction effects, the Dalitz plot probability density function would be where the dependence of A on the Dalitz plot position has been suppressed in the denominator for brevity. In a real experiment, the variation of the efficiency across the Dalitz plot and the contamination from background processes must be taken into account.
Since signal and background events tend to populate regions close to the kinematic boundaries of the conventional Dalitz plot, it is convenient to model the efficiencies and backgrounds using the so-called square Dalitz plot (SDP) defined by variables m and θ that have validity ranges between 0 and 1 and are given by where m max is the helicity angle of the D 0 K − system (the angle between the π and the D meson in the D 0 K − rest frame).
The primary results of a Dalitz plot analysis are the complex amplitudes given by c j in Eq. (2) that describe the relative contributions of each resonant component. However, the choice of normalisation, phase convention and amplitude formalism may not be the same for different implementations. Fit fractions and interference fit fractions provide a convenient convention-independent method to allow meaningful comparisons of results. The fit fraction is defined as the integral of a single decay amplitude squared divided by that of the coherent matrix element squared for the complete Dalitz plot, The sum of these fit fractions is not necessarily unity due to the potential presence of net constructive or destructive interference quantified by interference fit fractions defined for i < j only by where the dependence of F ( * ) i on the Dalitz plot position has been omitted. 6 Dalitz plot fit 6

.1 Background square Dalitz plot distributions
There are non-negligible background contributions in the signal region from combinatorial background and from B 0 → D ( * )0 π + π − and Λ 0 b → D ( * )0 pπ + decays. As shown in Table 3, these sources correspond to 7.4 %, 2.8 % and 2.3 % of the total number of candidates in the signal region, respectively, and therefore their Dalitz plot distributions need to be modelled. Small contributions from other sources of background are neglected. The shapes of all background sources in the SDP are described by histograms and are shown in Fig. 4.
The combinatorial background distribution is obtained from candidates in a high B 0 s mass sideband, in the range 5500-5900 MeV/c 2 . The result of the invariant mass fit described in Sec. 4 shows that this region contains only combinatorial background and a small amount of B 0 → D ( * )0 π + π − decays. The latter component is modelled using simulated decays as described below and subtracted from the sideband distribution. A sample of D 0 K ± π ± candidates is used to verify that the SDP distribution of combinatorial background does not depend significantly on the B 0 s candidate invariant mass, and therefore the sideband distribution can be considered a reliable description of the background in the signal region.
The SDP distributions of the Λ 0 b → D ( * )0 pπ + and B 0 → D ( * )0 π + π − backgrounds are derived from simulated events. In each shape, the components from the final states containing D 0 and D * 0 mesons are combined and the simulated samples reweighted as described in Sec. 4. The dominant contribution in the signal region comes, for both shapes, from the final state with a D 0 , not a D * 0 , meson.

Efficiency variation across the square Dalitz plot
Variation of the signal efficiency across the SDP is induced by the detector acceptance and by trigger, selection and PID requirements. The variation of the efficiency is studied using simulated samples of signal events generated uniformly over the SDP, with several data driven corrections. Statistical fluctuations from limited sample size are smoothed out by fitting the efficiency functions to a two-dimensional cubic spline across the SDP.
Corrections are applied for known differences between data and simulation in the tracking, trigger and PID efficiencies. A tracking correction is obtained for each of the four final state tracks as a function of η and p. The total correction is obtained from the product of the factors for each track.
The trigger efficiency correction is different for two mutually exclusive subsamples of the selected candidates. The first includes candidates that are triggered at hardware level by clusters in the hadronic calorimeter created by one or more of the final state particles, and the second contains those triggered only by particles in the rest of the event. For the first subsample, a correction is calculated from the probability of an energy deposit in the hadronic calorimeter to fire the trigger, evaluated from calibration data samples as a function of particle type, dipole magnet polarity, transverse energy and position in the calorimeter. In the second subsample, a smaller correction is applied to account for the requirement that the signal decay products did not fire the hadronic calorimeter hardware trigger. The efficiency is evaluated for each subsample as a function of SDP position, and these are combined into a single efficiency map according to their proportions in data.
The PID efficiency is evaluated using a calibration sample of D 0 → K + π − decays as described in Sec. 3. Efficiencies for background-subtracted samples of kaons and pions are obtained as functions of their p, p T and of the number of tracks in the event. The kinematic properties of the four final state signal particles are obtained from simulation while the distribution of the number of tracks in the event is taken from data. Efficiencies for each of the final state particles are evaluated and their product gives the efficiency for the candidate accounting for possible correlations between the kinematics of the four tracks.
Contributions from the various sources are then combined into a single efficiency map across the SDP that is used as an input to the Dalitz plot fit and is shown in Fig. 5. The largest source of variation arises due to the reconstruction, which causes a rapid drop of the efficiency at the smallest values of m , which corresponds to high m(D 0 K − ) and hence slow π + tracks. The largest source of efficiency variation induced by the selection arises due to the PID requirements, which lead to a maximum efficiency variation of about ±20 % across the SDP.

Amplitude model for
The Dalitz plot fit is performed using the Laura++ [90] package. The likelihood function that is optimised is given by where the indices i and k run over the N c selected candidates and the signal and background categories, respectively. The signal and background yields N k are given in Table 3. The signal probability density function P sig is a modified version of Eq. (18), where factors of |A m 2 (D 0 K − ), m 2 (K − π + ) | 2 in both numerator and in the integral in the denominator are multiplied by the efficiency function described in Sec. 6.2. The mass resolution is below 2 MeV/c 2 , much less than the width of the narrowest structures on the Dalitz plot, and therefore has negligible effect on the likelihood. The background SDP distributions are discussed in Sec. 6.1 and shown in Fig. 4. The free parameters of the fit are the real and imaginary parts of the complex coefficients, c j in Eq. (2), for each amplitude included in the fit model, except for the D * s2 (2573) − component for which the real and imaginary parts of the amplitude are fixed to 1 and 0, respectively, as a reference. Several parameters of the lineshapes are also determined from the fit, as described below. Results for the complex amplitudes are also presented in terms of their magnitudes and phases, and in addition the fit fractions and interference fit fractions are determined. Instead of being parameters of the fit, these are obtained from the fit result using large samples of simulated pseudoexperiments. This approach allows effects of non-trivial correlations between fit parameters to be appropriately treated.

Resonance
Spin Dalitz plot axis Model It is possible for the minimisation procedure to find a local minimum of the negative logarithm of the likelihood (NLL) function. Therefore to find the true global minimum the fit is repeated many times with randomised initial values of the complex amplitude.
The baseline amplitude model for B 0 s → D 0 K − π + decays is defined by considering many possible resonant, virtual or nonresonant contributions and removing those that do not significantly affect the fit. Resonances with unnatural spin-parity, that do not decay to two pseudoscalars, are not considered. The resulting signal fit model consists of the contributions shown in Table 4. There are a total of fourteen components: six K − π + resonances, four D 0 K − resonances, three virtual resonances and a D 0 K − nonresonant contribution. The majority are modelled with the RBW lineshape, the exceptions being: (i) the K − π + S-wave, including the K * 0 (1430) 0 resonance, which is modelled by the LASS lineshape with an additional contribution from the K * 0 (1950) 0 state; and (ii) the D 0 K − nonresonant component, which is modelled with an exponential form factor (EFF).
As discussed further in Sec. 8, a highly significant improvement in the likelihood is obtained when including two resonances, one spin-1 and another spin-3, both with m(D 0 K − ) ≈ 2.86 GeV/c 2 . Previous studies of the D * sJ (2860) − state [4,5], have assumed a single resonance in this region, and therefore values of the mass and width obtained from those analyses cannot be used in the fit. Instead, the parameters of these states are obtained from the data. The sensitivity of the data to the parameters of the D * s2 (2573) − resonance exceeds that of previous measurements [73], and therefore these parameters are also obtained from the fit.
The slope parameter, α, of the EFF model for the D 0 K − nonresonant contribution,  Table 5 for the fit fractions and complex coefficients, and in Table 6 for the masses and widths. Results for the interference fit fractions are presented in App. A. In Table 5, and for all results for fit fractions, values are given both for the nonresonant and K * 0 (1430) 0 parts of the LASS function separately and for the two combined taking into account their interference. The interference effects between the components of the K − π + S-wave explain most of the excess of the total fit fraction from unity. Other local minima of the NLL function are found to be separated from the global minimum by at least 10 units.
The fit quality is evaluated by determining a χ 2 value by comparing the data and the fit model in N bins = 576 SDP bins that are defined adaptively to ensure approximately equal population with a minimum bin content of 21 entries. The effective number of degrees of freedom of the χ 2 is bounded by N bins − N pars − 1 and N bins − 1, where N pars is the number of parameters determined by the data. The former choice gives a higher reduced χ 2 value of 1.21, where only statistical uncertainties are included in the calculation. The effects of systematic uncertainties on the χ 2 value are discussed at the end of Sec. 7. The distribution across the SDP of the pull, defined as the difference between the data and the fit model divided by the uncertainty, is shown in Fig. 6. Other unbinned tests [91] of the fit quality also show that the fit provides a good, but not perfect, model of the data.
Projections of the data and the baseline fit result onto m(K − π + ), m(D 0 K − ) and   m(D 0 π + ) are shown in Fig. 7. The dip visible in m(K − π + ) is due to the D 0 veto described in Sec. 3. Zooms around the main resonant contributions are shown in Fig. 8. Good, but not perfect, agreement between the data and the fit is seen. Further comparisons of regions of the data with the fit result are given in Figs. 9 and 10. These show projections of the cosine of the helicity angle of the K − π + and D 0 K − systems, respectively, and show that the spin content of the fit model matches well that of the data. In particular, Fig. 10(d) shows that the region around the D * sJ (2860) − states is well modelled by a combination of spin-1 and spin-3 states. This is confirmed by the χ 2 value of 56 that is found by comparing the data and the fit model in only the 70 SDP bins, defined with the adaptive binning scheme discussed above, that overlap or are contained in this region of phase-space (0.71 < m < 0.77). The distinctive angular distribution of the spin-3 state enables the comparatively precise determination of its properties (Table 6).
To test whether any other combination of resonances can provide a comparably good description of the data, the fit is repeated with different hypotheses. The results are shown in Table 7. The values of √ 2∆NLL are given as a crude indication of the significance but  are not otherwise used in the analysis -numerical values of the significance are instead obtained from pseudoexperiments as described in Sec. 8. Some of the results in Table 7 are labelled with * to indicate that the fit prefers to position one of the resonances in a different mass region. For spin-0 this is subthreshold, for spin-2 it is either very near to the D * s2 (2573) − mass or at higher mass. The spin of the D * s2 (2573) − state has not previously been determined experimentally [73]. As seen in Fig. 10(b), the helicity angle distribution in this region follows closely the expectation for a spin-2 state. No alternative spin hypothesis can give a reasonable description of the data -the closest is a fit assuming spin-0, which gives a value of √ 2∆NLL above 40. The helicity angle distributions for the best fits with spin-2 and spin-0 hypotheses are compared to the data in Fig. 11.
Another approach to assess the agreement between the data and the fit result is to compare their angular moments, obtained by weighting the events in each m(D 0 K − ) (m(K − π + )) bin by the Legendre polynomial of order L in cos θ(D 0 K − ) (cos θ(K − π + )), where θ(D 0 K − ) (θ(K − π + )) is the angle between the π + and the D 0 meson (the D 0 and the K − meson) in the D 0 K − (K − π + ) rest frame. This approach is very powerful in the case that resonances are only present in one invariant mass combination, since then structures are seen in moments up to 2 × J max , where J max is the highest spin of the contributing resonances. When resonances in other invariant mass combinations cause reflections, higher moments are introduced in a way that is hard to interpret.
The angular moments of the data and the fit model in m(D 0 K − ) and m(K − π + ) are compared in Figs. 12 and 13, respectively. Significant structures in the K * (892) 0 peak region are observed in moments up to order 2, as expected for a spin-1 resonance in the absence of reflections. The moments in the regions of other resonances are affected by reflections, as can be seen in the Dalitz plot (Fig. 3). Nonetheless, the large structures in the D * s (2573) − peak region in moments up to order 4 unambiguously determine that its spin is 2. At higher masses, interpretation of the moments becomes more difficult. Nonetheless, the reasonable agreement between data and the fit model provides confidence that the two-dimensional structures in the data are well described.

Systematic uncertainties
The considered sources of systematic uncertainty are divided into two main categories: experimental and model uncertainties. The experimental systematic uncertainties arise from imperfect knowledge of: the relative amount of signal and background in the selected events; the distributions of each of the background components across the SDP; the variation of the efficiency across the SDP; the possible bias induced by the fit procedure; the momentum calibration; the fixed masses of the B 0 s and D 0 mesons used to define the boundaries of the Dalitz plot. Model uncertainties occur due to: fixed parameters in the Dalitz plot model; the decision to include or exclude marginal components in the baseline fit model; the choice of models for the K − π + S-wave and the D 0 K − S-and P-waves. The systematic uncertainties from each source are combined in quadrature. The data are shown as black points, the result of the baseline fit with a spin-2 resonance is given as a solid blue curve, and the result of the fit from the best model with a spin-0 resonance is shown as a dashed red line.   the Dalitz plot fit are assigned as uncertainties.
The uncertainty due to the imperfect knowledge of the background distributions across the SDP is estimated by varying the histograms used to model the shape within their statistical uncertainties. In addition, the relative contributions from decays with D 0 and D * 0 mesons in the Λ 0 b → D ( * )0 pπ and B 0 → D ( * )0 π + π − backgrounds are varied. The effect on the results of not reweighting the SDP distribution of the D * 0 component in these backgrounds is also included as a source of systematic uncertainty.
The uncertainty arising due to the imperfect knowledge of the efficiency variation across the SDP is determined by varying the content of the histogram from which the spline function used in the fit is obtained. Since sources of systematic bias may affect the bins of this histogram in a correlated way, only the central bin in each cell of 3 × 3 bins is varied, and interpolation is used to obtain the values of the adjacent bins. The effects on the results of the Dalitz plot fit are assigned as uncertainties. In addition, the effect of binning the D 0 → K − π + control sample used to obtain the PID efficiencies is evaluated by varying the binning scheme.
An ensemble of pseudoexperiments is used to search for intrinsic bias in the fit procedure. The differences between the inputs and the mean values obtained from the ensemble are all found to be small. Systematic uncertainties are assigned as the sum in quadrature of the difference between the input and output values and the uncertainty on the mean from the fit to the ensemble of pseudoexperiments.
The uncertainty due to the momentum calibration is estimated by varying the calibration factor within its uncertainty [74,75]. The differences with respect to the default results are assigned as the corresponding systematic uncertainties.
The masses of the B 0 s and D 0 mesons are fixed to their known values [73] when the Dalitz plot coordinates are calculated. The analysis is repeated after varying the B 0 s and D 0 meson masses up and down by one standard deviation independently, and the changes in the fitted values are taken as the corresponding uncertainty.
The uncertainties due to fixed model parameters are evaluated by repeating the fit after varying these parameters within their uncertainties, with correlations taken into account where appropriate. The parameters that are modified are the masses and widths given in Table 4 and the Blatt-Weisskopf radius parameter, which is varied between 3 and 5 GeV −1 . As a cross-check, different Blatt-Weisskopf radius parameters are used for the K − π + and D 0 K − resonances, and the likelihood is minimised with respect to these parameters with results r BW (K − π + ) = 3. The models used to describe the K − π + S-wave and the D 0 K − S-and P-waves are known to be approximate forms, and therefore additional uncertainties are assigned due to the changes in the fitted values of the other parameters when these are replaced with alternative models. The LASS shape is replaced with a Flatté shape [92] for the K * 0 (1430) and a resonant term with a modified mass-dependent width for the κ (or K * 0 (800)) resonance at low m(K − π + ) [93]. The alternative model for the K − π + S-wave given in Ref. [94] is also used to fit the data, with the larger variation from the two alternative models assigned as systematic uncertainty. A K-matrix implementation of the K − π + S-wave [95] is also attempted but does not provide stable fit results. As an alternative model for the D 0 K − S-wave, the exponential form factor is replaced with a power-law dependence. To estimate the dependence of the results on the modelling of the D 0 K − P-wave, the two broad spin-1 D 0 K − resonances (D * s1 (2700) − and D * s1 (2860) − ) are described with a modified version of the Gounaris-Sakurai lineshape [96] instead of relativistic Breit-Wigner functions. In addition, the dependence of the results on the choice of description of the effective pole mass for virtual components [Eq. (14)] is evaluated by using a constant width instead of Eq. (13).
Summaries of the experimental systematic uncertainties on the fit fractions and complex amplitudes are given in Table 8. A breakdown is given in Table 9 for the fit fractions, and in Table 10 for the masses and widths. Similarly, summaries of the model uncertainties on the fit fractions and complex amplitudes are given in Table 11, with breakdowns for the fit fractions and masses and widths in Tables 12 and 13, respectively. The largest sources of experimental systematic uncertainties on the fit fractions are, in general, those due to the efficiency variation across the SDP, the signal and background fractions and the description of the background SDP distributions. The largest sources of model uncertainties on these parameters are, in general, from the description of the K − π + S-wave and from removing the K * (1680) 0 and B * + v components from the model. These are also the largest sources of uncertainty on the mass and width measurements. The magnitudes of the complex amplitudes are more robust against systematic uncertainties than the relative phases.
The reduced χ 2 value of 1.21 obtained by comparing the data and the default fit model in SDP bins, discussed in Sec. 6.3, corresponds to a tiny p-value, given the large number of degrees of freedom. Such a situation is not uncommon for high statistics Dalitz plot analyses, see e.g. Refs. [26,28]. Moreover, the χ 2 is evaluated accounting only for statistical uncertainties. Some disagreement between the data and the fit model is visible in the helicity angle projections in the regions of the peaks with the largest statistics, namely the K * (892) 0 ( Fig. 9(b)) and the D * s2 (2573) − (Fig. 10(b)) resonances. The latter is also visible in Fig. 8(a) as the reflection from one lobe of the D * s2 (2573) − structure overlaps with K * 0 resonances in the m(K − π + ) ≈ 1430 MeV/c 2 region. These regions correspond to bins with large pulls in Fig. 6. The small peak in Fig. 8(d) at m(D 0 K − ) ≈ 2.96 GeV/c 2 is not statistically significant.
As seen in this section, both experimental systematic and model uncertainties are comparable in size to the statistical uncertainties on the parameters associated with those resonances, suggesting that these uncertainties may significantly affect the χ 2 value. In addition, certain aspects of the modelling, such as the description of the K − π + and D 0 K −  S-waves, are known to be approximations. The default model gives the best agreement with the data among the alternatives considered. Nonetheless, the change in reduced χ 2  value when alternative models are used, which is typically in the range 0.05-0.10, gives an estimate of how much the approximations used may affect the goodness-of-fit. Therefore, the description of the data is considered to be acceptable.

Results
As discussed in Sec. 6.3, the data require both a spin-1 and a spin-3 resonance in the m(D 0 K − ) ≈ 2.86 GeV/c 2 region. Figure 14 shows the result of the baseline fit compared to alternative models containing only a single resonance, either spin-1 or spin-3, in this region. As shown in Table 7 These significances include only statistical uncertainties, so the effect of the largest systematic uncertainties is tested by repeating the procedure with the variations in the models discussed in Sec. 7 that give the largest effects on the fit fractions, masses and widths of the D * sJ (2860) − states. For the spin-1 only model, the effect of using the κ model to describe the K − π + S-wave is evaluated. For the spin-3 only model, the κ description of the K − π + S-wave, the addition of the K * 4 (2045) 0 state and the variation of the D 0 mass are considered. The conclusion is that two states are required in this region with significance of at least 10 standard deviations.
The masses and widths of these three states are determined to be where the first uncertainty is statistical, the second is due to experimental systematic effects and the third due to model variations. The results for the complex amplitudes, expressed both as real and imaginary parts and as magnitudes and phases, are given in Table 14. The results for the fit fractions are given in Table 15, while results for the interference fit fractions are given in App. A. For resonances without a significant signal, it is possible to set upper limits on their fit fractions, and therefore on their branching fractions. This is done for the K * (1680) 0 ,  Table 7). obtained, and converted into likelihood functions. The effect of systematic uncertainties is included by convolving the likelihood function with a Gaussian of width given by the systematic uncertainty. These are then used to set 90 % and 95 % confidence level (CL) upper limits by integrating the likelihood. The upper limits obtained with this procedure are included in Table 15.

Summary
The first amplitude analysis of the B 0 s → D 0 K − π + decay has been presented. The B 0 s → D 0 K − π + decay amplitude model contains a total of fourteen components: six K − π + resonances, four D 0 K − resonances, three virtual resonances and a nonresonant contribution. The complex amplitudes of these are determined, and fit fractions and interference fit fractions are reported in addition to enable convention-independent comparisons of the model. The fit fraction results are converted into branching fraction measurements. The result for B(B 0 s → D 0 K * (892) 0 ) is significantly more precise than the previous measurement [36], which was obtained from a much smaller and statistically independent data sample collected by LHCb during 2010. All other branching fraction results are first reported measurements.
A structure at m(D 0 K − ) ≈ 2.86 GeV/c 2 is found to be an admixture of a spin-1 and a spin-3 resonance with a significance of at least 10 standard deviations. Therefore the D * sJ (2860) − state previously observed by the BaBar collaboration in inclusive e + e − → D 0 K − X production [4] and by the LHCb collaboration in pp → D 0 K − X processes [5] consists of at least these two resonances. The properties of those states and of the D * s2 (2573) − resonance are measured. The spin of the D * s2 (2573) − resonance is experimentally determined for the first time, and is confirmed to be 2. The mass and width of this state are determined with significantly better precision than previous measurements [73]. The result for the width is consistent with the previous world average. The result for the mass, however, is somewhat below the previous average, which is dominated by a measurement by the BaBar collaboration [3] based on inclusive production in e + e − collisions. The Dalitz plot analysis technique used in this paper ensures that the background under the D * s2 (2573) − peak is small and does not contain contributions from decays of higher D − s resonances, resulting in a much lower systematic uncertainty compared to the inclusive approach.
The masses of the D * s1 (2860) − and D * s3 (2860) − states are found to be consistent within uncertainties, while a larger width of the spin-1 state than of the spin-3 state is preferred. These results appear to support an interpretation of these states being the J P = 1 − and 3 − members of the 1D family, though the 1 − state may be partially mixed with the vector member of the 2S family to give the physical D * s1 (2700) − and D * s1 (2860) − states. The discovery of the D * s3 (2860) − resonance represents the first observation of a heavy flavoured spin-3 particle, and the first time that a spin-3 state is seen to be produced in B decays. This discovery demonstrates that 1D charm resonances can be investigated experimentally, and therefore opens a new window for potential studies of the spectroscopy of heavy flavoured mesons. and GENCAT (Spain), Royal Society and Royal Commission for the Exhibition of 1851 (United Kingdom).  [40] M. Gronau, Improving bounds on γ in B ± → DK ± and B ±,0 → DX ±,0 s , Phys. Lett. B557 (2003) 198, arXiv:hep-ph/0211282.

A Results for interference fit fractions
The central values of the interference fit fractions are given in Table 17. The statistical, experimental systematic and model uncertainties on these quantities are given in Tables 18, 19 and 20, respectively.   Table 5.  Table 5.    Table 8.  Table 11.