Measurement of Singly Cabibbo-Suppressed Decays $D^{0}\to\pi^{0}\pi^{0}\pi^{0}$, $\pi^{0}\pi^{0}\eta$, $\pi^{0}\eta\eta$ and $\eta\eta\eta$

Using a data sample of $e^+e^-$ collision data corresponding to an integrated luminosity of 2.93 $fb^{-1}$ collected with the BESIII detector at a center-of-mass energy of $\sqrt{s}= 3.773~GeV$,we search for the singly Cabibbo-suppressed decays $D^{0}\to\pi^{0}\pi^{0}\pi^{0}$, $\pi^{0}\pi^{0}\eta$, $\pi^{0}\eta\eta$ and $\eta\eta\eta$ using the double tag method. The absolute branching fractions are measured to be $\mathcal{B}(D^{0}\to\pi^{0}\pi^{0}\pi^{0}) = (2.0 \pm 0.4 \pm 0.3)\times 10^{-4}$, $\mathcal{B}(D^{0}\to\pi^{0}\pi^{0}\eta) = (3.8 \pm 1.1 \pm 0.7)\times 10^{-4}$ and $\mathcal{B}(D^{0}\to\pi^{0}\eta\eta) = (7.3 \pm 1.6 \pm 1.5)\times 10^{-4}$ with the statistical significances of $4.8\sigma$, $3.8\sigma$ and $5.5\sigma$, respectively, where the first uncertainties are statistical and the second ones systematic. No significant signal of $D^{0}\to\eta\eta\eta$ is found, and the upper limit on its decay branching fraction is set to be $\mathcal{B}(D^{0}\to\eta\eta\eta)<1.3 \times 10^{-4}$ at the $90\%$ confidence level.


Introduction
The study of charmed meson decays, which involve both strong and weak interactions, is an interesting and challenging field in particle physics. Experimental measurements of charmed meson decays yield essential information for understanding the intrinsic decay mechanism and provide inputs to theoretical calculations and predictions. For example, Ref. [1] suggests that the measurement of the branching fraction (BF) of the hadronic decay D 0 → π 0 π 0 π 0 may shed light on the understanding of the role of isospin symmetry in D 0 decays to three-pion final states, and the isospin nature of the non-resonant contribution. Additionally, the study of the hadronic decays of charmed mesons provides important inputs for the studies of B physics [2].
The singly Cabibbo-suppressed (SCS) decays of the D 0 meson to three neutral pseudoscalar particles, D 0 → π 0 π 0 π 0 , π 0 π 0 η, π 0 ηη and ηηη, proceed dominantly through internal W -emission and W -exchange diagrams. Experimental studies of these decays are challenging due to the dominant presence of neutral particles (photons) in the final states, low BFs and high backgrounds. Until now, only a search for D 0 → π 0 π 0 π 0 decay has been performed by the CLEO Collaboration with a ψ(3770) data sample of 281 pb −1 in 2006 [3]. Using the "single tag" (ST) method, in which one D 0 orD 0 meson is found in each event, they obtained a BF upper limit of 3.5 × 10 −4 at the 90% confidence level (C.L.).
In this Letter, we present measurements of the BFs of the SCS decays D 0 → π 0 π 0 π 0 , π 0 π 0 η, π 0 ηη and ηηη with the "double tag" (DT) technique and a data sample corresponding to an integrated luminosity of 2.93 fb −1 [4], collected at a center-of-mass energy of √ s = 3.773 GeV with the BESIII detector at the BEPCII e + e − collider. Throughout the Letter, charge conjugate modes are always implied, unless explicitly mentioned.

BESIII Detector and Monte Carlo Simulation
BESIII [5] is a cylindrical spectrometer composed of a helium-gas-based main drift chamber (MDC), a plastic scintillator time-of-flight (TOF) system, a CsI(Tl) electromagnetic calorimeter (EMC), a superconducting solenoid providing a 1.0 T magnetic field, and a muon counter. The charged particle momentum resolution in the MDC is 0.5% at a transverse momentum of 1 GeV/c and the photon energy resolution in the EMC at 1 GeV, is 2.5% in the barrel region and 5.0% in the end-cap region. Particle identification (PID) combines the ionization energy loss (dE/dx) in the MDC with information from the TOF to identify particle types. More details about the design and performance of the detector are given in Ref. [5].
GEANT4-based [6] Monte Carlo (MC) simulation software is used to understand the backgrounds and to determine the detection efficiencies. The generator KKMC [7,8] is used to simulate the e + e − collision incorporating the effects of beam-energy spread and initial-state radiation (ISR). An inclusive MC sample including D 0D0 , D + D − and non-DD events, ISR production of ψ(3686) and J/ψ, and continuum processes e + e − → qq (q = u, d, s) is used to study the potential backgrounds. The known decay modes as specified in the Particle Data Group (PDG) [9] are generated by EVTGEN [10,11], while the remaining unknown decays of charmonium are modeled by LundCharm [12].

