Observation of a near-threshold structure in the $K^+$ recoil-mass spectra in $e^+e^-\to K^+ (D_s^- D^{*0} + D^{*-}_s D^0)$

We report a study of the processes of $e^+e^-\to K^+ (D_s^- D^{*0} + D^{*-}_s D^0)$ based on $e^+e^-$ annihilation samples collected with the BESIII detector operating at BEPCII at five center-of-mass energies ranging from 4.628 to 4.698 GeV with a total integrated luminosity of 3.7 fb$^{-1}$. An excess over the known contributions of the conventional charmed mesons is observed near the $D_s^- D^{*0}$ and $D^{*-}_s D^0$ mass thresholds in the $K^{+}$ recoil-mass spectrum for events collected at $\sqrt{s}=4.681$ GeV. The structure matches a mass-dependent-width Breit-Wigner line shape, whose pole mass and width are determined as $(3982.5^{+1.8}_{-2.6}\pm2.1)$ MeV/$c^2$ and $(12.8^{+5.3}_{-4.4}\pm3.0)$ MeV, respectively. The first uncertainties are statistical and the second are systematic. The significance of the resonance hypothesis is estimated to be 5.3 $\sigma$ over the pure contributions from the conventional charmed mesons. This is the first candidate of the charged hidden-charm tetraquark with strangeness, decaying into $D_s^- D^{*0}$ and $D^{*-}_s D^0$. However, the genuine properties of the excess need further exploration with more statistics.

