Measurement of the branching fraction for $\psi(3686)\to \omega K^0_SK^0_S$

Analyzing $(448.1\pm2.9)\times10^6$ $\psi(3686)$ events collected with the BESIII detector at the BEPCII collider, the $\psi(3686)\to \omega K_{S}^{0}K_{S}^{0}$ decay is observed for the first time. The branching fraction for this decay is determined to be $\mathcal{B}_{\psi(3686)\to \omega K_{S}^{0}K^{0}_{S}}$=$(7.04\pm0.39\pm0.36)$$\times10^{-5}$, where the first uncertainty is statistical and the second is systematic.


I. INTRODUCTION
Experimental studies of charmonium decays are an important tool to study Quantum Chromodynamics (QCD) [1]. In perturbative-QCD (pQCD), both J/ψ and ψ(3686) decays to light hadrons mainly proceed through annihilation of the cc pair into three gluons, with the decay width proportional to the square of the charmonium wave function at the origin [2]. One obtains the pQCD 12% rule for the ratio of branching fractions of ψ(3686) and J/ψ decays to light hadrons: A violation of this rule was first observed by Mark-II in the J/ψ → ρπ decay [3] and subsequently found in other decay channels. The inconsistency between pQCD predictions and experimental results is known as the ρπ puzzle. Various mechanisms have been proposed to explain the observed discrepancies. However, no single model provides an explanation of all available experimental results (for a review see Ref. [4]). Additional measurements of individual branching fractions of J/ψ and ψ(3686) decays provide complementary information to understand the ρπ puzzle and further investigate charmonium decay mechanisms [1]. In recent years, a lot of progress has been made in experimental studies of multi-body J/ψ and ψ(3686) decays. For example, the J/ψ → ωKK decays have been studied in Refs. [5,6] and the ψ(3686) → ωK + K − decay has been investigated in Refs. [7][8][9][10]. In Ref. [10], Q ωK + K − for the ψ → ωK + K − decays is determined to be (18.4±3.7)%, indicating no significant violation of the 12% rule. To date, however, no study of ψ(3686) → ωK 0 S K 0 S has been reported. A first measurement of the branching fraction of ψ(3686) → ωK 0 S K 0 S allows to test whether the 12% rule holds for ψ → ωK 0 S K 0 S decays. The branching fraction of ψ → φ(ω)K 0K 0 is expected to be equal to that of ψ → φ(ω)K + K − based on isospin symmetry (neglecting possible effects related to contin-uum production). Earlier measurements summarized in Ref. [11] roughly support this relation. Comparing the obtained branching fraction for ψ(3686) → ωK 0 S K 0 S to that for ψ(3686) → ωK + K − offers an opportunity to search for isospin symmetry violation in ψ(3686) decays.
In this paper, we present the first observation and branching fraction measurement of ψ(3686) → ωK 0 S K 0 S . This analysis is carried out using (448.1 ± 2.9) × 10 6 ψ(3686) events collected at the center-of-mass energy of 3.686 GeV with the BESIII detector in 2009 and 2012 [12]. An inclusive ψ(3686) MC sample, consisting of 506 × 10 6 events, is used to estimate potential backgrounds from ψ(3686) decays.
The detection efficiency is determined using 5×10 5 ψ(3686) → ωK 0 S K 0 S decays modeled to reproduce the kinematic distributions found in the data. The ω decays are described with a Dalitz decay model [13]. A data sample collected at the center-of-mass energy of 3.65 GeV with an integrated luminosity of 43.88 pb −1 is used to estimate the contamination from continuum processes.

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector [14] records symmetric e + e − collisions provided by the BEPCII storage ring [15], which operates with a peak luminosity of 1 × 10 33 cm −2 s −1 in the center-of-mass energy range from 2.0 to 4.95 GeV . BESIII has collected large data samples in this energy region [1]. The cylindrical core of the BESIII detector covers 93% of the full solid angle and 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 identification modules interleaved with steel. The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for 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 in the TOF barrel region is 68 ps, while that in the end cap region is 110 ps.
Simulated samples produced with the geant4based [16] 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 estimate the backgrounds.
The simulation includes the beam energy spread and initial state radiation (ISR) in the e + e − annihilations modeled with the generator kkmc [17].
The inclusive MC sample consists of the production of the charmonium resonances, and the continuum processes incorporated in kkmc [17]. The known decay modes are modeled with evtgen [18] using branching fractions taken from the Particle Data Group (PDG) [11], and the remaining unknown decays from the charmonium states are modeled with lundcharm [19]. Final state radiation (FSR) from charged final-state particles is incorporated with the photos package [20].

