Search for intermediate resonances and dark gauge bosons in J=ψ

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 , J. C. Chen, M. L. Chen, S. J. 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, Z. 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, Cheng Li, D. M. Li, F. Li, F. Y. Li, G. Li, H. B. Li, H. J. Li, J. C. Li, J. W. Li, Ke Li, L. K. Li, Lei Li, P. L. Li, P. R. Li, Q. Y. 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. 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, H. R. Qi, M. Qi, T. Y. Qi, S. Qian, C. F. Qiao, N. Qin, 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, Zongyuan Wang, T. Weber, D. H. Wei, P. Weidenkaff, 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, 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, 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, L. 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, PHYSICAL REVIEW D 102, 052005 (2020)

We report on an analysis of the decay J=ψ → γπ 0 η 0 using a sample of ð1310.6 AE 7.0Þ × 10 6 J=ψ events collected with the BESIII detector. We search for the CP-violating process η c → π 0 η 0 and a dark gauge boson U 0 in J=ψ → U 0 η 0 ; U 0 → γπ 0 ; π 0 → γγ. No evidence of an η c signal is observed in the π 0 η 0 invariantmass spectrum and the upper limit of the branching fraction is determined to be 5.6 × 10 −5 at the 90% confidence level. We also find no evidence of U 0 production and set upper limits at the 90% confidence level on the product branching fraction BðJ=ψ → U 0 η 0 Þ × BðU 0 → π 0 γÞ in the range between ð0.8 − 6.5Þ × 10 −7 for 0.2 ≤ m U 0 ≤ 2.1 GeV=c 2 . In addition, we study the process J=ψ → ωη 0 with ω → γπ 0 . The branching fraction of J=ψ → ωη 0 is found to be ð1.87 AE 0.09 AE 0.12Þ × 10 −4 , where the first uncertainty is statistical and the second is systematic, with a precision that is improved by a factor of 1.4 over the previously published BESIII measurement. DOI

