Observation of the decays $\chi_{cJ} \to \phi \phi \eta$

Using a data sample of $(448.1\pm2.9)\times10^{6}$ $\psi(3686)$ decays collected by the BESIII detector at the Beijing Electron Positron Collider (BEPCII), we observe the decays $\chi_{cJ}\to \phi\phi\eta~(J=0,~1,~2)$, where the $\chi_{cJ}$ are produced via the radiative processes $\psi(3686)\to\gamma\chi_{cJ}$. The branching fractions are measured to be $\mathcal B(\chi_{c0}\to\phi\phi\eta)=(8.41\pm0.74\pm0.62)\times10^{-4}$, $\mathcal B(\chi_{c1}\to\phi\phi\eta)=(2.96\pm0.43\pm0.22)\times 10^{-4}$, and $\mathcal B(\chi_{c2} \to \phi\phi\eta)=(5.33\pm0.52\pm0.39) \times 10^{-4}$, where the first uncertainties are statistical and the second are systematic. We also search for intermediate states in the $\phi\phi$ or $\eta\phi$ combinations, but no significant structure is seen due to the limited statistics.


I. INTRODUCTION
Studies of the properties of cc states play an important role in understanding the interplay between perturbative and non-perturbative effects in quantum chromodynamics (QCD). Besides J/ψ and ψ(3686) decays [1], the decays of the χ cJ (J = 0, 1, 2) [2,3] are also valuable to probe a wide variety of QCD phenomena.
To date, only a few measurements have been performed for decays of the form χ cJ → V V P , where V and P denote vector and pseudoscalar mesons, respectively [1], and no measurement of the branching fraction for χ cJ → φφη has previously been reported. The interest in these final states arises from the search for glueballs in the φφ invariant mass (M φφ ) spectrum. A previous partial wave analysis of the decay J/ψ → γφφ decay by the BESIII Collaboration [4] confirmed the existence of the η(2225) and observed the three tensor states f 2 (2010), f 2 (2300) and f 2 (2340), which were first observed in the process π − p → φφn [5]. Different experiments also searched for glueballs [6] decaying to φφ in B decays, but none have so far been observed [1]. Although there are no theoretical expectations, the decays χ cJ → φφη may contain contributions from intermediate states decaying to φφ and ηφ, and observations of the same resonances as those in J/ψ decays would provide supplementary and conclusive information regarding their existence.
Due to abundant χ cJ production in ψ(3686) radiative decays [1], the BESIII experiment provides an ideal place to search for new χ cJ decays based on the world's largest e + e − annihilation data sample of (448.1 ± 2.9) × 10 6 ψ(3686) events [7]. In this paper, we report the first measurements of the branching fractions of χ cJ decays to φφη. The φ meson can be reconstructed with φ → K + K − , φ → π + π − π 0 and φ → K 0 S K 0 L decays, and the η meson with η → γγ and η → π + π − π 0 decays. Compared to the φ → K + K − and η → γγ modes, other decay modes suffer from higher backgrounds and lower detection efficiencies. So in this analysis, the two φ mesons and the η meson are reconstructed with φ → K + K − and η → γγ processes.

II. DETECTOR AND MONTE CARLO SIMULATIONS
The BESIII detector is a magnetic spectrometer [8] located at the Beijing Electron Positron Collider (BEPCII) [9]. 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 for charged particles and photons is 93% over 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.
Large samples of simulated events are produced with a geant4-based [10] Monte Carlo (MC) package that includes the geometric description of the BESIII detector and the detector response. These samples 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 − annihilation modeled with the generator kkmc [11]. The 'inclusive' MC sample consists of the production of the ψ(3686) resonance, the ISR production of the J/ψ, and the continuum processes incorporated in kkmc [11]. The known decay modes are modeled with evtgen [12] using branching fractions taken from the Particle Data Group [1], and the remaining unknown decays of the charmonium states are modeled with lundcharm [13]. The final state radiation (FSR) from charged final state particles is simulated with the photos package [14].

