Observation of the decays $\chi_{cJ}\to\Sigma^{0}\bar{p}K^{+}+{\rm c.c.}~(J=0, 1, 2)$

The decays $\chi_{cJ}\to\Sigma^{0}\bar{p}K^{+}+{\rm c.c.}~(J = 0, 1, 2)$ are studied via the radiative transition $\psi(3686)\to\gamma\chi_{cJ}$ based on a data sample of $(448.1 \pm 2.9)\times10^{6}$ $\psi(3686)$ events collected with the BESIII detector. The branching fractions of $\chi_{cJ}\to\Sigma^{0}\bar{p}K^{+}+{\rm c.c.}~(J = 0, 1, 2)$ are measured to be $(3.03 \pm 0.12\pm 0.15)\times10^{-4}$, $(1.46 \pm 0.07\pm 0.07)\times10^{-4}$, and $(0.91 \pm 0.06\pm 0.05)\times10^{-4}$, respectively, where the first uncertainties are statistical and the second are systematic. In addition, no evident structure is found for excited baryon resonances on the two-body subsystems with the limited statistics.


I. INTRODUCTION
The P-wave charmonia χ cJ (J = 0, 1, 2) have been observed experimentally for a long time, however, most decay modes of them are still unknown. Though χ cJ can not be directly produced via electron-positron annihilation into a virtual photon, radiative decays of the ψ(3686) into χ cJ states make up about 10% of the total decay width of the ψ(3686) for each χ cJ [1]. Thus, the large ψ(3686) data sample containing (448.1 ± 2.9) × 10 6 events at BESIII can ideally be used to investigate χ cJ decays [2,3].
Many two-body decays of χ cJ → BB have been observed in experiments, but three-body decays of χ cJ → BBM are much less measured (B stands for a baryon, M stands for a meson), while the latter have advantages to search for and study excited baryons due to larger freedom of quantum numbers. For example, some experiments reported two excited Σ resonances around 1670 MeV/c 2 , which have the same mass and J P C quantum numbers but very different decay products and angular distributions [4][5][6][7]. Further experimental information will shed light on the understanding of these states.

II. BESIII DETECTOR
The BESIII detector [9] records symmetric e + e − collisions provided by the BEPCII storage ring [10], 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.7 GeV. 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 charged-particle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for the electrons of 1 GeV/c momentum. 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 [11] collected with the BESIII detector.
Simulated data samples produced with a geant4based [12] Monte Carlo (MC) package, which includes the geometric description of the BESIII detector and the detector response, are used to determine detection efficiencies and to estimate backgrounds. The simulation models the beam energy spread and initial state radiation (ISR) in the e + e − annihilations with the generator kkmc [14,15]. The inclusive MC sample includes 506 × 10 6 ψ(3686) events, the ISR production of the J/ψ, and the continuum processes incorporated in kkmc. The known decay modes are modelled with evtgen [16,17] using branching fractions taken from the Particle Data Group [1], and the remaining unknown charmonium decays are modelled with lundcharm [18]. Final state radiation (FSR) from charged particles is incorporated using the photos package [19].

