Search for new hadronic decays of hc and observation of hc → pp¯η\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p\overline{p}\eta $$\end{document}

A search for the hadronic decays of the hc meson to the final states pp¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p\overline{p} $$\end{document}π+π−π0, pp¯η\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p\overline{p}\eta $$\end{document}, and pp¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p\overline{p} $$\end{document}π0 via the process ψ(3686) → π0hc is performed using (4.48 ± 0.03) × 108ψ(3686) events collected with the BESIII detector. The decay channel hc → pp¯η\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p\overline{p}\eta $$\end{document} is observed for the first time with a significance greater than 5σ and a branching fraction of (6.41 ± 1.74 ± 0.53 ± 1.00) × 10−4, where the uncertainties are statistical, systematic, and that from the branching fraction of ψ(3686) → π0hc. Strong evidence for the decay hc → pp¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p\overline{p} $$\end{document}π+π−π0 is found with a significance of 4.9σ and a branching fraction of (3.84 ± 0.83 ± 0.69 ± 0.58) × 10−3. The significances include systematic uncertainties. No clear signal of the decay hc → pp¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p\overline{p} $$\end{document}π0 is found, and an upper limit of 6.59 × 10−4 on its branching fraction is set at the 90% confidence level.


BESIII detector and Monte Carlo simulation
The BESIII detector is a magnetic spectrometer [10] located at the Beijing Electron Positron Collider (BEPCII) [11]. The cylindrical core of the BESIII detector consists of a helium-based multi-layer 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 chamber muon identifier modules interleaved with steel.
The acceptance of charged particles and photons is 93% over the 4π solid angle. The charged-particle momentum resolution at 1 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 section is 68 ps, while that of the end cap section is 110 ps. The end cap TOF system was upgraded in 2015 with multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [12,13]. About 70% of the data sample used here was taken after this upgrade.
Simulated data samples produced with geant4-based [14] Monte Carlo (MC) software, which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate the background contributions. Inclusive MC samples are produced to estimate the contributions from possible background channels. The simulation includes the beam energy spread and initial-state radiation (ISR) in the e + e − annihilations modeled with the generator kkmc [15,16]. The ISR production of vector charmonium(like) states and the continuum processes are also incorporated in kkmc [15,16]. The known decay modes are modeled with evtgen [17,18], using branching fractions from the PDG [19], and the remaining unknown decays are generated with lundcharm [20]. Final state radiation from charged final state particles is incorporated with photos [21]. For the exclusive MC simulation samples, the three channels of interest are generated using the phase-space model (PHSP) for each signal mode.

Event selection and data analysis
Each charged track reconstructed in the MDC is required to originate from a region of 10 cm from the interaction point (IP) along the z axis, which is the symmetry axis of the MDC, and 1 cm in the plane perpendicular to z. The polar angle θ with respect to the z axis of the tracks must be within the fiducial volume of the MDC, |cos θ| < 0.93. The measurements of flight time in the TOF and dE/dx in the MDC 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, and that level is required to be greater than 0.001. Finally, a vertex fit constraining all charged tracks to come from a common IP is made.
Photon candidates are reconstructed from isolated electromagnetic showers produced in the crystals of the EMC. A shower is treated as a photon candidate if the deposited energy -2 -

JHEP05(2022)108
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.
The π 0 candidates are reconstructed from γγ combinations with invariant mass within (0.080, 0.200) GeV/c 2 . To improve the momentum resolution and suppress the wrong combination background, a one-constraint (1C) kinematic fit is performed to constrain the γγ invariant mass to the nominal π 0 mass [19], and the goodness-of-fit χ 2 1C (γγ) is required to be less than 20. The η candidates are reconstructed from γγ combinations with invariant mass within (0.450, 0.650) GeV/c 2 , and the invariant mass of the photon pair is constrained to the nominal η mass [19] with a 1C-kinematic fit requiring χ 2 1C (γγ) < 200. In order to reduce background events and to improve the mass resolution, a sixconstraint (6C) kinematic fit is performed constraining the final state energy-momentum to the total initial four-momentum of the colliding beams, and constraining the masses of the two π 0 s to the known π 0 mass in the h c → ppπ + π − π 0 and h c → ppπ 0 decays or constraining the masses of the π 0 and η mesons to their known values in the h c → ppη decay. The combination with the smallest value of the 6C-kinematic fit quality χ 2 6C is kept for further analysis. The χ 2 6C values for h c → ppπ + π − π 0 , h c → ppη and h c → ppπ 0 decays are required to be less than 45, 45, and 64, respectively. These values are obtained by optimizing the figure-of-merit (FOM) defined as S/ √ S + B, where S denotes the normalized number of signal events, obtained from MC simulation, while B is the number of background events, obtained from inclusive MC samples.
The J/ψ-related background is vetoed by requirements on the π 0 π 0 , π + π − , and η recoil masses. The bachelor π 0 from the decay ψ(3686) → π 0 h c is identified by its energy being closest to the expected energy, which is determined to be about 0.16 GeV according to ψ(3686) two-body decay. Furthermore, the bachelor π 0 when combined with other final state particles should not be consistent with coming from any resonance. Therefore, additional vetoes are applied to suppress background from ω → π + π − π 0 , η → π + π − π 0 , Σ + → pπ 0 andΣ − →pπ 0 , as given in table 1, where M and m denote the invariant mass and the known mass [19] of the indicated particle, respectively; RM denotes the recoiling mass, where p ψ(3686) and p X are the four momentums of ψ(3686) and X, respectively. The Λ/Λ-related background is also suppressed by requiring the invariant mass of pπ − /pπ + to be out of the Λ/Λ mass window, and the K 0 S background is rejected by requiring m π + π − to be out of the K 0 S mass window. All mass windows are obtained by optimizing the FOM and are listed in table 1. No significant intermediate process signal is observed in the study. (I) After applying all selection criteria, the invariant mass distributions for the three h c exclusive decay modes are shown in figure 1. Potential background channels from the inclusive MC sample are identified by TopoAna [22], which shows that the remaining background mainly originates from resonant production with the same final state particles as the signal. To investigate possible background from continuum processes, the same selection criteria are applied to a data sample of 44 pb −1 collected below the ψ(3686) resonance at √ s = 3.65 GeV. Only a few events survive in modes I and III, but they are outside the h c signal region.