III. EVENT SELECTION
The cascade decay of interest is ψ(3686) → γχ cJ , χ cJ → φφη, with φ → K + K − and η → γγ. Candidate events are required to have four charged tracks with zero net charge and at least three photons. Charged tracks in an event are required to have a polar angle θ with respect to the beam direction within the MDC acceptance | cos θ| < 0.93, and a distance of closest approach to the interaction point within 10 cm along the beam direction and 1 cm in the plane transverse to the beam direction. The TOF and dE/dx information are combined to evaluate particle identification (PID) confidence levels for the π and K hypotheses, and the particle type with the highest confidence level is assigned to each track. All charged tracks must be identified as kaons. Electromagnetic showers are reconstructed from clusters of energy deposited in the EMC. The energy deposited in nearby TOF counters is included to improve the reconstruction efficiency and energy resolution. Photon candidates must have a minimum energy of 25 MeV in the barrel region (| cos θ| < 0.80) or 50 MeV in the end cap region (0.86 < | cos θ| < 0.92). To exclude showers from charged particles, a photon must be separated by at least 10 • from the nearest charged track. The measured EMC time is required to be within 0 and 700 ns of the event start time to suppress electronic noise and energy deposits unrelated to the event of interest.
A four-constraint (4C) kinematic fit imposing overall energy-momentum conservation is performed with the γγγK + K − K + K − hypothesis, and the events with χ 2 4C < 40 are retained. The requirement is based on the optimization of the figure of merit (FOM), FOM ≡ N sig / √ N tot , where N sig and N tot are the number of signal events and total number of events estimated from the signal MC sample and data, respectively. For events with more than three photon candidates, the combination with the smallest χ 2 4C is retained. Further selection criteria are based on the four-momenta updated by the kinematic fit.
After the above requirements, the η candidate is reconstructed in its decay to γγ using the γγ pair with invariant mass M γγ closest to the nominal η mass [1]. The η signal region is defined as 0.52 ≤ M γγ ≤ 0.58 GeV/c 2 , with half-width approximately three times larger than the detector resolution (σ η = 10 MeV/c 2 ). Figure 1(a) shows a fit to the M γγ distribution. In the fit, the signal shape is modeled by the MC-simulated lineshape convolved with a Gaussian function with free width and the background is described by a linear function. The two signal φ candidates are chosen from the combination with the minimum value of [1], and i, j can be 0 or 1. MC studies show that the miscombination rates for both η and φ candidates are no more than 0.1%. The φ signal region is defined as 1.005 ≤ M K + K − ≤ 1.035 GeV/c 2 , with half-width about three times the sum of the detector resolution (σ φ = 1 MeV/c 2 ) and intrinsic width [1]. Figure 1(b) shows the fit to the M K + K − distribution obtained when one of the two φ candidates is randomly selected. In the fit, the signal shape is modeled as a P -wave Breit-Wigner convolved with a Gaussian function, and the background shape is represented by the function b( where m t is the K + K − mass threshold, and c and d are free parameters. The twodimensional (2-D) φ signal region is shown as the area "A" in Fig. 2, where M K + K − (1) and M K + K − (2) denote the invariant masses of the two φ candidates.
The mass recoiling against the η is required to be less than 3.05 GeV/c 2 to suppress background from the decay ψ(3686) → ηJ/ψ, J/ψ → γφφ. All combinations of M γγ are required to be outside the range [0.115, 0.150] GeV/c 2 to suppress background events with π 0 decays, and the invariant mass of γη must be outside the range [1.00, 1.04] GeV/c 2 to suppress background from the decay

IV. BACKGROUND AND SIGNAL YIELDS
According to a study of the inclusive MC sample, consisting of 5.06 × 10 8 ψ(3686) events, the background sources can be categorized into two classes. The class I background is from the decays with no η signal formed in the M γγ distribution, which is estimated by the events in   The quantum electrodynamics process under the ψ(3686) peak is studied based on the off-resonance sample of 48.8 pb −1 taken at the center-of-mass energy of 3.65 GeV [19]. With the same event selection criteria, no events survive, so this contribution is also negligible.
The signal yields are obtained from unbinned maximum likelihood fits to the M φK + K − γγ spectra, where at least one of the φ candidates has the invariant mass within the signal window. The fits are performed in the following three regions, "R1", "R2", and "R3", which correspond to the area "A" with the η in the signal region, the area "A" with the η in the sideband regions, and the areas "B" with the η in the signal region, respectively. In the fits, the signal shape is extracted from signal MC simulations, and the background shape is modeled as a constant. Figure 4 shows the fit results. The contribu- tion of the areas "B" with the η in the sideband region is negligible, since there are only 2 events. The signal yields for χ cJ → φφη decays are estimated by where N r obs is the number of observed events for the corresponding r region (r = R1, R2, or R3), and both the normalization factors f R2 and f R3 are 1.0 evaluated from the ratios of the background yields in the η and 2-D φ signal and sideband regions, respectively. The number of events obtained by the fits to M φK + K − γγ in different regions for χ cJ → φφη decays are summarized in Table I, along with their statistical significances.

V. BRANCHING FRACTIONS
The branching fractions of χ cJ → φφη decays are determined by where N ψ(3686) is the total number of ψ(3686) events, B 0 = B ψ(3686)→γχcJ B 2 φ→K + K − B η→γγ is the product of the branching fractions cited from the world average values [1], and ǫ is the detection efficiency.
In order to obtain the best possible estimate for the detection efficiencies, the signal MC samples are corrected in two aspects: (i) The track helix parameters [20] are corrected to reduce the difference of the kinematic fit χ 2 4C between the data and MC sample, where the correction factors are obtained from a clean control sample of ψ(3686) → K + K − π + π − decay.
(ii) Taking into account E1 transition effects on the lineshapes of χ cJ mesons generated with the Breit-Wigner functions, a weighting factor, ( Eγ1 Eγ10 ) 3 [21] is applied to the M φK + K − γγ spectra, where E γ1 is the radiative photon's energy in the rest frame of the ψ(3686) meson without detector reconstruction effects, and E γ10 is the most probable transition energy, Here m χcJ are the nominal masses of the χ cJ mesons [1] and E cms is the center-of-mass energy of 3.686 GeV.

VI. SYSTEMATIC UNCERTAINTIES
The sources of systematic uncertainty include the total number of ψ(3686) events, the MDC tracking efficiency, PID efficiency, photon detection efficiency, η and φ mass requirements, kinematic fit, fit procedure, peaking background estimation, and cited branching fractions.
The control samples of J/ψ → K 0 S K ± π ∓ , K 0 S → π + π − decays [22] has been used to investigate the MDC tracking efficiency, and the difference of 1% per K ± track between the data and MC simulation is assigned as the systematic uncertainty. By means of the same control sample, the uncertainty due to PID efficiency is estimated to be also 1% per K ± track. The systematic uncertainty from the photon detection efficiency is determined to be 1% per photon utilizing a control sample of J/ψ → ρ 0 π 0 with ρ 0 → π + π − and π 0 → γγ [23].
The systematic uncertainty arising from the η (φ) mass requirement is evaluated by changing the mass resolution and shifting the mass window. In the nominal fit, the η signal shape is described as the signal MC simulation convolved with a Gaussian function, and the φ signal shape is modeled as a P -wave Breit-Wigner function convolved with a Gaussian function. Alternative fits are performed by modeling the η signal shape with the signal MC simulation, changing the width of the Gaussian function for the φ signal shape to that obtained from the signal MC sample, and varying the η (φ) mass window by the respective mass resolution obtained in the fit to the signal shape. The difference of the efficiency of the η (φ) mass requirement between the data and MC sample is taken as the systematic uncertainty from the η (φ) mass requirement.
In the nominal analysis, the track helix parameters for charged tracks from signal MC samples are corrected to improve the agreement between the data and MC simulation. The alternative detection efficiency is obtained with no correction on the track helix parameters, and the difference is assigned as the systematic uncertainty associated with the kinematic fit.
The sources of systematic uncertainty from the fit procedure include the signal shape, background shape, and the fit range.
(i) In the nominal fit, the χ cJ signal shape is modeled with the MC simulation. An alternative fit is performed with the MC simulation convolved with a Gaussian function, and the difference of the signal yield is taken as the systematic uncertainty from the χ cJ signal shape.
(ii) Different order Chebychev functions instead of a constant are used in the alternative fits to describe the background. The largest difference of the signal yield is assigned as the systematic uncertainty from the background shape. The quadratic sum of the above three systematic uncertainties is taken as the systematic uncertainty from the fit procedure.
In the nominal fit, events in the η sideband and 2-D φ sideband regions are used to estimate contributions from peaking background sources with no η signal and only one φ signal, respectively. Alternative fits are performed by varying the η and φ sideband regions by the mass resolution, and the largest difference of the signal yield is taken as the corresponding systematic uncertainty. The quadratic sum of the two cases is taken as the systematic uncertainty from peaking backgrounds.
The uncertainties associated with the branching fractions of ψ(3686) → γχ cJ , φ → K + K − and η → γγ are estimated from the world average values [1]. The systematic uncertainty due to the trigger efficiency is negligible according to the studies in Ref. [24].
The total systematic uncertainty on the measured branching fractions for χ cJ → φφη decays is the quadratic sum of each individual contribution, as summarized in Table II.

VII. RESULTS AND DISCUSSION
The measured branching fractions of χ cJ → φφη decays are summarized in Table III, where the first uncertainties are statistical, and the second are systematic. Figure 5 shows the projections on the M φφ and M ηφ spectra. There are two combinations of M ηφ for each event. Compared with those from the signal MC samples, some excesses in data are observed. However, considering the limited statistics, it is hard to draw a conclusion that intermediate states appear in χ cJ → φφη decays. Perhaps in the future, utilizing more data samples, it would be worthwhile to combine other φ and η decay modes, such as φ → π + π − π 0 , φ → K 0 S K 0 L , and η → π + π − π 0 decays, to perform a partial wave analysis of χ cJ → φφη decays, so that we can make clear conclusions on the existences of intermediate states.