I. INTRODUCTION
The Standard Model (SM) has been successful in explaining a wide variety of experimental data; however it fails to explain several observations, such as dark matter, the baryon asymmetry in the Universe, the neutrino masses, and so on. Therefore, in recent years the search for new physics beyond the SM is one of the important activities of particle physicists worldwide. The BESIII (Beijing Electron Spectrometer) experiment is currently searching for beyond-the-SM physics using low-energy e þ e − collision data. This is complementary to experiments conducted at the Large Hadron Collider (LHC) at CERN, which use high-energy hadron collision data. Huge data samples accumulated by the BESIII detector and taken at centerof-mass energies corresponding to the masses of various charmonium resonances [J=ψ, ψð3686Þ and ψð3770Þ] offer a unique sensitivity to search for forbidden decays and dark matter particles in the low-energy region [1].
Charge conjugation and parity symmetry (CP) violation has only been observed in weak interactions, which in the SM, originates from a single complex phase in the Cabibbo-Kobayashi-Maskawa (CKM) quark-mixing matrix [2]. Therefore, searches for this phenomenon will provide new insights and will help to determine whether the phase in the CKM mixing matrix is the sole source of CP violation or whether there are other sources. The production of heavy pseudoscalar mesons, e.g., η, η 0 , and η c , in J=ψ decays offers an opportunity to test this fundamental symmetry. In the SM, the decays of η=η 0 → ππ can proceed only via the weak interactions and the expected branching fractions are at a level of 10 −29 -10 −27 [3], which are experimentally inaccessible. In the case of the CP violation taking place in an extended Higgs sector [3], the branching fraction of η → ππ may reach the level of 10 −12 , which is considerably larger than the expectation in the SM. The decay of an η c (J PC ¼ 0 −þ ) to two pseudoscalar mesons is forbidden due to CP conservation. The observation of these forbidden decays will be a clear indication of new physics beyond the SM. Using a sample of 225 million J=ψ events, BESIII reports the results of the search for η c → π þ π − and η c → π 0 π 0 and upper limits on the branching fractions are presented at the 90% confidence level (C.L.) [4]. In this paper, we present the first experimental search for η c → π 0 η 0 .
Except for gravitational effects, we still know very little about the constituents and interactions of dark matter. One possible model candidate for dark matter is an additional gauge boson [5,6]. If this additional boson corresponds to an extra Uð1Þ gauge symmetry, it is referred to as a "dark photon." A dark photon with a mass in the sub-GeV range can couple to the SM via kinetic mixing with the ordinary photon and parametrized by the mixing strength [5]. The dark photon occurs naturally in many proposed models and has been invoked to explain various experimental and observational anomalies [7]. This new gauge boson referred to as U 0 has the same quantum number, J PC ¼ 1 −− , as the ω meson. In the past, BESIII has reported on a search for the dark gauge photon in the initial-state radiation (ISR) reactions e þ e − → U 0 γ ISR → l þ l − γ ISR (l ¼ μ, e) [8] and electromagnetic Dalitz decays J=ψ → U 0 η=η 0 → e þ e − η=η 0 [9,10]. The same ISR method has been used by the BABAR experiment [11]. The BELLE and KLOE Collaborations report a search for a dark vector gauge boson decaying to π þ π − , where the dark vector gauge boson mass spans a range from 290 to 520 MeV=c 2 [12] and 519 to 973 MeV=c 2 [13], respectively.
In this paper, using a sample of 1.31 × 10 9 J=ψ events collected with the BESIII detector, we present the first study of J=ψ → γπ 0 η 0 , which allows us to search for the CP-violating decay of η c → π 0 η 0 and to search for a new gauge boson [5] by investigating the γπ 0 -mass spectrum. Additionally, we present the most accurate measurement of the J=ψ → ωη 0 branching fraction [current BESIII measurement value is (2.08 AE 0.30 AE 0.14Þ × 10 −4 [14] ].

II. THE BESIII EXPERIMENT AND MONTE CARLO SIMULATION
The BESIII detector is a cylindrical magnetic spectrometer [15] located at the Beijing Electron Positron Collider (BEPCII) [16], with an acceptance of charged particles and photons of 93% over 4π solid angle. 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 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. Particle identification (PID) for charged pions is performed by exploiting the TOF information and the specific ionization energy loss, dE=dx, measured by the MDC. The TOF and dE=dx information is combined to form PID probability for the pion, kaon, and proton hypotheses; each track is assigned to the particle type that corresponds to the hypothesis with the highest probability.
Simulated samples produced with the GEANT4-based [17] Monte Carlo (MC) package, which includes the geometric description [18,19] of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate the backgrounds. The inclusive MC sample consists of the production of the J=ψ resonance, and the continuum processes incorporated in KKMC [20].
The known decay modes are generated using the EVTGEN package [21] using branching fractions taken from the Particle Data Group (PDG) [22], and the remaining unknown decays from the charmonium states with the LUNDCHARM package [23]. The final-state radiations from charged final-state particles are incorporated with the PHOTOS package [24].
The three-body decay of J=ψ → γπ 0 η 0 without any intermediate states is simulated with a model based on a phase-space distribution of the final-state particles. The decays of J=ψ → γη c ; U 0 η 0 ; γη 0 , and ωη 0 are generated with an angular distribution of 1 þ cos 2 θ γ , where θ γ is the angle of radiative photon relative to the positron beam direction in the J=ψ-rest frame, while the subsequent η c ðη 0 Þ decays are generated with a phase-space model and the U 0 ðωÞ → γπ 0 decay is modeled by a P-wave [21].

III. EVENT SELECTION
Candidates of J=ψ → γπ 0 η 0 ; η 0 → π þ π − η; π 0 → γγ; η → γγ are required to have two oppositely charged tracks and at least five photon candidates. All charged tracks must originate from the interaction point with a distance of closest approach less than 10 cm in the beam direction and less than 1 cm in the transverse plane. Their polar angles, θ, with respect to the beam direction are required to satisfy j cos θj < 0.93.
Electromagnetic showers are reconstructed from clusters of firing EMC crystals. The energy deposited in nearby TOF counters is included to improve the reconstruction efficiency and energy resolution. The showers of the photon candidate must have a minimum energy of 25 MeV in the barrel region (j cos θ γ j < 0.80) and 50 MeV in the end-cap region (0.86 < j cos θ γ j < 0.92), where θ γ is the polar angle of the photon. To suppress showers originating from charged particles, a photon candidate must be separated by at least 10°from the nearest charged track. To suppress noise and energy deposits unrelated to the event, the time at which the photon is recorded in the EMC after the e þ e − collision is required to be within 0 ≤ t ≤ 700 ns.
After selecting the charged tracks and showers, a fourconstraint (4C) kinematic fit to the J=ψ → π þ π − 5γ hypothesis is performed using energy-momentum conservation. For events with more than five photon candidates, the combination with the smallest χ 2 4C is retained. To suppress background events with six photons in the final states, the χ 2 4C of the π þ π − 5γ hypothesis is required to be less than that for the π þ π − 6γ hypothesis.
To distinguish the photon from π 0 and η decays, we define the variable χ 2 This variable is used to choose from the five photon candidates two pairs of photons with two-photon invariant masses (M γγ ) closest to the nominal π 0 (m π 0 ) and η (m η ) masses. σ π 0 (σ η ) refers to the experimental mass resolution for a π 0 (η) decay. The four-photon combination with the smallest value for χ 2 π 0 η is chosen. To improve the mass resolution and to further suppress background events, we subsequently perform a fiveconstraint kinematic (5C) fit imposing energy-momentum conservation and an η-mass constraint under the hypothesis of π þ π − γγγη, where the η candidate is reconstructed with the selected pair of photons as described above. Events with a χ 2 5C less than 30 are accepted for further analysis.
To select π 0 candidates, the invariant mass of the two photons from π 0 decay, M γγ , must satisfy jM γγ − m π 0 j < 15 MeV=c 2 . To suppress background events with multi-π 0 in the final states, we require that the invariant mass of the radiative photon and any photon from the η decay is outside the π 0 -mass region of ½0.115; 0.155 GeV=c 2 . To select η 0 candidates, we calculate for each event the π þ π − η invariant mass, M π þ π − η , and require that jM π þ π − η −m η 0 j<15MeV=c 2 , where m η 0 is the nominal η 0 mass.
After applying the selection criteria, we obtain the π 0 η 0 invariant-mass distribution as shown in Fig. 1. No evident η c peak is seen. We found that the dominant background events are from decays with the η 0 as an intermediate state, such as J=ψ → γπ 0 η 0 , J=ψ → ωη 0 , and J=ψ → γη 0 , and the corresponding contributions are displayed in Fig. 1(a) as well. Other background contributions (non-η 0 background) are estimated from events for which the reconstructed η 0 mass falls within the η 0 -sideband regions (0.903<M π þ π − η <0.933GeV=c 2 and 0.983<M π þ π − η < 1.013GeV=c 2 ). The sum of the above contributions gives a reasonable description of the data.
An unbinned maximum-likelihood fit is performed to determine the signal yield as shown in Fig. 1(b). In the fit, the probability density function (PDF) of the signal is described by a MC simulated shape and the widths and masses are fixed to the world average values taken from the PDG [22]. The background shape from the J=ψ → γη 0 channel is described with a MC simulated shape, and the yield is fixed according to the published branching fractions [22]. The other nonpeaking background is described by a first-order Chebyshev polynomial function. The signal yield is N sig ¼ 7.2 AE 7.6 and the statistical significance of the η c signal is calculated to be 1.0σ using ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi −2 lnðL stat 0 =L stat max Þ p , where L stat max and L stat 0 are the maximum-likelihood values with the signal yield left free and fixed at zero, respectively. In addition, to account for the additive systematic uncertainties related to the fits, the fit range and the background shape are varied and the maximum signal yield among these cases is obtained as shown in Fig. 1(b).
Using the same selection criteria as used to search for η c → π 0 η 0 , we study the γπ 0 -mass (M γπ 0 ) distribution as shown in Fig. 2. A clear ω peak from J=ψ → ωη 0 decays can be observed. There is also a small background contribution from J=ψ → γη 0 decays which is smoothly distributed in the low-mass region of the M γπ 0 distribution. The contributions from non-η 0 backgrounds are described by events that are selected in the η 0 -sideband regions, 0.903 < M π þ π − η < 0.933 GeV=c 2 and 0.983 < M π þ π − η < 1.013 GeV=c 2 .
We search for the U 0 signal in steps of 10 MeV=c 2 in the M γπ 0 distribution ranging from 0.2 to 2.1 GeV=c 2 and excluding the mass region around the ω peak 1 GeV=c 2 , respectively. The U 0 signal and the tail of the ω signal are described by MC-simulated shapes, and the remaining background contribution is modeled with a linear Chebyshev polynomial. To take into account the additive systematic uncertainties related to the fits, alternative fits with different fit range and background shape are also performed, and the maximum upper limit among these cases has been selected. The number of extracted signal events, the significance, and the detection efficiency as a function of M U 0 are shown in Fig. 3. The largest local significance defined as before is computed to be 2.4σ at M U 0 ¼ 1.78 GeV=c 2 , the corresponding p-value is calculated to be 0.89. No significant signal for U 0 → γπ 0 is found.

VI. BRANCHING FRACTION MEASUREMENT
OF J=ψ → ωη 0 Figure 4 shows the mass distribution of M π þ π − η versus M γπ 0 . Events originating from the J=ψ → ωη 0 decay are  clearly visible. To extract the number of ωη 0 events, an unbinned extended maximum-likelihood fit using a twodimensional (2D) PDF including both variables, M π þ π − η and M γπ 0 , with the requirements of 0.6 < M γπ 0 < 1.0 GeV=c 2 and 0.908 < M π þ π − η < 1.008 GeV=c 2 is performed. Assuming zero correlation between the two discriminating variables M γπ 0 and M π þ π − η , the composite PDF in the 2D fit is constructed as follows: where the signal shapes for the ω (F ω sig ) and η 0 (F η 0 sig ) responses are modeled with a relativistic Breit-Wigner function convoluted with a Gaussian function. The widths and masses of the ω and η 0 are fixed in the fit. The parameters of the Gaussian function are free in the fit. N sig is the number of J=ψ → ωη 0 ; ω → γπ 0 ; η 0 → π þ π − η signal events. The backgrounds are divided into three categories, namely non-ω peaking background, non-η 0 peaking background, and non-ωη 0 background. The parameters N non−ω bkg , N non−η0 bkg , and N non−ωη0 bkg are the corresponding three background yields. The background shapes, F non−ω bkg and F non−η0 bkg , related to M γπ 0 and M π þ π − η , respectively, are described by first-order Chebyshev polynomials and all their corresponding parameters are free in the fit.
The fit results in N sig ¼ 506 AE 25 signal events. The projection plots of the fit on the M γπ 0 and M π þ π − η distributions are shown in Figs. 5(a) and 5(b), respectively.

VII. SYSTEMATIC UNCERTAINTY
The sources of systematic uncertainties and their corresponding contributions to the measurements of the upper limits and branching fraction are summarized in Table I.
The uncertainty of the number of J=ψ events is determined to be 0.54% by an analysis of inclusive hadronic events in J=ψ decays [25].
The uncertainty of the MDC tracking efficiency for each charged pion is studied by analyzing a nearly backgroundfree sample of J=ψ → ρπ events. The difference between the data and MC simulation is less than 1.0% for each charged track [26] whose value is taken as a systematic uncertainty. Similarly, the uncertainty related to the PID efficiencies of pions is also studied with the data sample, J=ψ → ρπ, and the average difference of the PID efficiencies between data and MC simulation is determined to be 1.0% for each charged pion, which is then taken as the corresponding systematic uncertainty. The photon detection efficiency is studied with the control sample J=ψ → π þ π − π 0 [27]. The difference in efficiency between the data and that predicted by MC simulations is found to be 0.5% per photon in the EMC barrel and 1.5% per photon in the end-cap part of the EMC. In our case, the uncertainty is on average 0.6% per photon whose value is obtained by weighting the uncertainties according to the angular distribution of the five photons found in our data sample. Thus, the uncertainty associated with the five reconstructed photons is 3.0%.  The uncertainty associated with the 5C kinematic fits comes from the inconsistency of the track helix parameters between the data and MC simulation. The helix parameters for the charged tracks of MC samples are corrected to eliminate part of the inconsistency, as described in Ref. [28]. We take half of the differences on the selection efficiencies with and without the correction as an estimate of the corresponding systematic uncertainties, which results in 0.4%.
Due to the difference in the mass resolution between the data and MC, the uncertainty related to the η 0 and π 0 masswindow requirements is investigated by smearing the MC simulation in accordance with the signal shape of the data.
The difference of the detection efficiency before and after smearing is assigned as the systematic uncertainty for the η 0 and π 0 mass-window requirements and found to be 0.2% and 1.1%, respectively.
The systematic uncertainty related to the finite statistics used by the MC simulation to obtain the overall reconstruction efficiency is calculated as ffiffiffiffiffiffiffiffiffi ffi ϵð1−ϵÞ n q , where ϵ is the detection efficiency and n is the number of generated MC events of the signal process. The corresponding systematic uncertainty is determined to be 1.0%.
The systematic uncertainties that affect the upper limits on the branching fraction of η c → π 0 η 0 and U 0 → γπ 0 are considered in two categories: additive and multiplicative. The additive systematic uncertainties on the fit range and background shapes are already accounted for in the analysis procedure that is applied to obtain the maximum upper limit of the signal yield. Therefore, here we only consider these uncertainties for the J=ψ → ωη 0 study. To study the uncertainty from the fit range, the fit is repeated with different fit ranges, and the largest difference in the signal yield, 1.8%, is taken as the systematic uncertainty. The uncertainty associated with the background shape in the fits to the M γπ 0 distribution is estimated using alternative fits by changing the linear Chebyshev polynomial to a second-order Chebyshev polynomial. The difference in signal yield (0.6%) is taken as the systematic uncertainty.
The uncertainty associated with the 2D fits of the J=ψ → ωη 0 channel is estimated by taking all parameters as free parameters in the fit. The change in signal yield (1.4%) is taken as the systematic uncertainty. The systematic uncertainty due to the π 0 veto is evaluated by varying the requirement on the mass window, and the difference in yield compared to the nominal choice (1.1%) is assigned as the systematic uncertainty.
The branching fractions of the intermediate processes of J=ψ → γη c , ω → γπ 0 , η 0 → π þ π − η, η → γγ, and π 0 → γγ  are taken from the PDG [22] and their errors are considered as a source of systematic uncertainty. For each case, the total systematic uncertainty is given by the quadratic sum of the individual contributions, assuming all sources to be independent.

VIII. RESULTS
Since no evident η c signal is seen in M π 0 η 0 , a Bayesian method is used to obtain the upper limit of the signal yield at the 90% C.L. To determine the upper limit on the η c signal, a series of unbinned maximum-likelihood fits are performed to the π 0 η 0 -mass spectrum with a varying number of expected η c signals. From this, we obtain the dependence of the likelihood on the number of signal events from which we extract the upper limit, taking into account the multiplicative systematic uncertainties as follows [29]: Here, L stat and L 0 are the likelihood curves before and after the inclusion of the multiplicative systematic uncertainty.
where Δ is the relative deviation of the estimated branching fraction from the nominal value, and σ syst is the multiplicative systematic uncertainties given in Table I. The branching fraction for a particular decay process is computed as where N sig is the number of extracted signal yield, ϵ is the signal selection efficiency, and B is the secondary branching fraction of the corresponding decay process. The normalized likelihood distribution for J=ψ → γη c ðη c → π 0 η 0 Þ candidates is shown in Fig. 6. The upper limit at the 90% C.L. of the signal yield (N UL ) and detection efficiency are determined to be 19.0% and 9.3% respectively, resulting in a branching fraction Bðη c → π 0 η 0 Þ of less than 5.6 × 10 −5 .
Due to no evident U 0 signal seen in M γπ 0 , we compute the upper limit on the product branching fraction BðJ=ψ → U 0 η 0 Þ × BðU 0 → π 0 γÞ at the 90% C.L. as a function of M U 0 using a Bayesian method after incorporating the systematic uncertainty by smearing the likelihood curve with a Gaussian function with a width of the systematic uncertainty as follows: where, L and L smear are the likelihood curves before and after the consideration of the systematic uncertainty. ϵ,ε, and σ ϵ are the detection efficiency, nominal efficiency, and the absolute total systematic uncertainty on the efficiency, respectively. As shown in Fig. 7, the combined limits on product branching fraction BðJ=ψ → U 0 η 0 Þ × BðU 0 → π 0 γÞ are established at the level of ð0.8-6.5Þ × 10 −7 for 0.2 ≤ M U 0 ≤ 2.1 GeV=c 2 .
With a detection efficiency of 14.9% obtained from a MC simulation, we obtain a branching fraction for the J=ψ → ωη 0 process of ð1.87 AE 0.09 AE 0.12Þ × 10 −4 , where the first uncertainty is statistical and the second systematic.

IX. SUMMARY
Using a sample of ð1310.6 AE 7.0Þ × 10 6 J=ψ events collected with the BESIII detector, the decay of  J=ψ → γη 0 π 0 is studied. We search for the CP-violating decay η c → π 0 η 0 and a dark gauge boson U 0 in J=ψ → U 0 η 0 ; U 0 → γπ 0 ; π 0 → γγ. No significant η c signal is observed in the π 0 η 0 invariant-mass spectrum, and the upper limit on the branching fraction is determined to be 5.6 × 10 −5 at the 90% C.L. Except for a clear ω peak in the γπ 0 mass spectrum, no significant excess is seen for any mass hypothesis in the range of 0.2 ≤ M U 0 ≤ 2.1 GeV=c 2 . The upper limits on the product branching fractions are calculated to be ð0.8-6.5Þ × 10 −7 at the 90% C.L. Due to lack of the theoretical predictions on the BðJ=ψ → U 0 η 0 Þ, we do not present the upper limit on the coupling of the dark vector gauge boson. In case of corresponding theoretical calculations in the future, we would like to present the detailed information, e.g., the detection efficiency, signal yield, and branching fraction, as shown in Tables II and III in the Appendix. The detection efficiencies increase first and then decrease, and the jumping of individual points is within the range of statistical error. In this case, it would be easy for readers or theorists to extract the coupling in case the corresponding prediction is available.
In addition, the branching fraction of J=ψ → ωη 0 is measured to be ð1.87 AE 0.09 AE 0.12Þ × 10 −4 , where the first uncertainty is statistical and the second systematic. This result is consistent with the previously published BESIII measurement but with an improvement in accuracy by a factor of 1.4.   The results of signal yield (N sig ), the upper limit at the 90% C.L. of the signal yield (N UL ), efficiency (ϵ), and branching fraction (B) as a function of M U 0 .