First observation of the decay $\chi_{cJ}\to \Sigma^{+}\bar{p}K_{S}^{0}+c.c. ~(J = 0, 1, 2)$

Using E1 radiative transitions $\psi(3686) \to \gamma\chi_{cJ}$ from a sample of $(448.1 \pm 2.9)\times10^{6}$ $\psi(3686)$ events collected with the BESIII detector, the decays $\chi_{cJ}\to \Sigma^{+}\bar{p}K_{S}^{0}+c.c.~(J = 0, 1, 2)$ are studied. The decay branching fractions are measured to be $\mathcal{B}(\chi_{c0}\to \Sigma^{+}\bar{p}K_{S}^{0}+c.c.) = (3.52 \pm 0.19\pm 0.21)\times10^{-4}$, $\mathcal{B}(\chi_{c1}\to \Sigma^{+}\bar{p}K_{S}^{0}+c.c.) = (1.53 \pm 0.10\pm 0.08)\times10^{-4}$, and $\mathcal{B}(\chi_{c2}\to \Sigma^{+}\bar{p}K_{S}^{0}+c.c.) = (8.25 \pm 0.83 \pm 0.49)\times10^{-5}$, where the first and second uncertainties are the statistical and systematic ones, respectively. No evident intermediate resonances are observed in the studied processes.


I. INTRODUCTION
The first charmonium states with J P C = J ++ discovered after the J/ψ and ψ(3686) were the χ cJ (J = 0, 1, 2) particles. Quarkonium systems, especially charm anti-charm states, are regarded as a unique laboratory to study the interplay between perturbative and nonperturbative effects in quantum chromodynamics (QCD). Experimental studies of charmonium decays can test QCD and QCD-based effective field theory calculations. The χ cJ states belong to the charmonium P -wave spin triplet, and therefore cannot be produced via a single virtual-photon exchange in electron-positron annihilations as are the J/ψ and ψ(3686). Until now the understanding of these states has been limited by the availability of experimental data. The world's largest data set of ψ(3686) events [1] collected with the BESIII [2] detector, provides a unique opportunity for detailed studies of χ cJ decays, since they are copiously produced in ψ(3686) radiative transitions with branching fractions of about 9% each [3].
Many excited baryon states have been discovered by BaBar, Belle, CLEO, BESIII, and other experiments in the past decades [3], but the overall picture of these states is still unclear. While many predicted states have not yet been observed, many states that do not agree with quark model predictions are observed (for a review see Ref. [4]). Therefore the search for new excited baryon states is important to improve knowledge of the baryon spectrum and the understanding of the underlying processes which describe confinement in the nonperturbative QCD regime. Experimentally, exclusive decays of χ cJ to baryon anti-baryon (BB) pairs, such as pp, ΣΣ, ΛΛ [5][6][7][8], have been investigated. How-ever, there are only a few experimental studies of χ cJ to BBM (M stands for meson). These channels are ideal to search for new excited baryons in intermediate states, which decay intoBM and BM . This paper reports the first measurements of the branching fractions of χ cJ → Σ +p K 0 S + c.c. via the radiative transition ψ(3686) → γχ cJ , where Σ + → pπ 0 , K 0 S → π + π − , and π 0 → γγ. The chargeconjugate state (c.c.) is included unless otherwise stated. We also report on a search for possible excited baryon states in the invariant-mass spectra of pK 0 S , and Σ + K 0 S .

II. BESIII DETECTOR
The BESIII detector is a magnetic spectrometer located at the Beijing Electron Positron Collider (BEPCII) [9]. The cylindrical core of the BE-SIII detector consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-offlight 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 4π solid angle. The charged-particle momentum resolution at 1 GeV is 0.5%, and the dE/dx resolution is 6% for the electrons from Bhabha scattering at 1 GeV. 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, while that of the end-cap part is 110 ps.

III. DATA SET AND MONTE CARLO SIMULATION
This analysis is based on a sample of (448.1 ± 2.9) × 10 6 ψ(3686) events [1] collected with the BE-SIII detector.
geant4-based [10] Monte Carlo (MC) simulation data are used to determine detector efficiency, optimize event selection, and estimate background contributions. Inclusive MC samples were produced to determine contributions from dominant background channels. The production of the initial ψ(3686) resonance is simulated by the MC event generator kkmc [11,12], and the known decay modes are modeled with evtgen [13,14] using the branching fractions summarized and averaged by the Particle Data Group (PDG) [3], while the remaining unknown decays are generated by lundcharm [15]. The final states are propagated through the detector system using geant4 software.

