Observation of the decays χ cJ → φφη

M. Ablikim, M. N. Achasov, P. Adlarson, S. Ahmed, M. Albrecht, M. Alekseev, A. Amoroso, F. F. An, Q. An, Y. Bai, O. Bakina, R. Baldini Ferroli, I. Balossino, Y. Ban, K. Begzsuren, J. V. Bennett, N. Berger, M. Bertani, D. Bettoni, F. Bianchi, J. Biernat, J. Bloms, I. Boyko, R. A. Briere, H. Cai, X. Cai, A. Calcaterra, G. F. Cao, N. Cao, S. A. Cetin, J. Chai, J. F. Chang, W. L. Chang, G. Chelkov, D. Y. Chen, G. Chen, H. S. Chen, J. Chen, M. L. Chen, S. J. Chen, X. R. Chen, Y. B. Chen, W. Cheng, G. Cibinetto, F. Cossio, X. F. Cui, H. L. Dai, J. P. Dai, X. C. Dai, A. Dbeyssi, D. Dedovich, Z. Y. Deng, A. Denig, I. Denysenko, M. Destefanis, F. De Mori, Y. Ding, C. Dong, J. Dong, L. Y. Dong, M. Y. Dong, Z. L. Dou, S. X. Du, J. Z. Fan, J. Fang, S. S. Fang, Y. Fang, R. Farinelli, L. Fava, F. Feldbauer, G. Felici, C. Q. Feng, M. Fritsch, C. D. Fu, Y. Fu, Q. Gao, X. L. Gao, Y. Gao, Y. Gao, Y. G. Gao, B. Garillon, I. Garzia, E. M. Gersabeck, A. Gilman, K. Goetzen, L. Gong, W. X. Gong, W. Gradl, M. Greco, L. M. Gu, M. H. Gu, S. Gu, Y. T. Gu, A. Q. Guo, L. B. Guo, R. P. Guo, Y. P. Guo, A. Guskov, S. Han, X. Q. Hao, F. A. Harris, K. L. He, F. H. Heinsius, T. Held, Y. K. Heng, M. Himmelreich, Y. R. Hou, Z. L. Hou, H. M. Hu, J. F. Hu, T. Hu, Y. Hu, G. S. Huang, J. S. Huang, X. T. Huang, X. Z. Huang, N. Huesken, T. Hussain, W. Ikegami Andersson, W. Imoehl, M. Irshad, Q. Ji, Q. P. Ji, X. B. Ji, X. L. Ji, H. L. Jiang, X. S. Jiang, X. Y. Jiang, J. B. Jiao, Z. Jiao, D. P. Jin, S. Jin, Y. Jin, T. Johansson, N. Kalantar-Nayestanaki, X. S. Kang, R. Kappert, M. Kavatsyuk, B. C. Ke, I. K. Keshk, A. Khoukaz, P. Kiese, R. Kiuchi, R. Kliemt, L. Koch, O. B. Kolcu, B. Kopf, M. Kuemmel, M. Kuessner, A. Kupsc, M. Kurth, M. G. Kurth, W. Kühn, J. S. Lange, P. Larin, L. Lavezzi, H. Leithoff, T. Lenz, C. Li, C. H. Li, Cheng Li, D. M. Li, F. Li, G. Li, H. B. Li, H. J. Li , J. C. Li, Ke Li, L. K. Li, Lei Li, P. L. Li, P. R. Li, W. D. Li, W. G. Li, X. H. Li, X. L. Li, X. N. Li, Z. B. Li, Z. Y. Li, H. Liang, H. Liang, Y. F. Liang, Y. T. Liang, G. R. Liao, L. Z. Liao, J. Libby, C. X. Lin, D. X. Lin, Y. J. Lin, B. Liu, B. J. Liu, C. X. Liu, D. Liu, D. Y. Liu, F. H. Liu, Fang Liu, Feng Liu, H. B. Liu, H. M. Liu, Huanhuan Liu, Huihui Liu, J. B. Liu, J. Y. Liu, K. Liu, K. Y. Liu, Ke Liu, L. Y. Liu, Q. Liu, S. B. Liu, T. Liu, X. Liu, X. Y. Liu, Y. B. Liu, Z. A. Liu, Zhiqing Liu, Y. F. Long, X. C. Lou, H. J. Lu, J. D. Lu, J. G. Lu, Y. Lu, Y. P. Lu, C. L. Luo, M. X. Luo, P. W. Luo, T. Luo, X. L. Luo, S. Lusso, X. R. Lyu, F. C. Ma, H. L. Ma, L. L. Ma, M.M. Ma, Q. M. Ma, X. N. Ma, X. X. Ma, X. Y. Ma, Y. M. Ma, F. E. Maas, M. Maggiora, S. Maldaner, S. Malde, Q. A. Malik, A. Mangoni, Y. J. Mao, Z. P. Mao, S. Marcello, Z. X. Meng, J. G. Messchendorp, G. Mezzadri, J. Min, T. J. Min, R. E. Mitchell, X. H. Mo, Y. J. Mo, C. Morales Morales, N. Yu. Muchnoi, H. Muramatsu, A. Mustafa, S. Nakhoul, Y. Nefedov, F. Nerling, I. B. Nikolaev, Z. Ning, S. Nisar, S. L. Niu, S. L. Olsen, Q. Ouyang, S. Pacetti, Y. Pan, M. Papenbrock, P. Patteri, M. Pelizaeus, H. P. Peng, K. Peters, J. Pettersson, J. L. Ping, R. G. Ping, A. Pitka, R. Poling, V. Prasad, M. Qi, S. Qian, C. F. Qiao, X. P. Qin, X. S. Qin, Z. H. Qin, J. F. Qiu, S. Q. Qu, K. H. Rashid, K. Ravindran, C. F. Redmer, M. Richter, A. Rivetti, V. Rodin, M. Rolo, G. Rong, Ch. Rosner, M. Rump, A. Sarantsev, M. Savrié, Y. Schelhaas, K. Schoenning, W. Shan, X. Y. Shan, M. Shao, C. P. Shen, P. X. Shen, X. Y. Shen, H. Y. Sheng, X. Shi, X. D. Shi, J. J. Song, Q. Q. Song, X. Y. Song, S. Sosio, C. Sowa, S. Spataro, F. F. Sui, G. X. Sun, J. F. Sun, L. Sun, S. S. Sun, X. H. Sun, Y. J. Sun, Y. K. Sun, Y. Z. Sun, Z. J. Sun, Z. T. Sun, Y. T. Tan, C. J. Tang, G. Y. Tang, X. Tang, V. Thoren, B. Tsednee, I. Uman, B. Wang, B. L. Wang, C.W. Wang, D. Y. Wang, K. Wang, L. L. Wang, L. S. Wang, M. Wang, M. Z. Wang, Meng Wang, P. L. Wang, R. M. Wang, W. P. Wang, X. Wang, X. F. Wang, X. L. Wang, Y. Wang, Y. Wang, Y. F. Wang, Y. Q. Wang, Z. Wang, Z. G. Wang, Z. Y. Wang, Z. Y. Wang, Zongyuan Wang, T. Weber, D. H. Wei, P. Weidenkaff, F. Weidner, H.W. Wen, S. P. Wen, U. Wiedner, G. Wilkinson, M. Wolke, L. H. Wu, L. J. Wu, Z. Wu, L. Xia, Y. Xia, S. Y. Xiao, Y. J. Xiao, Z. J. Xiao, Y. G. Xie, Y. H. Xie, T. Y. Xing, X. A. Xiong, Q. L. Xiu, G. F. Xu, J. J. Xu, L. Xu, Q. J. Xu, W. Xu, X. P. Xu, F. Yan, L. Yan, W. B. Yan, W. C. Yan, Y. H. Yan, H. J. Yang, H. X. Yang, L. Yang, R. X. Yang, S. L. Yang, Y. H. Yang, Y. X. Yang, Yifan Yang, Z. Q. Yang, Zhi Yang, M. Ye, M. H. Ye, J. H. Yin, Z. Y. You, B. X. Yu, C. X. Yu, J. S. Yu, T. Yu, C. Z. Yuan, X. Q. Yuan, Y. Yuan, C. X. Yue, A. Yuncu, A. A. Zafar, Y. Zeng, B. X. Zhang, B. Y. Zhang, C. C. Zhang, D. H. Zhang, H. H. Zhang, H. Y. Zhang, J. Zhang, J. L. Zhang, J. Q. Zhang, J. W. Zhang, J. Y. Zhang, J. Z. Zhang, K. Zhang, L. Zhang, Lei Zhang, S. F. Zhang, T. J. Zhang, X. Y. Zhang, Y. Zhang, Y. H. Zhang, Y. T. Zhang, Yang Zhang, Yao Zhang, Yi Zhang, Yu Zhang, Z. H. Zhang, Z. P. Zhang, Z. Y. Zhang, G. Zhao, J. Zhao, J. W. Zhao, J. Y. Zhao, J. Z. Zhao, PHYSICAL REVIEW D 101, 012012 (2020)


