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

Studying the hadronic decays of the $c\bar{c}$ states $J/\psi$, $\psi(3686)$, and $\chi_{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 $\chi_{cJ}$ states ref::Wong , may be applicable to further $\chi_{cJ}$ decays. Measurements of $\chi_{cJ}$ hadronic decays can provide new input on the color-octet mechanism and further assist in understanding $\chi_{cJ}$ decay mechanisms.

The BES Collaboration observed near-threshold structures in baryon-antibaryon invariant mass distributions in the radiative decay $J/\psi\to\gamma p\bar{p}$ ref::prl-91-022001 and the purely hadronic decay $J/\psi\to pK^{-}\bar{\Lambda}$ ref::prl-93-112002 . It was suggested theoretically that these near-threshold structures might be observation of baryonium ref::Datta1 ; ref::Datta2 ; ref::Datta3 or caused by final state interactions ref::Kerbikov1 ; ref::Kerbikov2 ; ref::Kerbikov3 . Excited $\Lambda$ and $N$ resonances were observed in the study of the decay $J/\psi\to nK_{S}^{0}\bar{\Lambda}$ ref::plb-659-789 . Studying the same decay modes in other charmonia can provide complementary information on these structures.

Anomalous enhancements near the threshold of $p\bar{\Lambda}+c.c.$ have been also observed in the decays of $\chi_{cJ}\to pK^{-}\bar{\Lambda}+c.c.$ by the BESIII Collaboration ref::paper-pkl . If isospin symmetry is conserved, the decay branching fraction (BF) ratio $\mathcal{B}(\chi_{cJ}\to pK^{-}\bar{\Lambda}+c.c.)/\mathcal{B}(\chi_{cJ}\to nK^{0}_{S}\bar{\Lambda}+c.c.)\sim 2$ should be satisfied, thereby implying the existence of the isospin conjugate decays $\chi_{cJ}\to nK^{0}_{S}\bar{\Lambda}+c.c.$.

In this analysis, we present the study of $\chi_{cJ}\to nK^{0}_{S}\bar{\Lambda}+c.c.$ using the $\psi(3686)$ data sample containing $(4.48\pm 0.03)\times 10^{8}$ events collected at BESIII ref::CPC42_023001 . The radiative decays $\psi(3686)\to\gamma\chi_{cJ}$, which have a BF of approximately 10% for each $\chi_{cJ}$ ref::pdg2014 , offer an ideal environment to investigate $\chi_{cJ}$ decays. Throughout this paper, charge conjugate modes are implied unless otherwise stated.

The BESIII detector is a magnetic spectrometer Ablikim:2009aa ; Ablikim:2019hff located at the Beijing Electron Positron Collider (BEPCII) Yu:IPAC2016-TUYA01 . 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\pi$ 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 GEANT4 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 KKMC . The inclusive MC samples consist of $5.06\times 10^{8}$ $\psi(3686)$ events, the ISR production of the $J/\psi$ state, and the continuum processes incorporated in kkmc. The known decay modes are modeled with eventgen EVTGEN1 using the BFs taken from the Particle Data Group (PDG) ref::pdg2014 , and the remaining unknown decays from the charmonium states with lundcharm LUNDCHARM1 . Radiation from charged final state particles is incorporated with the photos PHOTOS package.

The signal detection efficiencies are estimated with signal MC samples. The decays of $\psi(3686)\to\gamma\chi_{cJ}$ ($J=0$, 1, 2) are simulated following Ref. ref::paper-ggJp , in which the magnetic-quadrupole (M2) transition for $\psi(3686)\to\gamma\chi_{c1,2}$ and the electricoctupole (E3) transition for $\psi(3686)\to\gamma\chi_{c2}$ are considered in addition to the dominant electric-dipole (E1) transition. For the decay $\chi_{cJ}\to nK^{0}_{S}\bar{\Lambda}+c.c.$, a special generator based on results of Helicity Partial Wave Analysis (HelPWA) SM ; ref::chung ; ref::blatt ; ref::helb3 is used. The background MC samples are obtained from inclusive MC samples, in which the signal channels are removed with TopoAna TopoAna , and the samples are normalized to the data sample based on luminosity.

The signal process $\psi(3686)\to\gamma\chi_{cJ},\chi_{cJ}\to nK^{0}_{S}\bar{\Lambda}$ with $K^{0}_{S}\to\pi^{+}\pi^{-}$ and $\bar{\Lambda}\to\bar{p}\pi^{+}$ consists of the final state particles $\gamma n\bar{p}\pi^{+}\pi^{+}\pi^{-}$. Charged track candidates from the MDC must satisfy $|\!\cos\theta|<0.93$, where $\theta$ 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(\pi)$. Exactly four charged tracks are required in each candidate event.

Since $K_{S}^{0}$ and $\bar{\Lambda}$ have relatively long lifetimes, they are reconstructed by constraining the $\pi^{+}\pi^{-}$ pair and the $\bar{p}\pi^{+}$ pair to originate from secondary vertices, respectively. Charged track candidates, except the one used as a $\bar{p}$ in the $\bar{\Lambda}$ reconstruction, are assumed to be pions without applying particle identification. The decay lengths from the secondary vertex fits of $K_{S}^{0}$ and $\bar{\Lambda}$ divided by their corresponding uncertainties are required to be larger than two. The mass distributions of the reconstructed $K_{S}^{0}$ (denoted as $M_{\pi\pi}$) and $\bar{\Lambda}$ (denoted as $M_{\bar{p}\pi^{+}}$) 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_{S}^{0}$ and $[1.112,1.120]$ GeV/$c^{2}$ for $\bar{\Lambda}$.

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\theta|<0.80$ (barrel) are required to have deposited energies larger than $25$ MeV and those with $0.86<|\!\cos\theta|<0.92$ (end cap) must have deposited energies larger than $50$ MeV. The photon candidates must be at least $10\degree$ 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 $\psi(3686)\to\gamma nK^{0}_{S}\bar{\Lambda}$, where the neutron is treated as a missing particle. The value of $\chi^{2}_{\rm 1C}$ is required to be less than 200. If more than one combination survives in an event, the one with the smallest $\chi^{2}_{\rm 1C}$ is retained.

$M_{\rm miss}$ is defined as the invariant mass of the four momentum of the missing particle $p_{\rm miss}=p_{\rm CM}-p_{\gamma}-p_{K^{0}_{S}}-p_{\Lambda}$, where $p_{i}$ is the four momentum of the particle $i$ and $p_{\rm CM}$ is the four momentum of the initial $e^{+}e^{-}$ system. To further improve the purity, $M_{\rm miss}$ before the 1-C kinematic fit is required to satisfy $0.90<M_{\rm miss}<0.98$ GeV/$c^{2}$. The distribution of $M_{\rm miss}$ before the 1-C kinematic fit is shown in Fig. 1(c).

By analyzing the $\psi(3686)$ inclusive MC samples with TopoAna TopoAna , the only significant peaking background is found to be caused by the decays $\chi_{cJ}\to\Sigma^{\pm}\bar{\Lambda}\pi^{\mp}+c.c.$ with $\Sigma^{\pm}\to n\pi^{\pm}$. The two pions in this decay may accidentally fall into the $K_{S}^{0}$ mass window and fulfill the constraint for the secondary vertex fit, thus faking the $K_{S}^{0}$. This $\Sigma^{\pm}$ background peaks in the invariant mass spectrum of $nK^{0}_{S}\bar{\Lambda}$ (denoted as $M_{nK^{0}_{S}\bar{\Lambda}}$) but distributes uniformly in the $M_{\pi\pi}$ spectrum. Other backgrounds are smoothly distributed underneath the $\chi_{cJ}$ signal.

A simultaneous unbinned maximum-likelihood fit to the $M_{nK^{0}_{S}\bar{\Lambda}}$ spectra in both the $M_{\pi\pi}$ 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 $[0.430,0.466]$ and $[0.530,0.566]$ GeV/$c^{2}$, respectively. The fit model for the signal region is

and that for the sideband region is

The signal shape $f^{J}_{\rm signal}$ for each $\chi_{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 $BW(M_{nK^{0}_{S}\bar{\Lambda}})\times E_{\gamma}^{3}\times D(E_{\gamma})$, where $BW(M_{nK^{0}_{S}\bar{\Lambda}})=((M_{nK^{0}_{S}\bar{\Lambda}}-m_{\chi_{cJ}})^{2}+0.25\Gamma^{2}_{\chi_{cJ}})^{-1}$ is the nonrelativistic Breit-Wigner function with the width $\Gamma_{\chi_{cJ}}$ and the mass $m_{\chi_{cJ}}$ of the corresponding $\chi_{cJ}$ fixed to the PDG values ref::pdg2014 , $E_{\gamma}=(m_{\psi(3686)}^{2}-M_{nK^{0}_{S}\bar{\Lambda}}^{2})/2m_{\psi(3686)}$ is the energy of the transition photon in the rest frame of $\psi(3686)$, and $D(E_{\gamma})$ is a damping factor that suppresses the divergent tail due to $E_{\gamma}^{3}$. This damping factor is described by $D(E_{\gamma})={\rm exp}(-E_{\gamma}^{2}/8\beta^{2})$, where $\beta=(65.0\pm 2.5)$ MeV is measured by the CLEO collaboration ref::beta . 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}_{\rm peakbkg}$ are parameterized the same as the signal shapes. The yields $N_{2,J}$ in the signal region are normalized to $N^{\prime}_{2,J}$ in the $M_{\pi\pi}$ sideband regions according to the sizes of these two regions. The background shapes $f^{(\prime)}_{\rm flatbkg}$ in both regions are modeled as second-order Chebyshev polynomial functions. The numbers of fitted $\chi_{cJ}$ signal events, $N_{1,J}$, are listed in Table 1.

A special generator based on results of HelPWA ref::chung ; ref::blatt ; ref::helb3 for the decay $\chi_{cJ}\to nK^{0}\bar{\Lambda}+c.c.$ is developed to estimate the detection efficiencies. The description of HelPWA can be found in the supplemental material SM . For the $\chi_{cJ}$ data events used in HelPWA, the signal regions of $M_{nK_{S}^{0}\bar{\Lambda}}$ for $\chi_{c0}$, $\chi_{c1}$, and $\chi_{c2}$ are $[3.39,3.45]$, $[3.50,3.53]$, and $[3.54,3.57]$ GeV/$c^{2}$, respectively; the masses of $K^{0}$, $\bar{\Lambda}$, and $\chi_{cJ}$ are constrained to their known masses ref::pdg2014 . 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 $\chi_{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 $\chi_{cJ}\to nK_{S}^{0}\bar{\Lambda}$ and $\chi_{cJ}\to\bar{n}K_{S}^{0}\Lambda$ 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 $\chi_{cJ}$, $\epsilon_{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 $\chi_{cJ}$ channels and are consistent within the statistical uncertainties.

The BFs for $\chi_{cJ}\to nK^{0}_{S}\bar{\Lambda}$ are calculated using

where $N_{\psi(3686)}$ is the number of $\psi(3686)$ events ref::CPC42_023001 ; $\epsilon_{J}$ is the detection efficiency as listed in Table 1; $\prod_{i}\mathcal{B}_{i}=\mathcal{B}(\psi(3686)\to\gamma\chi_{cJ})\cdot\mathcal{B}(K^{0}_{S}\to\pi^{+}\pi^{-})\cdot\mathcal{B}(\bar{\Lambda}\to\bar{p}\pi^{+})$, where the BFs are taken from the PDG ref::pdg2014 . The results are summarized in Table 1.

The number of $\psi(3686)$ events is measured to be $(4.48\pm 0.03)\times 10^{8}$ based on inclusive hadronic events, as described in Ref. ref::CPC42_023001 , so the uncertainty is 0.6%. The systematic uncertainty due to the detection of $\gamma$ is studied with the well understood channel $J/\psi\to\rho^{0}\pi^{0}$ ref::gamma-recon . The efficiency difference between data and MC simulation is about 1% per photon. To estimate the uncertainties associated with $K^{0}_{S}$ and $\Lambda$ reconstruction, the decays $J/\psi\to K^{*\pm}(892)K^{\mp}$, $K^{*\pm}(892)\to K^{0}_{S}\pi^{\pm}$ and $J/\psi\to\Lambda\bar{\Lambda}$ are selected as the control samples. The uncertainties are determined to be $1.5\%$ per $K^{0}_{S}$ and $2.0\%$ per $\Lambda$. The systematic uncertainty caused by the 1-C kinematic fit is studied with the control sample $\psi(3686)\to\pi^{0}nK^{0}_{S}\bar{\Lambda}$ with purity about 98.5%, where $n$, $K^{0}_{S}$, and $\bar{\Lambda}$ 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_{S}^{0}$, $\Lambda$, and $M_{\rm miss}$ are estimated by repeating the analysis with alternative mass window requirements. We change the mass window of $K_{S}^{0}$ to $[0.473,0.521]$ and $[0.487,0.511]$ GeV/$c^{2}$, that of $\Lambda$ to $[1.110,1,123]$ and $[1.113,1.119]$ GeV/$c^{2}$, and that of $M_{\rm miss}$ to [0.880, 1.000] and [0.910, 0.970] GeV/$c^{2}$. The largest differences from the nominal BFs are assigned as the corresponding systematic uncertainties.

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 ref::damping , given by $D(E_{\gamma})=\frac{E_{\gamma 0}^{2}}{E_{\gamma}E_{\gamma 0}+(E_{\gamma}-E_{\gamma 0})^{2}}$, where $E_{\gamma 0}=(m_{\psi(3686)}^{2}-M_{\chi_{cJ}}^{2})/2m_{\psi(3686)}$. 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 $\chi_{c0}$, $\chi_{c1}$, and $\chi_{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 $\psi(3686)\to\gamma\chi_{c0}$ ($\chi_{c1}$, $\chi_{c2}$), $K^{0}_{S}\to\pi^{+}\pi^{-}$, and $\Lambda\to p\pi^{-}$ are 2.0% (2.5%, 2.1%), 0.1%, and 0.8%, respectively, according to the PDG ref::pdg2014 .

All systematic uncertainty contributions are summarized in Table 2. The total systematic uncertainty for each $\chi_{cJ}$ decay is obtained by adding all contributions in quadrature.

The decays of $\chi_{cJ}\to nK^{0}_{S}\bar{\Lambda}+c.c.$ ($J=0$, 1, 2) are observed for the first time using $(4.48\pm 0.03)\times 10^{8}$ $\psi(3686)$ events accumulated with the BESIII detector at the BEPCII collider. The BFs of $\chi_{cJ}\to nK^{0}_{S}\bar{\Lambda}+c.c.$ are determined to be $(6.67\pm 0.26_{\rm stat}\pm 0.41_{\rm syst})\times 10^{-4}$, $(1.71\pm 0.12_{\rm stat}\pm 0.12_{\rm syst})\times 10^{-4}$, and $(3.66\pm 0.17_{\rm stat}\pm 0.23_{\rm syst})\times 10^{-4}$ for $J=0$, 1, and 2, respectively. Isospin symmetry is examined by comparing our results with the isospin conjugate decays of $\chi_{cJ}\to pK^{-}\bar{\Lambda}+c.c.$ ref::paper-pkl . The ratios $\mathcal{B}(\chi_{cJ}\to pK^{-}\bar{\Lambda}+c.c.)/\mathcal{B}(\chi_{cJ}\to nK^{0}_{S}\bar{\Lambda}+c.c.)$ are measured to be $1.98\pm 0.09_{\rm stat}\pm 0.14_{\rm syst}$, $2.64\pm 0.23_{\rm stat}\pm 0.20_{\rm syst}$, and $2.29\pm 0.13_{\rm stat}\pm 0.16_{\rm syst}$ for $J=0$, 1, and 2, respectively, where common sources of systematic uncertainties are canceled. No obvious isospin violation is observed.

Enhancements are observed in the Dalitz plots shown in Fig. 3 and the mass distributions of two-body $n\bar{\Lambda}$ 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.