Analysis Strategy
At the ψ(3770) resonance, D 0D0 pairs are produced in a coherent 1 −− state without additional particles. A DT method, which was first developed by the MARK-III Collaboration [13,14], is used to measure the absolute BFs. We first select ST events in which ā D 0 meson is reconstructed in a specific hadronic decay mode. Then we search for D 0 decays in the remaining tracks, and DT events are those where D 0D0 pairs are fully reconstructed. The absolute BFs for D 0 decays are calculated by where the superscript 'sig' represents a specific D 0 signal decay, N α ST , ǫ α ST and ǫ sig,α DT are the yield of ST events, the ST detection efficiency and DT detection efficiency for a specific ST mode α, respectively, while N sig DT is the total yield for DT signal events, and B int is the product of the decay BFs for the intermediate states in the D 0 signal decay.

Data Analysis
Charged tracks are reconstructed from hits in the MDC and are required to have a polar angle θ satisfying | cos θ| < 0.93. The point of the closest approach of any charged track to the interaction point (IP) is required to be within 1 cm in the plane perpendicular to the beam and ±10 cm along the beam. Information from the TOF system and the dE/dx information in the MDC are combined to form PID C.L.s for the π and K hypotheses. Each track is assigned to the particle type with the highest PID C.L.
Photon candidates are reconstructed using clusters of energy deposited in the EMC crystals. The energy is required to be larger than 25 MeV in the barrel region (| cos θ| < 0.8) or 50 MeV in the end-cap region (0.86 < | cos θ| < 0.92). The energy deposited in nearby TOF counters is included to improve the reconstruction efficiency and energy resolution. The difference of the EMC time from the event start time is required to be within [0, 700] ns to suppress electronic noise and showers unrelated to the event.
The π 0 and η candidates are reconstructed from photon pairs by requiring the invariant masses M γγ to satisfy 115 < M γγ < 150 MeV/c 2 or 515 < M γγ < 570 MeV/c 2 , respectively. To improve the resolution, the photon pairs are fitted kinematically constraining their masses to the nominal π 0 or η masses [9], and the resulting energies and momenta of the two photons are used for subsequent analysis.
The ST candidates are selected by reconstruct-ingD 0 decays to K + π − , K + π − π 0 and K + π − π − π + . Two variables, the energy difference ∆E ≡ E D − E beam and the beam-energy-constrained mass M BC ≡ E 2 beam /c 4 − p 2 D /c 2 , are used to identify theD 0 candidates. Here, E beam is the beam energy, and E D (p D ) is the reconstructed energy (momentum) of theD 0 candidate in the e + e − center-of-mass system. ThoseD 0 candidates are accepted for further analysis that satisfy M BC > 1.83 GeV/c 2 and mode-dependent ∆E requirements, which are approximately three times the value of the resolution around theD 0 nominal mass [9], as summarized in Table 1. For each ST mode, if there is more than one candidate in the event, the one with the minimum |∆E| is selected.
The M BC distributions of the acceptedD 0 candidates are shown in Fig. 1, whereD 0 signals are observed with relatively low backgrounds. Binned maximum likelihood fits to the M BC distributions are performed to obtain the ST yields. In the fits, the signal shape is modeled by the MC simulated shape convolved with a Gaussian function representing the difference between data and MC simulation coming from the beamenergy spread, ISR, the ψ(3770) line shape, and resolution. The combinatorial background is modeled by an ARGUS function [15]. The ST yields are calculated by subtracting the integrated ARGUS background yields from the total events counted in the signal region 1.859 < M BC < 1.871 GeV/c 2 . The ST efficiency is studied using the same procedure on the inclusive MC sample. The resulting ST yields and the corresponding ST efficiencies are summarized in Table 1.
Points with an error bar are data, the blue solid lines are the total fit curves, the red dashed lines are the signal shapes, and the green long-dashed lines are the background shapes.
Candidates for the SCS decays, D 0 → π 0 π 0 π 0 , π 0 π 0 η, π 0 ηη and ηηη, are selected in the system recoiling against the taggedD 0 . Only events without any additional charged track are chosen. The D 0 signal decays are reconstructed with any combination of the selected π 0 and η candidates that have not been used in the ST side and do not share the same photon candidate. To distinguish the signal decay from combinatorial backgrounds, the energy difference ∆E and the beam-constrained mass M BC are also calculated for each accepted combination. A D 0 candidate is accepted if it satisfies a mode-dependent ∆E requirement, which corresponds to three times the value of the resolution around the ∆E peak based on MC simulation, as summarized in Table 2. The shift and asymmetry of the ∆E distributions are mainly due to the energy loss in the EMC for multi-photon final states. If there are multiple combinations for a given signal decay in an event, the one with the minimum |∆E| is selected.
Except for the decay D 0 → ηηη, MC studies indicate that the selected candidates have large backgrounds from D 0 → π 0 π 0 π 0 π 0 decay, which has a relatively large decay BF, and contain some background events and ǫ ηηη,α DT (in %)) efficiencies. The uncertainties are statistical only. BFs of π 0 and η decays to two photons are not included in the efficiencies.
However, MC studies indicate that backgrounds remain from photon mis-combinations in π 0 and η candidates. These are due to the matches of a good photon with noise in the EMC, which usually corresponds to a fake low energy photon. Furthermore, the MC indicates that this background can be reduced by requiring no other combination with the same final state and with χ 2 < 20. For instance for D 0 → π 0 π 0 π 0 , this requirement loses only 5% of signal events while it rejects 30% of mis-combination background.
With the above selection criteria, the M BC distributions of the accepted D 0 candidate events in data are shown in Fig. 2. The D 0 → π 0 π 0 π 0 , π 0 π 0 η and π 0 ηη signals are clear, but no obvious D 0 → ηηη signal is observed. The peaking backgrounds are dominated by the decay D 0 → π 0 π 0 π 0 π 0 , and the CF decays D 0 → K 0 S π 0 /η for D 0 → π 0 π 0 π 0 /η. The contributions from the cross feeds are small and will be considered in determining the signal yields. The miscombination background is negligible.
To determine the signal yields of the decays D 0 → π 0 π 0 π 0 , π 0 π 0 η, and π 0 ηη, unbinned maximum likelihood fits are performed to the M BC distributions. The probability density function (PDF) for signal is modeled with the MC simulated shape convolved with a Gaussian function representing the resolution difference and a potential mass shift between data and MC simulation. The peaking backgrounds from the CF decay D 0 → K 0 S π 0 /η (BKG I) and the decay D 0 → π 0 π 0 π 0 π 0 (BKG II) as well as the cross feeds (BKG III) are also included in the fit. The combinatorial background (BKG IV) is modeled by an ARGUS function [15]. The shapes of the various peaking backgrounds are modeled with those of MC simulations, and the corresponding magnitudes are fixed to the values estimated with a data driven method. We select a control sample of D 0 → π 0 π 0 π 0 π 0 from data with an approach similar to the signal selection, and obtain the yield N 4π 0 from a fit to the resulting M BC distribution. A mixed MC sample, which includes the possible resonant decays is generated with known BFs [9] and is subject to the selection criteria of D 0 → π 0 π 0 π 0 and D 0 → π 0 π 0 π 0 π 0 to evaluate the mis-identification rate ǫ 3π 0 and the detection efficiency ǫ 4π 0 , respectively. The magnitude of the background D 0 → π 0 π 0 π 0 π 0 in the selection of D 0 → π 0 π 0 π 0 is given by N 4π 0 ·ǫ 3π 0 /ǫ 4π 0 . Similar data driven approaches are applied to determine the magnitude of the peaking background D 0 → π 0 π 0 π 0 π 0 , the cross feed and the number of CF decays D 0 → K 0 S π 0 /η in each signal decay. The resulting fits for D 0 → π 0 π 0 π 0 , π 0 π 0 η and π 0 ηη are shown in Figs. 2 (a), (b) and (c), respectively. The signal yields and statistical significances, which are estimated from the likelihood difference between the fits with and without the signal included after considering the change in the number of degrees of freedom, are summarized in Table 2. Since no obvious D 0 → ηηη signal is observed, an upper limit on its decay BF is determined. We fit the M BC distribution of the D 0 → ηηη candidate events, where the signal is described by the MC simulated shape convoluted with a Gaussian function and the background by an ARGUS function. The parameters of the Gaussian function are fixed to those obtained in the fit of D 0 → π 0 ηη decay. The resultant best fit is shown in Fig. 2 (d). The PDF for the expected signal yield is taken to be the normalized likelihood L versus the BF in the fit, incorporating the systematic uncertainties as described below, and is shown as the inset plot in Fig. 2 (d). The upper limit on the BF at the 90% C.L., corresponding to The detection efficiencies for various decays of interest must take into account the effect of any intermediate states. The existence of intermediate states in the D 0 three-body decays is investigated by examining the corresponding Dalitz plots. Except for the decay D 0 → π 0 ηη, no obvious intermediate states are observed. Therefore, the detection efficiencies for the decays D 0 → π 0 π 0 π 0 , π 0 π 0 η and ηηη are obtained with MC samples of three-body phase space decay with uniform angular distributions.
For the decay D 0 → π 0 ηη, the a 0 (980) 0 is evident in the π 0 η invariant mass M π 0 η distribution. Figure 3 shows the M π 0 η spectrum of 23 events with two entries per event from the data sample with additional requirements −0.023 < ∆E < 0.020 GeV and 1.859 < M BC < 1.871 GeV/c 2 . An unbinned maximum likelihood fit is performed on the M π 0 η distribution to determine the a 0 (980) 0 signal yield.
In the fit, the shape of the a 0 (980) 0 is described with the shape from the MC sample of D 0 → a 0 (980) 0 η → π 0 ηη, which has two components: one with the π 0 combined with the correct η coming from the a 0 (980) 0 decay, and the other with the π 0 combined with the wrong η coming directly from the D 0 decay. The first peaks around the a 0 (980) 0 mass, while the second contributes a broad shape in the M π 0 η distribution. The MC shape is convolved with a Gaussian function to account for the mass resolution difference between data and MC simulation. In the MC simulation, the intermediate a 0 (980) 0 state is parameterized with the Flatté formula [16] with the central mass and the a 0 (980) 0 coupling constants coming from the Crystal Barrel experiment [17,18]. The component from the direct D 0 three-body decay is included in the fit, and its shape is the MC simulated shape, which is similar to that of the wrong η contribution in the a 0 (980) 0 shape. We also include the background in the fit, where its shape is determined from the inclusive MC sample. Both magnitudes for the D 0 three-body decay component and background are left free in the fit. The fit curves are shown in Fig. 3. The fit yields are 21 ± 5 events for the a 0 (980) 0 signal and 0 ± 4 events for the D 0 direct three-body decay, which  [9]. The first and second uncertainties are statistical and systematic, respectively. The upper limit is set at the 90% C.L..
implies the predominant process in the three-body decay of D 0 → π 0 ηη is D 0 → a 0 (980) 0 η. We also perform a fit without the a 0 (980) 0 signal included, and the statistical significance of the a 0 (980) 0 signal is calculated with the change of likelihood value with respect to that of the nominal fit taking into account the change of number of freedom in the fit. The significance for the a 0 (980) 0 signal is only 2.6σ, although it is the predominant component in the three-body decay. Therefore, in the decay of D 0 → π 0 ηη, the detection efficiency is estimated with the MC sample of D 0 → a 0 (980) 0 η → π 0 ηη as described above.
The resultant DT efficiencies for various decays are listed in Table 1. The BFs of these decays are calculated with Eq. 1, and summarized in Table 2.

