Search for New Hadronic Decays of $h_c$ and Observation of $h_c\rightarrow K^{+}K^{-}\pi^{+}\pi^{-}\pi^{0}$

Ten hadronic final states of the $h_c$ decays are investigated via the process $\psi(3686)\rightarrow \pi^0 h_c$, using a data sample of $(448.1 \pm 2.9) \times 10^6$ $\psi(3686)$ events collected with the BESIII detector. The decay channel $h_c\rightarrow K^{+}K^{-}\pi^{+}\pi^{-}\pi^{0}$ is observed for the first time with a significance of $6.0 \sigma$. The corresponding branching fraction is determined to be $\mathcal{B}(h_c\rightarrow K^{+}K^{-}\pi^{+}\pi^{-}\pi^{0}) =(3.3 \pm 0.6 \pm 0.6)\times 10^{-3}$ (the first uncertainty is statistical and the second systematical). Evidence for the decays $h_c\rightarrow \pi^{+} \pi^{-} \pi^{0} \eta$ and $h_c\rightarrow K^{0}_{S}K^{\pm}\pi^{\mp}\pi^{+}\pi^{-}$ is found with a significance of $3.6 \sigma$ and $3.8 \sigma$, respectively. The corresponding branching fractions (and upper limits) are obtained to be $\mathcal{B}(h_c\rightarrow \pi^{+} \pi^{-} \pi^{0} \eta ) =(7.2 \pm 1.8 \pm 1.3)\times 10^{-3}$ $(<1.8 \times 10^{-2})$ and $\mathcal{B}(h_c\rightarrow K^{0}_{S}K^{\pm}\pi^{\mp}\pi^{+}\pi^{-}) =(2.8 \pm 0.9 \pm 0.5)\times 10^{-3}$ $(<4.7\times 10^{-3})$. Upper limits on the branching fractions for the final states $h_c \rightarrow K^{+}K^{-}\pi^{0}$, $K^{+}K^{-}\eta$, $K^{+}K^{-}\pi^{+}\pi^{-}\eta$, $2(K^{+}K^{-})\pi^{0}$, $K^{+}K^{-}\pi^{0}\eta$, $K^{0}_{S}K^{\pm}\pi^{\mp}$, and $p\bar{p}\pi^{0}\pi^{0}$ are determined at a confidence level of 90\%.

monium states are located in the transition region of perturbative and non-perturbative quantum chromodynamics (QCD) and theoretical predictions suffer therefore from large uncertainties [1][2][3][4]. The study of charmonium states and their decays is therefore crucial for gaining a deeper understanding of the intermediate-energy regime of QCD, while QCD has been tested successfully at high energies [5]. For the h c , χ cJ and η c (2S) states, most of the decay channels are still unknown. After the discovery of the spin-singlet charmonium state h c ( 1 P 1 ) in 2005 [6,7], there were only few measurements of its decays. Contrary to the fact that h c → γη c is the prominent decay channel in every calculation, the predictions of the decay process h c → light hadrons range from 14 − 48% [1][2][3][4] depending on the theoretical model. Therefore experimental measurements are needed to test and improve the theoretical models.
Experimental challenges arise from the limited statistics since these non-vector states cannot be produced directly in e + e − annihilation. The best-measured decay mode is the radiative transition h c → γη c , occurring in 51% of all decays [8][9][10], while the sum of all other known branching fractions is less than 3% [11]. Among these measurements, the multi-pionic decay h c → 2(π + π − )π 0 has been confirmed recently by the BESIII collaboration [12] after the first evidence was reported by CLEO-c [13]. Furthermore, BESIII observed the decay mode h c → ppπ + π − and reported evidence for the decay h c → π + π − π 0 . Since the previous analyses mainly studied multi-pionic final states, this analysis focuses on hadronic final states containing kaons as they could lead to intermediate resonances such as e.g. φ and exited kaon states. After radiative decays of the h c to η ( ′ ) have been observed, the study of decays involving light vector-states different from the photon will be an extension of these observations. Finally, the observation of h c → ppπ + π − motivated us to study the decay with neutral pions, as it could give additional hints on baryonic intermediate states.
From these considerations, the following ten final states are chosen to search for undiscovered decay channels of the h c : These are referenced in this manuscript by roman numbers (i, ii, ..., x). In this analysis the h c meson is produced via ψ(3686) → π 0 h c using a data sample of (448.1±2.9)×10 6 ψ(3686) events [14] collected with the BESIII detector.

