Search for the rare decay η0 → π0π0π0π0 at BESIII

M. Ablikim, M. N. Achasov, P. Adlarson, S. Ahmed, M. Albrecht, M. Alekseev, A. Amoroso, Q. An, Y. Bai Anita, 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,28,§,∥ 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. Fang, S. S. Fang, Y. Fang, R. Farinelli, L. Fava, F. Feldbauer, G. Felici, C. Q. Feng, M. Fritsch, C. D. Fu, Y. Fu, X. L. Gao, Y. Gao, Y. Gao, Y. G. Gao, 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, T. T. Han, X. Q. Hao, F. A. Harris, K. L. He, F. H. Heinsius, T. Held, Y. K. Heng, M. Himmelreich, T. Holtmann, 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, S. Jaeger, Q. Ji, Q. P. Ji, X. B. Ji, X. L. Ji, H. B. 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. 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. Liu, L. Y. Liu, Q. Liu, S. B. Liu, T. Liu, X. Liu, X. Y. Liu, Y. B. Liu, Z. A. Liu, Z. Q. 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, R. T. 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. 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, H. Qi, M. Qi, T. Y. 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,28,∥∥ M. Savrié, Y. Schelhaas, C. Schnier, 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, Y. X. Song, S. Sosio, C. Sowa, S. Spataro, F. F. Sui, G. X. Sun, J. F. Sun, L. Sun, S. S. Sun, Y. J. Sun, Y. K. Sun, Y. Z. Sun, Z. J. Sun, Z. T. Sun, Y. X. 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, W. P. Wang, X. Wang, X. F. Wang, X. L. Wang, Y. Wang, Y. Wang, Y. D. Wang, Y. F. Wang, Y. Q. Wang, Z. Wang, Z. G. Wang, Z. Y. Wang, Ziyi 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. Wollenberg, L. H. Wu, L. J. Wu, Z. Wu, L. Xia, S. Y. Xiao, Y. J. Xiao, Z. J. Xiao, Y. G. Xie, Y. H. Xie, T. Y. Xing, X. A. Xiong, G. F. Xu, J. J. Xu, Q. J. Xu, W. Xu, X. P. Xu, F. Yan, L. Yan, W. B. Yan, W. C. Yan, W. C. Yan, H. J. Yang, H. X. Yang, L. Yang, R. X. Yang, S. L. Yang, Y. H. Yang, Y. X. Yang, Yifan 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. L. Zhang, J. Q. Zhang, J. W. Zhang, J. Y. Zhang, J. Z. Zhang, L. Zhang, Lei Zhang, S. F. Zhang, T. J. Zhang,41,∥ X. Y. Zhang, Y. H. Zhang, Y. T. Zhang, Yan 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, Lei Zhao, Ling Zhao, PHYSICAL REVIEW D 101, 032001 (2020)


I. INTRODUCTION
The η 0 meson has a special role in improving the understanding of low-energy quantum chromodynamics (QCD), and studies of its decays have attracted considerable theoretical [1] and experimental attention [2,3]. In addition to its important role in testing the fundamental discrete symmetries and searching for processes beyond the Standard Model (SM), η 0 decays offer unique opportunities to test the chiral perturbation theory (ChPT) [4] and the vector-meson dominance (VMD) [5] model.
In theory, η 0 → 4π 0 is a highly suppressed decay because of the S-wave CP-violation. In the light of an effective chiral Lagrangian approach, the S-wave CP-violation in η 0 → 4π 0 is induced by the so-called θ-term, which is an additional term in the QCD Lagrangian to account for the solution of the strong-CP problem. It was found that the S-wave CP-violation effect that contributed to this decay is at a level of 10 −23 [6,7]. While higher-order contributions, involving a D-wave pion loop or the production of two f 2 tensor mesons (see Fig. 1), provide a CP-conserving route through which the decay can occur. By ignoring the tiny contribution from the latter process, calculations based on ChPTand VMD model predict the branching fraction caused by D-wave CP-conserving to be at the level of 10 −8 [8]. However, the theoretical prediction is not strictly based on the effective field theory due to the lack of knowledge at such a high order in the chiral expansion and the use of a model to make an estimation. One does not know the reliability of that model a priori. Therefore, a search for the decay η 0 → 4π 0 is useful to check the reliability of it.
So far the decay η 0 → 4π 0 has not been observed. About three decades ago the first attempt to search for this mode was performed by the GAMS Collaboration, and the upper limit on the branching fraction was determined to be Bðη 0 → 4π 0 Þ < 5 × 10 −4 at the 90% confidence level (C.L.) [9]. The more recent upper limit of 3.2 × 10 −4 at 90% C.L. was obtained by the GAMS-4π Collaboration [10].
Although η 0 meson can not be produced directly in e þ e − annihilations, the decay J=ψ → γη 0 , with a branching fraction of ð5.13 AE 0.17Þ × 10 −3 [11], provides an abundant source of η 0 meson in this environment. The BESIII experiment has exploited this production mode to perform a series of studies of η 0 decays [12], based on a sample of ð1310.6 AE 7.0Þ × 10 6 J=ψ events taken at the center-ofmass energy of 3.097 GeV [13] with the BESIII detector, corresponding to 6.7 × 10 6 η 0 events. In this paper, using this same J=ψ sample, we perform a search for η 0 → 4π 0 via J=ψ → γη 0 .

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector [14] is a magnetic spectrometer located at the Beijing Electron Positron Collider (BEPCII) [15], which is a double-ring e þ e − collider with a design peak luminosity of 10 33 cm −2 s −1 at the center-of-mass energy of 3.773 GeV. 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 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 a GEANT4-based [16] Monte Carlo (MC) simulation framework, 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 e þ e − annihilation modeled with the generator KKMC [17]. The inclusive MC sample consists of the production of the J=ψ resonance, and the continuum processes incorporated in KKMC [17]. The known decay modes are modeled with [18,19] using branching fractions taken from the Particle Data Group (PDG) [20], and the remaining unknown decays of the charmonium states with LUNDCHARM [21]. The final state radiation (FSR) from charged final-state particles are incorporated with the PHOTOS package [22].

