Observation of the decays $\chi_{cJ}\to n K^{0}_{S}\bar\Lambda + c.c.$

By analyzing $4.48\times10^8$ $\psi(3686)$ events collected with the BESIII detector, we observe the decays $\chi_{cJ} \to n K^0_S\bar\Lambda + c.c.$ ($J=0$, 1, 2) for the first time, via the radiative transition $\psi(3686) \to \gamma \chi_{cJ}$. The branching fractions are determined to be $(6.67 \pm 0.26_{\rm stat} \pm 0.41_{\rm syst})\times10^{-4}$, $(1.71 \pm 0.12_{\rm stat} \pm 0.12_{\rm syst})\times10^{-4}$, and $(3.66 \pm 0.17_{\rm stat} \pm 0.23_{\rm syst})\times10^{-4}$ for $J=0$, 1, and 2, respectively.


Introduction
Studying the hadronic decays of the cc states J/ψ, ψ(3686), and χ cJ (J = 0, 1, 2) provides valuable information on perturbative QCD in the charmonium-mass regime and on the structure of charmonia. The color-octet mechanism, which successfully describes several decay patterns of P-wave χ cJ states [1], may be applicable to further χ cJ decays. Measurements of χ cJ hadronic decays can provide new input on the color-octet mechanism and further assist in understanding χ cJ decay mechanisms.
The BES Collaboration observed near-threshold structures in baryon-antibaryon invariant mass distributions in the radiative decay J/ψ → γpp [2] and the purely hadronic decay J/ψ → pK −Λ [3]. It was suggested theoretically that these near-threshold structures might be observation of baryonium [4][5][6] or caused by final state interactions [7][8][9]. Excited Λ and N resonances were observed in the study of the decay J/ψ → nK 0 SΛ [10]. Studying the same decay modes in other charmonia can provide complementary information on these structures.

Detector and data sets
The BESIII detector is a magnetic spectrometer [14,15] located at the Beijing Electron Positron Collider (BEPCII) [16]. 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 acceptance of charged particles and photons is 93% over 4π solid angle. The charged-particle momentum resolution at 1.0 GeV/c is 0.5%, and the specific energy loss (dE/dx) resolution is 6% for the 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 of the TOF barrel part is 68 ps, while that of the end cap part is 110 ps.
Simulated samples produced with the geant4-based [17] Monte-Carlo (MC) software, which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiencies and to estimate the background levels. The simulation takes into account the beam energy spread and initial state radiation (ISR) in the e + e − annihilation modeled with the generator kkmc [18]. The inclusive MC samples consist of 5.06×10 8 ψ(3686) events, the ISR production of the J/ψ state, and the continuum processes incorporated in kkmc. The known decay modes are modeled with eventgen [19] using the BFs taken from the Particle Data Group (PDG) [13], and the remaining unknown decays from the charmonium states with lundcharm [20]. Radiation from charged final state particles is incorporated with the photos [21] package.
The signal detection efficiencies are estimated with signal MC samples. The decays of ψ(3686) → γχ cJ (J = 0, 1, 2) are simulated following Ref. [22], in which the magneticquadrupole (M2) transition for ψ(3686) → γχ c1,2 and the electricoctupole (E3) transition for ψ(3686) → γχ c2 are considered in addition to the dominant electric-dipole (E1) transition. For the decay χ cJ → nK 0 SΛ + c.c., a special generator based on results of Helicity Partial Wave Analysis (HelPWA) [23-26] is used. The background MC samples are obtained from inclusive MC samples, in which the signal channels are removed with TopoAna [27], and the samples are normalized to the data sample based on luminosity.

