Observation of $\chi_{cJ}\to4K_{S}^{0}$

By analyzing $(448.1\pm2.9)\times10^6$ $\psi(3686)$ events collected with the BESIII detector operating at the BEPCII collider, the decays of $\chi_{cJ} \to 4K_S^0$ ($J=0,1,2$) are observed for the first time with statistical significances of 26.5$\sigma$, 5.9$\sigma$ and 11.4$\sigma$, respectively. The product branching fractions of $\psi(3686)\to\gamma\chi_{cJ}$, $\chi_{cJ}\to4K_{S}^{0}$ are presented and the branching fractions of $\chi_{cJ}\to4K_{S}^{0}$ decays are determined to be $\mathcal{B}_{\chi_{c0}\to 4K_S^0}$=$(5.76\pm0.34\pm0.38)$$\times10^{-4}$, $\mathcal{B}_{\chi_{c1}\to 4K_S^0}$=$(0.35\pm0.09\pm0.03)$$\times10^{-4}$ and $\mathcal{B}_{\chi_{c2}\to 4K_S^0}$=$(1.14\pm0.15\pm0.08)$$\times10^{-4}$, where the first uncertainties are statistical and the second are systematic, respectively.


I. INTRODUCTION
In the quark model, the χ cJ (J = 0, 1, 2) mesons are the 3 P J charmonium states. Since the χ cJ mesons cannot be directly produced in e + e − collisions, according to parity conservation, their decays are experimentally and theoretically not studied as extensively as the vector charmonium states J/ψ and ψ(3686). However, the χ cJ mesons can be produced in radiative decays of the ψ(3686) with branching fractions of about 9%, which provides a method to produce large χ cJ samples in order to study χ cJ decays.
Recent theoretical work indicates that the Color Octet Mechanism (COM) [1] could have large contributions to the decays of the P -wave charmonium states. However, many contradictions still exist between these theoretical calculations and experimental measurements. For instance, theoretical predictions of χ cJ decays to baryon anti-baryon pairs based on the COM [2][3] [4] were inconsistent with experimental measurements [5]. Thus more precise experimental results are mandatory to further understand χ cJ decay dynamics. Furthermore, the χ c0 and χ c2 states are expected to decay via two-gluon processes into light hadrons, giving access to the investiga- tion of glueball dynamics. Thus, comprehensive measurements of exclusive hadronic decays of χ cJ are valuable.
For the decay modes of χ cJ → 4K, the branching fractions of χ cJ decays into 2(K + K − ) and K + K − K 0 S K 0 S have been measured by Belle [6] and BES [7] with results summarized in Table I. In this paper, by analyzing (448.1 ± 2.9) × 10 6 ψ(3686) events [8] collected with the BESIII detector [9], we present the first measurements of the branching fractions of χ cJ decays to 4K 0 S .

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector is operated at the Beijing Electron Positron Collider II (BEPCII), which has reached a peak luminosity of 1.0 × 10 33 cm −2 s −1 at a center-of-mass energy of √ s = 3.773 GeV. The detector has a geometrical acceptance of 93% of the solid angle and is composed of four main components. A heliumgas based main drift chamber (MDC) is used to track charged particles. The single wire resolution is better than 130 µm, which, together with a magnetic field of 1 T, leads to a momentum resolution of 0.5% at 1 GeV/c. The energy loss per path length dE/dx is measured with a resolution of 6%. The MDC is surrounded by a time-offlight system built from plastic scintillators. It provides a 2σ K/π separation up to 1 GeV/c momentum with a time resolution of 80 (110) ps for the barrel (end-caps). Particle energies are measured in the CsI(Tl) electro-magnetic calorimeter (EMC), which achieves an energy resolution for electrons of 2.5% (5%) at 1 GeV/c momentum and a position resolution of 6 mm (9 mm) for the barrel (endcaps). Outside of the magnet coil, a muon counter composed of resistive plate chambers provides a spatial resolution of better than 2 cm. A more detailed description of the detector can be found in Ref. [9].
A GEANT4 [10] based Monte Carlo (MC) simulation package is used to optimize the event selections and estimate the signal efficiency and the background level. The event generator KKMC [11] simulates the electronpositron annihilation and the production of the ψ resonances. Particle decays are generated by EVTGEN [12] for the known decay modes with branching fractions from the Particle Data Group (PDG) [5] and Lundcharm [13] for the unknown ones. An inclusive MC sample contain-ing 506 × 10 6 generic ψ(3686) decays is used to study background. The signal MC samples of the χ cJ decays, generated according to a phase space model, are used to determine efficiencies.

