Evidence for a neutral near-threshold structure in the K 0S recoil-mass spectra in e + e − → K 0S D + s D ∗− and e + e − → K 0S D ∗ + s D −

We study the processes e + e − → K 0S D + s D ∗− and e + e − → K 0S D ∗ + s D − , as well as their charge conjugated processes, at ﬁve center-of-mass energies between 4.628 GeV and 4.699 GeV, using data samples corresponding to an integrated luminosity of 3 . 8 fb − 1 collected by the BESIII detector at the BEPCII storage ring. Based on a partial reconstruction technique, we ﬁnd evidence of a structure near the thresholds for D + s D ∗− and D ∗ + s D − production in the K 0S recoil-mass spectrum, which we refer to as the Z cs (3985) 0 . Fitting with a Breit-Wigner line shape, we ﬁnd the mass of the structure to be (3992 . 2 ± 1 . 7 ± 1 . 6

We study the processes e + e − → K 0 S D + s D * − and e + e − → K 0 S D * + s D − , as well as their charge conjugated processes, at five center-of-mass energies between 4.628 GeV and 4.699 GeV, using data samples corresponding to an integrated luminosity of 3.8 fb −1 collected by the BESIII detector at the BEPCII storage ring. Based on a partial reconstruction technique, we find evidence of a structure near the thresholds for D + s D * − and D * + s D − production in the K 0 S recoil-mass spectrum, which we refer to as the Zcs(3985) 0 . Fitting with a Breit-Wigner line shape, we find the mass of the structure to be (3992.2 ± 1.7 ± 1.6) MeV/c 2 and the width to be (7.7 +4.1 −3.8 ± 4.3) MeV, where the first uncertainties are statistical and the second are systematic. The significance of the Zcs(3985) 0 signal is found to be 4.6σ including both the statistical and systematic uncertainty. We report the Born cross section multiplied by the branching fraction at different energy points. The mass of the Zcs(3985) 0 is close to that of the Zcs(3985) + . Assuming SU(3) symmetry, the cross section of the neutral channel is consistent with that of the charged one. Hence, we conclude that the Zcs(3985) 0 is the isospin partner of the Zcs(3985) + .
The charged-tetraquark candidate Z cs (3985) + [1] was observed at BESIII in the D + sD * 0 and D * + sD 0 final states [22][23][24][25][26]. The mass of the Z cs (3985) + is close to the D + sD * 0 and D * + sD 0 thresholds, which is consistent with theoretical predictions [17][18][19][20]. Meanwhile, another charged-tetraquark candidate, Z cs (4000) + [27], was observed in the J/ψK + final states in an amplitude analysis of the decay B + → J/ψφK + at LHCb. However, the widths of these two Z + cs states are inconsistent with each other. The observation of these charged Z cs states motivates a search for a neutral isospin partner Z 0 cs . The mass of the Z 0 cs is expected to be heavier than that of the Z + cs by (0.05 ± 0.21) GeV/c 2 under the molecular hypothesis, or by (0.06 ± 0.12) GeV/c 2 under the tetraquark hypothesis [23]. A promising approach to this challenge at BESIII is to search for the process e + e − →K 0 Z cs (3985) 0 + c.c. and then compare its cross section to that of e + e − → K − Z cs (3985) + + c.c., which tests the isospin symmetry in the production and decay dynamics. A similar strategy was pursued in the analysis of the Z c charged and neutral states [10,11]. Observation and study of the Z 0 cs is crucial for understanding the nature of the Z cs states. In this letter, we study the processes e + e − → K 0 S D + s D * − and e + e − → K 0 S D * + s D − , which is denoted as e + e − → K 0 S (D + s D * − + D * + s D − ) in the context, as well as their charge conjugated modes, using e + e − collision data sets corresponding to an integrated luminosity of 3.8 fb −1 [28] at center-of-mass energies √ s = 4.628, 4.641, 4.661, 4.682 and 4.699 GeV [28]. These samples were collected by the BESIII detector at the Beijing Electron Positron Collider (BEPCII). Detailed information about BEPCII and BESIII can be found in Refs. [29][30][31]. We use a partial reconstruction technique to maximize the detection efficiency; only the K 0 S produced in association with the D + s D * − or D * + s D − (the bachelor K 0 S ) and one of the ground-state D mesons (here D subsequently denotes D + s or D − ) are detected, while the other final-state particles are not reconstructed. The Z 0 cs candidate is then searched for in the invariant mass distribution recoiling against the bachelor K 0 S candidate. Charge conjugation is implied throughout the discussion.
Simulated samples produced with a geant4-based [32] Monte Carlo (MC) package, which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to understand the backgrounds. The e + e − annihilations are simulated with the kkmc [33] generator, which includes the effects of the beam-energy spread and initial-state radiation (ISR). The inclusive MC sample consists of the production of open-charm hadronic systems, ISR production of vector charmonium(-like) states, and continuum processes incorporated in kkmc [33]. The known decay modes are modelled with evtgen [34] using branching fractions reported by the Particle Data Group (PDG) [2], and the remaining unknown decays from charmonium states are modelled with lundcharm [36]. The final-state radiation (FSR) from charged final-state particles is simulated with the photos package [37]. For the non-resonant three-body signal processes e + e − → K 0 S (D + s D * − + D * + s D − ), the momenta distributions of final-state particles are generated following phase space. For the resonant signal process e + e − → K 0 S Z 0 cs → K 0 S (D + s D * − + D * + s D − ), we assume that the Z 0 cs state has a spin-parity of 1 + , which corresponds to S-waves in both of the decays e + e − → K 0 S Z 0 cs and Z 0 cs → D + s D * − + D * + s D − , which we denote as (S, S). The corresponding angular distribution is taken into account in simulating the cascade decays. Other possibilities for the Z 0 cs spin-parity are tested to evaluate the systematic uncertainty related to this assumption. We carry out two types of partial reconstruction, which are referred as the D + s -tag and D − -tag methods, respectively. For D + s (D − )-tag method, only the bachelor K 0 S and D + s (D − ) candidates are reconstructed. We use the decay modes D + s → K + K − π + , K 0 S K + , K + K − π + π 0 , K 0 S K + π + π − and η ′ π + to form the D + s candidates; and the decay modes D − → K + π − π − , K 0 S π − and K 0 S π + π − π − to form the D − candidates. To ensure that each charged track, which is not associated to K 0 S detection, originates from the e + e − interaction point (IP), |V r | < 1 cm and |V z | < 10 cm are required. Here, |V r | is the distance between the charged track and the beam axis in the transverse plane, and |V z | is the closest distance of the charged track to the IP along the axis of beam. The polar angles of charged tracks are required to satisfy |cosθ| < 0.93. The flight time in the time-of-flight system and the energy deposited in the multilayer drift chamber for each charged track are used to identify particles by calculating the probabilities P(i), where i denotes K or π. We require P(K) (P(π)) to be greater than P(π) (P(K)) to classify a particle as a kaon (pion) candidate.
The K 0 S candidates are reconstructed through the π + π − decay mode without particle identification requirements. Both pions must satisfy |V z | < 20 cm, and |cosθ| < 0.93 and their trajectories are constrained to originate from a common vertex by applying a vertex fit, the χ 2 of which is required to be less than 100. The K 0 S candidate is then formed and the opposite direction of its momentum is constrained to point at the IP, with the corresponding χ 2 required to be less than 40. The decay length of K 0 S candidate must be greater than two standard deviations of the vertex resolution away from the IP. The invariant mass of π + π − pair, M (π + π − ), is required to be within (0.492, 0.503) GeV/c 2 .
The π 0 and η candidates are reconstructed through π 0 /η → γγ. The photon showers in the electromagnetic

Final state Requirement
calorimeter must have energies greater than 25 MeV in the barrel region (|cos θ| < 0.80) and greater than 50 MeV in the end-cap region (0.86 < | cos θ| < 0.92). Showers must have an associated time within 700 ns of the event start time. A kinematic fit is applied to constrain the invariant mass of the γγ pair to the known π 0 or η mass reported in the PDG [2], and the resultant χ 2 is required to be less than 10. The η ′ is reconstructed through η ′ → π + π − η. The mass of the π + π − η is required to be within 10 MeV/c 2 of the known η ′ [2] mass.
To improve the signal purity, the requirements listed in Table I are adopted to restrict the final states within the regions of the φ, K * and ρ resonances, which dominate the decays. In the selection of The recoil mass RM (K 0 S D) of the K 0 S D system is obtained according to RM (X) = ||p e + e − − p X ||, where p e + e − is the four-momentum of the initial e + e − system and p X is the four-momentum of the X system. The   We find the threshold enhancement can not be explained by these excited states, and hence, we consider its possible origin to be the neutral Z 0 cs state. A simultaneous unbinned maximum likelihood fit is applied to the distributions of RM (K 0 S ) at five energy points. We adopt two S−wave Breit-Wigner functions R 1 and R 2 to describe the Z 0 cs resonance where R 1 describes the decay Z 0 cs → D + s D * − , and R 2 describes Z 0 cs → D * + s D − , M equals RM (K 0 S ), m 0 is the mass of the Z 0 cs , and Γ 0 is the total width of the Z 0 cs . The momentum of the K 0 S in the initial e + e − system is q, the momentum of the D + s (D − ) in the rest frame of the D + s D * − (D * + s D − ) system is p 1 (2) , and the corresponding momentum at M = m 0 is p * 1 (2) . In the fit, under the assumption of the isospin symmetry, a Gaussian constraint is imposed to restrict the width of the Z 0 cs within the uncertainty of the Z cs (3985) + width, which is (13.8 +8.1 −5.2 ± 4.9) MeV [1]. The factor f denotes the ratio of the two signal channels The default value of f is chosen to be 0.5, with other possibilities considered as a systematic uncertainty.
The fit depends on the detector resolution and mass-dependent efficiency, which are derived from simulated samples. The detector resolution is determined using the Z 0 cs signal MC samples, in which the width of the Z 0 cs is set to be 0. The signal probability density function (PDF) is constructed as follows: where E 1(2) is the efficiency function and G is the Gaussian resolution function. The backgrounds in the fit include three components: the non-resonant process e + e − → K 0 S (D + s D * − + D * + s D − ), the excited D * * s D s backgrounds, and the combinatorial backgrounds. The first and second components are described using histogram PDFs extracted from MC samples, and the third component is described using the distribution from the D + s (D − ) sideband. In the fit, the yields of the excited D * * s D s backgrounds are estimated from isospin relations according to those calculated for the e + e − → K − Z cs (3985) + process, and the numbers are fixed in the fit [1], while the yields of the non-resonant process are free. The sizes of the combinatorial background are fixed to the values in Table II. The fitted mass and width of the Z 0 cs are given in Table III, where the Z cs (3985) + resonance parameters are included for comparison. The results are consistent with the theoretical predictions [17][18][19][20]23]. We sum up the RM (K 0 S ) distributions from all data sets, and superimpose the simultaneous fit curves in the last plot of Fig. 2. Comparing the fits with or without considering the contribution from the Z 0 cs , the number of degrees of freedom is changed by 7 (the mass and width of the Z 0 cs , together with the cross section of the Z 0 cs at the five center-of-mass energies). The value of 2 ln L, where L is the likelihood value, is changed by 42.0. This corresponds to a statistical significance of 5.0σ according to Wilks' theorem [40]. When also considering systematic uncertainties, which are described in the supplemental material [39], the significance of the Z 0 cs signal becomes 4.6σ. The reduced chi-squared of the fit in Fig. 2 is 0.9, indicating good compatibility between the model and the data. According to the fitted signal yields in Table IV  following equation where N obs are the signal yields,ǫ are the combined MC-determined reconstruction efficiencies in the two D-tag methods, L is the integrated luminosity, (1 + δ) is the radiative correction factor, and δ vac is the vacuum-polarization correction factor [41]; their values are given in Table IV. We assume B(Z 0 The factor of 2 in the denominator in Eq. (3) is necessary because of the equal transition rate of K 0 andK 0 to K 0 S . The cross section results at the five center-of-mass energies are listed in Table V. The χ 2 of each energy point is defined as the square of difference of the cross sections of two channels divided by the sum-of-squares of these uncertainties. The χ 2 total /ndf is the sum of the χ 2 divided by the number of energy points. The cross section results for the neutral channel are consistent with those for the charged one [1], which agree with the prediction based on isospin symmetry.
Systematic uncertainties on the measurement of the Z 0 cs resonance parameters and production cross sections are extensively investigated as detailed in Ref. [39]. An important contribution is associated with the background modelling in the fit and the Z 0 cs signal model. For the background modelling, we vary the size and shape of the combinatorial backgrounds according to the M (D) sideband control samples, as well as explore the additional contributions from the highly excited D * * (s) states. For the signal modelling, we test different J P assignments of the Z 0 cs by changing the matrix elements in the signal simulations. The total systematic uncertainties are, overall, similar to the statistical uncertainties on each measurement.
In summary, based on data sets with center-of-mass energies from 4.628 GeV to 4.699 GeV at BESIII, evidence of a neutral open-strange hidden-charm state, Z cs (3985) 0 , is found in the K 0 S recoil-mass spectrum of the e + e − → K 0 S (D + s D * − + D * + s D − ) + c.c. processes, with a resonance mass and width determined as (3992.2 ± 1.7 ± 1.6) MeV/c 2 and (7.7 +4.1 −3.8 ± 4.3) MeV, respectively. The significance of the state is determined to be 4.6σ. Since this state decays through D + s D * − and D * + s D − , it should contain at least four quarks, ccsd. The measured mass of the Z cs (3985) 0 is larger than that of the Z cs (3985) + , which is consistent with theoretical prediction [23]. In addition, the Born cross sections of e + e − →K 0 Z cs (3985) 0 + c.c. multiplied by the branching fraction of Z cs (3985) 0 → D + s D * − + D * + s D − at the five energy points are measured and found to be consistent with those of e + e − → K − Z cs (3985) + + c.c. [1], as is expected under isospin symmetry. Hence, we conclude that the Z cs (3985) 0 is the isospin partner of the Z cs (3985) + .
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 Key R&D Program of China under Contracts Nos. 2020YFA0406400, 2020YFA0406300; National Natural Science Foundation of China (NSFC) under Contracts Nos. 11521505, 11635010,  11735014, 11805086, 11822506, 11835012, 11935015, 11935016, 11935018, 11961141012, 12022510, 12025502 Supplemental Material for "Evidence for a neutral near-threshold structure in the K 0 S recoil-mass spectra in e + e − → K 0 S D + s D * − and e + e − → K 0 S D * + s D − " Appendix A: Fit results of Ds-tag and D-tag methods Figure 3 and Fig. 4 show the fits to the recoil mass distributions RM (K 0 S ) using only the D + s -tag or D − -tag method, respectively. Table VI lists   To avoid potential bias, the analysis strategy is firstly implemented using 1/3 of the data set at √ s = 4.682 GeV. The remaining 2/3 of data set at √ s = 4.682 GeV is then used for a consistency check.
The distributions of RQ(K 0 S D + s ) and RQ(K 0 S D − ) are shown in Fig. 5 and Fig. 6, respectively, for the two subsets. Using the method described in the paper, we determine the number of combinatorial background events in the D * − and D * + s signal regions that are listed in Table VII. The ratios between the numbers from two data subsets are all consistent with two, as expected.
The RM (K 0 S ) distributions from two data subsets are shown in Fig. 7. The yields of Z cs (3985) 0 are obtained from the fits to the distributions, which are listed in Table VIII. The ratio between the number of Z cs (3985) 0 signal events from two data subsets is also consistent with two.          The total systematic uncertainties on the Z cs (3985) 0 resonance parameters and cross sections are the quadrature sums of the assigned uncertainties arising from the sources discussed below. A summary of these contributions is listed in Table IX. When including all these sources of systematic uncertainty, the significance of the Z cs (3985) 0 signal becomes 4.6σ.
In the nominal fit to the RM (K 0 S ) spectra, we choose to constrain the width of the Z 0 cs with the uncertainty of the Z + cs width, to improve the precision of our measurement, according to the isospin symmetry of the Z 0 cs and Z + cs . If the constraint is removed, the fit width of the Z 0 cs becomes 4.1 +4.7 −3.9 (stat. only), which is consistent with the nominal result. Tracking, PID and reconstruction of intermediate states: The uncertainties on both the tracking and PID efficiencies for each charged track are assigned to be 1%. The uncertainties associated with K 0 S , π 0 and η reconstruction are assigned to be 2%. The uncertainties from tracking, PID and intermediate states reconstruction in different tag channels are weighted by the factor B l ε l . Here, "l" indicates each D + or D − s decay channel. D − /D + s signal window: The uncertaintiy associated with the defintion of the D − /D + s signal window is estimated by comparing the D − /D + s signal from data and MC. The widths of the D − /D + s peaks in data and MC are slightly different. We estimate that these differences in resolution lead to a relative 0.5% difference in efficiency, which is assigned as a systematic uncertainty.
Mass scale: A control sample of e + e − → K 0 S D − D + s events with √ s larger than 4.62 GeV is selected, in which the K 0 S and D + s are reconstructed. We fit the D − peak in the corrected recoil mass spectrum RM (K 0 S D + s ) + M (D + s ) − m(D + s ). The D − signal is modelled with a MC-determined signal shape convolved with a Gaussian function. The Gaussian parameters are determined to be µ=(0.07 ± 0.68) MeV/c 2 and σ=(0.60 ± 2.62) MeV. Since the corrected recoil mass RM (K 0 S D + s ) + M (D + s ) − m(D + s ) is largely insensitive to the resolution of the D + s mass, we attribute any mass shift to the bachelor K 0 S . Hence, cconsidering the central value and uncertainty of this study, we take a maximum mass shift of 0.8 MeV/c 2 as the systematic uncertainty.
Detector resolution: To understand the potential difference of detector resolution in data and MC simulations, the same control sample of e + e − → K 0 S D − D + s events is used. From the "Mass scale" study , the width of the smearing function is at most 3.2 MeV. We therefore smear the resolution function in the Z 0 cs fit by this amount and reperform the mass fit. The resultant differences on the final results are taken as systematic uncertainties.
f factor: In the default fit, the two signal processes Z 0 cs → D + s D * − and Z 0 cs → D * + s D − are combined and we assume their fraction factor is 0.5 in nominal calculation. To estimate the possible systematic bias arising from this source, we assume the probability distribution of f is uniform between 0 and 1 with no prior knowledge, we take the RMS value of 1/ √ 12 ( 0.3) as the uncertainty on f . Hence, we vary f to 0.2 and 0.8 and take the largest difference with respect to the nominal result as the systematic uncertainty from this source.