Event selection and background analysis
The signal process ψ(3686) → γχ cJ , χ cJ → nK 0 SΛ with K 0 S → π + π − andΛ →pπ + consists of the final state particles γnpπ + π + π − . Charged track candidates from the MDC must satisfy |cos θ| < 0.93, where θ is the polar angle with respect to the z axis, which is the axis of the MDC. The closest approach to the interaction point is required to be less than 20 cm along the z direction and less than 10 cm in the plane perpendicular to z. The TOF and dE/dx information are combined to calculate the particle identification probabilities (P ) for the hypotheses that a track is a pion, kaon, or proton. Proton candidates are required to satisfy P (p) > P (K) and P (p) > P (π). Exactly four charged tracks are required in each candidate event.
Since K 0 S andΛ have relatively long lifetimes, they are reconstructed by constraining the π + π − pair and thepπ + pair to originate from secondary vertices, respectively. Charged track candidates, except the one used as ap in theΛ reconstruction, are assumed to be pions without applying particle identification. The decay lengths from the secondary vertex fits of K 0 S andΛ divided by their corresponding uncertainties are required to be larger than two. The mass distributions of the reconstructed K 0 S (denoted as M ππ ) andΛ (denoted as Mp π + ) candidates are shown in Figs. 1(a) and 1(b), respectively, where the signal regions are defined as [0.480, 0.516] GeV/c 2 for K 0 S and [1.112, 1.120] GeV/c 2 forΛ. Photons are reconstructed as energy clusters in the EMC. The shower time is required to be within [0, 700] ns from the event start time. Photon candidates within |cos θ| < 0.80 (barrel) are required to have deposited energies larger than 25 MeV and those with 0.86 < |cos θ| < 0.92 (end cap) must have deposited energies larger than 50 MeV. The photon candidates must be at least 10°away from any charged track to suppress Bremsstrahlung photons or splitoffs. We require there is at least one photon candidate satisfying the above criteria. To select the best combination and to improve the resolution, a 1-C kinematic fit is applied under the hypothesis of ψ(3686) → γnK 0 SΛ , where the neutron is treated as a missing particle. The value of χ 2 1C is required to be less than 200. If more than one combination survives in an event, the one with the smallest χ 2 1C is retained. M miss is defined as the invariant mass of the four momentum of the missing particle where p i is the four momentum of the particle i and p CM is the four momentum of the initial e + e − system. To further improve the purity, M miss before the 1-C kinematic fit is required to satisfy 0.90 < M miss < 0.98 GeV/c 2 . The distribution of M miss before the 1-C kinematic fit is shown in Fig. 1(c).
By analyzing the ψ(3686) inclusive MC samples with TopoAna [27], the only significant peaking background is found to be caused by the decays χ cJ → Σ ±Λ π ∓ + c.c. with Σ ± → nπ ± . The two pions in this decay may accidentally fall into the K 0 S mass window and fulfill the constraint for the secondary vertex fit, thus faking the K 0 S . This Σ ± background peaks in the invariant mass spectrum of nK 0 SΛ (denoted as M nK 0 SΛ ) but distributes uniformly in the M ππ spectrum. Other backgrounds are smoothly distributed underneath the χ cJ signal.

Branching fraction measurement
A simultaneous unbinned maximum-likelihood fit to the M nK 0 SΛ spectra in both the M ππ signal and sideband regions, as shown in Fig. 2, is performed to determine the signal yields and peaking backgrounds. The lower and upper sideband regions are defined as and that for the sideband region is The signal shape f J signal for each χ cJ resonance is described by its line shape convolved with a double-Gaussian function to account for the mass resolution. Each signal line shape is modeled with is the nonrelativistic Breit-Wigner function with the width Γ χ cJ and the mass m χ cJ of the corresponding χ cJ fixed to the PDG values [13], is the energy of the transition photon in the rest frame of ψ(3686), and D(E γ ) is a damping factor that suppresses the divergent tail due to E 3 γ . This damping factor is described by D(E γ ) = exp(−E 2 γ /8β 2 ), where β = (65.0 ± 2.5) MeV is measured by the CLEO collaboration [28]. The two Gaussian functions in the convolution share the same mean value which is then floated in the fit. The relative width and size of the second Gaussian to the first Gaussian function are fixed to the results of MC studies, while the width of the first Gaussian function is floated. The peaking background shapes f J peakbkg are parameterized the same as the signal shapes. The yields N 2,J in the signal region are normalized to N 2,J in the M ππ sideband regions according to the sizes of these two regions.
The background shapes f ( ) flatbkg in both regions are modeled as second-order Chebyshev polynomial functions. The numbers of fitted χ cJ signal events, N 1,J , are listed in Table 1.
A special generator based on results of HelPWA [24][25][26] for the decay χ cJ → nK 0Λ +c.c. is developed to estimate the detection efficiencies. The description of HelPWA can be found in the supplemental material [23]. For the χ cJ data events used in HelPWA, the signal regions of M nK 0 SΛ for χ c0 , χ c1 , and χ c2 are [3.39, 3.45], [3.50, 3.53], and [3.54, 3.57] GeV/c 2 , respectively; the masses of K 0 ,Λ, and χ cJ are constrained to their known masses [13]. The inclusive background MC sample is used to calculate the background likelihood with negative weight.
The signal MC events are generated with the HelPWA model, in which the parameters of coupling constants are determined by fitting the model to the χ cJ data events. The Dalitz plots and the two-body invariant mass M ij distributions of the data sample are shown in Figs. 3 and 4, respectively, where i and j denote the final particles. The generated signal MC events based on HelPWA along with the simulated background events from the inclusive MC sample are represented as the solid lines in Fig. 4. The signal MC samples are generated with χ cJ → nK 0 SΛ and χ cJ →nK 0 S Λ separately. After applying the same event selection criteria to the signal MC samples, we fit the invariant mass spectra with the same methods used for the experimental data. The detection efficiencies of χ cJ , J , are averaged over both charge conjugate channels and listed in Table 1. The efficiency differences between the two charged conjugated modes are less than 0.5% for all three χ cJ channels and are consistent within the statistical uncertainties.    The BFs for χ cJ → nK 0 SΛ are calculated using where N ψ(3686) is the number of ψ(3686) events [12]; J is the detection efficiency as listed in Table 1; where the BFs are taken from the PDG [13]. The results are summarized in Table 1.

