Measurement of the decay 
η′→π0π0η
 at MAMI

An experimental study of the η0 → π0π0η → 6γ decay has been conducted with the best up-to-date 
statistical accuracy, by measuring η0 mesons produced in the γp → η0 
p reaction with the A2 tagged-photon 
facility at the Mainz Microtron, MAMI. The results obtained for the standard parametrization of the 
η0 → π0π0η matrix element are consistent with the most recent results for η0 → ππη decays, but have 
smaller uncertainties. The available statistics and experimental resolution allowed, for the first time, an 
observation of a structure below the πþπ− mass threshold, the magnitude and sign of which, checked within 
the framework of the nonrelativistic effective-field theory, demonstrated good agreement with the cusp that 
was predicted based on the ππ scattering length combination, a0 − a2, extracted from K → 3π decays.

From the η 0 → ππη decay, one could also learn about QCD-related models and, in particular, the low-energy effective-field theory of QCD, chiral perturbation theory (ChPT) [20]. In ChPT, the pseudoscalar singlet η 1 , associated with the physical state η 0 , is not included explicitly (but implicitly through low-energy constants) due to the Uð1Þ A anomaly, which also renders it massive in the chiral limit of massless quarks. To include the singlet, ChPT is extended by going to a large number of color charges (large-N C ) that extends the symmetry to Uð3Þ L × Uð3Þ R [20][21][22][23]. Here hadronic decays of η 0 mesons act as a probe for testing the validity and effectiveness of the ChPT extensions and other theoretical models. The η 0 → ππη decay is particularly suitable for these studies as final-state interactions between ππ and πη are mainly dominated by scalar contributions, whereas G-parity conservation suppresses vector resonances. Therefore, the properties of the lowest-lying scalar resonances, f 0 ð500Þ and a 0 ð980Þ, could, in principle, be studied through the η 0 → ππη decay [15].

A. Dalitz plot
The density of the η 0 → π 1 π 2 η Dalitz plot, with π 1 π 2 being either the π þ π − or the π 0 π 0 pair, is completely described by the matrix element of the corresponding three-body decay amplitude. The η 0 → π 1 π 2 η Dalitz plot is typically expressed in terms of variables X and Y, defined as The observables T π 1 , T π 2 , and T η denote the kinetic energies of the two final-state pions and η in the η 0 rest frame, and Q ¼ T η þ T π 1 þ T π 2 ¼ m η 0 − m η − 2m π . Typically, in recent measurements [24][25][26], the following parametrization: (with a, b, c, and d being real-valued parameters) was sufficient to describe the experimental Dalitz plots. From the theoretical point of view, the cX term should be zero for η 0 → π 0 π 0 η as a consequence of the Bose-Einstein symmetry of the π 0 π 0 wave function, and for η 0 → π þ π − η as a consequence of charge-parity conservation. For the η 0 → π 0 π 0 η Dalitz plot filled only with jXj values (as in Refs. [24,27]), the cX term makes no sense, as any linear dependence on X is canceled by adding points with the plot density of j1.0 − cXj and j1.0 þ cXj corresponding to coordinates −X and X, respectively. Another parametrization, which historically [28] was considered in previous experiments, is assuming a linear Y dependence of the η 0 decay amplitude with a complex-valued parameter α. Because Eq. (4) can be transformed into Eq. (3) via a ¼ 2ReðαÞ and b ¼ Re 2 ðαÞ þ Im 2 ðαÞ, the linear parametrization of Y constraints parameter b from being negative, which is in contradiction with all experimental results obtained for parameter b so far [24][25][26].

B. Previous measurements
The results of previous measurements for the standard parametrization of the η 0 → ππη matrix element with Eq. (3) are plotted in Fig. 1, along with the most recent theoretical calculations. For convenience, the results of this work are added in Fig. 1 as well. So far, the neutral decay mode η 0 → π 0 π 0 η was only measured with the GAMS spectrometers [24,27]. The first experimental data were taken with the GAMS-2000 spectrometer at the IHEP accelerator U-70 [27], where η 0 mesons were produced in the chargeexchange reaction π − p → η 0 n, with 5.4 × 10 3 η 0 → π 0 π 0 η → 6γ decays accumulated. A more comprehensive study of this decay was later made with the upgraded GAMS-2000 spectrometer, GAMS-4π, which accumulated 1.5 × 10 4 decays [24]. The decay with charged pions, η 0 → π þ π − η, was measured by the CLEO [29], VES [25], and BESIII [26] Collaborations. For their analyses, CLEO, VES, and BESIII accumulated 6.7 × 10 3 , 20.1 × 10 3 , and 43.8 × 10 3 η 0 → π þ π − η decays, respectively. Compared to the other experiments, the measurement by CLEO [29] determined only the linear parametrization in Eq. (4). Comparison of the recent results from BESIII [26] and GAMS-4π [24] indicates that the isospin limit is a good approximation. The latest, high-statistics results from BESIII [30] for both the decay modes, which are based on the same data sample that was used to analyze the η 0 → πππ decay [31], appeared after the submission of this work and will, therefore, not be considered in the present paper.