III. EVENT SELECTION
We reconstruct events from the decay chain of the charmonium transitions ψ(3686) → γχ cJ followed by the hadronic decays χ cJ → 4K 0 S and K 0 S → π + π − . A photon candidate is defined as a shower detected within the EMC exceeding an energy deposit of 25 MeV in the barrel region (covering the region | cos θ| < 0.8, where θ is the polar angle with respect to the positron beam direction) or of 50 MeV in the end-caps (0.86 < | cos θ| < 0.92). To suppress the electronics noise and beam background, the clusters are required to start within 700 ns after the event start time and fall outside a cone angle of 10 • around the nearest extrapolated charged track. All charged tracks are required to originate from the interaction region defined as |V z | < 20 cm and | cos θ| < 0.93, where V z denotes the distance of the closest approach of the reconstructed track to the interaction point (IP) in the z direction. Candidate events must have eight charged tracks with zero net charge and at least one good photon. The K 0 S candidates are reconstructed using vertex fits by looping over all oppositely charged track pairs in an event (assuming the tracks to be π ± without particle identification). To suppress the π + π − combinatorial background, the reconstructed decay lengths (L) of the K 0 S candidates are required to be more than twice their standard deviations (σ L ). The distribution of L/σ L for all K 0 S candidates is shown in Fig. 1. The invariant mass of π + π − (M π + π − ) must be within the K 0 S signal region, defined as 12 MeV/c 2 around the K 0 S nominal mass [5]. The M π + π − distribution for all K 0 S candidates is shown in  a four-momentum conservation constraint (4C) is applied to the events. The χ 2 4C of the kinematic fit is required to be less than 200. The spectrum of the invariant mass of the 4K 0 S (M 4K 0 S ) of the accepted candidate events is shown in Fig. 3. Clear χ c0 , χ c1 and χ c2 signals are observed.
The M 4K 0 S distribution is fitted using an unbinned maximum likelihood fit. In the fit, each χ cJ signal is described with the MC simulation shape convolved with a Gaussian function with free parameters. Since the background level is very low, as discussed in Sec. IV, the background shape is assumed to be flat. The signal yields of χ c0 , χ c1 and χ c2 are fitted to be 319.4 ± 19.0, 21.6 ± 5.2 and 68.0 ± 8.7, respectively. The statistical significances are estimated to be 26.5σ, 5.9σ and 11.4σ for χ c0 , χ c1 and χ c2 individually, which are determined by comparing the fit likelihood values with and without each χ cJ signal separately. S . The points with error bars are data, the blue curve is the overall fit, and the red curve is the fitted background.

IV. BACKGROUND STUDIES
The continuum data taken at √ s =3.65 GeV, corresponding to an integrated luminosity of 44.45 pb −1 [14], are used to estimate the QED background. No events within this sample satisfy the same selection criteria applied to the main data sample. In addition, the inclusive MC sample is used to study all potential backgrounds from ψ(3686) decays. Two background events are found to be from ψ(3686) →K * K 0 S f ′ 2 andK * K 0 S f 0 (1710). Further studies with large exclusive MC samples show that the two background sources make up a uniform distribution around the χ cJ signal regions. So, all peaking background components are negligible in this analysis.