III. EVENT SELECTION
In this analysis, the pseudoscalar mesons η 0 and π 0 are reconstructed in the modes η 0 → 4π 0 and π 0 → γγ. Candidate J=ψ → γ4π 0 decays are chosen by selecting ' η The photon candidate with the maximum energy deposited in the EMC is treated as the radiative photon directly originating from the J=ψ decay. For the two-body signal decay J=ψ → γη 0 , this photon carries a unique energy of 1.4 GeV. To reconstruct the π 0 candidate from the remaining photons, a one-constraint (1C) kinematic fit is performed to each photon pair with the invariant mass constrained to the π 0 mass, with the requirement that the goodness-of-fit χ 2 1C ðγγÞ < 25. And at least four π 0 candidates are required for further analysis. Then an eightconstraint (8C) kinematic fit is performed for the γπ 0 π 0 π 0 π 0 combination by enforcing energy-momentum conservation and constraining the invariant masses of each of the four photon pairs to the nominal π 0 mass. If more than one combination is found in an event, only the one with the smallest χ 2 8C is retained. The χ 2 8C distribution is shown in Fig. 2. Candidate events with χ 2 8C > 30 are rejected. The MC study shows that after the above selection, about 8% of events have a miscombination of photons, which mainly occurs between the different π 0 candidates. Since this is not the miscombination between the radiative photon and the photon from π 0 candidates, the invariant mass of 4π 0 is still in the η 0 mass region.

IV. FIT MODEL AND UPPER LIMIT ON BRANCHING FRACTION OF η 0 → 4π 0
An unbinned maximum likelihood fit is performed on the 4π 0 mass spectrum, allowing for background contributions and a possible signal component. In the fit, the line shape of the η 0 signal is determined by a MC simulation. The shape of peaking background J=ψ → γη 0 ; η 0 → π 0 π 0 η; η → π 0 π 0 π 0 is obtained from the dedicated MC simulation [23,24] with the number fixed to the normalized value while the nonpeaking background contribution is described by a third-order Chebychev polynomial function with the number and parameters of Chebychev function free.
A Bayesian approach is used to determine an upper limit on the branching fraction of η 0 → 4π 0 . A series of unbinned extended maximum likelihood fits are performed for different assumed values of the signal yield N, and for each fit the negative log-likelihood S is determined. For each value of N the branching fraction is where N J=ψ ¼ ð1310.6 AE 7.0Þ × 10 6 is the number of J=ψ events [13], ε is the detection efficiency, BðJ=ψ → γη 0 Þ and Bðπ 0 → γγÞ are the branching fractions of J=ψ → γη 0 and π 0 → γγ, respectively, which are taken from Ref. [11].
The distribution of normalized likelihood values, defined as LðBÞ ¼ expð−½SðBÞ − S min Þ, where S min is the lowest negative log-likelihood obtained from the ensemble of fits, is taken as the probability density function (PDF) for the expected branching fraction of η 0 → 4π 0 . The upper limit on the branching fraction at the 90% C.L., defined as B UL , corresponds to the branching fraction at 90% of the integral of the PDF, and is found to be 4.57 × 10 −5 , considering statistical uncertainties alone.

