Search for the decay $h_c\rightarrow\pi^0J/\psi$

A search for the decay $h_c\rightarrow\pi^0J/\psi$ is performed using a sample of $h_c$ produced in the reaction $e^+e^-\rightarrow\pi^+\pi^-h_c$. The data samples were collected with the BESIII detector at center-of-mass energies between 4.189 and 4.437 GeV, corresponding to a total integrated luminosity of 11 fb$^{-1}$. No significant signal is observed. Upper limits on the branching ratio $\mathcal{B}(h_c\rightarrow\pi^0J/\psi)/\mathcal{B}(h_c\rightarrow\gamma\eta_c\rightarrow\gamma K^+K^-\pi^0)$ and on the branching fraction $\mathcal{B}(h_c\rightarrow\pi^0J/\psi)$ are determined to be $7.5\times10^{-2}$ and $4.7\times10^{-4}$ at $90\%$ confidence level, respectively. The latter is derived from the former using the measured branching fraction of the normalization channel. This is the first determination of the upper limit of the decay $h_c\rightarrow\pi^0J/\psi$.


Introduction
Charmonium, the bound state of a charm quark and anticharm quark (cc), plays an important role in our understanding of quantum chromodynamics (QCD), which is the fundamental theory of the strong interaction.Low energy QCD remains a field of high interest both experimentally and theoretically.All charmonium states below the open-charm (D D) threshold have been observed experimentally and can be well described by QCD inspired models [1].However, knowledge is still sparse on the P -wave spin-singlet state, h c (1P ) [2].So far, only a few decay modes of the h c have been observed, such as h c → γη c (η ) [2].Searches for new decay modes of the h c can provide useful information to constrain theoretical models in the charmonium region.
In 1992, the E760 Collaboration reported an evidence of the h c in the π 0 J/ψ decay mode [3].The decay h c → π 0 J/ψ was not confirmed by the successor experiment E835 with higher statistics in 2005 [4], but E835 Collaboration found the evidence for h c in another decay mode γη c .More measurements are thus needed for clarification.In addition, several theoretical articles have addressed the decay h c → π 0 J/ψ, with predictions for the partial width around several keV [5][6][7][8].Until now, there is no experimental result to confirm this.
Many studies of the h c have been performed using the reaction ψ(3686) → π 0 h c from the ψ(3686) data sample at BESIII [9][10][11].However, it is hard to search for the decay h c → π 0 J/ψ using the ψ(3686) data sample due to the large background from ψ(3686) → π 0 π 0 J/ψ.BESIII has observed a sizeable cross section for the process e + e − → π + π − h c between 4.189 and 4.437 GeV [12].This process is advantageous in the search for the decay h c → π 0 J/ψ, because it avoids the dominant background of ψ(3686) → π 0 π 0 J/ψ present in the ψ(3686) decay.
In this paper, a search for h c → π 0 J/ψ is reported using e + e − → π + π − h c events from data samples collected at center-of-mass energies between 4.189 and 4.437 GeV with the BESIII detector [13] with a total integrated luminosity of 11 fb −1 .The J/ψ is reconstructed in its decay to a + − pair ( = e or µ), and the π 0 is reconstructed via the decay π 0 → γγ.Additionally, the decay h c → γη c with η c → K + K − π 0 is used as the normalization channel.

