Study of $D^{(*)+}_{sJ}$ mesons decaying to $D^{*+} K^0_{\rm S}$ and $D^{*0} K^+$ final states

A search is performed for $D^{(*)+}_{sJ}$ mesons in the reactions $pp \to D^{*+} K^0_{\rm S} X$ and $pp \to D^{*0} K^+ X$ using data collected at centre-of-mass energies of 7 and 8 TeV with the LHCb detector. For the $D^{*+} K^0_{\rm S}$ final state, the decays $D^{*+} \to D^0 \pi^+$ with $D^0 \to K^- \pi^+$ and $D^0 \to K^- \pi^+ \pi^+ \pi^-$ are used. For $D^{*0} K^+$, the decay $D^{*0} \to D^0 \pi^0$ with $D^0 \to K^- \pi^+$ is used. A prominent $D_{s1}(2536)^+$ signal is observed in both $D^{*+} K^0_{\rm S}$ and $D^{*0} K^+$ final states. The resonances $D^*_{s1}(2700)^+$ and $D^*_{s3}(2860)^+$ are also observed, yielding information on their properties, including spin-parity assignments. The decay $D^*_{s2}(2573)^+ \to D^{*+} K^0_{\rm S}$ is observed for the first time, at a significance of 6.9 $\sigma$, and its branching fraction relative to the $D^*_{s2}(2573)^+ \to D^+ K^0_{\rm S}$ decay mode is measured.

mesons in the reactions pp → D * + K 0 S X and pp → D * 0 K + X using data collected at centre-of-mass energies of 7 and 8 TeV with the LHCb detector. For the D * + K 0 S final state, the decays D * + → D 0 π + with D 0 → K − π + and D 0 → K − π + π + π − are used. For D * 0 K + , the decay D * 0 → D 0 π 0 with D 0 → K − π + is used. A prominent D s1 (2536) + signal is observed in both D * + K 0 S and D * 0 K + final states. The resonances D * s1 (2700) + and D * s3 (2860) + are also observed, yielding information on their properties, including spin-parity assignments. The decay D * s2 (2573) + → D * + K 0 S is observed for the first time, at a significance of 6.9 σ, and its branching fraction relative to the D * s2 (2573) + → D + K 0 S decay mode is measured.

Introduction
The discovery by the BaBar collaboration of a narrow state D * s0 (2317) + in the decay to D + s π 0 [1], and the subsequent discovery of a second narrow particle, D s1 (2460) + in the decay to D * + s π 0 [2][3][4], raised considerable interest in the spectroscopy of heavy mesons. 1 These discoveries were a surprise because quark model calculations based on heavy quark effective theory (HQET) [5] predicted the masses of these resonances to be above the DK and D * K thresholds, respectively. 2 Consequently their widths were expected to be very large, as for the corresponding J P = 0 + and J P = 1 + resonances in the D J spectrum.
The D + sJ mesons are expected to decay into the DK and D * K final states if they are above threshold. The BaBar collaboration has explored the DK and D * K mass spectra [6,7] observing two states, D * sJ (2700) + and D * sJ (2860) + , both decaying to DK and D * K with a natural parity (NP) assignment. 3 A third structure, D sJ (3040) + , is observed only in the D * K decay mode with a preferred unnatural parity (UP) assignment. The D * sJ (2700) + resonance was also observed by the Belle and BaBar collaborations in a study of B decays to DDK [8,9]. Both collaborations obtain a spin-parity assignment J P = 1 − for this state, and so it is labelled as D * s1 (2700) + . The LHCb experiment has performed studies of the DK final states in the inclusive process, pp → DKX [10], and in the Dalitz plot analysis of B 0 s → D 0 K − π + decays [11,12]. In the inclusive analysis, the D * s1 (2700) + and D * sJ (2860) + are observed with large statistical significance and their properties are found to be in agreement with previous measurements. In the exclusive Dalitz plot analysis of the B 0 s → D 0 K − π + decays, the D 0 K − mass spectrum shows a complex resonant structure in the 2860 MeV mass region. 4 This is described by a superposition of a broad J P = 1 − resonance and a narrow J P = 3 − resonance with no evidence for the production of D * s1 (2700) − . Since the narrow structure at 2860 MeV seen in inclusive DK and D * K analyses could contain contributions from various resonances with different spins, it is labelled as D * sJ (2860) + . In references [13][14][15][16][17][18] attempts are made to identify these states within the quark model and in Ref. [19] within molecular models. The expected spectrum for D + s mesons has recently been recomputed in Refs. [20,21]. In particular, Ref. [20] points out that six states are expected in the mass region between 2.7 and 3.0 GeV. To date, evidence has been found for three of the states; hence finding the rest would provide an important test of these models. In this paper we study the D * + K 0 S and D * 0 K + systems using pp collision data, collected at centre-of-mass energies of 7 and 8 TeV with the LHCb detector.

