Search for the rare decay $\eta'\rightarrow\pi^{0}\pi^{0}\pi^{0}\pi^{0}$ at BESIII

Based on a sample of 1.31 billion $J/\psi$ events collected with the BESIII detector, we perform a search for the rare decay $\eta'\rightarrow 4\pi^{0}$ via $J/\psi\rightarrow\gamma\eta'$. No significant $\eta'$ signal is observed in the invariant mass spectrum of 4$\pi^{0}$. With a Bayesian approach, the upper limit on the branching fraction of $\eta'\rightarrow 4\pi^{0}$ is determined to be $\mathcal{B}(\eta'\rightarrow 4\pi^{0})$ $<4.94\times10^{-5}$ at the 90\% confidence level, which is a factor of six smaller than the previous experimental limit.

Using a sample of 1.31 × 10 9 J/ψ events collected with the BESIII detector, we perform a search for the rare decay η ′ → 4π 0 via J/ψ → γη ′ . No significant η ′ signal is observed in the 4π 0 invariant mass spectrum. With a Bayesian approach, the upper limit on the branching fraction is determined to be B(η ′ → 4π 0 ) < 4.94 × 10 −5 at the 90% confidence level, which is a factor of six smaller than the previous experimental limit.

I. INTRODUCTION
The η ′ meson has a special role in improving the understanding of low-energy Quantum Chromodynamics (QCD), and studies of its decays have attracted considerable theoretical and experimental attention [1,2]. In addition to its important role in testing the fundamental discrete symmetries and searching for processes beyond the Standard Model (SM), η ′ decays offer unique opportunities to test chiral perturbation theory (ChPT) and the vector-meson dominance (VMD) model.
In theory, η ′ → 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 η ′ → 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 contributed to this decay is at a level of 10 −23 [3,4], which is far beyond current experimental sensitivity. However, 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 ChPT and VMD models predict the branching fraction caused by D-wave CP -conserving to be at the level of 10 −8 [5]. Therefore, an observation of η ′ → 4π 0 with a branching fraction at a level of 10 −8 would indicate ChPT and VMD models are reliable in calculating the decay η ′ (η) → 4π 0 .
So far the decay η ′ → 4π 0 has not been observed. About three decades ago the first attempt to search for this mode was performed by the joint CERN-IHEP experiment and the upper limit on the branching fraction was determined to be B(η ′ → 4π 0 ) < 5×10 −4 at the 90% confidence level (C.L.) [6]. The more recent upper limit of 3.2 × 10 −4 at 90% C.L. was obtained by the GAMS-4π experiment [7].

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector [11] is a magnetic spectrometer located at the Beijing Electron Positron Collider (BEPCII) [12]. 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 muonidentifier 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 [13] 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 [14]. The inclusive MC sample consists of the production of the J/ψ resonance, and the continuum processes incorporated in kkmc [14]. The known decay modes are modeled with evtgen [15] using branching fractions taken from the Particle Data Group (PDG) [8], and the remaining unknown decays of the charmonium states with lundcharm [16]. The final state radiation (FSR) from charged final-state particles are incorporated with the photos package [17].

III. EVENT SELECTION
In this analysis, the pseudoscalar mesons η ′ and π 0 are reconstructed in the modes η ′ → 4π 0 and π 0 → γγ. Candidate J/ψ → γ4π 0 decays are chosen by selecting events with at least nine isolated photons and no charged tracks. Photon candidates are reconstructed from clusters of energy deposited in the EMC. Photon candidates are required to have at least 25 MeV of energy for barrel showers (| cos θ| ≤ 0.8) or 50 MeV for end-cap showers (0.86 ≤ | cos θ| ≤ 0.92).
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/ψ → γη ′ , 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. Then an eight-constraint (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. To ensure a good description of data, a signal MC simulation is modeled with the decay amplitude in Ref. [5], which assumes the pion-loop contribution as shown in Fig. 1(a) (the contribution from f 2 mesons shown in Fig. 1(b) is considered to be much smaller). After the full event selection, the detection efficiency (ε) is determined to be 2.46%. Figure 3 shows the 4π 0 mass spectrum, M (4π 0 ), after selection. No significant η ′ signal is evident.
In order to investigate possible sources of contamination we apply the selection to an inclusive MC sample of 1.2 ×10 9 J/ψ events. All the significant background components are listed in Table I, showing their expected contribution to the data sample. Since the decay J/ψ → nπ 0 , where n is the number of π 0 mesons, is forbidden because of C conservation, we find that the background comes mainly from decays with the same final state as the signal, for example J/ψ → ηω, η → π 0 π 0 π 0 , ω → γπ 0 or from radiative decays with more than four π 0 s in the final states, of which the dominant mode is J/ψ → γη ′ , η ′ → π 0 π 0 η, η → π 0 π 0 π 0 . This latter channel contributes to a broad structure around 0.88 GeV/c 2 in the 4π 0 mass spectrum as indicated by the dashed line in Fig. 3. No other source of background peaks in the η ′ mass region.

IV. FIT MODEL AND UPPER LIMIT ON BRANCHING FRACTION OF η ′ → 4π 0
An unbinned maximum likelihood fit is performed on the 4π 0 mass spectrum, allowing for background contributions and a possible signal component. The peaking background J/ψ → γη ′ , η ′ → π 0 π 0 η, η → π 0 π 0 π 0 is described by a distribution obtained from a dedicated MC generator [18,19], and its normalisation fixed from 8C distribution of candidate events in the η ′ signal region showing data and the contribution from the considered background sources. The arrow indicates the selection requirement of χ 2 8C < 30.
knowledge of the number of J/ψ events and the branching fractions of the decays. The non-peaking background contribution is modeled with a third-order Chebychev polynomial function. The signal shape and resolution is taken from simulation. The fit result is shown superimposed in Fig. 3, where the signal contribution is negligible. The M(4π 0 ) distribution in data, together with the total fit result and the contributions from non-peaking background and the peaking background J/ψ → γη ′ , η ′ → π 0 π 0 η, η → π 0 π 0 π 0 . Also shown is the expected shape of the signal contribution, with arbitrary normalisation.
A Bayesian approach is used to determine an upper limit on the branching fraction of η ′ → 4π 0 . Many fits are performed for different assumed values of the signal yield N , which is a fixed parameter, 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 ± 7.0) × 10 6 is the number of J/ψ events [10], ε is the detection efficiency, B(J/ψ → γη ′ ) and B(π 0 → γγ) are the branching fractions of J/ψ → γη ′ and π 0 → γγ, respectively, which are taken from Ref. [8].
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 η ′ → 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, non-peaking background shape, peaking background shape and the number of the peaking background events.
The systematic uncertainty from the fit ranges is estimated by varying them by ±5 MeV/c 2 .
In the fit to the M(4π 0 ) distribution, signal shape is taken from 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 η, η → γγ.
The uncertainty from the non-peaking 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/ψ → γη ′ and its cascade decays, η ′ → π 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 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, summarised in Table II, has contributions from the knowledge the photon detection efficiency, the efficiency of the kinematic fit, signal model, the branching fractions of the sub-decays 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. [20], 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 resolution between data and MC simulation [21]. From the study of ψ(3686) → γχ c1 (χ c1 → 4π 0 ) decay [22], it is known that the energy resolution in data is 4% wider than in 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 η ′ → 4π 0 is described with the decay amplitude in Ref. [5]. A fit with an alternative signal model replacing η ′ → 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/ψ → γη ′ and π 0 → γγ [8] induces a corresponding uncertainty on the calculated upper limit on the branching fraction. The number of J/ψ events is determined from the measured number of hadronic decays, and is found to be (1310.6 ± 7.0) × 10 6 [10], 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 η ′ → 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 η ′ → 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 η ′ → 4π 0 is performed via J/ψ → γη ′ . No evidence for the rare decay η ′ → 4π 0 is found, and an upper limit of B(η ′ → 4π 0 ) < 4.94 × 10 −5 is set at the 90% confidence level. This limit is approximately a factor of six smaller than the previous most stringent result [7].
The current limit is still far above the level of 10 −8 , which is predicted in theory. Further studies of η ′ rare decays are still necessary to test ChPT and VMD model and look for the CP -violation (S-wave) η → 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.