Result
To determine the number of h c signal events N sig in each decay mode, unbinned maximum likelihood fits are performed to the corresponding mass spectra as shown in figure 1. In all fits, the signal distribution is described by a MC-simulated shape convolved with a Gaussian function accounting for the mass resolution difference between data and MC simulation. The background shape is described by an ARGUS function [23], where the threshold parameter of the ARGUS function is fixed to the kinematical threshold of 3.551 MeV/c 2 . The branching fractions of h c → ppX are determined by Here, i B i is the product of branching fractions of the decaying particles like B(π 0 → γγ) and B(η → γγ) taken from the PDG [19]. The number of ψ(3686) events is determined to be N ψ(3686) = (448.1 ± 2.9) × 10 6 [9]. The detection efficiencies ε are obtained from signal MC simulations, and are determined to be 6.0%, 18.1%, and 23.0% for the three decay modes, respectively.

JHEP05(2022)108
In case of mode II, there are two η decay modes, η → γγ and η → π + π − π 0 . A simultaneous unbinned maximum likelihood fit is performed to determine the branching fraction B(h c → ppη), which is taken as the common parameter among the different decay modes. The corresponding number of h c signal events in the two different final states is calculated by: For the η → γγ mode, an additional normalized peaking background component from h c → γη c , η c → ppπ 0 is included. For the η → π + π − π 0 mode, the accepted candidate events require the invariant mass of π + π − π 0 to be in the η signal region, i.e., 532 < M π + π − π 0 < 562 MeV/c 2 . The corresponding η side-band shows no obvious peaking background. The numerical results for B(h c → ppX) and the resulting branching fraction For mode III, no significant signal is observed, and an upper limit on the branching fraction is determined by a Bayesian approach [24]. To obtain the likelihood distribution, the signal yield is scanned using the fit function, eq. (4.1). Systematic 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 the 90% confidence level on the number of 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, (4. 3) The results are listed in table 2 and 3. Among the three h c decay modes, mode I is observed with a statistical significance of 5.1 standard deviations (σ). The significance for mode II is also determined to be 5.1σ by combining the two η decay modes, while the significance of mode III is 1σ. The statistical significance is estimated by the likelihood difference between the fits with and without signal component, taking the change in the degrees of freedom into account. To evaluate the effect of the systematic uncertainty on the signal significance, we repeat the fits with variations of the signal shape, background shape, and fit range, and find the statistical significance of mode II to be always larger than 5σ, and mode I to be larger than 4.9σ.

Systematic uncertainty
The sources of systematic uncertainties for the branching fractions include tracking, photon detection, π 0 reconstruction, PID, the kinematic fit, mass windows, fitting procedure, the branching fraction of the intermediate decay, the number of ψ(3686) events and the physics model describing the h c production and decay dynamics. All the systematic uncertainties are summarized in table 4 for modes I and III, and in table 5 for mode II. The overall systematic uncertainty for the product branching B(ψ(3686) → π 0 h c ) · B(h c → ppX) is obtained by summing all individual components in quadrature. The third uncertainty -5 -