Detector and simulation
The LHCb detector [22,23] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV. The minimum distance of a track to a primary vertex (PV), the impact parameter, is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors (RICH). Photons, electrons and hadrons 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. The online event selection is performed by a trigger, which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction.
In the simulation, pp collisions are generated using Pythia [24] with a specific LHCb configuration [25]. Decays of hadronic particles are described by EvtGen [26], in which final-state radiation is generated using Photos [27]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [28] as described in Ref. [29]. We also make use of simple generator-level simulations [30] to study kinematic effects.

Event selection
We search for D ( * )+ sJ mesons using the inclusive reactions and where X represents a system composed of any collection of charged and neutral particles. Use is made of both 7 and 8 TeV data for reaction (1), corresponding to an integrated luminosity of 3 fb −1 , and 8 TeV data only for reaction (2) which corresponds to an integrated luminosity of 2 fb −1 . The charmed mesons in the final state are reconstructed in the decay modes D * + → D 0 π + , with D 0 → K − π + and D 0 → K − π + π + π − , and D * 0 → D 0 π 0 , with D 0 → K − π + and π 0 → γγ. The K 0 S mesons are reconstructed in their K 0 S → π + π − decay mode. Because of their long lifetime, K 0 S mesons may decay inside or outside the vertex detector. Candidate K 0 S mesons that are reconstructed using vertex detector information are referred to as "long" while those reconstructed without vertex detector information are called "downstream". Those that decay within the vertex detector acceptance have a mass resolution about half as large as those that decay outside of its acceptance. Reaction (1) with D 0 → K − π + serves as the primary channel for studying the D ( * )+ sJ resonance structures and their parameters, while reaction (1) with D 0 → K − π + π + π − and reaction (2) are used for cross-checks and to confirm the observed signatures.
Charged tracks are required to have good track fit quality, momentum p > 3 GeV and p T > 250 MeV. These conditions are relaxed to p > 1 GeV and p T > 150 MeV for the "soft" pion originating directly from the D * + decay. In the reconstruction of the D 0 candidates we remove candidate tracks pointing to a PV, using an impact parameter requirement. All tracks used to reconstruct the D mesons are required to be consistent with forming a common vertex and the D meson candidate must be consistent with being produced at a PV. The D * + and K 0 S , and similarly the D 0 and K + candidates, are fitted to a common vertex, for which a good quality fit is required. The purity of the charmed meson sample is enhanced by requiring the decay products to be identified by the particle identification system, using the difference in the log-likelihood between the kaon and pion hypotheses ∆ ln L Kπ [31]. We impose a tight requirement of ∆ ln L Kπ > 3 for kaon tracks and a loose requirement of ∆ ln L Kπ < 10 for pions. The overlap region in the particle identification definition of a kaon and a pion is small and does not affect the measured yields, given the small number of multiple candidates per event.
Candidate D 0 mesons are required to be within ±2.5 σ of the fitted D 0 mass where the mass resolution σ is 8.3 MeV. The D 0 π + invariant mass is computed as where m D 0 is the world average value of the D 0 mass [32]. For the channel D 0 → K − π + π + π − , the invariant mass m(D 0 π + ) is defined similarly. Figure 1 shows the D 0 π + invariant mass spectrum for (a) D 0 → K − π + and (b) D 0 → K − π + π + π − . Clean D * + signals for both D 0 decay modes are observed. We fit the mass spectra using the sum of a Gaussian function for the signal and a second-order polynomial for the background. The signal regions are defined to be within ±2.5 σ of the peak values, where σ = 0.7 MeV for both channels.
The π + π − mass spectra for the two K 0 S types, the long K 0 S and downstream K 0 S , are shown in Figs. 1(c) and 1(d) and are fitted using the same model as for the D 0 π + invariant masses. The signal regions are similarly defined within ±2.5 σ of the peak, with σ = 4.1 MeV and 8.7 MeV for long and downstream K 0 S , respectively. The π 0 candidates are obtained by kinematically fitting to a π 0 hypothesis each pair of photon candidates with energy greater than 600 MeV, with the diphoton mass constrained to the nominal π 0 mass [32]. Candidate D * 0 mesons are formed by combining D 0 → K − π + decays with all π 0 candidates in the event that have p T > 450 MeV. The resulting D * 0 candidate is required to have p T > 6000 MeV. Figure 2 shows the ∆m(D 0 π 0 ) = m(K − π + π 0 ) − m(K − π + ) distribution, where a clear D * 0 signal can be seen. The mass  spectrum is fitted using for background the threshold function where in this case m = ∆m(D 0 π 0 ), m th is the ∆m(D 0 π 0 ) threshold mass and α, β and γ are free parameters. In Eq. (4) P (m) is the center of mass momentum of the two-body decay of a particle of mass m into two particles with masses m 1 and m 2 , The function B(m) gives the correct behaviour of the fit at threshold. The D * 0 signal is modelled using the sum of two Gaussian functions. We select the candidates in the ±2 σ window around the peak, where σ = 1.72 MeV is the width of the dominant Gaussian fitting function, and we form D * K pairings by combining D * + and K 0 S candidates for reaction (1), and D * 0 and K + candidates for reaction (2). To suppress the large combinatorial background, a set of additional criteria is applied. We define θ as the angle between the momentum direction of the kaon in the D * K rest frame and the momentum direction of the D * K system in the laboratory frame. Whereas the signal events are expected to be symmetrically distributed in the variable cos θ, after correcting for efficiency, more than 90% of the combinatorial background is found in the negative cos θ region. The cos θ requirements are optimized using the D * s1 (2700) + signal, an established resonance. We fit the D * K mass spectra (using the model described below) with different cos θ selections and obtain the yields for D * s1 (2700) + signal (N S ) and background events (N B ) in the D * s1 (2700) + signal region (defined in the window |m(D * K) − m(D * s1 (2700) + )| < Γ(D * s1 (2700) + )/2). We compute the signal significance S = N S / √ N S + N B and signal purity P = N S /(N S + N B ) and find that the requirements cos θ > 0 (for and cos θ > −0.1 (for D * 0 K + ) each provide a good compromise between significance and purity in the respective channel. With the same method it is also found that it is optimal to require p T > 4000 MeV for all three final states. Simulations show that the mass resolution is much smaller than the natural widths of the resonances.
The analysis of the D * K system, with D * → Dπ, is a three-body decay and therefore allows a spin analysis of the produced resonances and a separation of the different spinparity components. We define the helicity angle θ H as the angle between the K 0 S and the π + from the D * + decay, in the rest frame of the D * + K 0 S system. Simulated events are used to determine the efficiency as a function of cos θ H , which is found to be uniform only for the D * + K 0 S candidates formed from the downstream K 0 S sample. Therefore, for studying the angular distributions we do not use the long K 0 S sample, which removes approximately 30% of the data.

Mass spectra
In order to improve the mass resolution on the D * K mass spectra, we compute the D * + , K 0 S and D * 0 energies using the world average mass measurements [32]. The D * + K 0 S mass spectrum for D 0 → K − π + is shown in Fig. 3 and contains 5.72×10 5 combinations. We observe a strong D s1 (2536) + signal and weaker resonant contributions due to D * s2 (2573) + , D * s1 (2700) + , and D * sJ (2860) + states. The D * s2 (2573) + decay to D * + K 0 S is observed for the first time. A binned χ 2 fit to the mass spectrum is performed in which the narrow D s1 (2536) + is described by a Gaussian function with free parameters. Other resonances are described by relativistic Breit-Wigner (BW) functions (in D-, P -and F -wave for D * s2 (2573) + , D * s1 (2700) + , and D * s3 (2860) + respectively). Using the definition of the center-of-mass momentum P (m) given in Eq. (5), we parameterize the BW function for a resonance of mass M as where and are the Blatt-Weisskopf form factors [33]. No dependence of the resonance parameters on the Blatt-Weisskopf radius R is found and it is therefore fixed to 2.5 GeV −1 . The quantity L is the angular momentum between the two decay fragments: L = 1 for P -wave, L = 2 for D-wave and L = 3 for F -wave resonances. The D sJ (3040) + resonance is described by a nonrelativistic BW function multiplied by P (m). The D * s2 (2573) + parameters are fixed to the values obtained in the fit to the DK mass spectra [10]. The background is described by an empirical model [34], where P (m) is described in Eq. (5) and m 0 , a i=1,2 and b i=0,1,2 are free parameters. In Eq.(9) we impose continuity to B(m) and to its first derivative at the mass m 0 and therefore the number of free parameters is reduced by two. Resonances are included sequentially in order to test the χ 2 improvement when a new contribution is added. A better fit is obtained if a broad resonance in the 3000 MeV mass region is included. We find strong correlation between the parameters of this structure and the background and therefore we add the D sJ (3040) + resonance in the fit with parameters fixed to the values obtained by BaBar [7]. 5 We also study the D * + K 0 S in the D * + sideband region, defined as 2014.0 < m(D 0 π + ) < 2018.1 MeV. A smooth mass spectrum is obtained, well fitted by the above background model with no evidence for additional structures. Table 1(a) gives the resulting D * s1 (2700) + and D * sJ (2860) + fitted parameters. Statistical significances are computed as S = ∆χ 2 , where ∆χ 2 is the difference in χ 2 between fits with the resonance included and excluded from the fitting model. Large significances for D * s1 (2700) + and D * sJ (2860) + are obtained, especially for the D 0 → K − π + decay mode. The significance of the D sJ (3040) + enhancement is 2.4 σ.
A search is performed for the D * s1 (2860) + resonance previously observed in the B 0 s → D 0 K − π + Dalitz plot analysis [11,12]. We first introduce in the fit an incoherent BW function with parameters free to vary within their statistical uncertainties around the reported values in Ref. [11], but the fit returns a negligible contribution for this state. Since  Table 1: Results from the fits to the D * + K 0 S and D * 0 K + mass spectra. Resonances parameters are expressed in MeV. When two uncertainties are presented, the first is statistical and the second systematic. The symbol ndf indicates the number of degrees of freedom. two J P = 1 − overlapping resonances may be present in the mass spectrum, interference is allowed between the D * s1 (2860) + and the D * s1 (2700) + resonance by including the amplitude where c and φ are free parameters. In this fit we also add the D * s3 (2860) + resonance with parameters fixed to those from Refs. [11,12] and the D * s1 (2700) + with parameters fixed to those from the DK analysis [10]. The resulting fit quality is similar to that obtained without the presence of the D * s1 (2860) + resonance (χ 2 /ndf = 92/103). However it is found that the D * s1 (2860) + is accommodated by the fit with strong destructive interference. We conclude that the data are not sensitive to the D * s1 (2860) + resonance. Systematic uncertainties on the resonance parameters are computed as quadratic sums of the differences between the nominal fit and fits in which the following changes are made.
• The alternative background function Eq. (4) is used.
• The fit bias is evaluated by generating and fitting pseudoexperiments obtained using the parameters from the best fit. The deviations of the mean values of the distributions from the generated ones are taken as systematic uncertainties. • The parameters of the D sJ (3040) + state, fixed to the values of Ref. [7] in all the fits, have been varied according to their total uncertainties.
• From the study of high-statistics control samples, a systematic uncertainty of 0.0015 Q on the mass scale is added, where Q is the Q-value involved in the resonance decay.
• The fitting model that includes the D * s1 (2860) + resonance is tested with D * s3 (2860) + parameters fixed and the D * s1 (2700) + parameters left free. The different contributions to the systematic uncertainties are summed in quadrature and are summarized in Table 2. It can be noted that, combining statistical and systematic uncertainties, the resulting D * s1 (2700) + mass is about 3 σ higher than previous measurements while the D * sJ (2860) + parameters are consistent with those of the D * s3 (2860) + resonance [32].
The angular distributions are expected to be proportional to sin 2 θ H for NP resonances and proportional to 1 + h cos 2 θ H for UP resonances, where h is a free parameter. The D * K decay is forbidden for a J P = 0 + resonance. Therefore the selection of candidates in different ranges of cos θ H can enhance or suppress different spin-parity contributions. We separate the D * + K 0 S data into two different categories, the NP sample, obtained with the selection | cos θ H | < 0.5 and the UP sample, with the selection | cos θ H | > 0.5.
The D * + K 0 S mass spectra for the NP sample is shown in Fig. 4(a), while the corresponding mass spectrum for the UP sample is shown in Fig. 4(b). Most resonant structures are in the NP sample. An enhancement in the 2860 MeV mass region in Fig. 4(b) indicates the possible presence of additional UP contributions. The fitted parameters are given in Tables 1(b) and 1(c). Figure 5(a) shows the D * + K 0 S mass spectrum for D 0 → K − π + π + π − , which contains 3.92×10 4 combinations. Similar resonant structures to those seen for the D * + K 0 S final state with D 0 → K − π + are observed, albeit at lower significance. Table 1(d) provides the fitted resonance parameters. Due to the limited data samples, some parameters have been fixed to the values obtained from the fit to the D * + K 0 S sample with D 0 → K − π + . The mass values are found to be consistent with the results from the other measurements.
The D * 0 K + mass spectrum is affected by a high level of combinatorial background, mostly due to the D * 0 reconstruction (see Fig. 2). As observed previously, the D * K mass spectra are dominated by NP resonances and therefore in Fig. 5(b) we show the D * 0 K + mass spectrum for the NP sample. The mass spectrum contains 2.53×10 4 combinations. We observe similar resonant structures as seen in the study of the D * + K 0 S mass spectra. The fitted resonance parameters are given in Table 1(e); mass values are consistent with the results from the fits to the other mass spectra. We do not have the sensitivity to the parameters of the D * s1 (2700) + and D * s3 (2860) + resonances in the fits to the D * + K 0 S , D 0 → K − π + π + π − , and D * 0 K + mass spectra due to the low statistical significance of the signals.
5 Measurement of the branching fraction of the decay D * s2 (2573) + → D * + K 0 S We measure the branching fraction of the decay D * s2 (2573) + → D * + K 0 S , D 0 → K − π + relative to that of the decay D * s2 (2573) + → D + K 0 S . For this purpose the D + K 0 S mass spectrum from Ref. [10], collected at 7 TeV with an integrated luminosity of 1 fb −1 , is re-fitted. In this study both long and downstream K 0 S candidate types are used. The final states D * + K 0 S , with D 0 → K − π + π + π − and D * 0 K + are used as cross checks and to aid in determining the significance of the signal. Figure 6 shows the D + K 0 S mass spectrum from Ref. [10] along with the results of the fit described below. A narrow structure is seen near threshold, due to the cross-feed from the decay D s1 (2536) where the π 0 /γ are not reconstructed. In the higher mass region, a strong D * s2 (2573) +  signal and a weak signal due to the D * s1 (2700) + resonance are observed. Due to the difficulty of controlling the systematic uncertainties related to the determination of the relative efficiencies of the D * + K 0 S and D + K 0 S final states, we normalize the two mass spectra using the D s1 (2536) + signal which is observed as a peak in the D * + K 0 S and as cross-feed in the D + K 0 S final states. The D * s2 (2573) + resonance is a well known NP J P = 2 + state. To enhance the signal to background ratio, we plot in Fig. 7 the D * K mass spectra for the NP sample of the three final states. All three distributions show a strong D s1 (2536) + signal and an enhancement at the D * s2 (2573) + mass. The D + K 0 S mass spectrum and the three D * K mass spectra are fitted using the background function where P (m) is given in Eq. (5) and β and γ are free parameters. The D s1 (2536) + crossfeed into D + K 0 S is modelled using the sum of two Gaussian functions with the same mean, and the D * s2 (2573) + resonance is modelled as a relativistic BW function convolved with a Gaussian function describing the experimental resolution (σ = 3.5 MeV). Since the intrinsic width of the D s1 (2536) + state in the D * K spectra is much smaller than the experimental resolution, the D s1 (2536) + is modelled using the sum of two Gaussian functions with the same mean. We obtain m(D s1 (2536) + ) = 2535.00 ± 0.01 MeV, in good agreement with the PDG average. The D * s2 (2573) + resonance is modelled as a relativistic BW function convolved with the experimental resolution (σ = 2.5 MeV for D * s2 (2573) + → D * + K 0 S , D 0 → K − π + ) taking the mass value as a free parameter and with the full width constrained to the value obtained from the fit to the D + K 0 S mass spectrum (Γ = 17.5 ± 0.4 MeV). Table 3 summarizes the fit results. We note the large statistical significance of the Table 3: Results from the fits to the D + K 0 S and D * + K 0 S mass spectra for the evaluation of the D * s2 (2573) + → D * + K 0 S relative branching fraction.  Figure 7: Mass spectra, in the D s1 (2536) + mass region, for the NP sample of (a) D * + K 0 S with D 0 → K − π + , (b) D * + K 0 S with D 0 → K − π + π + π − , and (c) D * 0 K + final states. The full (red) lines describe the fitting function. The dashed lines show the fitted background and the dotted lines the D * s2 (2573) + contributions. The insets display the D * K mass spectra after subtracting the fitted background. D * s2 (2573) + in the D * + K 0 S final states, especially for the sample with D 0 → K − π + . Consistency is found, within the uncertainties, in the D * s2 (2573) + mass measurements for the different final states. We therefore identify the observed structure as the first observation of the D * s2 (2573) + → D * + K 0 S decay. The relative branching fraction

Final state
is determined using the results of fits to the D * s2 (2573) + → D * + K 0 S , D 0 → Kπ data shown in Fig. 7(a) and the D * s2 (2573) + → D + K 0 S data shown in Fig. 6, summarized in Table 3. Using the D * + K 0 S final state, we verify that the D s1 (2536) + cross-feed into the D 0 K 0 S mass spectrum, when the pion from the D * + decay is ignored, contains all the D s1 (2536) + signal. Similarly, using the D * 0 K + data, we ignore the π 0 from the D * 0 decay and plot the D 0 K + mass spectrum. Also in this case, it is found that the D s1 (2536) + cross-feed contains all the decays in the D s1 (2536) + signal region. It is assumed that the D s1 (2536) + meson decay to D * K is dominant. We test this hypothesis by studying the D 0 π 0 K + mass spectrum and find that no D s1 (2536) + signal is present outside the D * 0 → D 0 π 0 signal region. Indicating explicitly in brackets the D * + decay modes, we define and where N indicates the yields and D s1 (2536) + → (D + K 0 S ) f indicates the cross-feed from D s1 (2536) + → D * + K 0 S where D * + → D + (π 0 /γ) and the π 0 /γ are undetected. We measure the D * s2 (2573) + relative branching ratio as where indicates the efficiency for each final state. The ratio B D , defined below, is taken from Ref. [32], where D + (π 0 /γ) indicates both D + π 0 and D + γ decays and f NP is defined below.
In the evaluation of the D * s2 (2573) + relative branching fraction, we make use of the D * + K 0 S NP sample. This selection is used to improve the signal to background ratio for the D * s2 (2573) + resonance in the D * + K 0 S final state. We also fit the D * + K 0 S mass spectrum using the full dataset and we report the D s1 (2536) + yield indicated as Total in Table 3. In Eq. (16) the total D s1 (2536) + yield is used because of the unnatural parity of this state, and this requires a correction to the D * s2 (2573) + yield for the effects of the NP sample selection. The angular distribution for a NP resonance is expected to be proportional to sin 2 θ H and therefore the requirement | cos θ H | < 0.5 selects 69% of the candidates. This correction in Eq. (16) is included through the factor f NP = 1.45.
In Eq. (16) it can be noted that the efficiencies (D * s2 (2573) + → (D 0 π + )K 0 S ) and (D s1 (2536) + → (D 0 π + )K 0 S ) involve the same final state. They are determined from simulation and are found to be the same within uncertainties. Similarly, the efficiencies (D s1 (2536) + → (D + K 0 S ) f ) and (D * s2 (2573) + → (D + K 0 S )) are also found to be the same within uncertainties. Therefore, the efficiency ratios are set to unity. Table 4 summarizes the measurements used to estimate the D * s2 (2573) + relative branching fraction. We obtain The systematic uncertainty on the D * s2 (2573) + relative branching fraction is computed as the quadratic sum of the differences between the reference values and those obtained when the following changes are made.
• The D + K 0 S data are collected at 7 TeV, while the D * K 0 S data include 7 TeV and 8 TeV data samples. We compute systematic uncertainties on the R 1 and R 2 ratios using the D * K 0 S at 7 TeV only and include the deviation in the systematic uncertainty.
• The uncertainty on the B D parameter is propagated as a systematic uncertainty.
• Using simulation, we compute efficiency distributions as functions of m(D * + K 0 S ) and m(D + K 0 S ) and observe that they have weak variations in the regions used to evaluate the relative branching fraction. We assign a 10% systematic uncertainty to cover the assumptions that the efficiencies as functions of m(D * + K 0 S ) and m(D + K 0 S ) in Eq. (16) are the same.
• We vary the shape of the background function using Eq. (4) in the fits to the D * + K 0 S and D + K 0 S mass spectra and obtain new estimates for the resonance yields. We also remove the convolution with the resolution function or replace the relativistic BW functions with simple BW functions and include an additional Gaussian function to describe the D s1 (2536) + signal.
• We vary the D * s2 (2573) + width by its statistical uncertainty (0.4 MeV) simultaneously in the fits to the D + K 0 S and D * K 0 S mass spectra. The contributions to the systematic uncertainty are summarized in Table 5 with the dominant component arising from the use of different datasets collected at different centre-of-mass energies.
We also perform a new estimate of the D * s2 (2573) + significance in the D * + K 0 S final state by combining in quadrature the statistical and systematic uncertainties on the yield   (see Table 4) and obtain S = N signal /σ tot = 6.9, where σ tot is the total error. This estimate is in good agreement with that reported in Table 3. 6 Spin-parity analysis of the D * + K 0 S system We obtain information on the spin-parity of the states observed in the D * + K 0 S mass spectrum. The data for D 0 → K − π + are first divided into five equally spaced bins in cos θ H . The five mass spectra in the D * + K 0 S threshold region (m(D * + K 0 S ) < 2650 MeV) are fitted using the model described in Sec. 5 with fixed D s1 (2536) + and D * s2 (2573) + resonance parameters, to obtain the signal yields as functions of cos θ H for each resonance.
As stated previously, we determine from simulations that the efficiency as a function of cos θ H is consistent with being uniform; therefore we plot uncorrected angular distributions.

Resonance
J P Function The resulting distributions for D s1 (2536) + and D * s2 (2573) + are shown in Fig. 8(a) and Fig. 8(b), and are fitted using the functions described in Table 6. A good description of the data is obtained in terms of the expected angular distributions for J P = 1 + and J P = 2 + resonances. We note that the shape of the D s1 (2536) + angular distribution is in agreement with that measured in Ref. [35].
The D * + K 0 S data, with D 0 → K − π + , are then divided into eight equally spaced bins in cos θ H . The mass spectra are fitted (for m(D * + K 0 S ) < 3400 MeV) with the model described in Sec. 4 with fixed resonance parameters, to obtain the yields as functions of cos θ H for each resonance. The resulting distributions are shown in Fig. 9 and details of the fit results are given in Table 6.
We observe that the D * s1 (2700) + state is reasonably well described by the expected NP function (χ 2 /ndf = 11.4/7 with p-value 12.2%). The fit to the D * s3 (2860) + angular distribution has a slightly lower p-value (6.3%). Reference [18] suggests the possibility of the presence of UP state contributions in this mass range, which cannot be excluded in this fit: there is evidence for the presence of a small signal in the 2860 MeV mass region for the UP sample shown in Fig. 4(b). The consistency with the NP assignment confirms the presence of the decay D * s3 (2860) + → D * + K 0 S . We also show in Fig. 9(c) the cos θ H distribution for the enhancement at the D sJ (3040) + position and find it consistent with a UP assignment.

Summary
A study of the resonant structures in the D * + K 0 S and D * 0 K + systems is performed using pp collision data, collected at centre-of-mass energies of 7 and 8 TeV with the LHCb detector. For the D * + K 0 S final state, the decay chains D * + → D 0 π + with D 0 → K − π + and D 0 → K − π + π + π − are used, with an integrated luminosity of 3.0 fb −1 . For D * 0 K + , the decay chain D * 0 → D 0 π 0 , D 0 → K − π + is used, with an integrated luminosity of 2.0 fb −1 .
A prominent D s1 (2536) + resonance is observed in both D * + K 0 S and D * 0 K + final states. Resonances D * s1 (2700) + and D * s3 Study of the angular distributions supports natural parity assignments for both resonances, although the presence of an additional unnatural parity contribution in the 2860 MeV mass range cannot be excluded. The data are not sensitive to the presence of an additional D * s1 (2860) + resonance. The D * s2 (2573) + decay to D * + K 0 S is also observed for the first time, at a significance of  6.9 σ, with a branching fraction relative to the D + K 0 S decay mode of B(D * s2 (2573) + → D * + K 0 S ) B(D * s2 (2573) + → D + K 0 S ) = 0.044 ± 0.005 (stat) ± 0.011 (syst). (19) This measurement is in agreement with expectations from recent calculations of the charm and charm-strange mesons spectra [21] which predict a value of 0.058 for this ratio. A spinparity analysis of the decay D * s2 (2573) + → D * + K 0 S supports the natural parity assignment. The data also show weak evidence for further structure in the region around 3040 MeV consistent with contributions from unnatural parity states.