We report a study of the processes of e + e − → K + (D − s D * 0 + D * − s D 0 ) based on e + e − annihilation samples collected with the BESIII detector operating at BEPCII at five center-of-mass energies ranging from 4.628 to 4.698 GeV with a total integrated luminosity of 3.7 fb −1 . An excess over the known contributions of the conventional charmed mesons is observed near the D − s D * 0 and D * − s D 0 mass thresholds in the K + recoil-mass spectrum for events collected at √ s = 4.681 GeV. The structure matches a mass-dependent-width Breit-Wigner line shape, whose pole mass and width are determined as (3982.5 +1.8 −2.6 ± 2.1) MeV/c 2 and (12.8 +5.3 −4.4 ± 3.0) MeV, respectively. The first uncertainties are statistical and the second are systematic. The significance of the resonance hypothesis is estimated to be 5.3 σ over the pure contributions from the conventional charmed mesons. This is the first candidate of the charged hidden-charm tetraquark with strangeness, decaying into D − s D * 0 and D * − s D 0 . However, the genuine properties of the excess need further exploration with more statistics.
The existence of a Z cs state with a mass lying around the D − s D * 0 and D * − s D 0 thresholds has been predicted in several theoretical models, including tetraquark scenarios [18,19], the D sD * molecular model [20], the hadroquarkonium model [19], and in the initial-single-chiralparticle-emission mechanism [21]. Like the Z c states, the decay rate of the Z cs to open-charm final states is expected to be larger than the decay rate to charmonium final states [5]. Hence, one promising method to search for the Z cs state is through its decays to D − s D * 0 and D * − s D 0 . In this Letter, we report on a study of the process e + e − → K + (D − s D * 0 + D * − s D 0 ) at center-of-mass energies √ s = 4.628, 4.641, 4.661, 4.681, and 4.698 GeV. The data samples have a total integrated luminosity of 3.7 fb −1 and were accumulated by the BESIII detector at the BEPCII collider. Details about BEPCII and BESIII can be found in Refs. [22][23][24]. To improve the signalselection efficiency, a partial-reconstruction technique is implemented in which only the charged K + (the bachelor K + ) and the D − s are reconstructed. Here and elsewhere, charge-conjugated modes are always implied, unless explicitly stated otherwise. To improve the signal purity, we only reconstruct the decays D − s → K + K − π − and K 0 S K − , which have large branching fractions (BF). By reconstructing the D − s meson, the flavors of the missing D 0 and the bachelor K + are fixed. We observe an enhancement near the D − s D * 0 and D * − s D 0 mass thresholds in the K + recoil-mass spectrum for events collected at √ s = 4.681 GeV and carry out a fit to the en-hancement with a possible new Z cs candidate, denoted as Z cs (3985) − , in the K + recoil-mass spectra at different energy points. Simulated data samples produced with a geant4based [25] Monte Carlo (MC) simulation package, which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate background. The simulation includes the beam-energy spread and initialstate radiation (ISR) in the e + e − annihilations modeled with the generator kkmc [26]. Inclusive MC samples consist of the production of open-charm processes, the ISR production of vector charmonium(-like) states, and the continuum processes incorporated in kkmc. The known decay modes are modeled with evtgen [27] using branching fractions taken from the Particle Data Group (PDG) [1], and the remaining unknown decays from the charmonium states are simulated with lundcharm [29]. The final-state radiation from charged finalstate particles is incorporated with the photos package [30]. For the three-body non-resonant (NR) signal process, e + e − → K + (D − s D * 0 + D * − s D 0 ), the final-state particles are simulated assuming non-resonant production [31]. For the simulation of the Z cs (3985) − signal process, e + e − → K + Z cs (3985) − , we let the Z cs (3985) − decay into the D − s D * 0 and D * − s D 0 final states with equal rates. The Z cs (3985) − state is assigned a spin-parity of 1 + , as the corresponding production and subsequent decay processes are both in the most favored S-wave. However, other spin-parity assignments are allowed, and these are tested as systematic variations.
To identify the processes e + e − → K + (D − s D * 0 + D * − s D 0 ), we reconstruct combinations of the bachelor K + and the decays D − s → K + K − π − or K 0 S K − . Data taken at all five center-of-mass energy points are analyzed using the same procedure, but two-thirds of the data set at √ s = 4.681 GeV was kept blinded until after the analysis strategy was established and validated [32]. We first select events with at least four charged tracks. Charged tracks detected in the multilayer drift chamber (MDC) are required to be within a polar-angle range of |cosθ| < 0.93, where θ is defined with respect to the zaxis of the MDC. The distance of closest approach to the interaction point (IP) must be |V z | < 10 cm along the zaxis and |V xy | < 1 cm in the perpendicular plane, except for those pions used for the reconstruction of K 0 S decays. Particle identification (PID) for charged tracks combines measurements of the ionization energy loss in the MDC and the flight time recorded by the time-of-flight counter to form likelihoods L(h) (h = K, π) for each hadron h hypothesis. The charged K(π) particle is identified by requiring L(K) > L(π) (L(π) > L(K)). For the decay D − s → K + K − π − , to improve the signal purity, we only retain the D − s candidates within the Dalitz plot regions consistent with D − s → φπ − or D − s → K * (892) 0 K − decays by requiring that the invariant masses satisfy either M (K + K − ) < 1.05 GeV/c 2 or 0.850 < M (K + π − ) < 0.930 GeV/c 2 . For the decay D − s → K 0 S K − , K 0 S can- didates are reconstructed from two oppositely charged tracks satisfying |V z | < 20 cm. The two charged tracks are assigned as π + π − without imposing further PID criteria. They are constrained to originate from a common vertex and are required to have an invariant mass within 0.485 < M (π + π − ) < 0.511 GeV/c 2 . Furthermore, the decay length with respect to the IP of the K 0 S candidate is required to be greater than twice its resolution. s candidate, thus providing improved resolution compared to RM (K + D − s ). A clear peak is seen in this distribution at the nominal D * 0 mass, which corresponds to the final state K + D − s D * 0 . There is also a contribution from K + D * − s D 0 , which appears as a broader structure beneath the K + D − s D * 0 signal. Therefore, we require RM to be in the interval (1.990, 2.027) GeV/c 2 to isolate the signal candidates of both signal processes.
To estimate the shape of combinatorial background, we use wrong-sign (WS) combinations of D − s and K − candidates, rather than the right-sign D − s and K + candidates. The WS K − D − s recoil-mass distribution, scaled by a factor of 1.18, agrees with the data distribution in the sideband regions, (1.91, 1.95) GeV/c 2 and (2.08, 2.11) GeV/c 2 , as shown in Fig. 2. The number of background events within the signal region is estimated to be 282.6 ± 12.0 by a fit to the sideband data with a linear function, whose slope is determined from the WS data. In addition, the WS events are used to represent the combinatorial-background distribution of the recoil mass of the bachelor K + . This technique has been used previously in the observation of the Z c (4025) + at BESIII [10].
We validate the use of the WS data-driven background modeling of both the RM (K + D − s ) and RM (K + ) spectra by comparing the corresponding distributions between WS combinations and background-only contributions. Furthermore, the RM (K + ) distribution of the events in the sideband regions in Fig. 2 agrees well with that of the corresponding WS data. Figure 3(a) shows the RM (K + ) distribution for events at √ s = 4.681 GeV; an enhancement is evident in the region RM (K + ) < 4 GeV/c 2 compared to the expectation from the WS events. This is clearly illustrated in the RM (K + ) distribution in data with subtraction of the WS component in Fig. 4. The enhancement cannot be attributed to the NR signal processes e + e − → K + (D − s D * 0 +D * − s D 0 ). To understand potential contributions from the processes e + e − → D , we examine all known D * * (s) excited states [1,33] using MC simulation samples. Dedicated exclusive MC studies show that none of these processes, including possible interference effects, exhibit a narrow structure below 4.0 GeV/c 2 [32].
The following three processes that contain excited D * * + s background have potential contributions to the RM (K + ) spectrum: We estimate their production cross sections by studying several control samples. The yields for channel (1) are estimated by analyzing the D * s1 (2536) + peak in the D * 0 K + mass spectra using two separate partially reconstructed samples: K + D − s (with D * 0 missing) and K + D * 0 (with D − s missing). For channel (2), control samples are selected by reconstructing D 0 K + γ (with missing D − s ) or K + D * − s (with missing D 0 ). The D * s2 (2573) + yield is obtained from combined fits to the D 0 K + mass spectra. From this, the contribution from channel (2) to the signal candidates in Fig. 3 [34] to estimate the strength of this background contribution. The shapes in RM (K + ) of these three channels are extracted from MC samples, whereas the normalization is derived from the control samples. The estimated background contributions of the channels (1), (2) and (3) in the RM (K + ) spectrum at √ s = 4.681 GeV are 54.4 ± 8.0, 19.1 ± 7.6 and 15.0 ± 13.3 events, respectively. For the other energy points, the estimated yields of the three channels are given in Ref. [32].
Two processes with excited non-strangeD * * 0 states that produce potential enhancements around 4 GeV/c 2 . In these processes, the RM (K + ) spectrum is distorted due to limited production phase space. The first process is studied using an amplitude analysis of the control sample e + e − → D * 0D * 1 (2600) 0 (→ D − π + ) at all five energy points. Since the ratio B(D * is unknown, it is difficult to project the results of the amplitude analysis into our signal channel. Instead, we determine the ratio in our nominal fit, providing a constraint on the size of the D * 0D * 1 (2600) 0 (→ D − s K + ) component at the different energy points. For the second process, no significant signal is observed in the control [35], the contribution of the D 0 D * 3 (2750) 0 channel to Fig. 3 is estimated to be 0.0 ± 0.4 events, and the corresponding upper limit is taken into account as a source of systematic uncertainty.
As no known processes explain the observed enhancement in the RM (K + ) spectrum, which is very close to the threshold of D − s D * 0 (3975.2 MeV/c 2 ) and D * − s D 0 (3977.0 MeV/c 2 ), we consider the possibility of describing the structure as a D − s D * 0 and D * − s D 0 resonance with a mass-dependent-width Breit-Wigner line shape, denoted as Z cs (3985) − . A simultaneous unbinned maximum likelihood fit is performed to the RM (K + ) spectra at all five energy points, as shown in Fig. 3. The Z cs (3985) − component is modeled by the product of an S-wave Breit-Wigner shape with a mass-dependent width of the following form Here, M is the reconstructed mass; m 0 is the resonance mass; Γ 0 is the width; q is the K + momentum in the initial e + e − system; where B j is the BF of the j th decay. We assume f = 0.5 in the nominal fit and take variations of f into account in the studies of systematic uncertainty.
The Z cs (3985) − signal shape, which is used in the fit depicted in Fig. 3, is the f -dependent sum of the efficiency-weighted F j functions convolved with a resolution function, which is obtained from MC simulation. The resolution is about 5 MeV/c 2 and is asymmetric due to the contribution from ISR. The parametrization of the combinatorial-background shape is derived from the kernel estimate [36] of the WS distribution, whose normalization is fixed to the number of the fitted background events within the decorrelated RM (K + D − s ) signal window. The shapes of the NR and D * 0D * 1 (2600) 0 (→ D − s K + ) signals are taken from the MC simulation. The size of the NR component at each energy point and the ratio B(D * are free parameters in the fit. In addition, a component that describes the total contributions of the excited D * * + s processes is included, whose shape is taken from MC simulation and its size is fixed according to the yields estimated from the control-sample studies. From the fit, the parameters m 0 and Γ 0 are determined to be MeV. The significance of the signal is calculated taking into account the look-elsewhere effect (LEE) [37], where 5000 pseudo-datasets are produced with the sum of null-Z cs (3985) − models and fitted with the same strategy as the nominal fit to obtain the distribution of −2 ln(L 0 /L max ), where L 0 and L max are fitted likelihood values under the null-Z cs (3985) − hypothesis and alternative hypothesis, respectively. In the generation of the pseudo-data, all systematic uncertainties relevant to the mass and width measurements are considered. The resulting distribution is found to be well described by a χ 2 distribution with 13.8 degrees of freedom. With an observed value of −2 ln(L 0 /L max ) = 59.14, we obtain a significance of 5.3 σ. The number of Z cs (3985) − events observed at √ s = 4.681 GeV is the most prominent compared to the other four energy points. If we fit only to data at √ s = 4.681 GeV, we obtain consistent Z cs (3985) − resonance parameters.
The Born cross section σ B (e + e − → K + Z cs (3985) − + c.c.) times the sum of BFs of the decays Z cs (3985) − → D − s D * 0 + D * − s D 0 is equal to nsig Lintfcorrε , where n sig is the number of the observed signal events, L int is the integrated luminosity, andε is the BF-weighted detection efficiency. We define f corr ≡ (1 + δ ISR ) 1 |1−Π| 2 , where (1 + δ ISR ) is the radiative-correction factor and 1 |1−Π| 2 is the vacuumpolarization factor [38]. The numerical results are listed in Table I.
Sources of systematic uncertainties on the measurement of the Z cs (3985) − resonance parameters and the cross section are studied, in which the main sources include the mass scaling, detector resolution, the signal model, background models and the input cross-section line shape for σ B (e + e − → K + Z cs (3985) − ).
We select a control sample of e + e − → D s1 (2536) + D * − s → K + D * 0 D * − s at √ s = 4.681 GeV by detecting K + D * 0 with D * 0 → π 0 D 0 , D 0 → K − π + , K − π + π 0 as well as K − π + π + π − with a missing D * − s in the final state to study the mass scaling of the recoil mass of the low-momentum bachelor K + . We fit the D * − s peak in the spectra of the recoil mass of K + D * 0 , where the D * − s signal is modeled with a MC-determined signal shape convolved with a Gaussian function to represent a potential difference between data and MC simulation. The fitted Gaussian parameters are determined to be µ = −0.2 ± 0.5 MeV/c 2 and σ upper < 1.43 MeV (68% C.L.), which are used to determine the systematic effects due to mass scaling and detection resolution. After incorporating the evaluated detection resolution difference up to the upper uncertainty, we find the maximum change on the result of the fitted width to be 1.0 MeV.
In this work the two Z cs signal processes are difficult to distinguish due to the partial-reconstruction method and the limited sample size. Hence, without any a priori knowledge, we vary the BF ratio f in the range from 0.2 to 0.8, corresponding to the standard deviation of a uniform distribution from 0 to 1. We find the resulting changes on the mass and width to be 0.2 MeV/c 2 and 1.0 MeV, respectively. In the nominal fit, we assume that the spin-parity of the Z cs (3985) − is 1 + and that the relative momentum between K + and Z cs (3985) − in the rest frame of the e + e − system and the relative momentum between D − s (D * − s ) and D * 0 (D 0 ) in Z cs (3985) − system are both in an S-wave state, denoted as 1 + (S, S). This hypothesis can only be verified by an amplitude analysis of the signal final states, which is not feasible with the current statistics. Therefore, as systematic variations, we test the assumptions of spin-parity and angular momentum with 1 + (D, S), 0 − (P , P ), 1 − (P , P ) and 2 − (P , P ) configurations. These tests give maximum changes of 1.0 MeV/c 2 in the mass and 2.6 MeV in the width. The systematic uncertainty related to the combinatorial background is estimated by varying both the sideband yield within its uncertainties and the background parametrization; the quadrature sums of each largest difference from the nominal fit are 0.5 MeV/c 2 and 0.5 MeV for the mass and width, respectively, which are taken as the systematic uncertainties. The efficiency curves adopted in the resonance fit are varied within the uncertainties of their parametrizations, and the differences of 0.1 MeV/c 2 in mass and 0.2 MeV in width to the nominal fit are taken as the related systematic uncertainty.
Any potential effects of the known D * * (s) states (as described further in Ref. [32]) on the measurements are evaluated. We vary the size of the D * * + s andD * 3 (2750) 0 background components within their uncertainties in the fit and take the variations as systematic uncertainties. For the knownD * * 0 states, which have RM (K + ) distributions similar to that of the NR signal, the fit is repeated with each state as an additional component with its shape taken from MC simulation and the yield as a free parameter. To further check theD * 1 (2600) 0 component, we remove the NR component from the simultaneous fit. The ratio B(D * 1 (2600) 0 → D − s K + )/B(D * 1 (2600) 0 → D − π + ) then increases from 0.00±0.02 to 0.12±0.02. We evaluate the quadrature sum of the mass and width differences be-tween each of the results from these alternative fits with respect to the nominal fit and assign the quadrature sums as related systematic uncertainties of 1.0 MeV/c 2 for the mass and 3.4 MeV for the width. We vary the input Born cross section σ B (e + e − → K + Z cs (3985) − ) within the uncertainties and repeat the signal extraction, which gives a maximum change of 0.6 MeV/c 2 for the mass and 1.7 MeV for the width.
Other systematic effects mostly influence the measurement of the cross section. Average uncertainties associated with the tracking, PID and K 0 S reconstruction efficiencies are estimated to be 3.6%, 3.6% and 0.4%, respectively. The efficiency of the RM (K + D − s ) requirement is re-estimated by changing the MC-simulated resolution according to the observed difference with respect to data and the resulting change is taken as the systematic uncertainty on the cross section. The integrated-luminosity uncertainty, measured with large-angle Bhabha events, is estimated to be 1%. The uncertainties on the quoted BFs for the involved decays [1] are included as part of the systematic uncertainty.
In summary, we study the reactions e + e − → The first uncertainties are statistical and the second are systematic. The significance of this resonance hypothesis is estimated to be 5.3 σ over the pure contributions from the conventional charmed mesons. If Z cs (3985) − is established as it claims, it couples to at least one of D − s D * 0 and D * − s D 0 , and has unit charge, the quark composition is most likely ccsū. Hence, it would become the first Z cs tetraquark candidate observed in experiment. The measured mass is close to the mass threshold of D sD * and D * sD , which is consistent with the theoretical calculations in Ref. [18,20,21]. In addition, the Born cross sections σ B (e + e − → K + Z cs (3985) − + c.c.) times the sum of the branching fractions for Z cs (3985) − → D − s D * 0 + D * − s D 0 decays are measured at the five energy points. Due to the limited size of the statistics, only a one-dimensional fit is implemented and the potential interference effects are neglected. This means that the properties of the observed excess might not be fully explored and there exists other possibilities of explaining the near-threshold enhancement. To further improve studies of the excess, more statistics are necessary in order to carry out a sophisticated amplitude analysis.
The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center for their strong support. This work is supported in part by National Supplemental Material for "Observation of a near-threshold excess in the K + recoil mass spectra in e + e − → K + (D − s D * 0 + D * − s D 0 )" Figure 5 shows the distribution of the K + D − s recoil-mass in data and MC simulation samples at √ s = 4.628, 4.641, 4.661 and 4.698 GeV, after the same selection criteria as those imposed for the data shown in Fig. 2 of the main letter. Table II         To understand the potential backgrounds from excited D * * (s) states, all reported states in the PDG [1] whose production and decay is allowed within the available phase-space at √ s = 4.681 GeV are investigated. The corresponding RM (K + ) distributions of the MC simulations are plotted in Figs. 7 and 8. Furthermore, possible interferences among those excited D * * (s) states are systematically scanned, and the choices with the largest interferences around RM (K + ) = 4.0 GeV/c 2 are compared with the distributions in data, shown in Fig. 9 and Fig. 10. It is evident that none of the states can explain the narrow peaking structure below 4.0 GeV/c 2 .   . K + recoil-mass spectra in data with the WS background contributions subtracted, and MC simulations of two possible background processes for the K + D * − s D 0 final state, whose interferences are taken into account. The interference effect is tuned to be largest around 4.0 GeV/c 2 . In the non-resonant (NR) process, the angular momentum (L K + X , L D − s D * 0 ) denotes the angular momentum between K + and X D − s D * 0 , and D − s and D * 0 in the e + e − (X D − s D * 0 ) rest frame, respectively. Individual contributions are scaled according to the observed yields in the control samples.  To avoid potential bias, the analysis strategy is firstly implemented and validated using the first one-third of data set at √ s = 4.681 GeV, where the fit result is shown in Fig. 11(left) and given in Table IV. Afterward, we split the two-thirds of data into two parts for consistency check by implementing the same fit procedures, the results of which are depicted in Fig. 11(middle) and (right). The corresponding numerical results are listed in Table IV. The fitted resonance parameters between the 1st and 2nd one-third of data set are consistent within statistical uncertainty, while the comparison between the 1st and 3rd one-third of data set shows that the fitted masses and widths are in agreement within 1.5σ and 1.0σ, respectively. Overall, the three sets of fit results are compatible and we can assume they are due to the same source. Hence, the three parts of data at √ s = 4.681 GeV are combined to obtain the nominal fit results listed in Table IV.