V. SYSTEMATIC UNCERTAINTIES
Two categories of systematic uncertainty are considered: those associated with the fit model and procedure, and those which enter when using Eq. (1) to express the signal yield as a branching fraction.
The fit-related uncertainties come mainly from the fitting ranges, signal shape, nonpeaking background shape, peaking background shape and the number of the peaking background events.
The uncertainties due to the fit range are considered by varying the fit ranges of 5 MeV=c 2 or 10 MeV=c 2 . And the maximum change on the results is taken as the systematic uncertainty from fit range. In the fit to the Mð4π 0 Þ distribution, signal shape is taken from the MC simulation. To assess the uncertainty due to the signal shape, an alternative fit is performed by convolving a Gaussian function with a fixed resolution of 2.6 MeV and mean of 2.1 MeV which are obtained from a high purity control sample of J=ψ → γη 0 ; η 0 → π 0 π 0 η; η → γγ. The uncertainty from the nonpeaking background shape is determined by using a fourth-order Chebychev polynomial in place of the third-order Chebychev polynomial. To assess the uncertainty associated with the number of the peaking background events, its contribution is recalculated after varying the branching fractions of J=ψ → γη 0 and its cascade decays, η 0 → π 0 π 0 η and η → π 0 π 0 π 0 , within their uncertainties, and new fits are performed. The systematic uncertainty associated with peaking background shape is evaluated by convolving a Gaussian function with the resolution and mean value left free. Among these cases, the dominant fit-model uncertainty arises from fitting range [0.705, 1.095] GeV=c 2 , and it changes the upper limit at the 90% C.L. to B UL ¼ 4.88 × 10 −5 .
The other category of systematic uncertainties, summarized in Table II, has contributions from the knowledge of the photon detection efficiency, the efficiency of the kinematic fit, signal model, MC statistical uncertainty for the detection efficiency, the branching fractions of the subdecays involved in the signal process, and the total number of J=ψ events.
The uncertainty from the photon detection is investigated with a high purity control sample of J=ψ → π þ π − π 0 . It is found that the differences between data and MC simulation are 0.5% and 1.5% for each photon deposited in the barrel and end cap of the EMC, respectively. With the same approach as used in Ref. [25], the uncertainty on the detection efficiency for each photon in the signal decay is estimated to be 0.53%, and thus the nine photons in the final state induce an overall uncertainty of 4.8%. The uncertainty associated with the kinematic fit is estimated by adjusting the components of the photonenergy error matrix in the signal MC sample to reflect the known difference in the resolution between data and the MC simulation [26]. From the study of the ψð3686Þ → γχ c1 ðχ c1 → 4π 0 Þ decay [27], it is known that the energy resolution in data is 4% wider than in a MC simulation. The relative difference in efficiency, 4.1%, is taken as the systematic uncertainty from the kinematic fit.
In the normal fit, the cascade decay η 0 → 4π 0 is described with the decay amplitude in Ref. [8]. A fit with an alternative signal model replacing η 0 → 4π 0 decay amplitude with a phase space (PHSP) distribution is performed. The change of the efficiency, 2.0%, is taken as the uncertainty due to the signal model.
The relative uncertainty in the knowledge of the branching fractions of J=ψ → γη 0 and π 0 → γγ [11] induces a corresponding uncertainty on the calculated upper limit on the branching fraction. The uncertainty of the detection efficiency, 0.4%, caused by MC statistics is also taken as a source of systematic uncertainty. The number of J=ψ events is determined from the measured number of hadronic decays and is found to be ð1310.6 AE 7.0Þ × 10 6 [13], which corresponds to a relative uncertainty of 0.54%.
Assuming all systematic uncertainties presented in Table II are independent, the total relative uncertainty is obtained to be 7.3%, by adding all individual uncertainties in quadrature.

VI. RESULT
The final upper limit on the branching fraction is determined by convolving the likelihood distribution L with the systematic uncertainties to obtain the smeared likelihood L smear , In this exercise all components listed in Table II, whatever their nature, can be considered as an uncertainty on the detection efficiency ε. The nominal efficiency value isε, σ ε is the absolute total systematic uncertainty on the efficiency, and B is the branching fraction of η 0 → 4π 0 . Figure 4 shows the normalized likelihood distribution after taking all systematic uncertainties into account. The corresponding upper limit of the branching fraction of η 0 → 4π 0 at the 90% C.L. is determined to be 4.94 × 10 −5 .

VII. SUMMARY
Using a sample of 1.31 × 10 9 J=ψ events collected with the BESIII detector, a search for the decay η 0 → 4π 0 is performed via J=ψ → γη 0 . No evidence for the rare decay η 0 → 4π 0 is found, and an upper limit of Bðη 0 → 4π 0 Þ < 4.94 × 10 −5 is set at the 90% confidence level. This limit is approximately a factor of 6 smaller than the previous most stringent result [10].
The current limit is still far to reach the theoretical predication with a level of 10 −8 [8]. Further studies of η 0 rare decays are still necessary to test the ChPT and VMD model and look for the CP-violation (S-wave) η 0 → 4π 0 decay. A sample of 10 10 J=ψ events has now been collected at BESIII, which will allow for even more sensitive searches to be performed for this important decay mode. University of Groningen (RuG) and the Helmholtzzentrum fuer Schwerionenforschung GmbH (GSI), Darmstadt.