Search for the $X(2370)$ and observation of $\eta_{c}\to\eta\eta\eta^\prime$ in $J/\psi\to\gamma\eta\eta\eta^{\prime}$

Using a sample of $1.31\times10^{9} ~J/\psi$ events collected with the BESIII detector, we perform a study of $J/\psi\to\gamma\eta\eta\eta^{\prime}$ to search for the $X(2370)$ and $\eta_{c}$ in the $\eta\eta\eta^{\prime}$ invariant mass distribution. No significant signal for the $X(2370)$ is observed, and we set an upper limit for the product branching fraction of ${\cal B}(J/\psi\to\gamma X(2370)\cdot{\cal B}(X(2370)\to\eta\eta\eta^{\prime})<9.2\times10^{-6}$ at the 90% confidence level. A clear $\eta_{c}$ signal is observed for the first time, yielding a product branching fraction of ${\cal B}(J/\psi\to \gamma\eta_{c})\cdot{\cal B}(\eta_{c}\to \eta\eta\eta^{\prime}) = (4.86\pm0.62~({\rm stat.})\pm0.45~({\rm sys.}))\times10^{-5}$.


I. INTRODUCTION
The non-Abelian property of quantum chromodynamics (QCD) permits the existence of glueballs formed by gluons, the gauge bosons of the strong force [1][2][3]. The search for glueballs is an important field of research in hadron physics. However, the identification of glueballs is difficult in both experiment and theory due to the possible mixing of the pure glueball states with nearby qq nonet mesons. Lattice QCD (in the quenched approximation) predicts the lowest-lying glueballs are scalar (mass 1.5−1.7 GeV/c 2 ), tensor (mass 2.3−2.4 GeV/c 2 ), and pseudoscalar (mass 2.3−2.6 GeV/c 2 ) [4][5][6][7][8].
The radiative decay J/ψ → γgg is a gluon-rich process and is therefore regarded as one of the most promising hunting grounds for glueballs [9,10]. A possible pseudoscalar glueball candidate, the X(2370), is observed in the π + π − η invariant mass distribution through the decays of J/ψ → γπ + π − η [11] and in the KKη invariant mass distribution in the decays of J/ψ → γKKη [12] with statistical significances of 6.4σ and 8.3σ, respectively. The measured mass is consistent with the LQCD prediction for the pseudoscalar glueball [6]. In a calculation using an effective Lagrangian that couples the pseudoscalar glueball to scalar and pseudoscalar mesons, the ratios of the branching fractions of the pseudoscalar glueball decays Γ G→ηηη /Γ tot G , Γ G→KKη /Γ tot G and Γ G→ππη /Γ tot G are predicted to be 0.00082, 0.011 and 0.090 [13], respectively, for an assumed glueball mass of 2.370 GeV/c 2 . An observation of the X(2370) in J/ψ → γηηη would contribute to our understanding of this state. In parallel, we search for the η c since this charmonium state has never been observed decaying to ηηη [14].

II. DETECTOR AND MONTE CARLO SIMULATIONS
The BESIII detector is a magnetic spectrometer [16] located at the Beijing Electron Positron Collider (BEPCII) [17]. 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 (0.9 T in 2012) 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 charged-particle 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.
Simulated samples produced with the geant4based [18] Monte Carlo (MC) package which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate the backgrounds. The simulation includes the beam energy spread and initial state radiation (ISR) in the e + e − annihilations modeled with the generator kkmc [19]. The inclusive MC sample includes production of the J/ψ resonance as well as continuum processes incorporated with kkmc [19]. The known decay modes are modeled with evtgen [20] using branching fractions taken from the Particle Data Group (PDG) [14], and the remaining unknown decays from the charmonium states with lundcharm [21]. Final state radiation (FSR) from charged final state particles is incorporated with the photos package [22]. To estimate the selection efficiency and to optimize the selection criteria, signal MC events are generated for J/ψ → γX(2370), γη c → γηηη . The polar angle of the photon in the J/ψ center of mass system, θ γ , follows a 1 + cos 2 θ γ distribution. The decay of X(2370)/η c → ηηη is simulated using phase-space (PHSP) generator. So does the process η → ηπ + π − . To obtain the efficiency curves, MC events are generated for J/ψ → γX, X → ηηη , where X means 0 −+ non-resonant state. For the process η → γπ + π − , a generator taking into account both the ρ − ω interference and the box anomaly is used [23]. The analysis is performed in the framework of the BESIII offline software system (BOSS) [24] incorporating the detector calibration, event reconstruction and data storage.

III. EVENT SELECTION
Charged tracks in the polar angle range |cosθ| < 0.93 are reconstructed from hits in the MDC. Tracks must extrapolate to within 10 cm of the interaction point in the beam direction and 1 cm in the plane perpendicular to the beam. Each track is assumed to be a pion and no particle identification is applied. Candidate events are required to have two charged tracks and zero net charge.
Photon candidates are required to have an energy deposition above 25 MeV in the barrel region (| cos θ| < 0.80) and 50 MeV in the end cap (0.86 < | cos θ| < 0.92). To exclude showers from charged tracks, the angle between the shower position and the charged tracks extrapolated to the EMC must be greater than 10 • . A timing requirement in the EMC is used to suppress electronic noise and energy deposits unrelated to the event. At least six (seven) photons are required for the η → γπ + π − (η → π + π − η) mode.
For the J/ψ → γηηη , η → γπ + π − channel, a sixconstraint (6C) kinematic fit is performed to the hypothesis of J/ψ → γγηηπ + π − . This includes a 4C fit to the J/ψ initial four-momentum and 1C fit of each pair of photons to have an invariant mass equal to that of an η. For events with more than six photon candidates, the combination with the minimum χ 2 6C is selected, and χ 2 6C < 30 is required. Events with |M γγ −m π 0 | < 0.02 GeV/c 2 are rejected to suppress background containing a π 0 , where the m π 0 is the nominal mass of the π 0 [14]. In order to reduce the background due to misreconstruction of the event, events with |Mγγ − m η | < 0.02 GeV/c 2 are rejected, where the Mγγ is the invariant mass of all photon pairs except the pairs from the constrained η candidates and m η is the nominal mass of η [14]. A clear η signal is observed in the invariant mass distribution of γπ + π − (M γπ + π − ), as shown in Fig. 1(a). The π + π − invariant mass is required to be near the ρ mass region, M π + π − > 0.5 GeV/c 2 . Candidate η is reconstructed from the γπ + π − pair with |M γπ + π − −m η | < 0.015 GeV/c 2 , where the m η is the nominal mass of the η [14]. If there is more than one combination, we select the one with M γπ + π − closest to m η .
After applying the requirements above, we obtain the invariant mass distribution of ηηη (M ηηη ), in which clear η c signal is observed, as shown in Fig. 1 For the J/ψ → γηηη , η → π + π − η channel, a sevenconstraint (7C) kinematic fit is performed to the hypothesis of J/ψ → γηηηπ + π − in order to improve the η mass resolution. If there are more than seven photon candidates, the combination with the minimum χ 2 7C is retained, and χ 2 7C < 50 is required. To suppress background from π 0 → γγ, |M γγ − m π 0 | > 0.02 GeV/c 2 is required for all photon pairs. In order to reduce the background due to wrong reconstruction of the event, events with |M γrγη −m η | < 0.02 GeV/c 2 are rejected, where the M γrγη is the invariant mass of the radiative photon (γ r ) directly from J/ψ decays paired with any photon from an η candidate decay (γ η ). The η candidates are formed from π + π − η combination satisfying |M π + π − η − m η | < 0.015 GeV/c 2 and the combination with M π + π − η closest to m η is selected, where M π + π − η is the invariant mass of π + π − η, as shown in Fig. 1(c). Finally, the invariant mass distribution of ηηη (M ηηη ), with a clear signal of η c , is shown in Fig. 1(d).
FIG. 1. Invariant mass distributions for the selected candidates of J/ψ → γηηη . Plots (a) and (b) are the invariant mass distributions of γπ + π − and ηηη for η → γπ + π − , respectively; (c) and (d) are the invariant mass distributions of π + π − η and ηηη for η → π + π − η, respectively. The dots with error bars are data and the histograms are for the signal MC samples (arbitrary normalization).

IV. SIGNAL EXTRACTION
Potential backgrounds are studied using an inclusive MC sample of 1.2 × 10 9 J/ψ decays. No significant peaking background is observed in the invariant mass distribution of ηηη . Non-η processes are studied using the η mass sidebands, which are [0.890, 0.920] and [0.995, 1.025] GeV/c 2 . No clear peak is observed in X(2370) and η c mass region from sideband study.
There is no obvious signal for the X(2370) in ηηη invariant mass distributions in Figs. 1(b) Plots (a) and (b) are efficiency curves for the decays of η → γπ + π − and η → π + π − η obtained from J/ψ → γX → γηηη MC simulation, where X means 0 −+ non-resonant state. Plots (c) and (d) are the simultaneous fit results for the X(2370) in the invariant mass distribution of ηηη for the decays of η → γπ + π − and η → π + π − η, respectively. Plots (e) and (f) are the fit results for ηc in the invariant mass distribution of ηηη for the decays of η → γπ + π − and η → π + π − η, respectively. The dots with error bars represent the data; the red solid curves show the fit results; the hatched areas represent the signal of the X(2370) scaled to the upper limit or the signal of the ηc; the brown dashed lines show the events from η sideband; the green hyphenated lines represent the Chebychev polynomial function or the ARGUS function.

Breit-Wigner (BW) function convolved with a double
Gaussian function to account for the mass resolution. Due to low statistics, the mass and width of the BW function are fixed to previously published BESIII results [11] while the parameters of the double Gaussian function are fixed to the results obtained from the fit of signal MC samples generated with zero width of the X(2370). Interference between the X(2370) and other components is ignored. The non-η background events are described using η mass sidebands and the yields are fixed in the fit; the remaining background is described by a second order Chebychev polynomial function with free parameters. In the simultaneous fit, the signal ratio for the two η decay modes are fixed with a factor calculated by their branching fractions and efficiencies.
To obtain an upper limit of the signal yield, the distribution of normalized likelihood values for a series of expected signal event yields is taken as the probability density function (PDF). The 90% confidence level (C.L.) yield, N U L , is set such that 90% of the PDF area above zero yield is contained between 0 and N U L . We repeat this procedure with different X(2370) signal shape parameters, fit ranges, η sideband regions and background shapes, and the maximum upper limit among these cases is selected. The obtained upper limits of the signal yields are listed in Table I. The MC detection efficiencies of J/ψ → γX(2370) → γηηη for the two η decay modes are determined to be 2.95% (η → γπ + π − ) and 2.32% (η → π + π − η). The upper limit of the product branching fraction is B (J/ψ → γX(2370) · B(X(2370) → ηηη ) < 8.70 × 10 −6 .  E0Eγ +(E0−Eγ ) 2 to damp the divergent tail at low mass arising from the E 3 is the nominal energy of the transition photon [25]. The mass and width of the η c are fixed to PDG values [14]. Interference between η c and other components is ignored. Backgrounds are modeled with similar components as for the fit of the X(2370) discussed above, while the Chebychev polynomial is replaced with an ARGUS function [26]. The obtained signal yields, which have correlated uncertainties due to the constrianed fit, for J/ψ → γη c → γηηη are listed in Table I. The detection efficiencies of J/ψ → γη c → γηηη for two η decay modes are determined to be 2.94% (η → γπ + π − ) and 2.35% (η → π + π − η). We observe some disagreements in the data vs. MC simulated ηη, ηη , and ηηη invariant mass spectra. We employ a machine learning (ML) method [27] to re-weight the signal MC events based on the meson candidate's four-momenta. This reduces the inconsistency between data and sig-nal MC, providing an accurate efficiency. The product branching fraction of J/ψ → γη c → γηηη is then determined to be (4.86 ± 0.62(stat.)) × 10 −5 . The statistical significance of η c is determined to be 8.1σ.

V. SYSTEMATIC UNCERTAINTIES
Several sources of systematic uncertainties are considered, including the data-MC efficiency differences in the MDC tracking and the photon detection efficiency, the kinematic fit, and the mass window requirements for the π 0 , η, ρ and η . Uncertainties associated with the fit ranges, the background shapes, the sideband regions, quantum number of X(2370), the signal shape parameters of η c , damping factor, efficiency calculation, intermediate resonance decay branching fractions and the total number of J/ψ events are considered. TABLE II. Systematic uncertainties for determination of the upper limit of branching fraction of J/ψ → γX(2370) → γηηη (in %). The items with * are common uncertainties of both η decay modes.

A. Efficiency estimation
The MDC tracking efficiencies of charged pions are investigated using a clean control sample of J/ψ → ppπ + π − [28]. The difference in tracking efficiencies between data and MC simulation is 1.0% for each charged pion. The photon detection efficiency is studied with a clean sample of J/ψ → ρ 0 π 0 [30]. The result shows that the data-MC efficiency difference is 1.0% per photon.
The systematic uncertainties associated with the kinematic fit are studied with the track helix parameter correction method, as described in Ref. [31]. The differences with respect to those without corrections are taken as the systematic uncertainties.
Due to the difference in the mass resolution between data and MC simulation, the uncertainties related to the M π + π − and η mass window requirements are investigated by smearing the MC simulation to improve the consistency between data and MC simulation. The differences of the detection efficiency before and after smearing are assigned as systematic uncertainties for the M π + π − and η mass window requirements. The uncertainties from the π 0 and η mass window requirements are estimated by varying those mass windows. The changes in the resultant branching fractions are assigned as the systematic uncertainties from these items.
To study uncertainties related to the efficiency calculation with the ML method, we generate a generic MC sample with J/ψ → γη c , η c → f 2 (1810)η (f 2 (1810) → ηη) process to represent the signal and J/ψ → γηηη as the non-η c background. The numbers of signal and background events are fixed to fitting results. The efficiency difference between the generic MC sample and the ML method are taken as systematic uncertainty from this item. Furthermore, we consider the effects arising from different quantum numbers of the X(2370). We generate J/ψ → γX(2370) decays under the assumption of a sin 2 θ γ angular distribution. The resulting difference of efficiency with respect to the nominal value is taken as systematic uncertainty.

B. Fit to the signal
Systematic uncertainties related to the X(2370) signal treatment are already accounted for in the upper limit yield, as discussed previously; here, we discuss the treatment of the η c signal. To study the uncertainties from the fit range, the fits are repeated with different fit ranges, and the largest difference among these signal yields is taken as systematic uncertainty. The uncertainties from the η sideband region are estimated by using alternative sideband regions. The maximum difference among signal yields with respect to the nominal value is taken as the uncertainty. To estimate the uncertainty associated with the background shape, alternative fit with a truncated second order polynomial for the background is performed. The maximum difference in signal yields with respect to the nominal value is taken as systematic uncertainty. To study the uncertainty associated with the parameters of η c , we change these values by ±1σ and repeat the fit. The largest difference from our nominal result among these alternative fits is taken as the uncertainty. The uncertainty due to damping factor is estimated by using an alternative form of the damping factor, which was used by the CLEO collaboration [32], where E γ is the energy of the transition photon and β = 0.065 GeV. The difference between the results with different damping factor forms is taken as the systematic uncertainty.
A summary of all the uncertainties are shown in Tables II and III. The total systematic uncertainties are obtained by adding all individual uncertainties in quadrature, assuming all sources to be independent.
In this paper, J/ψ → γηηη is studied with two η decay modes. The measurements from the two η decay modes are, therefore, combined by considering the difference of uncertainties for these two measurements. The combination of common and independent systematic uncertainties for the two η decay modes are calculated with weighted least squares method [33]. The total systematic uncertainties are 12.8% and 9.2% for B(J/ψ → γX(2370)) · B(X(2370) → ηηη ) and B(J/ψ → γη c ) · B(η c → ηηη ), respectively.
No evident signal for the X(2370) is observed in the ηηη invariant mass distribution. The final upper limit of the product branching fraction of J/ψ → γX(2370) → ηηη incorporates the 12.8% relative systematic uncer-tainty by convolving the likelihood distribution with a Gaussian function: (1) where L(N ) is the likelihood distribution, σ sys = 0.128N , and N is the input signal yield. The resulting upper limit of B(J/ψ → γX(2370) → γηηη ) is estimated to be 9.2 × 10 −6 , which is not in contradiction with the value predicted in Ref. [13] where X(2370) is assumed as a pseudoscalar glueball. To understand the nature of X(2370), it is mandatory to measure its spin and parity and to search for it in more decay modes with higher statistics.