Systematic Uncertainties
With the DT technique, the BF measurements are insensitive to systematics coming from the ST side since they mostly cancel. For the signal side, systematic uncertainties come mainly from the π 0 (η) reconstruction efficiency, ∆E resolution, CF background veto, χ 2 requirement, M BC fit, MC model, MC statistics, BFs of π 0 and η decays, and strong phase correction.
The π 0 reconstruction efficiency, including the photon detection efficiency, is studied as a function of π 0 momentum using a control sample of D 0 → K − π + π 0 events. The difference of the π 0 reconstruction efficiencies between data and MC simulation is regarded as the uncertainty related to π 0 reconstruction. We assume that the uncertainty due to reconstruction of the η is the same as that for the π 0 . The momentum weighted uncertainties of π 0 (η) reconstruction efficiencies are taken as the associated systematic uncertainties and are listed in Table 3 for each decay.
Uncertainty in the ∆E resolution is studied by widening the ∆E requirement from 3 to 3.5 times the resolution around the ∆E peak. For each decay, the resultant change of the BF is taken as the systematic uncertainty.
To estimate the uncertainty due to the K 0 S veto for D 0 → π 0 π 0 π 0 and π 0 π 0 η decays, the measurement is repeated with an alternative K 0 S mass window rejection region of 450 < M π 0 π 0 < 530 MeV/c 2 , which is enlarged from 4.5σ to 5.0σ of the resolution. The change of the BF for each decay is taken as the relevant systematic uncertainty.
The uncertainty arising from the χ 2 4π(ABC) requirements is investigated by repeating the measurement with an alternative requirement χ 2 4π(ABC) < 25. The resultant difference of the BF is taken as the corresponding systematic uncertainty for each decay.
Several aspects are considered to estimate the uncertainty related to the M BC fit. To examine the uncertainty in the fit range, a fit with an alternative range of (1.835, 1.890) GeV/c 2 is performed. The uncertainty of the signal shape is examined with an alternative fit, in which a Crystal Ball function is used to model the D 0 signal. Due to the long lifetime of K 0 S , the photons from π 0 (which are from K 0 S decay) decay do not originate from the IP. To study the uncertainty due to the imperfect simulation of the photon production vertex and its abnormal incidence into the EMC, an alternative MC sample, in which the K 0 S lifetime is set to zero, is used to determine the magnitude of BKG I. The uncertainty in BKG II is investigated with an alternative MC sample of D 0 → π 0 π 0 π 0 π 0 generated as phase space decay. The uncertainty from BKG III is checked by varying its magnitude by one standard deviation in the fit. The uncertainty from BKG IV is investigated by replacing the ARGUS function with the inclusive MC simulated background shape. For each of these sources, the resultant difference of the signal yield is treated as the corresponding systematic uncertainty for each decay. The total uncertainty associated with the M BC fit is the quadratic sum of the above individual values.
The uncertainty in the MC model is examined by analyzing the alternative MC events with and without involving the resonances f 0 (980) and a 0 (980) 0 . The maximum change in the detection efficiency is taken as the systematic uncertainty. For the decay D 0 → π 0 π 0 π 0 , the MC sample with f 0 (980) intermediate state, D 0 → f 0 (980)π 0 → π 0 π 0 π 0 , is selected. For the decay D 0 → π 0 π 0 η, the MC samples with f 0 (980) or a 0 (980) 0 intermediate states, D 0 → f 0 (980)η → π 0 π 0 η or D 0 → a 0 (980) 0 π 0 → π 0 π 0 η, are chosen. For the decay D 0 → π 0 ηη, the MC sample of the direct phase space decay is used. As for the decay D 0 → ηηη, no uncertainty in the MC model is assigned due to the relatively small phase space.
The uncertainty on the efficiency due to limited MC statistics is determined by ǫ (1 − ǫ)/N . Here, ǫ is the detection efficiency, and N is the number of the generated MC events. The uncertainties of the BFs for π 0 and η decays to two photons are taken from the PDG [9].
The uncertainty due to the quantum-correlation of the D 0D0 pair is considered via the strong phase factor. The absolute BF is calculated by B sig CP± = 1 1∓C f B sig , where B sig is calculated from Eq. 1, C f is the strong phase factor [19], which is (−12.4 ± 1.8)%, (−8.7 ± 1.6)% and (−7.0 ± 1.3)% for the ST mode ofD 0 → K + π − , K + π − π 0 and K + π − π − π + , respectively. The value of CP+ or CP− that determine the largest difference in BFs is used to give the systematic uncertainty.
Assuming all uncertainties, summarized in Table 3, are independent, the total uncertainties in the BF measurements are obtained by adding the individual uncertainties in quadrature.