I. INTRODUCTION
Studies of the properties of cc states play an important role in understanding the interplay between perturbative and nonperturbative 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 → VVP, 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 decaying to ϕϕ in B decays [6], but none have so far been observed [7]. 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. The BESIII detector has collected ð448.1 AE 2.9Þ × 10 6 ψð3686Þ decays [8], which is the world's largest data sample of ψð3686Þ decays produced in e þ e − annihilation. 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 [9] located at the Beijing Electron Positron Collider (BEPCII) [10]. 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 a 4π solid angle. The momentum resolution for charged particles 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 [11] 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 [12]. 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 [12]. The known decay modes are modeled with EvtGen [13] using branching fractions taken from the Particle Data Group [1], and the remaining unknown decays of the charmonium states are modeled with LUNDCHARM [14]. The final state radiation (FSR) from charged final state particles is simulated with the PHOTOS package [15].

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 j cos θj < 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 higher 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 (j cos θj < 0.80) or 50 MeV in the end cap region (0.86 < j cos θj < 0.92). To exclude showers from charged particles, a photon must be separated by at least 10°f rom the nearest charged track. The measured EMC time is required to be within 0 and 700 ns of the collision 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 = ffiffiffiffiffiffiffi ffi N tot p , 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-momentum 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 its 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 line shape 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. In each event, one of the ϕ candidates is randomly chosen to be ϕ 1 with the other ϕ 2 . 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 its 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 combinations 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 [1], and c and d are free parameters. The two-dimensional (2D) ϕ 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 randomly assigned K þ K − ð1Þ and K þ K − ð2Þ, respectively.
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 ψð3686Þ → ϕϕϕ, where one ϕ decays to γη.
A total of 495 candidate events survive in the R1 region, which corresponds to the area A with M γγ in the signal region, as shown in Fig. 3(a). The distributions