BESIII detector and Monte Carlo simulation
The BESIII detector [13] records symmetric e + e − collisions provided by the BEPCII storage ring [14], 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.95 GeV.BESIII has collected large data samples in this energy region [15].The cylindrical core of the BESIII detector covers 93% of the full solid angle and consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator timeof-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 identification modules interleaved with steel.The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for 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 in the TOF barrel region is 68 ps, while that in the end cap region is 110 ps.The end cap TOF system was upgraded in 2015 using multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [16].About 70% of the data sample used here was taken after this upgrade.Simulated data samples produced with a geant4-based [17] 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 [18].The inclusive MC sample includes the production of open charm processes, the ISR production of vector charmonium(-like) states, and the continuum processes incorporated in kkmc [18].The known decay modes are modelled with evtgen [19] using branching fractions taken from the Particle Data Group (PDG) [2], and the remaining unknown charmonium decays are modelled with lundcharm [20].Final state radiation (FSR) from charged final state particles is incorporated using photos [21].Signal MC samples for e + e − → π + π − h c are generated using isotropic phase space populations at each center-of-mass energy point, assuming that the input line-shape follows a coherent sum of Y (4220) and Y (4390) Breit-Wigner functions, whose parameters are fixed to the measured values in Ref. [12].The subsequent h c → π 0 J/ψ decay is generated uniformly in phase space, and the transition h c → γη c is generated with an angular distribution of 1 + cos 2 θ * , where θ * is the angle of the photon with respect to the h c helicity direction in the h c rest frame.Furthermore, the three-body process η c → K + K − π 0 is simulated according to the measured two-body invariant-mass distributions [10].

Event selection
For each charged track, the distance of closest approach to the interaction point is required to be within ±10 cm in the beam direction and within 1 cm in the plane perpendicular to the beam direction.The polar angle (θ) of the charged tracks must be within the fiducial volume of the MDC (|cos θ| < 0.93).Photons are reconstructed from isolated showers in the EMC, which are at least 10 • away from the nearest charged track.The photon energy is required to be at least 25 MeV in the barrel region (|cos θ| < 0.80) or 50 MeV in the end cap region (0.86 < |cos θ| < 0.92).To suppress electronic noise and energy depositions unrelated to the event, the time at which the photon is recorded in the EMC is required to be within 700 ns of the event start time.
Since the decays h c → π 0 J/ψ and h c → γη c result in the final states γγπ + π − + − and γγγπ + π − K + K − , respectively, candidate events are required to have four charged tracks with zero net charge, at least two photons for h c → π 0 J/ψ, and at least three for h c → γη c .For the h c → π 0 J/ψ decay, tracks with momenta larger than 1.0 GeV/c are assigned to be leptons from J/ψ decay.Otherwise, they are considered as pions.Leptons from the J/ψ decay with an energy deposited in the EMC larger than 1.0 GeV are identified as electrons, and those with less than 0.4 GeV as muons.For the h c → γη c decay, information from TOF and dE/dx measurements is combined to form particle identification (PID) likelihoods for the π, K and p hypotheses.Each track is assigned with a particle type corresponding to the hypothesis of the highest PID likelihood.Exactly two oppositely charged π and K each are required in each event.To reduce the background contributions and to improve the mass resolution, a five-constraint kinematic fit is performed for the two channels.The total fourmomentum is constrained to the initial four-momentum of the e + e − system.Additionally, the invariant mass of the two photons from the π 0 decay is constrained to the π 0 mass taken from the PDG [2].If there is more than one candidate in an event, the one with the smallest χ 2 of the kinematic fit is selected.The χ 2 is required to be less than 30 for h c → π 0 J/ψ or 40 for h c → γη c .

