Measurements of $e^+e^- \to K_{S}^{0}K^{\pm}\pi^{\mp}\pi^0$ and $K_{S}^{0}K^{\pm}\pi^{\mp}\eta$ at center-of-mass energies from $3.90$ to $4.60~\mathrm{GeV}$

Using $5.2 \ \mathrm{fb}^{-1}$ $e^+ e^-$ annihilation data samples collected with the BESIII detector, we measure the cross sections of $e^+e^- \to K_S^0 K^\pm \pi^\mp \pi^0$ and $K_{S}^{0}K^{\pm}\pi^{\mp}\eta$ at center-of-mass energies from $3.90$ to $4.60$ GeV. In addition, we search for the charmonium-like resonance $Y(4260)$ decays into $K_{S}^{0}K^{\pm}\pi^{\mp}\pi^0$ and $K_{S}^{0}K^{\pm}\pi^{\mp}\eta$, and $Z_c^{0,\pm}(3900)$ decays into $K_{S}^{0}K^{\pm}\pi^{\mp,0}$ and $K_{S}^{0}K^{\pm}\eta$. Corresponding upper limits are provided since no clear signal is observed.

Using 5.2 fb −1 e + e − annihilation data samples collected with the BESIII detector, we measure the cross sections of e + e − → K 0 S K ± π ∓ π 0 and K 0 S K ± π ∓ η at center-of-mass energies from 3.90 to 4.60 GeV. In addition, we search for the charmonium-like resonance Y (4260) decays into K 0 S K ± π ∓ π 0 and K 0 S K ± π ∓ η, and Z 0,± c (3900) decays into K 0 S K ± π ∓,0 and K 0 S K ± η. Corresponding upper limits are provided since no clear signal is observed.

I. INTRODUCTION
With the experimental progress in the past decade, many charmonium-like (XY Z) states were observed, which can not be accommodated within the naive quark model and are proposed as the candidate of the hiddencharm exotic mesons [1,2]. In this paper, we focus on the exotic states Y (4260) and Z c (3900).
Despite many interpretations, the natures of the Y (4260) and Z c (3900) are still unclear. To comprehend these states, it is necessary to study more decay modes. All of the observed decay modes of the Y (4260) and Z c (3900) are associated with the charm sector and no light hadron decay modes have been found yet [51][52][53]. A search for light hadron decay modes of the Y (4260) and Z c (3900) is complementary to previous studies and may help to distinguish between different theoretical models and to understand strong interaction effects in this energy region. Among the large number of potential light hadron decay modes, the branching fractions (BFs) of charmonium states decaying into K S K ± π ∓ π 0 and K S K ± π ∓ η are usually large [54]. In four-body final state, there should be abundant intermediate states, which may supply more possible decay channels for searching Y (4260) and Z c (3900). Furthermore, the existence of charged and neutral pions in the final states enable a study of isospin multiplets. In this paper, we present a measurement of the Born cross sections (σ B ) for e + e − → K 0 S K ± π ± π 0 and K 0 S K ± π ∓ η. We also report upper limits of