V. BRANCHING FRACTIONS
For each decay of ψ(3686) → γχ cJ , χ cJ → 4K 0 S , K 0 S → π + π − , 5 × 10 5 signal MC events are generated using a 1 + λ cos 2 θ distribution, where θ is the angle between the direction of the radiative photon and the beam, and λ = 1, −1/3, 1/13 for J = 0, 1, 2 in accordance with expectations for electric dipole transitions. Since no obvious substructures are found in the M 2K 0 S and M 3K 0 S distributions of the accepted χ cJ → 4K 0 S candidate events, as shown in Fig. 4, the χ cJ decay products are generated using phase space (PHSP). Intrinsic width and mass values as given in Ref. [5] are used for the χ cJ states in the simulation. To reduce the difference of the distributions of χ 2 of the 4C kinematic fit (χ 2 4C ) between data and MC simulation, we correct the track helix parameters of MC simulation in the 4C kinematic fit. The χ 2 4C distribution after corrections is shown in Fig. 5, in which the consistency between data and MC simulation is reasonable. The obtained corrected efficiencies for χ cJ → 4K 0 S are (5.51 ± 0.03)%, (6.19 ± 0.04)% and (6.08 ± 0.04)%, respectively, including detector acceptance as well as reconstruction and selection efficiencies.
The signal yields N J obs are obtained by fitting to the M 4K 0 S distribution. The branching fraction is calculated with where ǫ is the efficiency, N ψ(3686) is the number of ψ(3686) events, B ψ(3686)→γχcJ and B K 0 S →π + π − are the branching fractions of the PDG fit of ψ(3686) → γχ cJ decays and K 0 S → π + π − decay [5].

VI. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties in the measurements of B χcJ →4K 0 S originate from several sources, as summarized in Table II. They are estimated and described below.  The number of ψ(3686) events has been measured to be N ψ(3686) = (448.1 ± 2.9)× 10 6 with the inclusive hadronic data sample, as described in Ref. [8]. The uncertainty of the total number is 0.6%. The systematic uncertainty due to the photon detection is assumed to be 1.0% per photon with the control sample J/ψ → ρ 0 π 0 [15].
The systematic uncertainty associated with K 0 S reconstruction is determined to be 1.5% per K 0 S with the control samples of J/ψ → K * ± (892)K ∓ , K * ± (892) → K 0 S π ± and J/ψ → φK 0 S K ∓ π ± in Ref. [16]. To estimate the systematic uncertainties of the MC model for the χ cJ → 4K 0 S decay, we compare our nominal efficiency with that determined from the signal MC events after mixing some possible subresonant decays, including χ cJ → f 0 (1500)f 0 (1500), The systematic uncertainties are estimated as the relative changes of efficiencies, which are 0.4%, 0.2% and 0.2% for χ c0 , χ c1 and χ c2 decays, respectively.
We correct the track helix parameters for MC simulation in the 4C kinematic fit. The change in detection efficiency is not more than 1.0% when varying the correction factors within one standard deviation around the nominal value. We therefore assume 1.0% as the systematic uncertainty of the 4C kinematic fit.
To estimate the systematic uncertainties in the angular distribution, we use a reweighting method. New signal MC events are obtained by reweighting the angular distribution of the K 0 S in the signal MC events to data. The changes to the detection efficiencies are taken as the systematic uncertainties, which are 0.7%, 0.5% and 0.7% for χ c0 , χ c1 and χ c2 decays, respectively.
The systematic uncertainties due to the fit range are estimated by a series of fits with alternative intervals. The standard deviations of the resulting branching fractions are assigned as the systematic uncertainties, which are 0.6%, 1.5% and 0.9% for χ c0 , χ c1 and χ c2 decays, respectively.
To estimate the systematic uncertainties due to the signal shape, we use alternative signal shapes, a Breit Wigner function smeared with a double Gaussian and a MC shape (including E 3 dependence) convolved with a Gaussian function, to describe each χ cJ signal. The maximum deviations of the resulting branching fractions are assigned as the relevant systematic uncertainties, which are 0.4%, 2.8% and 1.7% for χ c0 , χ c1 and χ c2 decays, respectively.
The systematic uncertainties from the branching fractions of ψ(3686) → γχ cJ and K 0 S → π + π − decays quoted from the PDG [5] are 2.0%, 2.5% and 2.1% for χ c0 , χ c1 and χ c2 decays and 0.07% for K 0 S , respectively. We assume that all systematic uncertainties are independent and add them in quadrature to obtain the total systematic uncertainty for each decay.

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