IV. DATA ANALYSIS
For the reaction channel ψ(3686) → γχ cJ , with χ cJ → Σ +p K 0 S , Σ + → pπ 0 , π 0 → γγ, and K 0 S → π + π − , the final-state particles are ppπ + π − γγγ. Charged tracks must be in the active region of the MDC, corresponding to | cos θ| < 0.93, where θ is the polar angle of the charged track with respect to the beam direction. For the anti-proton (p), the point of closest approach to the interaction point (IP) must be within ±1 cm in the plane perpendicular to the beam (R xy ) and ±10 cm along the beam direction (V z ). Due to the long lifetime of the K 0 S and Σ + , there is no requirement on R xy or V z for the track candidates used to form the K 0 S or Σ + candidates. Photon candidates are reconstructed by summing the energy deposition in the EMC crystals produced by the electromagnetic showers. The minimum energy necessary for counting a photon as a photon candidate is 25 MeV for barrel showers (| cos θ| < 0.8) and 50 MeV for end-cap showers (0.86 < | cos θ| < 0.92). To eliminate showers originating from charged particles, a photon cluster must be separated by at least 10 • from any charged track. The timing of the shower is required to be within 700 ns from the reconstructed event start time to suppress noise and energy deposits unrelated to the event. Events with two positively charged tracks, two negatively charged tracks, and at least three good photons are selected for further analysis. The TOF (both end-cap and barrel) and dE/dx measurements for each charged track are used to calculate the p-value based on the χ 2 PID values for the hypotheses that a track is a pion, kaon, or proton. Two oppositely charged tracks are identified as a proton/anti-proton pair if their proton hypothesis pvalues are greater than their kaon or pion hypothesis p-values. The remaining charged tracks are considered as pions by default. The numbers of protons and anti-protons as well as the negatively and positively charged pions should be equal to one.
The K 0 S candidate is reconstructed with a pair of oppositely charged pions. To suppress events from combinatorial background contributions, we require that the π + π − pair is produced at a common vertex [18].
Next a four-constraint (4C) kinematic fit imposing energy-momentum conservation is performed under the ppπ + π − γγγ hypothesis. If there are more than three photon candidates in an event, the combination with the smallest χ 2 4C is retained, and its χ 2 4C is required to be less than those for the ppπ + π − γγ and ppπ + π − γγγγ hypotheses. The value of χ 2 4C is required to be less than 50. For the selected signal candidates, the γγ combination (γ 1 γ 2 ) with an invariant mass closest to the π 0 mass is reconstructed as π 0 candidate, and the remaining one (γ 3 ) is considered to be the radiative photon from the ψ(3686) decay. The γγ invariant mass is required to satisfy |M γγ −m π 0 | < 15 MeV/c 2 . Here and throughout the text, M i represents a measured invariant mass and m i represents the nominal mass of the particle(s) i [3]. To reduce background events withΛ →pπ + , |Mp π + − m Λ | > 6 MeV/c 2 is required. Figure 1 shows the scatter plot of the π + π − invariant mass versus the pπ 0 invariant mass of data. To select events which contain both a K 0 S and a Σ + candidate, Fig. 1). The widths of the mass intervals are chosen to be 3 times the invariant-mass resolution.
A hint of a structure in the invariant-mass distribution of thepK 0 S subsystem in the χ c0 signal region can be seen in Fig. 3(a). Considering the width and mass, it is most likely theΣ(1940) − with M = 1940 MeV/c 2 , Γ = 220 MeV, and I(J P ) = 1( 3 2 − ) [3]. Other excited Σ * states are most likely excluded because their widths are much larger. For the fit to the invariant-mass distribution Mp K 0 S , several contributions are considered, namely the line shape from the phase-space model, the normalized K 0 S and Σ + mass sidebands in the χ c0 signal region (described in detail in the background analysis), and theΣ(1940) − signal from the MC simulation, where the mass and width ofΣ(1940) − are fixed to the world average values [3]. To estimate the statistical signal significance of theΣ(1940) − contribution, we use the quantity −2 ln(L 0 /L max ), where L 0 and L max are the likelihoods of the fits without and withΣ(1940) − signal, respectively. The statistical significance of thē Σ(1940) − signal is obtained to be 3.2σ. The signal significance is reduced to 2.3σ if the width of Σ(1940) − is taken as the lower value of 150 MeV [3]. For all other invariant-mass distributions of the twobody subsystems, the description using the phasespace model is in good agreement with data. For example, thepK 0 S mass distributions from data and MC simulations in the χ c1 and χ c2 signal regions are shown in Figs. 3(b) and 3(c). Possible background contributions are studied with the inclusive MC sample of 5.06 × 10 8 simulated ψ(3686) decays. Peaking background contributions in the χ cJ mass regions are dominated by the channels χ cJ →∆ − π + ∆ 0 (∆ − → pπ 0 , ∆ 0 → pπ − ) and χ cJ → ppρ + π − (ρ + → π + π 0 ). Other background events, mainly from the channels and ψ(3686) → J/ψπ 0 π 0 (J/ψ → p∆ 0 π − ,∆ 0 →pπ + ) are not peaking in the χ cJ mass regions. The amount of background events is estimated by using the normalized K 0 S and Σ + mass sideband events, as shown in Fig. 1. The blue long dashed boxes are the selected K 0 S mass sidebands (1.1694 < M pπ 0 < 1.2094 GeV/c 2 , 0.466 < M π + π − < 0.482 GeV/c 2 and 0.514 < M π + π − < 0.530 GeV/c 2 ) and the Σ + mass sidebands (0.49 < M K 0 S < 0.506 GeV/c 2 , 1.1094 < M pπ 0 < 1.1494 GeV/c 2 and 1.2294 < M pπ 0 < 1.2694 GeV/c 2 ), and the green dashed boxes are those from non-K 0 S and non-Σ + sidebands (1.1094 < M pπ 0 < 1.1494 GeV/c 2 and 1.2294 < M pπ 0 < 1.2694 GeV/c 2 , 0.466 < M π + π − < 0.482 GeV/c 2 and 0.514 < M π + π − < 0.530 GeV/c 2 ). The normalized background contribution in the χ cJ mass regions is estimated as half of the total number of events in the four blue sideband regions minus one quarter of the total number of events in the four green sideband regions of Fig. 1, and shown as a green-shaded histogram in Fig. 4.  Figure 4. Fit to the Σ +p K 0 S invariant-mass distribution in the χcJ mass region of [3.3, 3.6] GeV/c 2 . Dots with error bars are data, the red solid curve shows the result of the unbinned maximum-likelihood fit, the green-shaded histograms are the events from the normalized K 0 S and Σ + mass sidebands, the blue solid line is sum of the peaking and flat background components, and the violet long dashed curve is the contribution of the peaking background normalized according to the sideband events.
An unbinned maximum-likelihood fit to the Σ +p K 0 S invariant-mass distribution is performed for the total selected signal candidates, as shown in Fig. 4. The complete formula for the fit is P DF total = N 1 × P DF signal + N 2 × P DF peakingbkg + N 3 × P DF flatbkg . The parameters N 1 and N 3 are free, and N 2 is fixed to the number of the events determined from the K 0 S and Σ + mass sidebands. Here, P DF signal is the sum of the signal line shapes of the three χ cJ resonances each convolved with a Gaussian function related to the χ cJ mass resolution, where the width of the Guassian function is fixed to each of the MC simulated value. The line shape of each resonance is described by: with m χcJ and Γ χcJ the mass and width of the cor- is the energy of the transition photon in the rest frame of ψ(3686) and D(E γ ) is the damping factor which suppresses the divergent tail due to the E 3 γ dependence of P DF signal . It is described by exp(−E 2 γ /8β 2 ) where β is one of the free parameters in the fit. For all three resonances the same β value is required. The fit result β = (68.7 ± 13.0) MeV is consistent with the value measured by the CLEO experiment [19].
The peaking background component P DF peakingbkg is the same as the signal distribution. It is used to describe the distribution of the normalized events from the K 0 S and Σ + mass sidebands where clearly the three χ cJ resonances can be identified. The P DF flatbkg is described by a first-order polynomial.
For the unbinned maximum-likelihood fit, β, the masses and widths of the χ cJ resonances, and the two coefficients of the polynomial are taken as free parameters. The event yields of the fitted χ cJ → Σ +p K 0 S signals are listed in Table I. The branching fractions for χ cJ → Σ +p K 0 S are calculated by where N ψ(3686) is the total number of ψ(3686) events, ǫ is the corresponding detection efficiency as listed in Table I, which is obtained by weighting the simulated Dalitz plot distribution with the distribution from data, and j B j = B(ψ(3686) → γχ cJ ) × B(Σ + → pπ 0 ) × B(K 0 S → π + π − ) × B(π 0 → γγ), where the branching fractions are taken from the PDG [3]. The results of the branching-fraction calculation for the decays χ cJ → Σ +p K 0 S are also listed Table I. Number of signal events (N χ cJ obs ), detection efficiency (ǫ), and branching fractions B(χcJ → Σ +p K 0 S ), where the first uncertainty is statistical and the second is systematic.

V. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties on the χ cJ → Σ +p K 0 S branching-fraction measurements are listed in Table II.
The systematic uncertainty of the photondetection efficiency is studied by considering the decay J/ψ → π + π − π 0 [20] and is about 1% for each photon, so 3% is assigned for the three photons in the final states.
The uncertainty related to the particle identification (PID) and tracking of the proton and antiproton is studied with the control samples of J/ψ and ψ(3686) → ppπ + π − [21]. The average differences of efficiencies between MC simulations and data are 0.4%, 0.4%, and 0.3% for the proton from χ c0 , χ c1 , and χ c2 decays, respectively, with the transverse momentum and angle region of our signal channel considered. Similarly forp, they are 0.4%, 0.3%, and 0.3%, respectively, so the uncertainties on the proton and anti-proton pair PID and tracking are 0.6%, 0.5%, and 0.4% for χ c0 , χ c1 , and χ c2 decays, respectively.
The uncertainty associated with the 4C kinematic fit comes from the inconsistency between data and MC simulation, as described in detail in Ref. [22]. In this analysis, we take the efficiency with the correction as the nominal value, and the differences between the efficiencies with and without correction, 0.4%, 0.4%, and 0.3% for χ c0 , χ c1 , and χ c2 , respectively, as the systematic uncertainties from the kinematic fit.
The uncertainty related with the π 0 (K 0 S , Σ + ) mass window requirement is studied by fitting the π 0 (K 0 S , Σ + ) mass distributions of data and signal MC simulation with a free Crystal Ball (Gaussian, Gaussian) function and a first-order Chebyshev polynomial function. We obtained the selection efficiency of the π 0 (K 0 S , Σ + ) mass region, which is the ratio of the numbers of π 0 (K 0 S , Σ + ) events with and without the π 0 (K 0 S , Σ + ) mass window, determined by integrating the fitted signal shape. The difference in efficiency between data and MC simulation, 0.3% (0.3%, 0.1%), is assigned as the systematic uncertainty. The systematic uncertainty from the veto of the Λ mass window is negligible due to the high detection efficiency.
The uncertainty of the detection efficiency is studied by changing the number of bins in the Dalitz plot. The maximum differences of the signal detection efficiency, 1.0%, 0.5% and 0.4% , are taken as uncertainties for χ c0 , χ c1 , and χ c2 decays, respectively. The uncertainty of assuming ψ(3686) → γχ c1 (χ c2 ) as pure E1 transition is studied by considering the contribution from higher order multiple amplitudes [24] in the MC simulation, the differences of the efficiency, 0.8% for χ c1 and 0.2% for χ c2 , are taken as the systematic uncertainties. For χ c0 → Σ +p K 0 S , there is a possible structure in thepK 0 S invariant distribution. The corresponding systematic uncertainty is estimated by mixing χ c0 → Σ +Σ (1940) − MC sample and the PHSP signal MC sample in a proportion, which is obtained from fitting Mp K 0 S distribution. The difference between the efficiencies before and after mixing, 0.1%, is considered to be the systematic uncertainty. The total uncertainties associated with the efficiency for χ c0 , χ c1 , and χ c2 are 1.0%, 0.9%, and 0.4%, respectively.
The systematic uncertainty due to the signal line shape is considered by changing the damping factor from exp(−E 2 γ /8β 2 ) to E0Eγ +(E0−Eγ ) 2 used by KEDR [25], is the peak energy of the transition photon, the differences in the fit results for χ c0 , χ c1 , and χ c2 , 1.4%, 1.9%, and 0.4% are assigned as the systematic uncertainties.
The uncertainty associated with the detector resolution is studied by making the width of the Gaussion function to be free, no changes are found for the χ c0 , χ c1 , and χ c2 signal yields, thus these uncertainties are neglected.
The systematic uncertainties due to the χ c0 , χ c1 , and χ c2 mass and width in the fit are studied by changing them from free to the world average values [3]. The differences of the χ c0 , χ c1 , and χ c2 signal yields, 3.0%, 0.4% and 3.9% are taken as the systematic uncertainties.
The uncertainty from the determination of χ cJ signal events due to the fit range is obtained from the maximum difference in the fit results by changing the fit range from [3.30, 3.60] GeV/c 2 to [3.30, 3.65] GeV/c 2 or [3.25, 3.60] GeV/c 2 . The maximum differences in the fitted yields for χ c0 , χ c1 , and χ c2 are 0.9%, 1.4%, and 0.8%, respectively.
The uncertainty due to the estimation of the background contribution using the K 0 S and Σ + mass sidebands can be estimated by changing the sideband ranges. Changing the mass range of K S , non-Σ + mass region accordingly, the differences of χ c0 , χ c1 , and χ c2 signal yields are 0.3%, 0.1%, and 0.5%, respectively. The uncertainty from the shape of the non-χ cJ background is estimated by changing the polynomial degree from the first to the second in fitting the Σ +p K 0 S invariant mass, and the differences in the fit results are 2.8%, 1.4%, and 1.4%, respectively. The total uncertainties associated with the background shape are 2.8%, 1.4%, and 1.5% for χ c0 , χ c1 , and χ c2 decays, respectively.
The total systematic uncertainty is the sum in quadrature of all uncertainties added for each χ cJ decay.