IV. BACKGROUND AND SIGNAL YIELDS
According to a study of the inclusive MC sample, consisting of 5.06 × 10 8 ψð3686Þ decays, the background sources can be categorized into two classes. Class I background is from events that do not form an η signal in the M γγ distribution. This background is estimated using events in the η sideband regions of M γγ ∈ ½0.47; 0.50 ∪ ½0.60; 0.63 GeV=c 2 . Class II background arises from decays having only one ϕ present, which is described using events in the 2D ϕ sideband regions (the areas B of Fig. 2), where one M K þ K − lies in the ϕ signal region and the other is in the ϕ sideband region of M K þ K − ∈ ½1.045; 1.075 GeV=c 2 . Since there are no events observed in the area C of Fig. 2, in which both M K þ K − lie in the ϕ sideband region, we ignore this contribution.
The possible quantum electrodynamics (QED) processes under the ψð3686Þ peak consist of the peaking and nonpeaking parts. According to the Born cross sections measured by the BESIII [20] and Belle experiments [21] and the nominal solution below, the upper limits on the numbers of events at 90% confidence level are estimated to be less than 0.1 for the QED peaking e þ e − → γχ cJ processes, and thus their contributions can be negligible. The QED nonpeaking processes that do not have the χ cJ signals, like e þ e − → γϕϕη, can contribute to the flat distributions in Fig. 3. The QED processes are also studied using a 48.8 pb −1 off-resonance data sample taken at a center-of-mass energy of 3.65 GeV [22], where the scaling factor of the integrated luminosity is 0.07 compared with the ψð3686Þ data sample [8]. But no events survive after applying the same event selection criteria.
The signal yields are obtained from unbinned maximum likelihood fits to the M ϕK þ K − γγ spectra, where at least one of the ϕ candidates has an invariant mass within the signal window. The fits are performed in the R1, R2, and R3 regions, where the latter two regions correspond to the area A with M γγ in the sideband regions, and the areas B with M γγ in the signal region, respectively. In the fits, the signal shape in the M ϕK þ K − γγ distribution is extracted from signal MC simulations, and the background shape is modeled as a constant, where the background is from the events that do not form the χ cJ signals, such as other ψð3686Þ decays and the possible QED nonpeaking processes. Figure 3 shows the fit results. The contribution of the areas B with M γγ in the sideband region is negligible, since there are only two 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, which have  been evaluated from the ratios of the background yields in the η and 2D ϕ signal and sideband regions, respectively. The numbers 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 for χ cJ → ϕϕη decays are determined by where N ψð3686Þ is the total number of ψð3686Þ decays, 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 [23] are modified to reduce the difference of the kinematic fit χ 2

4C
between the data and MC sample, where the correction factors for their resolutions are obtained from a clean control sample of the ψð3686Þ → K þ K − π þ π − decay. (ii) Taking into account E1 transition effects on the line shapes of χ cJ mesons generated with the Breit-Wigner functions, a weighting factor ð E γ1 E γ10 Þ 3 [24] 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. The detection efficiencies are determined to be ð5.30 AE 0.02Þ%, ð6.77 AE 0.03Þ%, and ð6.62 AE 0.03Þ% for χ c0 , χ c1 , and χ c2 → ϕϕη decays, respectively.

VI. SYSTEMATIC UNCERTAINTIES
The sources of systematic uncertainty include the total number of ψð3686Þ decays, 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 AE π ∓ ; K 0 S → π þ π − decays [25] have been used to investigate the MDC tracking efficiency, and the difference of 1% per K AE 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 AE 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 → γγ [26].
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 shape derived from 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 shape from 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 modified to improve the agreement between the data and MC simulation. An alternative detection efficiency is obtained with no modification to 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 with free width, and the difference of the signal yield is taken as the systematic uncertainty from the χ cJ signal shape.  The largest difference of the signal yield is taken as the systematic uncertainty associated with the fit range. 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 2D ϕ 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 width of the η (ϕ) sideband regions within one standard deviation of the corresponding 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 extracted from the world average values [1]. The systematic uncertainty due to the trigger efficiency is negligible according to the studies in Ref. [27].
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 existence of intermediate states.