II. BESIII Detector and Monte Carlo Simulation
The BESIII detector is a magnetic spectrometer [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 calorime-ter (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 the 4π solid angle. The chargedparticle momentum resolution at 1 GeV/c is 0.5%, and the 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.
geant4-based [17,18] Monte Carlo (MC) simulations are used to study the detector response and to estimate background contributions. Inclusive MC samples are produced to estimate the contributions from possible background channels. The production of the initial ψ(3686) resonance in e + e − annihilation is simulated using the MC event generator kkmc [19,20]. Its known decay modes are modeled with evtgen [17,18] using the world average branching fraction values [21], while the remaining unknown decays are generated using lundcharm [22]. geant4 is used to simulate the particle propagation through the detector system. The simulation includes the beam-energy spread and initial-state radiation (ISR) in the e + e − annihilations. In addition, exclusive MC samples containing one million events are generated using the phase-space model (PHSP) for each signal mode to optimize the selection criteria and to study the efficiency.

III. Data Analysis
Each charged track reconstructed in the MDC is required to originate from a region of 10 cm of the interaction point along the beam direction and 1 cm in the plane perpendicular to the beam. The polar angle θ of the tracks must be within the fiducial volume of the MDC |cosθ| < 0.93. Tracks used in reconstructing K 0 S mesons are exempted from these requirements and |cosθ| < 0.93 is required only for the daughter pions. The TOF and dE/dx measurements for each charged track are combined to compute particle identification (PID) confidence levels for pion, kaon and proton hypotheses. The track is assigned to the particle type with the highest confidence level which has to be larger than 0.001.
Photon candidates are reconstructed from electromagnetic showers produced in the crystals of the EMC. A shower is treated as a photon candidate if the deposited energy is larger than 25 MeV in the barrel region (|cosθ| < 0.8) or 50 MeV in the end cap region (0.86 < |cosθ| < 0.92). 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. To remove Bremsstrahlung photons, the angle between the photon and the extrapolated impact point in the EMC of the nearest charged track must be larger than 10 • for charged pions and kaons and 20 • for protons, respectively.
Following the application of a vertex fit that con- straints all charged tracks to arise from a common interaction point, a kinematic fit, constraining the total energy and momentum to the initial four momentum, is performed to further improve the momentum resolution and to suppress background. Final states containing a K 0 S undergo a secondary vertex fit, which ensures the daughter pions being produced in a common vertex. A K 0 S candidate is accepted if 487 < M (π + π − ) < 511 MeV/c 2 and λ/∆λ > 2, where ∆λ is the uncertainty on the decay length λ obtained from the secondary vertex fit.
A pair of photons is treated as a π 0 or η candidate if it satisfies |M (γγ) − M π 0 | < 30 MeV/c 2 or |M (γγ) − M η | < 50 MeV/c 2 , which corresponds to an interval of about ±3 times the mass resolution. In all final states containing π 0 or η mesons, the kinematic fit also constraints γγ pairs to have the expected nominal masses. The combination with the best χ 2 value is kept for the further analysis in case there are more than one γγ combinations.
To suppress contamination from decays with different numbers of photons, such as the dominant decay ψ(3686) → γχ c2 , where the χ c2 decays to the same final states as the h c , the following procedure is applied. The χ 2 4C,nγ value is obtained from a four-constraint fit including the expected number of photons n for a given signal hypothesis with respect to the initial four momentum. The value χ 2 4C,(n−1)γ is determined from an additional 4C fit with one missing photon compared to the desired signal process. An event is rejected if χ 2 4C,nγ > χ 2 4C,(n−1)γ . To further suppress background, the χ 2 (4+N )C value of the total kinematic fit, including additional mass constraints for π 0 , η and K 0 S candidates Table I. Applied requirements on the χ 2 (4+N)C and mass windows used as vetoes in each exclusive mode. The lower case m denotes the nominal particle mass [11]. Mode (denoted by N), is limited depending on the final state (see Table I). Additional vetoes in the π 0 π 0 , π + π − and η recoil masses, as listed in Table I, are applied to suppress background from ψ(3686) → (π 0 π 0 , π + π − , η) J/ψ. Since the π 0 from the decay ψ(3686) → π 0 h c (denoted by π 0 1 and identified among all π 0 candidates by its energy being closest to the expected) should not create any structure together with other final state particles. Therefore, additional vetoes are applied to suppress background from ω → π + π − π 0 , η → π + π − π 0 , η → π 0 π 0 π 0 , K * ± → π 0 K ± , ∆(1232) + → pπ 0 and Σ + → pπ 0 as given in Table I. The mass windows for these vetoes and for the χ 2 criterion are optimized simultaneously for each channel using the figure of merit of S/ √ S + B and are listed in Table I. Here, S denotes the number of signal events, obtained from signal MC, which is scaled to the branching fraction as determined in this analysis. Therefore, the unoptimized selection criteria were used in a first iteration to obtain a preliminary branching fraction (or upper limit). This preliminary result is then fed into the next iteration of optimization until the procedure converges within the uncertainties. The number of background events B is obtained from the ψ(3686) inclusive MC and scaled to the expected number of events. Figure 1 shows the obtained invariant mass distributions of the different decay modes. After applying all selection criteria, the remaining background originates mostly from the nonresonant production of the same final state particles as the signal and thus cannot be suppressed further.

IV. Determination of Branching Fractions
To determine the number of signal events N hc , an unbinned maximum likelihood fit to the invariant mass spectra of the particles to reconstruct the h c is performed as shown in Fig. 1. In each fit, the signal contribution is described by a Breit-Wigner function convoluted with a detector resolution function as given in Ref. [23]. Here, the mass and width of h c in the Breit-Wigner function are fixed to their world average values [11], and the parameters in the resolution function are determined with the signal MC simulation. The background shape is described by an ARGUS function [24], where the threshold parameter of the ARGUS function is fixed to the kinematical threshold of 3551 MeV/c 2 . In case of the modes h c → K + K − π + π − π 0 and h c → π + π − π 0 η, additional peaking background from the processes h c → γη c , η c → K + K − π + π − and η c → π + π − η is included and scaled to the expected number of events based on the world average values. The resulting branching fractions are determined by: (1) Here B(h c → X) denotes the branching fraction of the h c meson decaying to final state X, the branching fraction of ψ(3686) → π 0 h c is given by B(ψ(3686) → π 0 h c ) = (8.6 ± 1.3) × 10 −4 [11]. The number of ψ(3686) events is given by N ψ(3686) = (448.1 ± 2.9) × 10 6 [14]. And i B i is the product of branching fractions of the decaying particles like B(π 0 → γγ), B(η → γγ) and B(K 0 S → π + π − ) taken from [11]. The efficiency ε is obtained from signal MC simulations and N hc is the number of signal events obtained by the fit. In order to determine B(h c → X), the product branching fraction is divided by B(ψ(3686) → π 0 h c ). Both values, the product branching fraction and the h c decay branching fraction are listed in Table II.
In case no significant signal contribution is observed for the modes (iv)-(x), upper limits on the branching fractions are determined by the Bayesian approach [25]. To obtain the likelihood distribution, the signal yield is scanned using the fit function described earlier. Sys- Table II. Overview of the branching fractions and upper limits obtained in this analysis for decay processes of the hc meson. The first uncertainty shown is the statistical and the second the systematical uncertainty of the measurement method which includes the uncertainty that arises due to the use of external branching fractions.
ppπ 0 π 0 < 11.8 8.7 < 4.4 × 10 −7 < 5.2 × 10 −4 tematic uncertainties are considered by smearing the obtained likelihood curve with a Gaussian function with the width of the systematic uncertainty of the respective decay mode. The upper limit at a confidence level of 90% is finally obtained by: The upper limit on the number of observed events N up hc is determined by integrating the smeared likelihood function L(N ) up to the value N up hc , which corresponds to 90% of the integral. The results are listed in Table II.
Among the ten final states, the decay h c → K + K − π + π − π 0 is observed with a statistical significance of 6σ and evidences for the decays h c → π + π − π 0 η and h c → K 0 S K ± π ∓ π + π − are found with statistical significances of 3.6σ and 3.8σ, respectively. The combined significance of the modes (iv)-(x) is determined to be 3.5σ. The statistical significance is determined by the likelihood ratio between a fit with and without signal component by taking the change in the number of fit parameters into account.
For the final state K + K − π + π − π 0 in the h c decay, a search for intermediate resonances is performed to obtain information about underlying sub-processes. Despite the large background contamination, the signal content is determined by an unbinned maximum likelihood fit to the invariant K + K − π + π − π 0 mass in slices of masses of possible sub-systems. The resulting distributions are shown in the Fig. 2.
No firm conclusions about contributions of intermediate resonances can be drawn based only on the extracted projections with the present statistics. The Kπ distribution shows a possible structure in the K * (892) region, which may signal the production of this resonance. In the invariant K ± π ∓ π 0 mass distribution there may be a hint in the mass region of 1.9 − 2.0 GeV/c 2 for the production of an excited kaon, such as the K * 2 (1980) or K 2 (1820). Conceivable sub-processes would then be h c → K * (892)/K * 0,2 (1430) (K 2 (1820)/K * 2 (1980)). As shown in this analysis, evidence for the decay h c → K 0 S K ± π ∓ π + π − has been found.

V. Systematic Uncertainties
The sources of systematic uncertainties for the branching fractions include tracking, PID, selection, uncertainties caused by the signal fitting procedure and efficiency determination, and they are explained in the following. All the systematic uncertainties are summarized in Table III. The overall systematic uncertainty for the product branching fraction B(ψ(3686) → π 0 h c ) × B(h c → X) is obtained by summing all individual components (from column two to six of Table III) in quadrature. An additional systematic uncertainty of ∆ ext =15.1% is added in quadrature due to the branching fraction of ψ(3686) → π 0 h c used in the calculation of branching fractions of B(h c → X).
kaons and four photons reconstructed as π 0 are involved. This gives the following contributions to the systematic uncertainty: 4% (1% per track) due to PID, 4% (1% per track) due to track reconstruction, 4% due to photon reconstruction (1% per photon) and additional 2% for π 0 reconstruction (1% for each π 0 ). The total uncertainty due to PID and event reconstruction is given for each final state in the second column of Table III. The systematic uncertainty due to the selection criteria is determined by varying the nominal selection criteria. For each mass window requirement, the nominal value of the criterion is varied by ±10 MeV/c 2 in increments of ±0.5 MeV/c 2 . The maximum deviation from the nominal branching fraction is quoted as a systematic uncertainty as given in column three of Table III. The uncertainties associated with the kinematic fit are determined by comparing the efficiencies with and without the helix parameter correction. For charged particles, differences in the χ 2 distributions of the kinematic fit between data and MC have been studied by using a control sample J/ψ → φf 0 (980), φ → K + K − , f 0 (980) → π + π − , which ensures high statistics and purity [31]. The helix parameters of the corresponding tracks are corrected accordingly [31]. The difference between the result determined with and without this correction applied are assigned as systematic uncertainty of the kinematic fit.
The systematic uncertainty due to the physics model and the efficiency used to simulate signal MC arises from the limited knowledge of intermediate states in the h c decay. Therefore have MC samples been generated including additional intermediate states and the results are compared with those of the nominal phase space sample. The systematic uncertainty of the fit results from the choice of background parametrization is determined by using a second order Chebychev polynomial. In case of a peaking background, the uncertainty of the branching fraction (e.g. B(η c → π + π − K + K − ) = (6.9 ± 1.1) × 10 −3 ) [11] has been used to determine the uncertainty of the peaking background. This uncertainty contributes 1.2% in case of h c → K + K − π + π − π 0 and 2% in h c → π + π − π 0 η, respectively. Another contribution to the systematic uncertainty of the fit is caused by the limited range of the invariant mass in which the fit is applied. Therefore, the fit range is extended from 3.50 − 3.55 GeV/c 2 to 3.40 − 3.65 GeV/c 2 , and the difference of the branching fraction result is used as a systematic uncertainty. Further uncertainties arise from the parametrization of the resolution distributions. Instead of using the default parametrization, a Crystal Ball distribution has been used. The uncertainty arising from the number of ψ(3686) events is 0.7% [14].

VI. Summary
In this analysis, ten final states of the h c decays have been searched for using a data sample of (448.1±2.9)×10 6 ψ(3686) events collected at BESIII. The decay h c → K + K − π + π − π 0 is observed for the first time. Furthermore, evidence for the decays h c → π + π − π 0 η and h c → K 0 S K ± π ∓ π + π − are found with statistical significances of 3.6σ and 3.8σ, respectively. The combined significance of the modes (iv)-(x) is determined to be 3.5σ. Upper limits are determined in case there was no signal observed. The measured branching fractions and upper limits at 90% confidence level are listed in Table II. Summing up the branching fractions obtained in this analysis, shows that these decays contribute at a level of ∼ 1.3% to all decays and contribute at the same level Table III. Summary of systematic uncertainties. Here ∆BB = Σi∆ 2 i is the systematic uncertainty for B(ψ(3686) → π 0 hc) × B(hc → X) and ∆B = Σi∆ 2 i + ∆ 2 ext is that for B(hc → X) . ∆i represents the individual uncertainties given in columns two to six and ∆ext=15.1% is an additional uncertainty for B(hc → X) due to the external uncertainty of B(ψ(3686) → π 0 hc) [11]. All values are given in %.
tion Fit (i)