III. EVENT SELECTION AND DATA ANALYSIS
The two charged pion candidates from the ω → π + π − π 0 decay are required to have a polar angle θ with respect to the symmetry axis of the detector of |cosθ| < 0.93. Furthermore, they are required to originate from the interaction point (IP). The distance of closest approach to the IP must be less than 10 cm along the z-axis, and less than 1 cm in the transverse plane. Particle identification (PID) for charged pion and kaon candidates is performed using the dE/dx and TOF information. The confidence levels for pion and kaon hypotheses (CL π and CL K ) are calculated. Pion candidates are required to satisfy CL π > CL K .
The K 0 S candidates are reconstructed from two oppositely charged tracks, which are assumed to be pions without imposing PID criteria. The two charged tracks from the K 0 S → π + π − decay must satisfy |cosθ| < 0.93. Their distance of closest approach to the IP along the MDC axis direction is required to be less than 20 cm. A secondary vertex fit is performed for both K 0 S candidates. The chosen π + π − pairs are constrained to originate from a common secondary vertex. Their invariant mass is required to be in the range (0.486, 0.510) GeV/c 2 , which corresponds to an interval of about three times the fitted resolution around the K 0 S nominal mass. The decay length of the K 0 S candidate is required to be at least two standard deviations of the vertex resolution away from the IP.
The π 0 candidates are reconstructed from photon pairs reconstructed from EMC showers. Each EMC shower must start within 700 ns of the event start time. Its energy is required to be greater than 25 MeV in the barrel region (|cosθ| < 0.8) and 50 MeV in the end cap region (0.86 < |cosθ| < 0.92). The opening angle between the photon candidate and the extrapolated impact point of the nearest charged track in the EMC must be greater than 10 • . Each photon pair with an invariant mass within the range (0.115, 0.150) GeV/c 2 is used as a π 0 candidate.
A five-constraint kinematic fit for the e + e − → π + π − π + π − π + π − π 0 hypothesis is performed including energy-momentum conservation (4C) and with the invariant mass of the two photons constrained to the mass of the π 0 (1C). The helix parameters of charged tracks of the MC events are corrected to improve consistency between data and MC simulation (see Ref. [21] for details). Events satisfying χ 2 5C < 80 are kept for further analysis. If there are multiple candidates in an event, the one with the smallest value of χ 2 5C + χ 2 S 2 represent the fit-quality of the secondary vertex fits for the two K 0 S candidates. Studies of the inclusive ψ(3686) MC sample using a generic event type analysis tool, TopoAna [22], show that those background events surviving all selection criteria contain contributions including an ω → π + π − π 0 decay. These events, however, do not contain a K 0 S K 0 S pair and can thus be described by K 0 S K 0 S sideband events. The one-dimensional K 0 S signal and sideband regions are defined as (0.486, 0.510) and (0.454, 0.478) or (0.518, 0.542) GeV/c 2 , respectively.
The twodimensional (2D) K 0 S K 0 S signal region (called 2D signal region) is defined as the rectangle with both π + π − combinations lying in the K 0 S signal regions. The 2D sideband regions (sideband 1 and sideband 2) are defined based on the one-dimensional sideband as shown in Fig. 1(b).  Fig. 1. (a) The M π + π − distribution of the K 0 S candidates in data and the scaled signal and background contributions from the inclusive MC sample, where the pairs of red and pink arrows show the 1D K 0 S signal and sideband regions, respectively. (b) The M π + π − (1) versus M π + π − (2) distribution of the K 0 S candidates in data, where the red solid, pink dashed, and green long-dashed rectangles denote the 2D signal region, sideband 1 region, and sideband 2 region, respectively. Besides all requirements mentioned in text, the events displayed here have been required to be within the region (0.777, 0.807) GeV/c 2 in the invariant mass of the π + π − π 0 system. Figure 2 shows the M π + π − π 0 distributions of the candidate events for e + e − → π + π − π 0 K 0 S K 0 S for the signal and sideband regions. In these figures, a clear ω signal is observed. Unbinned maximum likelihood fits are performed to these events. In the fits, the ω signal shape is described by a Breit-Wigner function convolved with a Gaussian function and the background shape is described by a first order Chebychev polynomial. The width of the Breit-Wigner function is fixed to the world average value of the ω meson [11]. The mass of the ω meson as well as the mean and width of the Gaussian are free parameters. In the fits to the events in the 2D sideband 1 and sideband 2 regions, the parameters of the Gaussian function are fixed to the values obtained from the fit to the 2D signal region. The results of the best fits to the data are shown in Fig. 2.
From these fits, we obtain the numbers of e + e − → ωK 0 S K 0 S in the 2D K 0 S signal, sideband 1 and sideband 2 region, denoted as N sig , N side1 , and N side2 , respectively, where N side2 represents the yield of background events with two mis-formed K 0 S s, N side1 indicates that with one or two mis-formed K 0 S s, and N sig is the number of events with two successfully formed K 0 S s. Assuming that the combinatorial backgrounds distribute evenly in the 2D spectrum and in the projected distribution within K 0 S signal region, the total number of e + e − → ωK 0 S K 0 S candidate events in the ψ(3686) data is calculate by The obtained values of N sig , N side1 , N side2 , and N tot in data are summarized in the first row of Table 1.
A similar analysis and fit are performed for the continuum data at √ s = 3.65 GeV. The parameters of the Gaussian function are fixed to the values obtained from the fit to the ψ(3686) data. The results of the best fit to the data are shown in Fig. 3. The obtained values of N sig , N side1 , N side2 , and N tot in the continuum data are summarized in the second row of Table 1.  Figure 4 shows the distributions of cos θ, where θ is the polar angle in the lab frame, and the momenta of final state particles as well as the invariant mass spectra of all two-particle combinations for those candidates surviving all selection criteria in data and MC simulation. Figure 5 shows the Dalitz plots of the accepted candiates in data and BODY3 signal MC events. The K 1 (1270) resonance is observed in the ωK 0 S mass spectrum in the data. To determine the efficiency ( ψ(3686)→ωK 0 S K 0 S ) for the signal process, the ψ(3686) → ωK 0 S K 0 S decays are simulated with the modified data-driven generator BODY3 [18], which was developed to simulate different intermediate states in data for a given three-body final state. The Dalitz plot of M ωK 0 S (1) vs. M ωK 0 S (2) found in data, including a bin-wise correction for backgrounds and efficiencies, is taken as an input to the BODY3 generator. The consistency between data and BODY3 signal MC events is much better than that between data and phase space (PHSP) signal MC events. The BODY3 signal MC events are used to determine the signal selection efficiency in the further analysis. With the same analysis procedure as data analysis, the detection efficiency for ψ(3686) → ωK 0 S K 0 S is determined to be (11.23 ± 0.05)%, where the branching fractions of ω → π + π − π 0 and K 0 S → π + π − are not included. The reliability of the detection efficiency is validated due to good consistency between data and MC simulation, as shown in Fig. 4.

