Measurements of Born Cross Sections of $e^+e^-\to D_s^{*+} D_{sJ}^{-} +c.c.$

The Born cross sections are measured for the first time for the processes $e^+e^-\to D_s^{*+}D_{s0}^*(2317)^- +c.c.$ and $e^+e^-\to D_s^{*+}D_{s1}(2460)^- +c.c.$ at the center-of-mass energy $\sqrt{s}=$ 4.600~GeV, 4.612~GeV, 4.626~GeV, 4.640~GeV, 4.660~GeV, 4.68~GeV, and 4.700~GeV, and for $e^+e^-\to D_s^{*+}D_{s1}(2536)^- +c.c.$ at $\sqrt{s}=$ 4.660~GeV, 4.680~GeV, and 4.700~GeV, using data samples collected with the BESIII detector at the BEPCII collider. No structures are observed in cross-section distributions for any of the processes.


II. DETECTOR, DATA SAMPLES AND MONTE CARLO SIMULATIONS
The BESIII detector located at the Beijing Electron Positron Collider (BEPCII) [35] is a major upgrade of the BESII experiment at the BEPC accelerator [36]. The cylindrical core of the BESIII detector consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon-identifier modules interleaved with steel. The acceptance of charged particles and photons is 93% over the 4π solid angle. The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for the electrons from Bhabha scattering. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end-cap) region. The time resolution of the TOF barrel part is 68 ps. The time resolution of the end cap part was 110 ps before 2015; at that time the end cap TOF system was upgraded with multi-gap resistive plate chambers, and the time resolution improved to 60 ps.
The results reported in this article for the processes e + e − → D * + s D * s0 (2317) − and e + e − → D * + s D s1 (2460) − are determined from seven energy points from 4.600 to 4.700 GeV, where the data at 4.600 GeV were accumulated in 2014, and the remainder in 2020. At energies of 4.660 GeV and above, the data are also used to study e + e − → D * + s D s1 (2536) − . Table I lists the individual c.m. energies and their integrated luminosities.
The GEANT4-based [37] Monte Carlo (MC) simulation framework BOOST [38], which includes the description of the detector geometry and response, is used to produce large simulated event samples. These samples are exploited to optimize the event selection criteria, to determine the detection efficiency, and to evaluate the initial-state radiation (ISR) correction factor (1 + δ). Exclusive phase space (PHSP) MC samples of e + e − → D * + s D * s0 (2317) − , e + e − → D * + s D s1 (2460) − , and e + e − → D * + s D s1 (2536) − are generated using KKMC [39][40][41], where the effects of ISR and beam-energy spread are taken into account. Generic MC samples of open-charm processes are used to estimate background contribu-tions. The known decay modes are modeled with BE-SEVTGEN [42,43], using branching fractions taken from PDG [10], while the unknown decays of charmonium states are modeled using LUNDCHARM [44]. Final-state radiation effects are simulated by the PHO-TOS [45] package. In the signal MC samples, D * + s are simulated to decay into γD + s , with the subsequent decay of D + s into K + K − π + , while the D * s0 (2317) − , D s1 (2460) − , and D s1 (2536) − mesons decay inclusively. A P -wave model and a Dalitz plot decay model [46][47][48] are used to simulate D * + s → γD + s and D + s → K + K − π + , respectively.