Systematic uncertainty
The number of ψ(3686) events is measured to be (4.48 ± 0.03) × 10 8 based on inclusive hadronic events, as described in Ref. [12], so the uncertainty is 0.6%. The systematic uncertainty due to the detection of γ is studied with the well understood channel J/ψ → ρ 0 π 0 [29]. The efficiency difference between data and MC simulation is about 1% per photon. To estimate the uncertainties associated with K 0 S and Λ reconstruction, the decays J/ψ → K * ± (892)K ∓ , K * ± (892) → K 0 S π ± and J/ψ → ΛΛ are selected as the control samples. The uncertainties are determined to be 1.5% per K 0 S and 2.0% per Λ. The systematic uncertainty caused by the 1-C kinematic fit is studied with the control sample ψ(3686) → π 0 nK 0 SΛ with purity about 98.5%, where n, K 0 S , andΛ are selected using the same criteria as the nominal analysis but the number of good photons is required to be at least two. The difference of the selection efficiencies between data and MC simulation is determined to be 4.1% and assigned as the corresponding systematic uncertainty. The uncertainties associated with the mass windows of K 0 S , Λ, and M miss are estimated by repeating the analysis with alternative mass window requirements. We change the mass window of K The systematic uncertainty related to the fitting procedure includes multiple sources. For the signal line shape, the parameterization of the damping factor may introduce a systematic uncertainty. The nominal damping factor is changed to another popular form used by KEDR [30], given by D(E γ ) = . The resulting differences in the fit are assigned as the related systematic uncertainties. In addition, the background function is changed from a second to a third order Chebyshev function, and the differences in signal yields are taken as the systematic uncertainties. The systematic uncertainty due to the fit range is evaluated by changing the fit range from [3.35, 3.60] to [3.35, 3.65] and [3.30, 3.65] GeV/c 2 , and the maximum differences in the fitted yields are considered as the associated systematic uncertainties. The total uncertainties of the fitting procedure are estimated to be 2.8%, 4.1%, and 2.3% for χ c0 , χ c1 , and χ c2 , respectively.
The systematic uncertainties arising from MC modeling are estimated by using different model parameters and some unimportant intermediate processes in HelPWA. The differences of efficiencies based on the new HelPWA results and the nominal ones are taken as the uncertainties. The systematic uncertainties due to the input BFs of ψ(3686) → γχ c0 (χ c1 , χ c2 ), K 0 S → π + π − , and Λ → pπ − are 2.0% (2.5%, 2.1%), 0.1%, and 0.8%, respectively, according to the PDG [13].
All systematic uncertainty contributions are summarized in Table 2. The total systematic uncertainty for each χ cJ decay is obtained by adding all contributions in quadrature.
Enhancements are observed in the Dalitz plots shown in Fig. 3 and the mass distributions of two-body nΛ subsystems shown in Fig. 4. We perform a HelPWA with the main goal to produce MC samples that describe the data well enough to obtain a good estimate of the efficiency. While the HelPWA does describe the data nicely, the complexity of the model we used here does not allow to draw any firm conclusions on the relative contributions of individual resonances.    Table 1. Helicity angles and amplitudes of the sequential decays (a), (b) and (c). λ i denotes the value of helicity for the corresponding particle, and m denotes the spin z projection of virtual photon χ cJ .
Helicity amplitude is expressed with partial waves in terms of LS-coupling scheme [4][5][6]. For two-body decays with spin J → s + σ, it follows where L and S are the orbital angular momentum and spin quantum numbers for the two-body decay, r is the magnitude of the relative momentum of the final states, r 0 is the corresponding momentum at the R nominal mass, g ls is a coupling constant to be determined, and B L (r) is the Blatt-Weisskopf barrier factor [7]. The non-resonant decay is described with the helicity amplitude constructed in the event system, in which the X-Y plane is chosen as the three-body decay plane, and the Z-axis is defined as the normal to the decay plane. To describe the n, K 0 S and Λ orientations in the X-Y -Z system, one needs to introduce three rotations to carry the momenta of final states from the χ cJ rest frame, x-y-z to the three-body helicity system, X-Y -Z by three Euler angles, α, β,and γ as shown in Fig. 3. The calculations of Euler angles can be found in [8].
The amplitude is written as [9], where m denotes the z-component of the χ cJ spin in its rest frame. The decay amplitude T λ 1 ,λ 2 ,λ 3 is a fit constant depending on the helicity values λ i of the final state particles, and µ is the helicity of the χ cJ in the helicity system fixed to the three-body reference, which is summed over.

Fit method
To determine the decay coupling constants g ls and decay amplitude T , we fit the model to the data events by minimizing the function defined as S = −(ln L data − ln L bg ), (1.10) X Y Z γ x y z α β N Figure 3. Illustration of rotations to carry the n, K 0 S and Λ orientations from the χ cJ rest frame xyz to the three body helicity system XY Z by the three Euler angles α, β and γ.
where the contributions from background events are subtracted from the likelihood of observed data events L data , and the likelihood function L for the N data events is defined with the amplitudes as  [10] in the fit. The possible interferences between the non-resonant process and the processes with intermediate states are also included, and parity conservation is required to reduce the number of independent helicity amplitudes.