JHEP05(2022)108
Significance(σ) 4.9 - Table 2. The number of observed signal events N hc , the absolute branching fraction B(h c → ppX), the product branching fraction B(ψ(3686) → π 0 h c ) × B(h c → ppX), and the statistical significance, including systematic uncertainties. Here, the first uncertainty is statistical, the second is systematic, and the third one arises from the branching fraction of ψ(3686) → π 0 h c [19].
Mode II for B(h c → ppX) of ∆ ext = 15.1% is due to the uncertainty of the branching fraction of ψ(3686) → π 0 h c [19].
• π 0 and η reconstruction efficiencies. There are two π 0 candidates with different momentum distributions for mode I and III, while there is only one π 0 candidate in mode II. The uncertainties of the two π 0 reconstructions differ due to the different momentum distributions. The corresponding uncertainties for the three decay modes are determined to be 0.6%, 0.3%, 1.8%, respectively. The uncertainty due to η reconstruction from γγ final states is determined by using a high purity control sample of J/ψ → ηpp decays [26]. The difference of η reconstruction efficiency between data and MC simulation gives an uncertainty of 1% per η.
• PID. The uncertainty due to PID is determined to be 1.0% per pion [27], 1.3% per proton, and 1.6% per antiproton, based on the same samples used to estimate tracking (d) Figure 1. Fits to the invariant mass spectra of (a) ppπ + π − π 0 and (b) ppπ 0 and simultaneous fits to the invariant mass spectra of (c) ppη with η → π + π − π 0 and (d) η → γγ for data. Data are shown as points with error bars, the total fit result is shown by the red solid line, the background contribution is denoted by the green dashed line (including the peaking background contribution shown barely discernible in pink in (d)), and the signal contribution is illustrated by the gold dashed line. The background obtained from inclusive MC samples is shown by the gray shaded histogram. efficiencies [6]. Tables 4 and 5 list the relative systematic uncertainties due to PID for the different decay modes.
• Kinematic fit. The uncertainties associated with the kinematic fit are studied with the track helix parameter correction method, as described in ref. [28]. In the standard analysis, these corrections are applied. The difference of the MC signal efficiencies without corrected track parameters are taken as the corresponding systematic uncertainties.
• Mass windows. The uncertainties associated with the mass windows are estimated by repeating the analysis with alternative mass window requirements. The largest differences from the nominal branching fractions are assigned as the corresponding systematic uncertainties. In addition, the systematic uncertainties due to the Λ/Λ, Σ + /Σ − mass window requirements are estimated by the control samples, ψ(3686) → ΛΛπ 0 π 0 and ψ(3686) → ΣΣπ + π − , respectively. The differences of the selection efficiencies between data and MC simulation from control samples are taken as the corresponding systematic uncertainties.
-7 -JHEP05(2022)108 • Fit range. The uncertainty due to the fitting range is obtained by changing the range by ±0.01 GeV/c 2 , and the largest difference in the branching fraction is taken as the systematic uncertainty.
• Signal shape. The uncertainty due to signal shape is estimated by replacing the MC-simulated shape convolved with a Gaussian function by only the MC-simulated shape. The difference in the measured branching fraction is taken as the systematic uncertainty.
• Background shape. The uncertainties caused by the background shape are estimated by alternatively using different shapes. For mode I, we use an MC-simulated shape which is obtained by an MC sample of ψ(3686) → ph 1 (1170)∆ + + c.c. and ψ(3686) → ∆ ++ ∆ −− π 0 π 0 , to replace the ARGUS function. For mode II, we use a 2 nd -order Chebychev polynomial to replace the nominal ARGUS function. The differences in the measured branching fractions are taken as the systematic uncertainties.
The uncertainty due to η → γγ peaking background in mode II is taken from the uncertainty on the branching fraction. For mode III, we use the ARGUS function to replace the nominal MC-simulated shape.
• Intermediate decays. The maximum deviations from the known branching fractions are taken as systematic uncertainties as listed in tables 4 and 5.
• N ψ(3686) . The uncertainties due to the number of ψ(3686) events (N ψ(3686) ) are determined with inclusive hadronic ψ(3686) decays and estimated to be 0.7% [9]. → π 0 h c , we take the efficiency differences between the signal MC generated from PHSP and the HELAMP model [17,18] as the systematic uncertainties. The uncertainties for the three decay modes are determined to be 2.4%, 0.7%, 0.4%, respectively.
Tables 4 and 5 summarize all the systematic uncertainties of the different decay modes. The overall systematic uncertainties are obtained by adding all systematic uncertainties in quadrature assuming they are independent. For mode II, there are two η decay channels. Therefore, the uncommon items are determined using the weighted average of the detection efficiency and the branching fraction of the subsequent decays in individual decay modes [29].

Summary
Using a data sample of (448.1 ± 2.9) × 10 6 ψ(3686) events collected with the BESIII detector, three decay modes of the h c have been searched for. The decay channel h c → ppη is -8 -
-9 -observed for the first time with a 5.1σ statistical significance, and evidence for the decay h c → ppπ + π − π 0 is found with a statistical significance of 4.9σ. No obvious signal for h c → ppπ 0 is seen. The product branching fractions, B(ψ(3686) → π 0 h c ) × B(h c → ppX), and the absolute branching fractions, B(h c → ppX), are listed in table 2. The branching fractions obtained in this analysis are at the level of ∼ 10 −3 , which is the same level as the previously observed decays of h c → 2(π + π − )π 0 , ppπ + π − [6] and h c → K + K − π + π − π 0 , π + π − π 0 η, K 0 S K ± π ∓ π + π − [7]. These measurements are essential to test the theoretical prediction [2]. Finally, it is still unclear whether the hadronic decay width of the h c is of the same order as the radiative decay width predicted in [30]. Future experimental measurements searching for more decay modes based on the larger data set of ψ(3686) events [31], together with improved theoretical calculations can help us to answer this question.