III. SELECTION CRITERIA
Thus, there are three charged tracks from D + s decays and one additional photon candidate from D * + s → γD + s . Each track must originate from the interaction point (IP), which means that the distance of the closest approach to the IP of each track is required to be within 10 cm in the beam direction and within 1 cm in the plane perpendicular to the beam direction. Additionally, each track must reside within the active region of the MDC, meaning that its polar angle θ must satisfy | cos θ| < 0.93. The dE/dx and TOF information are used to perform particle identification (PID). Pion candidates are required to satisfy Prob(π) > Prob(K) and Prob(π) > 0.001, where Prob(π) and Prob(K) are the PID probabilities for a track to be a pion and kaon, respectively. Kaon candidates are required to satisfy Prob(K) > Prob(π) and Prob(K) > 0.001. With these PID requirements, the probability of misidentifying a K + as π + or a π + as K + is below 3%.
The photon candidates are selected from EMC showers that are not associated with charged tracks. The deposited energy is required to be larger than 25 MeV in the barrel EMC (|cos θ| < 0.8), or larger than 50 MeV in the end-cap EMC (0.86 < |cos θ| < 0.92). To eliminate the showers from charged particles, the angle between the photon and the extrapolated impact point of any good charged track at the EMC front face must be larger than 20 • . The timing of the shower is required to be within [50,700] ns after the reconstructed event start time to suppress noise and energy deposits unrelated to the event.
To reconstruct D * + s mesons, two kaons, one pion, and one photon candidates are required. All the K + K − π + in all combinations of γK + K − π + , which pass a vertex fit are kept. Mass-constrained fits to the nominal masses of D + s and D * + s (2C) are performed, and the χ 2 of this fit is required to be less than 10 to suppress the backgrounds.
Here P c.m. , P γ , P K + , P K − , and P π + are the fourmomenta of the initial e + e − system, the selected γ, K + , K − , and π + , respectively, M (γK + K − π + ) is the invariant mass of the γK + K − π + system, and m D * +  [49], indicate that there are no peaking backgrounds in the signal region. We considered possible contributions of processes such as e + e − → D + s γD − sJ and e + e − → D 0 (→ K − π + )K + γD * s0 (2317) − by plotting the recoil mass distributions for events in the D * + s mass sidebands; no peaking structures were observed. An unbinned maximum-likelihood fit is performed to the D * + s recoil-mass distributions to determine the signal yields of D − sJ mesons. The signal probability density function is modeled according to a MC-derived signal shape, while the background is modeled with an ARGUS function [50]. The fitted results are summarized in Fig. 1 and Table I Table I. The Born cross sections of e + e − → D * + s D − sJ are calculated using the formula where N fit is the fitted D − sJ signal yield, 1 + δ is the radiative correction factor obtained from a QED calculation with 1% accuracy [51] using the KKMC generator, 1 + δ vp is the vacuum polarization factor, which is taken from Ref. [52] (δ vp is around 0.055 for all studied energy points), and L int is the integrated luminosity at each energy point. The D * + s reconstruction efficiency ǫ D * + s is the product of the detection efficiency ǫ and the branching fractions for D * + s → γD + s and D + s → K + K − π + (93.5% and 5.39%, respectively [10]). The calculation of the upper limits on Born cross sections at 90% C.L. is performed analogously, replacing N fit with N U.L. .  Table I. The systematic uncertainties and the method to take them into account in the upper limits are discussed in Sec. V. The Born cross sections, with statistical uncertainties only, are shown in Fig. 2.

V. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties of the measured cross sections of e + e − → D * + s D * s0 (2317) − , e + e − → D * + s D s1 (2460) − , and e + e − → D * + s D s1 (2536) − are divided into two categories: multiplicative systematic uncertainties and additive systematic uncertainties. The multiplicative systematic uncertainties are associated with tracking and PID efficiencies, photon-detection efficiency, MC sample size, ISR and vacuum-polarization corrections, measurement of luminosity, the branching fractions of intermediate states, and the kinematic fit. The additive systematic uncertainties are associated with the D − sJ mass, the background shape, and the fit range. The uncertainties related to PID and tracking are determined to be 3.0% respectively [53]. The uncertainty of the photon reconstruction efficiency is 1.0%, which is derived from the study of J/ψ → ρ 0 (→ π + π − )π 0 (→ γγ) [54] events. The uncertainties due to finite sizes of the MC samples are determined to be at most 0.9% at each energy point, which arises from the statistical uncertainty of the selection efficiency measured from these samples. From an MC study, we find that the systematic uncertainty due to the choice of D + s → K + K − π + decay model is negligible. The shapes of the cross section of the processes e + e − → D * + s D − sJ affect the radiative correction factor and the detection efficiency. Due to the limited number of events in the data sample, a detailed determination of the energy dependence ("line shape"), which would allow for an iterative determination of radiative correction factors, is not possible. Therefore, the input line shapes are changed to a first order polynomial multiplied by √ E m − E 0 , a function that reasonably describes the shape of available data, where E m stands for the c.m. energy and E 0 stands for the threshold energy for each process, and the differences in the selection efficiency ε(1 + δ) are taken as the systematic uncertainties. The uncertainty from vacuum polarization is less than 0.1% [52], which is negligible compared to other sources of uncertainties. The integrated luminosities of the data samples are measured using large-angle Bhabha scattering events with an uncertainty less than 1.0%. The uncertainties in the branching fractions B(D + s → K + K − π + ) and B(D * + s → γD + s ) are 2.8% and 0.7% [10], respectively. The uncertainty of the 2C kinematic fit is estimated using control samples of e + e − → D * + s D * − s events at √ s = 4.42 GeV and 4.6 GeV. The difference between the data and MC efficiencies due to the 2C kinematic fit is 1.7%, which is taken as the systematic uncertainty. Using a control sample of e + e − → D * + s D * − s , the mass resolutions of D * − s candidates between signal MC simulation and data are consistent. Thus, the systematic uncertainty due to mass resolution is negligible.
The uncertainties due to D − sJ mass are estimated by varying its value by the measured uncertainty [10]. The differences in the fitted D − sJ signal yields are taken as the systematic uncertainties. Using a control sample of e + e − → D * + s D * − s decays, we find the resolution in missing mass to be essentially the same in data as in the MC, at the current level of statistics. The systematic uncertainties due to fitting itself, the background shape, and the fit range are estimated with a toy MC method. We simulate an ensemble of experiments, generating D * + s recoil mass distributions based on the nominal fitted results. We generate 1000 distributions and subsequently fit them to obtain D − sJ signal yields. We plot these results in Gaussian distributions. The difference between the mean values of these distributions and the input signal yields represents the systematic bias or uncertainty due to the fitting procedure. The uncertainty is found to be negligible for all D − sJ signals. A similar method is used to estimate the systematic uncertainties due to the fit range and the background shape. For the fit range, the lower bound is changed from 2.20 GeV/c 2 to 2.18 GeV/c 2 and to 2.22 GeV/c 2 . For the background shape, instead of an ARGUS function [50], For those energy points with a D − sJ statistical significance larger than 3σ, the central values of the cross section with statistical and systematic uncertainties are reported, and the systematic uncertainties are summarized in Table II. For the other energy points, the upper limits on the cross section at 90% C.L. are reported and the systematic uncertainties are taken into account in two steps. Firstly, among the additive systematic uncertainties described above, the most conservative upper limit at 90% C.L. is kept. Then, to take into account the multiplicative systematic uncertainty, the likelihood curve is convolved with a Gaussian function with a width parameter equal to the corresponding total multiplicative systematic uncertainty. All of the multiplicative systematic uncertainties for the energy points with a D − sJ signal significance less than 3σ are summarized in Table III. Assuming that all the sources are independent, the total systematic uncertainty is obtained by adding them in quadrature. The final results of the Born cross section with systematic uncertainties considered are listed in Table I Table I, and shown in Fig. 2 with statistical er- sJ signal events N fit , the upper limit at 90% C.L. on the number of fitted D − sJ signal yields NU.L., the ISR radiative correction factor (1 + δ), the statistical signal significance, and the measured Born cross section (σB) and its upper limit (σ U.L.

VII. ACKNOWLEDGEMENT
The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center for their strong sup-