II. DETECTORS AND DATA SAMPLES
The BESIII detector [55], operating at the BEPCII collider [56], is a general purpose spectrometer with a geometrical acceptance of 93% of 4π solid angle. It has four main components: (1) a small-cell, helium-based (60% He, 40% C 3 H 8 ) multi-layer drift chamber (MDC) with 43 layers providing an average single-hit resolution of 135 µm, a momentum resolution of 0.5% at 1.0 GeV/c in a 1.0 T magnetic field, and a specific ionization energy loss (dE/dx) resolution better than 6%, (2) a time-offlight (TOF) detector constructed of 5 cm thick plastic scintillators, with 176 strips of 2.4 m length in two layers in the barrel and 96 fans of the end-caps with time resolutions of 80 and 110 ps, respectively, which provide a 2σ K/π separation for momenta up to ∼ 1.0 GeV/c, (3) an electromagnetic calorimeter (EMC) consisting of 6240 CsI (Tl) crystals in a cylindrical barrel structure and two end-caps with an energy resolution of 2.5% (5%) at 1.0 GeV and a position resolution of 6 mm (9 mm) in the barrel (end-caps), and (4) a muon counter (MUC) consisting of resistive plate chambers (RPCs) in nine barrel and eight end-cap layers, which provide a 2 cm position resolution. More details of the BESIII detector can be found in Ref. [55].
This analysis is based on 5.2 fb −1 e + e − annihilation data samples [57] collected with the BESIII detector at center-of-mass energies ( √ s) from 3.90 to 4.60 GeV [58], which are listed in Table I. Monte-Carlo (MC) simulations are used to optimize the event selection criteria, to study the detector acceptance and to understand the potential backgrounds. The geant4-based [59] simulation software boost [60] is implemented to simulate the detector response, describe geometry and material, realize digitization, and incorporate time-dependent beam backgrounds. Six generic MC samples, equivalent to the integrated luminosity of the data at the energy points 4.009, 4.230, 4.260, 4.360, 4.420 and 4.600 GeV are generated to study the backgrounds. The primary known decay channels are generated using evtgen [61] with the BFs set to the world average values [54] while the unknown decay modes are generated with lundcharm [62]. Continuum hadronic events are generated with kkmc [63] and QED processes such as Bhabha scattering, dimuon, and digamma events are generated with kkmc and babayaga [64].
To study the efficiency of each final state, a sample of 1 × 10 5 signal events is generated at each energy point using kkmc, which simulates e + e − annihilation, including beam energy spread and ISR effects.
Candidate events for e + e − → K 0 S K ± π ∓ π 0 /η, with K 0 S → π + π − and π 0 /η → γγ are selected according to the following steps. First, K 0 S candidates are selected by looping over all pairs of oppositely charged tracks, which are assumed to be pions. Next, primary and secondary vertex fits [65] are performed and the decay length of the secondary vertex fit is required to be greater than twice its uncertainty. Furthermore, the invariant mass of the pion pair is required to satisfy 54]. If there are multiple K 0 S candidates in one event, the one with the smallest χ 2 from the secondary vertex fit is selected.
In addition to the two charged tracks that make up the K 0 S , two oppositely charged tracks are required. For the latter two charged tracks, the polar angle θ must satisfy | cos θ| < 0.93 and the distance of closest approach to the interaction point must be less than 10.0 cm and 1.0 cm along the beam direction and in the plane perpendicular to the beam direction, respectively. The particle type for each charged track is determined by selecting the hypotheses with the highest probability, which is calculated with the combined information from TOF and dE/dx measurements for different particle hypotheses. One charged track must be identified as a kaon and the other as a pion.
Photons are reconstructed from clusters deposited in the EMC, with the energy measured in the TOF included to improve reconstruction efficiency and energy resolution. At least two photons are required per event. The energy of a photon candidate is required to be larger than 25 MeV in the barrel region (| cos θ| < 0.80) or 50 MeV in the end-cap region (0.86 < | cos θ| < 0.92). The cluster timing is required to be between 0 and 700 ns to suppress electronic noise and energy depositions unrelated to the event of interest. To eliminate showers associated with charged particles, the opening angle between a photon candidate and the extrapolated position of the closest charged track should be larger than 20 degrees.
Finally, a four-constraint (4C) kinematic fit imposing energy-momentum conservation is performed to the final states. Only events with χ 2 4C < 60 are accepted. For events with more than two photon candidates, the photon pair with the smallest χ 2 4C from the kinematic fit is accepted. After the 4C kinematic fit, no peaking background is observed in the generic MC samples. The invariant mass distributions of π + π − versus γγ at 4.258 GeV are shown in Fig. 1 as an example in which obvious K 0 S and π 0 /η peaks are observed. The signal regions are defined as M (π + π − ) ∈ (0.488, 0.508) GeV/c 2 , M (γγ) ∈ (0.12, 0.15) GeV/c 2 (for the π 0 mode) and M (γγ) ∈ (0.52, 0.58) GeV/c 2 (for the η mode). The sideband regions are defined as M (π + π − ) ∈ (0.463, 0.483) ∪ (0.513, 0.533) GeV/c 2 , M (γγ) ∈ (0.08, 0.11) ∪ (0.16, 0.19) GeV/c 2 (for the π 0 mode) and M (γγ) ∈ (0.44, 0.50) ∪ (0.60, 0.66) GeV/c 2 (for the η mode). The signal yields at each energy, presented in Tables I and II, are obtained according to N sig = N A − N B /2 + N C /4, where N is the number of events and the subscript A denotes the signal region, and the subscripts B and C denote the sideband regions.
The Born cross section is calculated from where L is the integrated luminosity, ǫ is the detection efficiency, B is the product of the BF of K 0 S → π + π − and that of π 0 /η → γγ [54], 1 |1−Π(s)| 2 is the vacuum polarization correction factor [66], and (1 + δ ISR ) is the ISR correction factor [67] which is determined by the MC simulation programmer kkmc. The ISR factors are set to 1.0 to get the initial cross section lineshape as input to kkmc. From kkmc, the updated ISR factors are obtained, then the cross section lineshape is updated too. We repeat this process till both ISR factors and cross section converge.
The  couplings tuned to approximately match those appear in the data sample and is weighted according to the momentum distributions observed in the data sample. As illustrated in Figs. 2 and 3, the mixing MC sample gives a much better description of the data than a phase space (PHSP) MC sample. The observed cross sections are presented in Tables I and II, and illustrated in Fig. 4.
Since there is no obvious structure in the line shapes of the Born cross sections for e + e − → K 0 S K ± π ∓ π 0 and e + e − → K 0 S K ± π ∓ η, as shown in Fig. 4, the upper limits of Y (4260) → K 0 S K ± π ∓ π 0 and Y (4260) → K 0 S K ± π ∓ TABLE I. Data sets and results of the Born cross section measurement for e + e − → K 0 S K ± π ∓ π 0 . The table includes the integrated luminosity L, the number of observed signals events Nsig, the efficiency ǫ, the ISR correction factor (1 + δ ISR ), the vacuum polarization correction factor 1 |1−Π(s)| 2 , and the Born cross section σB. The first errors are statistical and the second ones are systematic. The details of systematic uncertainties are described in Sec. III D.
describes the continuum process e + e − → K 0 S K ± π ∓ π 0 /η, the parameters p 0 and p 1 are determined by fitting the line shapes with only the continuum process. BW( √ s) given in Eq.   is obtained. Here σ B and σ fit B are the measured and fitted Born cross sections, δ i is the energy dependent part of the total uncertainty, which includes the statistical uncertainty and the energy dependent part of systematic uncertainty, the δ c is the energy independent part of the systematic uncertainty (the systematic uncertainties are described in detail in Sec. III D), h is a free parameter introduced to take into account the correlation of different energy points, and the subscript i indicates the index of each energy point [68] . The Q 2 is used to calculate the likelihood L = e −0.5Q 2 , whose normalized distribution is used to get the upper limits of Γ e + e − B at the 90% confidence level (C.L.), which is determined to be 0.050 eV and 0.19 eV for the π 0 mode and the η mode, respectively. p 0 ( √ s) p 1 and parameters p0 = (6.14 ± 1.54) × 10 6 (pb) · (GeV) p 1 and p1 = 6.68 ± 0.17 in the π 0 mode and p0 = (1.86 ± 0.97) × 10 5 (pb) · (GeV) p 1 and p1 = 6.56 ± 0.36 in the η mode. The pink dash-dotted lines are the MC shape of the Y (4260) with an arbitrary scale factor.
C. Upper limits on σB e + e − → π 0,∓ Z 0,± c (3900), Since there is no obvious Z c (3900) signal in the invariant mass distributions of K 0 S K ± π ∓,0 (π 0 mode) and K 0 S K ± η (η mode), as shown in Figs. 2 and 3, the upper limits at the 90% C.L. for the production cross section σ B e + e − → π Z c (3900) , with Z c (3900) → K 0 S K π/η are determined with an unbinned maximum likelihood fit to the invariant mass of K 0 S Kπ/η in the range (3.7, 4.1) GeV/c 2 , at the five energy points 4.226, 4.258, 4.358, 4.416, and 4.600 GeV. The contribution of non-K 0 S or non-π 0 /η backgrounds is negligible. In the fit, the Z c (3900) signal is described by the MC simulated shape, and the mass and width of the Z c (3900) are set to theirs world average value 3886.6 ± 2.4 MeV/c 2 and 28.2 ± 2.6 MeV/c 2 [54], respectively. The background is described by a second order polynomial function. The normalized likelihood distribution of the Born cross section L(σ B ) is determined by changing the number of signal events from 0 to 150 with a step size of 1. The upper limit (UL) at the 90% C.L. is calculated by solving the equation The final upper limits are shown in Table III, where all of the systematic uncertainties have been considered, the details of which are explained in Sec. III D. The ratio is also given in Table III, where the cross sections for e + e − → π ∓ Z c (3900) ± → π + π − J/ψ and e + e − → π 0 Z c (3900) 0 → π 0 π 0 J/ψ are from Ref. [45] and Ref. [69], respectively.