IV. SYSTEMATIC UNCERTAINTY
In this section, the systematic uncertainties of the measurement of the ψ(3686) → ωK 0 S K 0 S branching fraction are discussed.
The total number of ψ(3686) events and its relative uncertainty are taken from Ref. [12]. Its uncertainty of 0.65% is considered as a systematic uncertainty. The statistical uncertainty is negligible.
The π ± tracking efficiency was studied using a control sample of ψ(3686) → π + π − J/ψ with J/ψ → + − . The difference in the tracking efficiencies between data and MC simulation is 1% per π ± . The π ± PID efficiency was investigated with a control sample of J/ψ → ρπ [23]. The difference in the PID efficiencies between data and MC simulation is 1% per π ± .
The systematic uncertainty originating from the π 0 mass window is examined by comparing the branching fractions measured with and without this requirement. The change of the obtained branching fraction, 0.7%, is taken as an uncertainty.
The K 0 S reconstruction efficiency, including the efficiencies of tracking for π + π − , the M π + π − mass requirement, and the decay length requirement was studied using control samples of J/ψ → K * (892) ± K ± , K * (892) ± → K 0 S π ± and J/ψ → φK 0 S K ∓ π ± [24]. The difference in the K 0 S reconstruction efficiencies between data and MC simulation is 1.6% per K 0 S . To estimate the systematic uncertainty related to the K 0 S sideband, we examine the branching fraction using a K 0 S sideband region shifted by ±4 MeV/c 2 . The maximum change of  the branching fraction of 0.3% is taken as the systematic uncertainty.
For the fit to the M π + π − π 0 spectrum, three potential sources of systematic uncertainties are considered. First, the uncertainty from the background shape is estimated using an alternative background shape determined from the inclusive MC sample. The change of the fitted signal yield, 0.1%, is assigned as an uncertainty. Second, the uncertainty from the signal shape is estimated by varying the fixed width of the ω meson by ±1σ, where σ is the uncertainty of the world average value. The change of the fitted signal yield is negligible. Third, the uncertainty due to the fit range is estimated using fit ranges of (0.55, 0.95), (0.60, 1.00), (0.65, 0.95), and (0.65, 0.90) GeV/c 2 instead. The maximum change of the fitted signal yield, 0.8%, is taken as an uncertainty. Adding these three uncertainties in quadrature leads to the systematic uncertainty due to the fit to the M π + π − π 0 spectrum of 0.8%.
The systematic uncertainty of the signal modeling with the BODY3 generator is estimated by varying the bin size of the input Dalitz plot by ±20%. The maximum change of the signal efficiency, 0.8%, is taken as the corresponding uncertainty.
The systematic uncertainty due to the 5C kinematic fit is evaluated by varying the helix parameters by ±1 standard deviations. The change of signal efficiency, 0.9%, is taken as the relevant systematic uncertainty. The systematic uncertainty related to the scale factor f c for the continuum production of e + e − → ωK 0 S K 0 S is estimated by comparing the difference between 1/s and 1/s 3 dependencies of the cross section of the process e + e − → ωK 0 S K 0 S . The effect on the branching fraction, 0.6%, is taken as the corresponding systematic uncertainty.
The uncertainties of the quoted branching fractions for ω → π + π − π 0 and K 0 S → π + π − are 0.8% and 0.1%, respectively. The uncertainty due to the limited MC statistics is considered as a source of systematic uncertainty.
Assuming that all sources are independent, the total systematic uncertainty is determined to be 5.2%. The systematic uncertainties are summarized in Table 2.

VI. ACKNOWLEDGEMENT
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 Research and Development Program of China under Contracts Nos. 2020YFA0406300, 2020YFA0406400; National Natural Science Foundation of China (NSFC) under Contracts Nos. 11875170, 11775230, 11625523,  11635010, 11735014, 11822506, 11835012, 11935015,