C. Theoretical calculations
There are a few recent theoretical studies of η 0 → ππη decays. One of the studies was made within the framework of Uð3Þ chiral effective-field theory in combination with a relativistic coupled-channels approach [18,33,34], using the VES results and parameter a from GAMS-4π in their fitting procedure. The results of Ref. [18] are shown in Fig. 1 as Uð3Þ n and Uð3Þ c for the η 0 → π 0 π 0 η and η 0 → π þ π − η decays, respectively.
Another study was conducted within the large-N C ChPT, at lowest and next-to-leading orders, and resonance chiral theory (RChT) in the leading 1=N C approximation, with higher-order effects, such as ππ final-state interactions, being taken into account through a detailed unitarization procedure [15]. In the isospin limit, the intrinsic calculations for the charged and neutral decay amplitudes are the same. For the unitarization procedure that works within large-N C ChPT and RChT, experimental results are needed to fix particular theoretical constants. If one takes parameter a as its average from Refs. [24,25] and the uncertainty in a as the absolute difference in the two results, with their own uncertainties added, then the large-N C ChPT predicts b ¼ −0.050ð1Þ and d ¼ −0.092ð8Þ, and RChT gives b ¼ −0.033ð1Þ and d ¼ −0.072ð8Þ. The latter predictions are in slightly better agreement with the BESIII and GAMS-4π results, compared to the calculations with the Uð3Þ chiral effective-field theory [18]. In addition to the standard parametrization of η 0 → ππη, an expanded parametrization that includes two additional higher-order terms was also checked in Ref. [15]: where the notation of parameters κ 21 and κ 40 is left as in Ref. [15]. Such a parametrization was motivated by the fact that, in the chiral limit, expansions of Y should be of similar magnitude to those of order X 2 , with the expected hierarchy of the parameters to be jaj ∼ jdj ≫ jbj ∼ jκ 21 j ∼ jκ 40 j, which was not confirmed by existing data, giving jaj ∼ jbj ∼ jdj. For these additional terms, large-N C ChPT predicts κ 21 ¼ 0.003ð2Þ and κ 40 ¼ 0.002ð1Þ, and RChT gives κ 21 ¼ −0.009ð2Þ and κ 40 ¼ 0.001ð1Þ.
One more calculation of the η 0 → ππη decay, which includes mixing between scalar mesons, was studied within a generalized linear sigma model containing two nonets of scalar mesons and two nonets of pseudoscalar mesons [32]. The four nonets are two quark-antiquark and two fourquark states, respectively. For this calculation, denoted as UGLσ in Fig. 1, agreement with the experimental data was obtained only for parameter a, but not for b or d.
The main difference between the η 0 → π 0 π 0 η and η 0 → π þ π − η decays, which is caused by strong final-state interactions, could especially be visible in the η 0 → π 0 π 0 η spectrum below the π þ π − threshold, where a FIG. 1. Experimental results and theoretical calculations for the η 0 → ππη matrix-element parameters, with statistical and systematic uncertainties for experimental data points plotted by the blue error bars and boxes, respectively. The η 0 → π 0 π 0 η results from this work and by GAMS-4π [24] are shown with red boxes and the η 0 → π þ π − η results by VES [25] and BESIII [26] with green. The theoretical calculations are shown by black crosses. For Ref. [18], they are denoted as Uð3Þ n and Uð3Þ c for the neutral and charged decay mode, respectively. The predictions from Ref. [15] are denoted as L − N C and RChT for large-N C ChPT and resonance chiral theory, respectively. The calculation from Ref. [32] is denoted as UGLσ.
pronounced cusp is expected. The contributions from the final-state interactions for the η 0 → ππη decays were calculated within the framework of nonrelativistic effective-field theory (NREFT) [1], which was developed and successfully used for extracting ππ scattering lengths from K → 3π decays [5,7,8]. Based on the previously determined scattering lengths, the expected magnitude of the cusp is 6% [35] (it was 8% in the original work [1]). To check the reliability of such a prediction, a new measurement of η 0 → π 0 π 0 η with good statistical accuracy and experimental resolution is needed.
The most recent dispersive analysis of the η 0 → ππη decay amplitude, which is based on the fundamental principles of analyticity and unitarity, is presented in Ref. [19]. In the analysis framework, final-state interactions are fully taken into account, and the dispersive representation relies only on input for the ππ and πη scattering phase shifts. Because the dispersion relation contains subtraction constants that cannot be fixed by unitarity, these parameters were determined by fitting existing Dalitz-plot data. The prediction of a low-energy theorem was studied and the dispersive fit was compared to variants of ChPT.
This work presents an experimental study of the η 0 → π 0 π 0 η amplitude based on an analysis of ∼1.2 × 10 5 η 0 → π 0 π 0 η decays. In addition to the determination of the standard matrix-element parameters a, b, and d, the sensitivity to higher-order terms κ 21 YX 2 and κ 40 X 4 , suggested in Ref. [15], is tested. The magnitude of the cusp effect at the π þ π − threshold is checked within the NREFT framework [1] by using ππ scattering lengths extracted from K → 3π decays. Acceptance-corrected Dalitz plots and ratios of the X, Y, mðπ 0 π 0 Þ, and mðπ 0 ηÞ distributions to phase space are also provided, which could be used in theoretical approaches that determine some of their parameters from fitting to experimental data.
The CB detector is a sphere consisting of 672 optically isolated NaI(Tl) crystals, shaped as truncated triangular pyramids, pointing toward the center of the sphere. The crystals are arranged in two hemispheres that cover 93% of 4π sr, sitting outside a central spherical cavity with a radius of 25 cm, which is designed to hold the target and inner detectors. In the A2 experiments from 2007, TAPS was initially arranged in a plane consisting of 384 BaF 2 counters of hexagonal cross section. It was installed 1.45 m downstream from the CB center covering the full azimuthal range for polar angles from 1°to 20°. Later on, 18 BaF 2 crystals, covering polar angles from 1°to 5°, were replaced with 72 PbWO 4 crystals. This allowed running with a high MAMI electron current without decreasing the TAPS efficiency due to the high count rate in the crystals close to the photon-beam line. More details on the calorimeters and their resolutions are given in Ref. [14] and references therein.
The target is surrounded by a particle identification (PID) detector [42] used to distinguish between charged and neutral particles. It is made of 24 scintillator bars (50 cm long, 4 mm thick) arranged as a cylinder with a radius of 12 cm. A general sketch of the A2 setup is shown in Fig. 2. A charged particle tracker, MWPC, also shown in this figure (consisting of two cylindrical multiwire proportional chambers inside each other), was not used in the present measurement.
The present measurement was conducted in 2014 with a 1604 MeV unpolarized electron beam from the Mainz Microtron, MAMI C [39,40]. Bremsstrahlung photons, produced by the beam electrons in a 10 μm Cu radiator and collimated by a 4 mm diameter Pb collimator, impinged on a 10 cm long liquid hydrogen (LH 2 ) target located in the center of the CB. The energies of the incident photons were analyzed from 1426 MeV up to 1577 MeV by detecting the postbremsstrahlung electrons in 47 focalplane detectors of the EPT magnetic spectrometer. It was especially built to conduct η 0 measurements by covering the low-energy range of postbremsstrahlung electrons. The uncertainty of the EPT in E γ due to the width of its focal-plane detectors was about AE1.6 MeV, with a similar value (i.e., ∼1.6 MeV) in the systematic uncertainty in E γ due to the EPT energy calibration. The energy calibration of the EPT is based only on the simulation of electron tracing, using measured magnetic-field maps. The correctness of this calibration, as well as its uncertainty, was checked by measuring the position of the η 0 threshold, E γ ≈ 1447 MeV.
Because the EPT experiments were mainly dedicated to studying η 0 physics, the experimental trigger required the total energy deposit in the CB to be greater than ∼540 MeV to suppress triggering on reactions below the η 0 mass threshold.