D. Systematic Uncertainties
Various sources of systematic uncertainty are investigated in the K 0 S K ± π ∓ π 0 and K 0 S K ± π ∓ η lineshape measurement. We assume that the systematic uncertainties associated with the physics model used in the MC simulation, the luminosity, tracking, PID, γ reconstruction efficiency, K 0 S reconstruction efficiency, ISR correction factor, vacuum polarization factor and quoted BFs are energy independent, while the other systematic effects are energy dependent.
For the π 0 mode, a data-driven MC method is developed to obtain the efficiency. To estimate the uncertainty of this method, one thousand testing samples of e + e − → K 0 S K ± π ∓ π 0 are generated with eighteen different physics processes with random ratios, the ratio of each process is generated using uniform distribution between 0 to 1 and then normalized by the summation of these eighteen ratios. The difference between the estimated and the real efficiencies is fitted with a Gaussian function. The fit results give a mean of 0.4% which is neglected, and a width of 0.9% which is taken as the systematic uncertainty from the data-driven MC method. For the η mode, which has much lower statistics than the π 0 mode, alternative mixing ratios are used to generate a new MC sample and the efficiency difference between the two MC samples is adopted as the systematic uncertainty.
The uncertainty on the integrated luminosity is estimated to be 1.0% using Bhabha events [57].
The ISR correction factor introduces a 1.0% uncertainty since the termination condition of the recursion method used to get the correction factor is 1.0% between the last two iterations.
The uncertainty due to the vacuum polarization factor is found to be negligible [66]. The uncertainties of the quoted BFs are also considered.
The energy dependent ones include the systematic uncertainties from the choosing about mass window and sideband regions of K 0 S , π 0 , and η and the kinematic fit. The uncertainties associated with the K 0 S , π 0 , and η invariant mass regions are determined by changing them from (0.488, 0.508) to (0.483, 0.513) GeV/c 2 , (0.12, 0.15) to (0.115, 0.155) GeV/c 2 and (0.52, 0.58) to (0.51, 0.59) GeV/c 2 for the K 0 S , π 0 and η, respectively. The differences in the efficiencies are taken as the corresponding systematic uncertainties.
The uncertainty associated with the kinematic fit is determined by comparing the efficiencies with and without corrections to the track helix parameters [73].
Assuming all sources of systematic uncertainties are independent, the total uncertainties are the sums of the individual values in quadrature (Table IV).
The systematic uncertainties that affect the upper limits on σ B e + e − → πZ c (3900), Z c (3900) → K 0 S Kπ/η are considered in two categories: multiplicative and nonmultiplicative. The non-multiplicative systematic uncertainties on the signal shape and the background shape are considered by changing the signal shape to a Breit-Wigner function and varying the fit range, the parameters of the Z c (3900), and the order of the polynomial functions in the fit. The maximum upper limits are adopted for all combinations of these variations. The intermediate states in the Z c (3900) decay are considered by generating signal MC samples with alternative processes Z c (3900) → K * (892)K, K * (892) → K(K 0 S )π (π 0 mode), and Z c (3900) → a 0 (980)η, a 0 (980) → K 0 S K (η mode). The efficiency difference is considered as a multiplicative systematic uncertainty. All of the systematic uncertainties, which are listed in Table IV, excluding the side-band item and mixing MC item, are considered as the multiplicative systematic uncertainties. The effects of multiplicative systematic uncertainties are taken into account by convolving the distribution of L(σ B ) with a probability distribution function of sensitivity (S), which is assumed to be a Gaussian function with central valuê S and standard deviation δ S [74]: Here S is the sensitivity that refers to the denominator of Eq. (1) and δ s is the total multiplicative systematic uncertainty. L ′ (σ B ) is the likelihood distribution of the Born cross section after the multiplicative systematic uncertainties are incorporated.