Branching fraction measurement
Distributions of the invariant mass of + − (or versus the recoil mass of π + π − , RM (π + π − ) for h c → π 0 J/ψ (or γη c ) are shown in Fig. 1.Here, RM (π + π − ) = (P e + e − − P π + − P π − ) 2 , where P e + e − and P π ± are the four-momenta of the initial e + e − system and the π ± , respectively.No significant signal is observed for  the h c → π 0 J/ψ decay.As expected, a high-density area can be observed originating from the h c → γη c decay.
Figure 2 shows the distributions of RM (π + π − ) for the decays h c → π 0 J/ψ and h c → γη c for data in the J/ψ and η c signal regions.No significant signal is seen for the h c → π 0 J/ψ decay, while a clear peak is present for the h c → γη c decay.The green shaded histograms correspond to the normalized events from the J/ψ and η c sideband regions.No significant peaks are found in sideband events, and the sideband is not used in the following fit.A detailed study of the inclusive MC sample indicates that there are no peaking background contributions in the h c signal region [22].In order to extract the h c signal yield, a simultaneous unbinned maximum-likelihood fit to the RM (π + π − ) in two decay channels is performed.The signal shape of RM (π + π − ) is modeled by a shape obtained from the simulation convolved with a Gaussian function.The mean value and width of the Gaussian function are allowed to float, but are constrained to be the same for the two channels in the simultaneous fit.The background is described by a first-order polynomial function.The solid curves in Fig. 2 show the fit results.The h c signal yields are 0.5 ± 1.6 and 451.6±26.7 for the decays h c → π 0 J/ψ and h c → γη c , respectively.Since no significant signal is observed for the h c → π 0 J/ψ decay, an upper limit at 90% confidence level (C.L.) using the Bayesian method is given.With the fit function described above, the h c → π 0 J/ψ signal yield is scanned to obtain the likelihood distribution, which is then convolved with the systematic uncertainty.The upper limit on the h c → π 0 J/ψ signal yield N up π 0 J/ψ at 90% C.L. is obtained by solving the equation , where x is h c → π 0 J/ψ signal yield and F (x) is the probability density function of the likelihood distribution.The upper limit N up π 0 J/ψ is determined to be 4.8.Results of the simultaneous fit to the RM (π + π − ) distributions from the decays h c → π 0 J/ψ (left panel) and h c → γη c (right panel).The red solid lines are the total fit results and the blue dashed lines are the background components.The green shaded histograms correspond to the normalized events from the J/ψ and η c sideband regions.

The branching ratio B(h
where N is the yield of signal events, L is the integrated luminosity [23], σ is the cross section of e + e − → π + π − h c [12], 1 + δ is the radiative correction factor [18,24], is the efficiency, B(J/ψ → + − ) is the branching fraction of J/ψ → + − from the PDG [2], and i denotes each energy point.Table 1 shows the luminosity, cross section, radiative correction factor and efficiency at each center-of-mass energy.So the ratio is 0.847, and the upper limit on the branching ratio With the world averages of B(h c → γη c ) and B(η c → K + K − π 0 ) from the PDG [2], which is B(h c → γη c ) = (51 ± 6)% and B(η c → K + K − π 0 ) = B(η c → KKπ)/6 = (1.217± 0.067)% according to the isospin symmetry, an upper limit on the branching fraction B(h c → π 0 J/ψ) at 90% C.L. can be determined to be B(h c → π 0 J/ψ) < 4.7 × 10 −4 .The systematic uncertainties have been considered in these upper limit calculations.

Systematic uncertainty
The sources of systematic uncertainty related to the branching ratio B(h c → π 0 J/ψ)/B(h c → γη c → γK + K − π 0 ) and branching fraction B(h c → π 0 J/ψ) are summarized in Table 2, where the uncertainties associated with the charged track PID efficiencies and normalization channel yields, as well as the branching fraction of the normalization channel used to obtain the absolute branching fraction of h c → π 0 J/ψ.
The uncertainty in the photon efficiency and charged track PID efficiency is 1% per photon or per track [25][26][27].The uncertainty due to the kinematic fit is estimated by correcting the helix parameters of charged tracks, and the difference between the results with and without this correction is taken as the uncertainty [28].To estimate the uncertainty related to the input line-shape of the process e + e − → π + π − h c , the input line-shape is changed to a flat line-shape and the difference is taken as the uncertainty.The process e + e − → π + π − h c is generated by a three-body phase-space model.The uncertainty related to the MC decay model is obtained by substitution with a model e + e − → π ± Z c (4020) ∓ → π + π − h c .The angular distribution of h c → π 0 J/ψ is varied from phase space to 1 ± cos 2 θ * * , where θ * * is the angle of the π 0 with respect to the h c helicity direction in the h c rest frame.These two items are combined as the total uncertainty from the MC decay model.
The uncertainty from the J/ψ mass window requirement is estimated using e + e − → γ ISR ψ(3686) → γ ISR π + π − J/ψ events [29].The efficiency difference between data and MC simulation is taken as the systematic uncertainty, which arises from the different mass resolutions in the data and the simulation.The uncertainty from the η c mass window requirement is obtained by shifting the η c mass window by ±10 MeV/c 2 , and we take the difference of B(h c → γη c → γK + K − π 0 ) to the nominal one as the systematic uncertainty.These two items are combined as the total uncertainty of the mass window.
The number of events for the normalization channel h c → γη c is determined to be N γηc = 451.6 ± 26.7.The uncertainty from the yield of the normalization channel is thus 5.9%.The uncertainties from the branching fractions B(J/ψ → l + l − ), B(h c → γη c ) and B(η c → K + K − π 0 ) are taken from the PDG [2].
The overall systematic uncertainties except for the fit procedure are obtained by adding all sources of systematic uncertainties in quadrature, assuming they are uncorrelated, and are summarized in Table 2.The sources of uncertainty in the fit procedure include the fit range and background shape.The limit of the fit range is varied by ±5 MeV/c 2 , and the background shape is replaced from the first-order polynomial to a constant function or a second-order polynomial function.For the uncertainty from the fit procedure for the upper limit measurement, different combinations of fit range and background shape are used to get the upper limits, and the largest upper limit is chosen as the nominal one.Then the likelihood distribution of the nominal upper limit is convolved with the above overall systematic uncertainties to get the final upper limit.

Summary and discussion
The decay h c → π 0 J/ψ is searched for using the process e + e − → π + π − h c with data samples collected at center-of-mass energies between 4.189 and 4.437 GeV with the BESIII detector corresponding to a total integrated luminosity of 11 fb −1 .No significant signal is observed for the decay channel h c → π 0 J/ψ.The upper limits on the branching ratio B(h c → π 0 J/ψ)/B(h c → γη c → γK + K − π 0 ) and the branching fraction B(h c → π 0 J/ψ) at 90% confidence level are determined to be 7.5 × 10 −2 and 4.7 × 10 −4 , respectively.The latter is derived from the former using the measured branching fraction of the normalization channel.This is the first upper limit measurement of the branching fraction for the decay h c → π 0 J/ψ.The measured results are not consistent with the measurements by the E760 Collaboration [3], while in agreement with the E835 Collaboration [4].The h c total width is 0.7 ± 0.4 MeV from the PDG [2], and we take 1.1 MeV as the h c total width conservatively.Therefore, the upper limit on the partial width for h c → π 0 J/ψ is conservatively determined to be Γ(h c → π 0 J/ψ) < 0.52 keV, which is one order-of-magnitude lower than the current theoretical predictions (several keV) [5][6][7][8].
Figure 2.Results of the simultaneous fit to the RM (π + π − ) distributions from the decays h c → π 0 J/ψ (left panel) and h c → γη c (right panel).The red solid lines are the total fit results and the blue dashed lines are the background components.The green shaded histograms correspond to the normalized events from the J/ψ and η c sideband regions.

Table 1 .
The luminosity L, cross section σ, radiative correction factor 1 + δ and efficiency at each center-of-mass energy √ s.
100871, People's Republic of China i Also at School of Physics and Electronics, Hunan University, Changsha 410082, China j Also at Guangdong Provincial Key Laboratory of Nuclear Science, Institute of Quantum Matter, South China Normal University, Guangzhou 510006, China k Also at Frontiers Science Center for Rare Isotopes, Lanzhou University, Lanzhou 730000, People's Republic of China l Also at Lanzhou Center for Theoretical Physics, Lanzhou University, Lanzhou 730000, People's Republic of China m Henan University of Technology, Zhengzhou 450001, People's Republic of China