IV. EVENT SELECTION AND BACKGROUND ANALYSIS
For ψ(3686) → γχ cJ , χ cJ → Σ 0p K + with Σ 0 → γΛ and Λ → pπ − , the final state consists of ppK + π − γγ. 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 symmetry axis of the detector. For the two charged tracks from the Λ decay, the distance between their point of closest approach and the primary vertex is required to be less than 20 cm along the beam direction, and less than 10 cm in the plane perpendicular to the beam direction. For the remaining charged tracks, the same distance is required to be less than 10 cm along the beam direction and less than 1 cm in the plane perpendicular to the beam direction. The total number of charged tracks needs to be equal to or greater than four.
The TOF and dE/dx information is used to calculate a particle identification (PID) likelihood (P ) for the hy-potheses that a track is a pion, kaon, or proton. Tracks from the primary vertex are required to be identified as either an anti-proton (P (p) > P (K) and P (p) > P (π)) or a kaon (P (K) > P (p) and P (K) > P (π)). In case of daughter particles of a Λ candidate, the track with the larger momentum is identified as the proton, and the other is identified as the pion. For each candidate event, exactly onep, K + , and p, π − from the Λ decay are required. For all combinations of positively and negatively charged tracks, secondary vertex fits are performed [21], and the combination with the smallest χ 2 Λ is retained as the Λ candidate. In addition, the ratio of the decay length (L) to its resolution (σ L ) is required to be larger than zero. The mass distribution of the reconstructed Λ candidates is shown in Fig. 1(a). A mass window of |M pπ − − m Λ | < 0.004 GeV/c 2 is required to select the Λ signal events, where M pπ − is the invariant mass of selected proton-pion pairs and m Λ is the nominal mass of Λ taken from the PDG [1].
Photon candidates are reconstructed from the energy deposition in the EMC crystals produced by electromagnetic showers. The minimum energy requirement for a photon candidate is 25 MeV in the barrel region (| cos θ| < 0.80) and 50 MeV in the end-cap region (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 tracks. The time-information 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. The total number of photons is required to be at least two. To reduce background events from π 0 → γγ, we require |M γγ − m π 0 | > 0.015 GeV/c 2 .
A four-constraint (4C) kinematic fit imposing four-momentum conservation is performed using the ppK + π − γγ hypothesis. If there are more than two photon candidates in one event, the combination with the smallest χ 2 4C is retained, and its χ 2 4C is required to be smaller than those for the alternative ppK + π − γ and ppK + π − γγγ hypotheses. In addition, the value of χ 2 4C is required to be less than 40. For the selected signal candidates, the γΛ combination with the invariant mass closest to the nominal Σ 0 mass according to the PDG [1] is taken as the Σ 0 candidate. The distribution of the γΛ invariant mass is shown in Fig. 1 Fig. 1(b).
The Σ 0p K + invariant mass distribution after application of all selection conditions is shown in Fig. 2, where clear χ c0 , χ c1 , and χ c2 signals are observed. The signal MC simulation also shown in Fig. 2 agrees with the data very well. The ψ(3686) inclusive MC sample is used to study possible peaking backgrounds. Applying the same re-quirements as the data, the two main remaining background channels involve either ψ(3686) → K * +p Λ with K * + → K + π 0 (π 0 → γγ) decays or belong to the peaking background channel ψ(3686) → γχ cJ → γK +p Λ (Λ → pπ − ) that is missing the intermediate Σ 0 decay. Other small backgrounds are smoothly distributed below the χ cJ signal region. These backgrounds can be estimated by the Σ 0 sideband events normalized to the background level below the Σ 0 signal peak. The normalized sideband events are shown as the histogram with the dot line in Fig. 2.
The result of an unbinned maximum likelihood fit to the M Σ 0p K + distribution is shown in Fig. 3. Here, we signal is the probability density function describing the χ cJ resonances, f J peakbkg is the normalized shape of the Σ 0 sidebands, and f flatbkg is given by a secondorder polynomial. The line shape of each resonance f signal is modeled with the same formula BW 2 is the Breit-Wigner function, Γ χcJ is the width of the corresponding is the energy of the transition photon in the rest frame of the ψ(3686), and D(E γ ) is the damping factor which suppresses the divergent tail due to the E 3 γ dependence of f J signal . It is described by exp(−E 2 γ /8β 2 ), where β = (65.0 ± 2.5) MeV was measured by the CLEO experiment [22]. The signal shapes are convolved with Gaussian functions to account for the mass resolution. The parameters N 1,J , N 3 and two coefficients of the polynomial are taken as the free parameters in the fit, while N 2,J is fixed to the number of the normalized Σ 0 sideband events. In the description of f J signal , the masses and widths of the χ cJ states are fixed to the PDG values. The Gaussian resolution parameters in the region of the three χ cJ states are also free parameters, and are found to be 5.7, 5.1, and 4.1 MeV/c 2 for χ c0 , χ c1 , and χ c2 , respectively. The yields of signal events of all three χ cJ → Σ 0p K + decays are listed in Table I.
Dalitz plots and the one dimensional projections of χ cJ → Σ 0p K + events are shown in the left, middle and right columns of Fig. 4 for the χ c0 , χ c1 , and χ c2 signal regions, respectively, together with the distributions of MC simulated signal events based on a pure phase-space decay model.
ForpK + mass spectra of the data, it seems there are two structures around 1.7 and 2.0 GeV/c 2 for χ c0 decays, they are likelyΣ(1750) 0 andΣ(1940) 0 . There seems to be two structures around 1.9 GeV/c 2 for χ c1 decays and around 1.8 GeV/c 2 for χ c2 decays. For Σ 0 K + mass spectra, it seems there is a jump around 1.8 GeV/c 2 and a dip around 2.0 GeV/c 2 for χ c0 decays, the jump may be N (1880) with J P = 1 There is an indication around 1.95 GeV/c 2 for χ c1 decays, which may be N (1900) with J P = 3 2 + . There is no evident structure for χ c2 decays. For Σ 0p mass spectra, the data are consistent with the phase-space MC shapes, there is no evident structure for χ c0 , χ c1 and χ c2 decays. The mass distributions of two-body subsystems of the data are not completely consistent with the phase-space MC simulations, but it is difficult to draw any conclusions to them due to present limited statistics.
The differences between data and MC simulation indicate that these signal MC events cannot be used to calculate the selection efficiency directly. Instead, the detection efficiency is obtained by weighting the simulated Dalitz plot distribution with the distribution from data. We divide the Dalitz plots of M 2 pK + versus M 2 Σ 0 K + into 12 × 12, 8 × 7, and 6 × 8 bins in the χ c0 , χ c1 , and χ c2 regions, respectively. First, we obtain the weight factor ω i in each bin as the ratio between the Dalitz plot distribution of data and the normalized signal MC sample. In a second step, ω i is used to correct the Dalitz distributions of both the generated and reconstructed MC simulations. Finally, we determine the corrected detection efficiency as the ratio between the sum of event weights in reconstructed and generated MC. The results are listed in Table I. The branching fractions for χ cJ → Σ 0p K + + c.c. are calculated using where N obs is the number of signal events obtained from the fit, N ψ(3686) is the total number of ψ(3686) events, ǫ is the corresponding detection efficiency as listed in Table I, and j B j = B(ψ(3686) → γχ cJ ) × B(Σ 0 → γΛ)×B(Λ → pπ − ) is the product branching fraction with individual values taken from the PDG [1]. The results for χ cJ → Σ 0p K + + c.c. (J = 0, 1, 2) are listed in Table I.

VI. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties on the measurement of the branching fractions of χ cJ → Σ 0p K + + c.c. are discussed below.
Using the control samples of J/ψ → ppπ + π − and J/ψ → K * K , the difference of tracking efficiencies between MC simulation and data is within 1% forp and K + . Therefore, 2% is taken as the tracking systematic uncertainty.
Thep/K + PID efficiency is studied using J/ψ → ppπ + π − and J/ψ → K 0 S K ± π ± control samples [23,24], with the result being that the PID efficiency for data agrees with that of the MC simulation within 1% per p/K + . So 2% is taken as the systematic uncertainty associated with the PID efficiency.
The photon detection efficiency is studied from a J/ψ → π + π − π 0 control sample [25]. The efficiency difference between data and MC simulation is about 1% per photon, so that 2% is assigned as the systematic uncertainty from the two photons.
In order to determine the uncertainty associated with the secondary vertex fit and the decay length requirement, we determine the efficiency of these selection criteria by comparing the Λ → pπ − signal yields with and without those selections for both data and signal MC. From a fit to the pπ − invariant mass distributions, we find a data-MC difference of 0.7% that is assigned as the systematic uncertainty. For each track stemming from Λ → pπ − decays, the systematic uncertainty from the tracking efficiency is 1.0% according to an analysis of J/ψ →pK + Λ [26]. The total uncertainty of the Λ reconstruction is 2.1%.
The uncertainty associated with the 4C kinematic fit comes from a potential inconsistency between data and MC simulation; this difference is reduced by correcting the track helix parameters in the MC simulation, as described in detail in Ref. [27]. The difference of the efficiency with and without the helix correction is considered as the systematic uncertainty from the kinematic fit.
The uncertainty related to the Λ and Σ 0 mass windows is studied by determining the yield of Λ (Σ 0 ) inside the mass windows for both data and signal MC simulation. The difference between data and MC simulation is found to be negligible for Λ, and to be 0.2% for Σ 0 .
In the weighting procedure, the Dalitz plots were divided into 12 × 12, 8 × 7 and 6 × 8 bins in order to calculate the event-weights used in the efficiency determination. We repeat this procedure with different bin configurations. The maximum difference between the nominal binning and the alternate configuration is taken as the weighting related uncertainty listed in Table II. The sta- Table I. Summary of the number of fitted signal events (N obs ), detection efficiency (ǫ), and branching fraction B(χcJ → Σ 0p K + + c.c.), where the first uncertainty is statistical and the second one is systematic.     tistical uncertainty of the efficiency is determined directly from MC simulations and amounts to less than 0.5%.
The systematic uncertainty related to the fitting procedure includes multiple sources. Concerning the signal line shape, the damping factor is changed from exp(−E 2 γ /8β 2 ) as used by CLEO to E0Eγ +(E0−Eγ ) 2 as used by KEDR [28]. The resulting differences in the fit are assigned as the systematic uncertainties. In addition, the fit range is varied from [3.30, 3.60] GeV/c 2 to [3.30, 3.65] GeV/c 2 and [3.25, 3.60] GeV/c 2 and the maximum differences in the fitted yields are considered as the associated systematic uncertainties. Regarding the peaking background contributions, the Σ 0 sideband ranges were changed from [ With regard to non-χ cJ backgrounds, the fit function is changed from a second to a third order polynomial in the fit to the Σ 0p K + invariant mass distribution and the difference between the two fits is taken as the systematic uncertainty.
Although there is no evident intermediate resonances on two-body subsystems of χ cJ decays, the mass distributions of two-body subsystems are not completely consistent with the phase-space MC simulations. This implies the existence of intermediate baryon resonances. With the present statistics, it is difficult to study them in detail and draw any conclusions to them. More ψ(3686) events in the future in combination with advanced analysis technique, such as partial wave analysis, may shed light on the intermediate structures.