IV. SUMMARY
The Born cross sections for e + e − → K 0 S K ± π ∓ π 0 and K 0 S K ± π ∓ η are measured with data samples collected at center-of-mass energies from 3.90 to 4.60 GeV. Since no clear structure is observed, the upper limits of the product Γ e + e − B Y (4260) →K 0 S K ± π ∓ π 0 at 90% C.L. is estimated to be less than 0.05 eV and that of Γ e + e − B Y (4260) →K 0 S K ± π ∓ η is estimated to be smaller than 0.19 eV. Ref. [75] reported four solutions of the product Γ e + e − B Y (4260) → π + π − J/ψ ,  in which the maximum is 13.3 ± 1.4 eV and the minimum is 1.5 ± 0.3 eV. Comparing them with our results, the branching fraction of the Y (4260) decaying into K 0 S K ± π ∓ π 0 and K 0 S K ± π ∓ η is much smaller, which indicates a much smaller coupling of the Y (4260) to the light hadrons K 0 S K ± π ∓ π 0 and K 0 S K ± π ∓ η. We also search for e + e − → πZ c (3900), Z c (3900) → K 0 S Kπ/η and no obvious Z c (3900) signal is observed in the charged nor neutral mode. The 90% C.L. upper limits on the cross sections are given at √ s = 4.226, 4.258, 4.358, 4.416, and 4.600 GeV. The absence of a signal suggests that the cross sections for light hadron decay modes are small and that the annihilation of cc in the Y (4260) and Z c (3900) is suppressed. Additional exploration of light hadron decay modes is needed to confirm the hypotheses.