III. EVENT SELECTION
The process γp → η 0 p → π 0 π 0 ηp → 6γp was searched for in events with seven energy-deposit clusters detected in both the CB and TAPS within the trigger time window, assuming that one of the seven clusters was from the recoil proton. The off-line cluster algorithm is optimized for finding a group of adjacent crystals in which the energy is deposited by a single-photon electromagnetic (e=m) shower, but it also works well for recoil protons. In the present beam-energy range, the recoil protons from the reaction γp → η 0 p were produced within the polar-angle range covered solely by TAPS, which required only events with at least one cluster detected in TAPS to be considered. In case of more than one cluster in TAPS, the recoil proton could be identified by a time-of-flight (TOF) method, based on the information from TAPS time-to-digital converters (TDCs). In general, the selection of event candidates and the reconstruction of the reaction kinematics were based on the kinematic-fit technique, which uses a constrained leastsquares fit with constraints based on energy and momentum conservation [43]. Typically, such a technique is sufficient to identify which cluster is from the recoil proton. Because the LH 2 target was 10 cm long, the vertex coordinate z (along the beam line) was used as another kinematic-fit variable, which can be used as either its measured or its free parameter. This improves angular resolution as the kinematic fit then determines the recoil-particle angles by varying z instead of using the cluster angles, which are calculated by assuming z ¼ 0. In addition to seven-cluster events, the cases with eight clusters were also checked when the least energetic cluster in TAPS had energy below 40 MeV. For the process under the study, there are e=m showers leaking from the CB through its downstream tunnel to TAPS, and the events with small leakage still can be reconstructed by the kinematic fit, though with poorer resolution. The recoil proton cannot be eliminated with such an approach as their kinetic energy is much larger than 40 MeV. In the end, including eight-cluster events in the analysis gained the experimental statistic by just 3%.
To search for candidates to η 0 → π 0 π 0 η → 6γ decays, the γp → π 0 π 0 ηp → 6γp kinematic-fit hypothesis was checked by using the four main constraints (conservation of energy and three momentum projections) and three additional constraints on the invariant masses of two neutral pions and η decaying into two photons. The number of 45 possible combinations to pair six photons to form the three final-state mesons could be decreased by checking invariant masses of cluster pairs before testing the corresponding hypothesis. The number of 7 possible combinations to pick the proton cluster can be decreased to the number of the clusters detected in TAPS. The events for which at least one pairing combination satisfied the tested hypothesis at the 1% confidence level (C.L.) (i.e., with a probability greater than 1%) were selected for further analysis. The pairing combination with the largest C.L. was used to reconstruct the reaction kinematics. The combinatorial background from mispairing six photons to three mesons was found to be quite small, with further decreasing if tightening a selection criterion on the kinematic-fit C.L. or adding a constraint on the η 0 mass.
The determination of the experimental acceptance for the process γp → η 0 p → π 0 π 0 ηp → 6γp and the study of its possible background reactions were based on their Monte Carlo (MC) simulation. All MC events were propagated through a GEANT simulation (version 4.9.6) of the experimental setup. To reproduce the experimental resolutions, the GEANT output was additionally smeared, allowing the simulated and experimental data to be analyzed in the same way. The MC smearing was adjusted to match the experimental invariant-mass resolutions and kinematic-fit confidence level (C.L.) distributions. The simulated events were also checked whether they passed the trigger requirements.
The reaction γp → π 0 π 0 π 0 p → 6γp, in which π 0 mesons are produced via intermediate baryon states, was used for the energy calibration of the calorimeters and the determination of their energy resolution. This reaction was also used to determine the additional smearing that was needed to reproduce the experimental resolution in the analysis of the MC simulations. The final adjustment for the energy calibration and the experimental resolution was made by comparing the η 0 → π 0 π 0 η → 6γ peak in the analysis of the experimental data and the corresponding MC simulation.
Because the experimental acceptance and the invariantmass resolution depend on the η 0 production angle, systematic differences can appear if the simulated data do not reflect the actual production distributions. To diminish such systematic effects, the reaction γp → η 0 p was generated according to its actual yield and angular spectra measured as a function of energy in the same experiment [44]. In that measurement, the γp → η 0 p differential cross sections obtained from the analysis of η 0 → π 0 π 0 η → 6γ and η 0 → γγ decays were found to be in good agreement, which confirmed the data quality and the analysis reliability. To diminish the impact of the experimental resolution on the shape of the acceptance-corrected Dalitz plot, the detection efficiency for it was determined from the MC simulation in which the η 0 → π 0 π 0 η decay was generated close to the actual density distribution of the experimental plot, initially obtained by using the MC simulation with the η 0 → π 0 π 0 η decay generated as phase space.
The study of possible background reactions via their MC simulation showed that the process γp → 3π 0 p → 6γp can mimic γp → π 0 π 0 ηp → 6γp via the combinatorial mismatching of six photons to three π 0 , which can occur when the invariant mass mð6γÞ > 820 MeV. Similarly, the process γp → π 0 ηp → 4π 0 p → 8γp can mimic the signal channel when two of the eight final-state photons are not detected. The suppression of the γp → 3π 0 p → 6γp background was done by testing the corresponding kinematic-fit hypothesis and rejecting events based on its C.L., which also resulted in some losses of the actual γp → π 0 π 0 ηp → 6γp events. The suppression of the γp → π 0 ηp → 4π 0 p → 8γp background can be achieved solely by tightening the cut on the C.L. of the γp → π 0 π 0 ηp → 6γp hypothesis itself.
In addition to background from other reactions, there are two more background sources. The first comes from interactions of the beam photons with the target cell windows. This background can be studied by analyzing data samples with an empty target cell, i.e., containing no liquid hydrogen. The contribution from such events was found to be very small and was neglected in the further analysis. A second background is caused by random coincidences of the tagger counts with the experimental trigger. This background was subtracted directly from the experimental spectra by using event samples with only random tagger coincidences, which is a standard technique for experiments using a beam of tagged bremsstrahlung photons [14,45,46].
The event-selection procedure included tests of several kinematic-fit hypotheses. First, all events had to pass the γp → 6γp hypothesis with its C:L: > 0.01. As mentioned above, to test the γp → π 0 π 0 ηp → 6γp hypothesis, 45 combinations of pairing the six photons to two π 0 and one η meson were checked, and the combination with the best C.L. was selected for further analysis. Similarly, for the γp → π 0 π 0 π 0 p → 6γp hypothesis, the combination with the best C.L., from 15 possible pairings of the six photons to three π 0 mesons, was selected. The probability (or C.L.) distribution for experimental events selected as γp → π 0 π 0 ηp → 6γp with C:L: > 0.01 is shown in Fig. 3(a), after partially suppressing the γp → π 0 π 0 π 0 p → 6γp background by requiring C:L:ðπ 0 π 0 π 0 Þ < 0.08. This experimental distribution is described well by the sum of similar probability distributions for the MC simulations of γp → η 0 p → π 0 π 0 ηp → 6γp, γp → π 0 π 0 π 0 p → 6γp, and γp → π 0 ηp → 4π 0 p → 8γp, the events of which were selected with the same criteria. For further suppression of the background reactions, only events with C:L:ðπ 0 π 0 ηÞ > 0.07, shown by the vertical line in Fig. 3(a), were selected for the final analysis. The invariant-mass mðπ 0 π 0 ηÞ distributions for the selected experimental and MC events are depicted in Fig. 3(b), which shows the level of the background contributions remaining under the η 0 peak. Because fitting the γp → η 0 p → π 0 π 0 ηp → 6γp hypothesis with an additional constraint on the η 0 mass does not eliminate the background under the η 0 peak, the experimental number of η 0 mesons in every individual bin of various distributions was then determined by fitting the mðπ 0 π 0 ηÞ spectra corresponding to those bins. ) distributions from testing the γp → π 0 π 0 ηp → 6γp hypothesis for the experimental data (black points) and the sum of MC simulations for γp → η 0 p → π 0 π 0 ηp → 6γp, γp → π 0 π 0 π 0 p → 6γp, and γp → π 0 ηp → 4π 0 p → 8γp, shown by the red, blue, and green areas, respectively. All events also passed the requirement C:L:ðγp → π 0 π 0 π 0 p → 6γpÞ < 0.08. (b) Invariantmass mðπ 0 π 0 ηÞ distributions for the same samples after passing an additional requirement C:L:ðγp → π 0 π 0 ηp → 6γpÞ > 0.07, shown by the vertical line in (a). This fitting procedure is illustrated in Fig. 4 for two different bins of the η 0 → π 0 π 0 η Dalitz plot. To determine the number of measured η 0 mesons, the experimental mðπ 0 π 0 ηÞ distribution was fitted with the sum of the signal line shape, fixed from a high-statistics γp → η 0 p → π 0 π 0 ηp → 6γp MC sample, and a polynomial of order 4 for describing the background. The measured number of η 0 mesons in a given bin was determined from the number of events in the corresponding MC distribution for η 0 → π 0 π 0 η decays times the weight factor obtained from the fit. The uncertainty in the number of measured η 0 in a given bin was based on the uncertainties in the signal MC distribution, the uncertainty in its weight factor from the fit, and the polynomial-fit uncertainty. The total number of measured η 0 decays integrated over the entire Dalitz plot was obtained as ∼1.241ð4Þð8Þ × 10 5 , with the first uncertainty being determined by the experimental statistics and the second being the systematic uncertainty due to the fitting method. The latter uncertainty was estimated in two ways. First, the same fitting procedure was applied to the sum of MC simulations for the signal and background contributions, with the number of signal events taken as experimentally observed. Second, the number of η 0 decays was obtained by subtracting the background from the experimental spectra, based on the fit results (as shown in Fig. 4) for a polynomial describing the background distribution.
The acceptance-corrected Dalitz plot can be obtained by correcting the measured number of η 0 decays in each bin by the ratio of the reconstructed and generated events from the γp → η 0 p → π 0 π 0 ηp → 6γp MC simulation in the given bin. Figure 6(a) shows such a plot for the bin width 0.1 in both X and Y, which allows sufficient statistics for a reliable determination of the number of η 0 decays in each bin. Based on the analysis of the MC simulation for η 0 → π 0 π 0 η decays, the Dalitz-plot acceptance was found to be almost uniform (with the averaged efficiency of 24.5%); but with its boundaries smeared by the experimental resolution, the average of which for X and Y was determined as δX ¼ 0.09 and δY ¼ 0.076. Because the η 0 mass constraint was not used to obtain the background-free Dalitz plot, the reconstructed events can migrate outside the physical boundaries of the plot determined by the η 0 mass. To determine the η 0 → π 0 π 0 η matrix-element parameters, only bins that are fully inside the Dalitz plot's physical boundaries were considered. These bins contained, for the given bin size, 1.13 × 10 5 measured η 0 decays.
As a cross-check of the analysis discussed above (Analysis I), an independent analysis (Analysis II) of the same data was performed to measure the η 0 → π 0 π 0 η Dalitz plot by using the η 0 mass as a constraint in the kinematic fit. The main advantage of using the η 0 mass constraint is in the possibility of creating a sample of η 0 → π 0 π 0 η decays, which allows fitting on an event-by-event basis, and in improving resolution in measured observables. Also, all events are reconstructed within the η 0 → π 0 π 0 η physical boundaries. The main disadvantage of this approach is the background contributions from γp → π 0 π 0 π 0 p → 6γp and γp → π 0 ηp → 4π 0 p → 8γp remaining in the experimental data sample. Although such a background could be suppressed by tightening the selection criteria, this would also decrease the number of η 0 → π 0 π 0 η decays available for the analysis.
The main experimental details of Analysis II are given in Ref. [44]. That work included the measurement of the γp → η 0 p differential cross sections by using η 0 → π 0 π 0 η → 6γ and η 0 → γγ decays with the method similar to Analysis I, i.e., by fitting the η 0 signal above the background in bins of differential cross sections. Such an analysis tests the γp → π 0 π 0 ηp → 6γp kinematic-fit hypothesis, involving three additional constraints on the invariant masses of two neutral pions and η decaying into two photons. To measure the η 0 → π 0 π 0 η Dalitz plot, the γp → η 0 p → π 0 π 0 ηp → 6γp kinematic-fit hypothesis (with the fourth additional constraint on the invariant mass of the two neutral pions and η to equal the η 0 mass) was also tested. The experimental sample of ∼1.23 × 10 5 η 0 → π 0 π 0 η → 6γ decays, reconstructed by the fit with the four additional constraints, was then selected by requiring C:L:ðγp → η 0 p → π 0 π 0 ηp → 6γpÞ > 0.04 for the main hypothesis and C:L:ðγp → π 0 π 0 π 0 pÞ < 0.0075 for the background hypothesis (involving three additional constraints on the invariant masses of the three neutral pions). The resolution in X and Y for such selection criteria was determined as 0.07 and 0.06, respectively. The backgroundestimation procedure of Analysis II is illustrated in Fig. 5. Invariant-mass distributions, mð6γÞ, for events reconstructed by testing the γp → 6γp hypothesis, without any additional constraints on invariant masses, and selected by requiring the corresponding C:L: > 1%, are shown in Fig. 5(a) for the experimental data, the MC simulations of the signal channel γp → η 0 p → π 0 π 0 ηp → 6γp, and the background channel γp → 3π 0 p → 6γp combined with γp → π 0 ηp → 4π 0 p → 8γp. As seen, in the six-photon data sample, the background level under the η 0 peak is significantly larger than the η 0 signal. The events reconstructed by testing the γp → π 0 π 0 ηp → 6γp hypothesis (with the three additional invariant-mass constraints) and selected with the corresponding C:L: > 1% are shown in Fig. 5(b). The background level is much smaller here, but it is not sufficient yet for reliable analysis of η 0 → π 0 π 0 η → 6γ decays. The result of suppressing background contributions with C:L:ðγp → π 0 π 0 ηp → 6γpÞ > 0.04 for the main hypothesis and C:L:ðγp → π 0 π 0 π 0 pÞ < 0.0075 for the background hypothesis is shown in Fig. 5(c). As seen, the background level under the η 0 peak became sufficiently small not to cause much impact on the analysis of η 0 → π 0 π 0 η decays. However, the background suppression resulted also in discarding 20% of good η 0 → π 0 π 0 η → 6γ decays, preventing from further tightening background cuts. The same figure also includes the invariant-mass distributions for the events that did not pass the hypothesis also involving the fourth additional constraint on the η 0 mass and with the corresponding C:L:ðγp → η 0 p → π 0 π 0 ηp → 6γpÞ < 4% for them. Because of such a cut, the latter distributions have a dip in the region of the η 0 mass. Then the events from Fig. 5(c) for which C:L:ðγp → η 0 p → π 0 π 0 ηp → 6γpÞ > 4% are shown in Fig. 5(d), illustrating good agreement in the energy calibration and resolution for the experimental and MC η 0 → π 0 π 0 η → 6γ decays. The events shown in Fig. 5(d) also represent the final experimental and MC samples of η 0 → π 0 π 0 η → 6γ decays, but with the kinematics reconstructed from the results of testing the γp → η 0 p → π 0 π 0 ηp → 6γp hypothesis, involving the fourth additional constraint on the η 0 mass. The estimation of background remaining in the final experimental sample is demonstrated in Fig. 5(d) by the solid magenta line. The normalization of this background is based on the spectra shown in Fig. 5(c), where the combined MC simulation of the background reactions is normalized on the parts of the experimental spectrum that are away from the η 0 signal. For the case demonstrated in Fig. 5(d), the fraction of the remaining background is estimated as 5.5%. With a more conservative background  5. Experimental mð6γÞ and mðπ 0 π 0 η → 6γÞ invariant-mass distributions (black points) from Analysis II compared to various MC simulations, with the combined background from γp → π 0 π 0 π 0 p → 6γp and γp → π 0 ηp → 4π 0 p → 8γp shown by the solid magenta line and the γp → η 0 p → π 0 π 0 ηp → 6γp MC simulation shown in (a), (b), and (c) by the solid black line: (a) events reconstructed by testing the γp → 6γp hypothesis (no additional constraints on invariant masses of the final-state photons) and selected with the corresponding C:L: > 1%; (b) events reconstructed by testing the γp → π 0 π 0 ηp → 6γp hypothesis (with the three additional constraints on the photons' invariant masses) and selected with the corresponding C:L: > 1%; (c) events reconstructed as in (b), but selected by requiring C:L:ðγp → π 0 π 0 ηp → 6γpÞ > 4% and C:L:ðγp → 3π 0 p → 6γpÞ < 0.75%; events for which C:L:ðγp → η 0 p → π 0 π 0 ηp → 6γpÞ < 4% from the fit with the fourth additional constraint on the η 0 mass are shown by red points for experimental data and by the cyan line for the combined background; (d) events reconstructed and selected as in (c), but for which C:L:ðγp → η 0 p → π 0 π 0 ηp → 6γpÞ > 4%, with the γp → η 0 p → π 0 π 0 ηp → 6γp MC simulation shown by red points. normalization in Fig. 5(c), the background fraction in the final sample could be larger, reaching up to 7.5%.
Based on the MC simulation, the background events were found to be distributed randomly over the Dalitz plot, looking similar to the distribution from the phase-space MC simulation of η 0 → π 0 π 0 η → 6γ decays. Thus, such a behavior of the background events cannot result in mimicking any narrow structure as cusp, but the Dalitz-plot slopes, as well as a cusp structure, could look more shallow, introducing corresponding systematic effects in the  extracted matrix-element parameters. The acceptance-corrected η 0 → π 0 π 0 η Dalitz plot from Analysis II is shown in Fig. 6(b). Compared to Analysis I, the latter plot has smaller binning and includes bins that also cover the physical boundaries of the η 0 → π 0 π 0 η Dalitz plot. To include the boundary bins in a fit, the density function should be corrected for phase space available in those bins. Other informative distributions that illustrate the deviation of the actual η 0 → π 0 π 0 η decay from phase space are ratios of the experimental X, Y, mðπ 0 π 0 Þ, and mðπ 0 ηÞ spectra to their phase-space MC simulation. In Analysis I, those backgroundfree experimental spectra were obtained with a procedure similar to that used to measure the Dalitz plot itself and with the same selection criteria. The ratios obtained for each observable are shown in Fig. 7 by blue open squares. The results for larger masses in the mðπ 0 π 0 Þ and mðπ 0 ηÞ spectra are not included in Fig. 7, as the fitting procedure used to measure the experimental signal in those bins was giving large uncertainties in the results. The ratios from Analysis II are depicted in Fig. 7 by red open circles. Because the experimental sample includes the remaining background events, the systematic uncertainties that reflect the level of this background were added linearly to the statistical uncertainties. The magnitudes of these systematic uncertainties were determined from the change in the ratios depending on the kinematic-fit C.L. used for the event selection, under the assumption that such a change occurs solely due to a different level of the remaining background. The MC simulation of the two background reactions demonstrated that the spectra with the remaining background are close to the η 0 → π 0 π 0 η phasespace behavior and cannot mimic any narrow structures when divided by the corresponding phase-space spectra from the η 0 MC simulation. The normalization of the ratio distributions is based on the ratio in the number of events in the experimental and MC spectra. A smaller binning in Y and mðπ 0 π 0 Þ was chosen such that the expected cusp structure could be visible. Figure 7 shows that the data points from Analyses I and II are in good agreement within their uncertainties, but demonstrate a significant deviation from phase space. These points could be used, together with the Dalitz plots, for testing different models; therefore, all data points from the Dalitz plots and the four ratio distributions are provided as Supplemental Material [47] to this paper.

IV. RESULTS AND DISCUSSION
In Analysis I, to determine the η 0 → π 0 π 0 η matrixelement parameters, the Dalitz-plot fitting procedure was based on the minimization of where, for bin i, N i exp is the measured number of η 0 mesons, ϵ i is the corresponding detection efficiency, and f i ðX; YÞ ∼ jM i ðX; YÞj 2 is the theoretical function used in the fit and integrated over the given bin. The uncertainty σ i includes the uncertainties both in N i exp and in ϵ i . In Analysis II, the χ 2 calculation was based on the differences between the bin contents of the measured (i.e., uncorrected for the acceptance) Dalitz plot and the corresponding plot with the η 0 → π 0 π 0 η MC events weighted with the fit function. For every fit iteration, a weight of each MC event was calculated with the fit function taken from the generated values for X and Y and with its current parameters taken from the fit. Then the entry in the MC Dalitz plot was based on the reconstructed X and Y, thus taking the experimental acceptance and resolution into account.
The main results from fitting experimental Dalitz plots of Analyses I and II with different functions are listed in Table I. The first uncertainty in parameter values represents the errors from the fits. The systematic uncertainties in the matrix-element parameters were evaluated only for its standard parametrization with Eq. (3). The study of systematic effects was more scrupulous in Analysis I, which provides a near background-free Dalitz plot, but whose results could be sensitive to the fitting procedure of measuring the η 0 signal in every bin and to the experimental resolution, which smeared η 0 events out of the physical region. To estimate systematic uncertainties (given as the second uncertainty) in the matrix-element parameters, their sensitivity was tested to the changes in the procedure measuring the signal, the C.L. selection criteria, the Dalitzplot bin width, and the period of data taking.
Cross-checks of systematic effects in Analysis I, which are divided into two categories, are summarized in Table II. All tests that were made without changing the selection criteria for events in the final data sample were included in the first category. Test No. 1 accumulates the tests made to check the procedure that measures the η 0 signal. Those tests included changes in the order of the polynomial used to describe the background in each Dalitz-plot bin, by varying it from 2 to 5. Also, instead of using the signal line shape from its MC simulation, the number of η 0 decays was obtained by subtracting the background from the experimental spectra, based on the fit results for the polynomial of order 4, describing the background contribution. In test No. 2, the bin width in the mðπ 0 π 0 ηÞ spectra, used to measure the η 0 signal, varied from 2 to 8 MeV=c 2 . In test No. 3, the Dalitz-plot bin width varied from 0.05 to 0.15 in both X and Y. Because none of these tests reveals a deviation in parameter values more than the corresponding fit errors, they did not contribute to the systematic uncertainties of the matrix-element parameters.
The second category included all tests that were made by changing experimental statistics. The consistency of the results was checked by comparing the uncorrelated differences [48] between the parameter values obtained from fitting to the final, x f AE σ f , and test, x t AE σ t , Dalitz plots: with a systematic effect revealed if Δx uncor > 2.
Because the total data set was collected during three different periods of data taking, the consistency of these three subsets with each other was checked for the Dalitz plots with the main bin width, which was 0.1 in both X and Y (test No. 4), and with the bins enlarged to 0.15 (test No. 5), which decreased the corresponding statistical uncertainties closer to the level of the final Dalitz plot. The largest deviation, Δx uncor ¼ 2.0, was found for parameter a with the bin width 0.1, but it became smaller with the enlarged bin width. Based on the results of these tests, it was concluded that the data from the different periods of data taking are consistent with each other, and there are no systematic differences between them.
Test No. 6 involved changes in the selection criteria based on the kinematic-fit C.L. for the γp → π 0 π 0 ηp → 6γp and γp → π 0 π 0 π 0 p → 6γp hypotheses. Variation of C.L. values results in both a different level of background events in the experimental mðπ 0 π 0 ηÞ spectra and a change in the resolution of selected experimental and MC events. As Table II shows, the value of Δx uncor for each parameter exceeded the magnitude of 2, which was chosen to expose a systematic effect. The corresponding systematic uncertainties were then taken as half of the difference between the maximum and minimum parameter values obtained in test No. 6. Such an evaluation could be quite conservative as the tests giving parameters between their maximum and minimum values, as well as the parameter errors in individual fits, are neglected in this evaluation. Taking these systematic uncertainties into account resulted in the following values for the η 0 → π 0 π 0 η matrix-element parameters: For convenience, these results are also plotted in Fig. 1, which compares the existing results and calculations for the TABLE I. Results of fitting the η 0 → π 0 π 0 η Dalitz plots with different functions describing its density. The first uncertainty in parameter values represents their errors from the fits. The second error, which represents the systematic uncertainty, is evaluated for the fit with Eq. (3) only. The parameters held fixed during the fits with the NREFT amplitude [1] are marked as (fix).
Fit # jMðX; YÞj 2 ∼ 1 þ aY þ bY 2 þ dX 2 χ 2 =dof a b d In Analysis II, it was found that the systematic uncertainties in the parameter values are mostly caused by the background remaining in the selected η 0 → π 0 π 0 η experimental decays. Similar to Analysis I, their magnitudes were determined by comparing the results from fitting the experimental Dalitz plots that were obtained with looser and tighter selection criteria on the kinematic-fit C.L. for both the γp → η 0 p → π 0 π 0 ηp → 6γ and the γp → π 0 π 0 π 0 p → 6γp hypotheses. This assumes that the changes in the parameter values are solely caused by different fractions of the remaining background, rather than rejecting η 0 decays with poorer resolution by tightening the C:L:ðγp → η 0 p → π 0 π 0 ηp → 6γÞ criterion, in combination with lowering experimental statistics. The results of Analysis II for the standard parametrization of the η 0 → π 0 π 0 η Dalitz plot are listed in Table I (No. 2), demonstrating good agreement with Analysis I (No. 1) within the uncertainties. Similar to the results for the X, Y, mðπ 0 π 0 Þ, and mðπ 0 ηÞ spectra of Analysis II, the systematic uncertainties in the results for the standard parameters should be added linearly to the fit uncertainties. Because the results of Analysis II were obtained by fitting the measured plot, it was also checked that fitting the corresponding acceptance-corrected Dalitz plot, shown in Fig. 6(b), gives similar results. For the standard matrix-element parametrization, such a fit resulted in a ¼ −0.071ð7Þ, b ¼ −0.069ð11Þ, and d ¼ −0.061ð7Þ, with χ 2 =dof ¼ 1.085, demonstrating good agreement with fit No. 2 and, in such a way, confirming the reliability of the procedure used for the acceptance correction, which minimizes the smearing effect from the experimental resolution. For comparison, the fit to the acceptance-corrected Dalitz plot obtained with just a phase-space MC simulation results in a ¼ −0.073ð7Þ, b ¼ −0.061ð11Þ, and d ¼ −0.056ð8Þ, with χ 2 =dof ¼ 1.123, which deviates more from fit No. 2, but still is in good agreement with it and fit No. 1 within the uncertainties.
The present results of Analysis I (shown also in Fig. 1) and Analysis II are consistent with the previous η 0 → π 0 π 0 η measurement by GAMS-4π, a ¼ −0.067ð16Þð4Þ, b ¼ −0.064ð29Þð5Þ, and d ¼ −0.067ð20Þð3Þ [24], but improve upon their uncertainties. The agreement with the η 0 → π þ π − η results is better with the BESIII data [26], compared to the VES data [25]. At the same time, the present results do not improve much the situation existing between the experimental data and the calculations [15,18,32]. For instance, the results obtained for parameter a deviate from the prediction of the Uð3Þ chiral effective-field theory [18], but are consistent with the value from the generalized linear-sigma model [32]. The situation is opposite for parameters b and d, for which only the latter model cannot reproduce the experimental results. It is expected, however, that the agreement with the linear-sigma model could be improved if higher-order corrections were included in the calculations [49].
The quality of fitting the experimental Dalitz plots with its standard parametrization, based on Eq. (3), can be seen from both the χ 2 =degrees-of-freedom (χ 2 =dof) values listed in Table I and the fit results compared to the X, Y, mðπ 0 π 0 Þ, and mðπ 0 ηÞ spectra shown in Fig. 7. As seen in this figure, the Dalitz-plot fit results, depicted by the magenta solid line for Eq. (3), are in good agreement with the experimental data points, except for the region in the Y and mðπ 0 π 0 Þ spectra where the cusp structure was expected.
More fits were made to check the sensitivity of the measured Dalitz plots to the higher-order terms κ 21 YX 2 and κ 40 X 4 , the magnitudes of which were evaluated in Ref. [15]. The fit results after adding only the YX 2 term are listed in Table I (Nos. 3 and 4). They show practically no improvement in χ 2 =dof. The magnitude found for parameter κ 21 is somewhat larger than values predicted in Ref. [15] [with the closest prediction κ 21 ¼ −0.009ð2Þ], but its large uncertainty and correlation with other parameters cannot justify the results obtained. A similar situation was observed after adding the κ 40 X 4 term, which resulted in no improvement in χ 2 =dof and which had a strong correlation with the κ 21 X 2 term. It appears that η 0 → ππη data with much higher statistical accuracy are needed for a more reliable estimation of possible YX 2 and X 4 contributions.
As discussed in Sec. I, a negative value of parameter b excludes a linear Y dependence of the η 0 → ππη decay amplitude, with the matrix element parametrized according to Eq. (4). To allow a comparison to earlier measurements [24,26,29], this parametrization was also tested with the present data. The numerical results for Analyses I and II are listed in Table I (Nos. 5 and 6), and the fit from Analysis I is shown in Fig. 7 by the green dashed lines. As seen, the χ 2 =dof value becomes worse compared to the standard parametrization, and both the Y distribution and the related mðπ 0 π 0 Þ spectrum cannot be described with such a linear dependence. Nevertheless, both Analyses I and II give very similar results for ReðαÞ and ImðαÞ, which are also consistent with the most recent η 0 → π 0 π 0 η results from GAMS-4π: ReðαÞ ¼ −0.042ð8Þ and ImðαÞ ¼ 0.00ð7Þ [24].
As also seen in Fig. 7, both Analyses I and II indicate a cusp at the π þ π − threshold in the Y and mðπ 0 π 0 Þ spectra (marked by the vertical dashed lines), and none of the parametrizations tested so far were able to describe this region. To investigate this effect, the experimental Dalitz plots were fitted with the η 0 → ππη decay amplitude parametrized within the NREFT framework [1]. In this model, the decay amplitude is decomposed up to two loops, Aðη 0 → ππηÞ ¼ A tree þ A 1−loop þ A 2−loop , with the tree amplitude complemented by final-state interactions of one and two loops. The one-and two-loop amplitudes that describe the final-state interactions are calculated in Ref. [1] based on the magnitudes of ππ scattering lengths, which were previously extracted from the analysis of K → 3π decays [5,7,8]. In the same work, there are no calculations for the tree amplitude. To obtain couplings for the η 0 → ππη Lagrangian by matching the standard Dalitz-plot parametrization, jMðX;YÞj 2 ∼1þa 0 Y þb 0 Y 2 þ d 0 X 2 þÁÁÁ, the tree amplitude can be parametrized as where a 0 , b 0 , and d 0 are the tree-amplitude parameters that describe the η 0 → ππη dynamics before the contributions from final-state interactions. Those parameters are also involved in the calculation of the one-and two-loop amplitudes. Then the total amplitude determines the final dependence on Y and mðπ 0 π 0 Þ. At the same time, the dependence on X and mðπ 0 ηÞ is still mostly defined by the tree amplitude and its parameter d 0 , as πη final-state interactions turned out to be very small and were neglected in the NREFTapproach. Because the tree amplitude is the same for both the η 0 → π 0 π 0 η and η 0 → π þ π − η decays, the difference in them is determined solely by the final-state interactions, which also produce the cusp structure in the spectra from the neutral decay mode. The magnitude and the sign of this cusp structure is mostly determined by the scattering length combination a 0 − a 2 ¼ 0.2644 AE 0.0051, where a 0 ¼ 0.220 AE 0.005 and a 2 ¼ −0.0444 AE 0.0010 [1,50]. Because, in Analysis I with the Dalitz-plot bin width 0.1, the cusp region was mostly located in the boundary bins that were rejected from the fits, the experimental plot was remeasured with a narrower bin width, 0.05, in both X and Y, which allowed better access to the cusp. The results of fitting this Dalitz plot with the NREFT amplitude, having only three free parameters from the tree term, are listed in Table I (No. 7), and the fit is shown in Fig. 7 by the black dash-dotted lines. The results obtained here for the treeamplitude parameters are strongly model dependent and can be used only for qualitative purposes. From the comparison with the fit results based on the standard parametrization (No. 1), the χ 2 =dof value is of the same magnitude, but the NREFT fit describes better the Y and mðπ 0 π 0 Þ data points below the π þ π − threshold. The magnitude of parameters a 0 and b 0 , reflecting the decay dynamics before final-state interaction, now look closer to the linear Y dependence, which is typically expected for the tree amplitude [15,27]. As expected, the magnitude of parameter d 0 is very close to d from the fit with the standard parametrization. The parameters obtained for the tree amplitude in the fits with the NREFTamplitude (Nos. 7-11 in Table I) could be compared to the corresponding calculations for the leading-order terms of the η 0 → π 0 π 0 η decay amplitude, which are, e.g., provided in Ref. [15].
To test more reliably if the structure seen below the π þ π − threshold is in agreement with the magnitude of the cusp predicted by NREFT, the difference a 0 − a 2 was implemented in the NREFT code as its free parameter, and the fit repeated with the four parameters. The results of this fit are listed in Table I (No. 8). As seen, the four-parameter fit does not improve the χ 2 =dof value. Furthermore, such a fit is not well justified as the overall Y and mðπ 0 π 0 Þ distributions in NREFT depend on the individual values of a 0 and a 2 as well. A fairer procedure would be to fix the three treeamplitude parameters according to the fit made with fixed a 0 − a 2 (No. 7), and then to release only this difference to test solely the cusp region. The results for the latter fit are listed in Table I (No. 9). It demonstrates good consistency with the known value a 0 − a 2 ¼ 0.2644 AE 0.0051 [50], though with much larger uncertainties, compared to it.
The results of fitting the η 0 → π 0 π 0 η Dalitz plot from Analysis II with the NREFT amplitude, which are also listed in Table I (Nos. 10 and 11), confirm the corresponding results obtained from Analysis I. All fit results obtained with the NREFT amplitude speak for the quality of this model, together with the previously extracted scatteringlength combinations.
In Fig. 8, the data points for the Y dependence and the corresponding NREFT fit No. 7 are also compared to the dispersive analysis of the η 0 → ππη decay amplitude from Ref. [19]. Namely, their prediction based on fitting the BESIII data [26] with three subtraction constants is shown, including two error bands representing different uncertainties: one from fitting to the data and the other from the variation of the phase input. As the figure shows, the present data points are in good  Fig. 7(b), including the NREFT fit No. 7, but compared to the dispersive analysis of the η 0 → ππη decay amplitude from Ref. [19], the prediction of which is shown for the case of fitting to the BESIII data [26] with three subtraction constants. The light-gray error band represents the uncertainties obtained from fitting to the data, and the dark-gray band from the variation of the phase input. agreement within the uncertainties with this most recent calculation.

V. SUMMARY AND CONCLUSIONS
An experimental study of the η 0 → π 0 π 0 η → 6γ decay has been conducted with the best up-to-date statistical accuracy, by measuring η 0 mesons produced in the γp → η 0 p reaction with the A2 tagged-photon facility at the Mainz Microtron, MAMI. The results of this work obtained for the standard parametrization of the η 0 → π 0 π 0 η matrix element, a ¼ −0.074ð8Þ stat ð6Þ syst , b ¼ −0.063ð14Þ stat ð5Þ syst , and d ¼ −0.050ð9Þ stat ð5Þ syst are consistent with the most recent results for η 0 → ππη decays, but have smaller uncertainties. It was tested that including higher-order terms does not improve the description of the η 0 → π 0 π 0 η Dalitz plot. The available statistics and experimental resolution allowed, for the first time, an observation of a structure below the π þ π − mass threshold, the magnitude and sign of which, checked within the framework of the nonrelativistic effective-field theory, demonstrated good agreement (within a one-σ level) with the cusp that was predicted based on the ππ scattering length combination, a 0 − a 2 , extracted from K → 3π decays. The data points from the experimental Dalitz plots and ratios of the X, Y, mðπ 0 π 0 Þ, and mðπ 0 ηÞ distributions to phase space are provided as Supplemental Material [47] to the paper.