Measurement of the decay $\eta^{\prime}\to\pi^{0}\pi^{0}\eta$ at MAMI

An experimental study of the $\eta'\to \pi^0\pi^0\eta \to 6\gamma$ decay has been conducted with the best up-to-date statistical accuracy, by measuring $\eta'$ mesons produced in the $\gamma p \to \eta' p$ reaction with the A2 tagged-photon facility at the Mainz Microtron, MAMI. The results obtained for the standard parametrization of the $\eta'\to \pi^0\pi^0\eta$ matrix element are consistent with the most recent results for $\eta'\to\pi\pi\eta$ decays, but have smaller uncertainties. The available statistics and experimental resolution allowed, for the first time, an observation of a structure below the $\pi^+\pi^-$ 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 $\pi\pi$ scattering length combination, $a_0-a_2$, extracted from $K \to 3\pi$ decays.

Another η → π 0 π 0 η test of QCD could also be made by comparing it to the isospin-violating η → 3π 0 decay, which gives access to the light quark masses and the mixing properties of π 0 and η mesons [15,16]. As an initial hypothesis [17], the amplitudes of η → π 0 π 0 π 0 and η → π 0 π 0 η can be related as is the π 0 -η mixing arXiv:1709.04230v2 [hep-ex] 16 Jun 2018 angle, m q is the mass of quark q, andm is the averaged mass of quarks u and d. Equation (1) assumes that the η → π 0 π 0 π 0 decay occurs entirely through π 0 -η mixing in the η → ηπ 0 π 0 decay. According to Ref. [18], such an assumption is too strong, but there should still be a nonnegligible contribution from ηπ → ππ rescattering, which can, e.g., be described through dispersion relations [19]. From the η → ππη 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 η , 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 η mesons act as a probe for testing the validity and effectiveness of the ChPT extensions and other theoretical models. The η → ππη 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 η → ππη decay [15].

A. Dalitz plot
The density of the η → π 1 π 2 η Dalitz plot, with π 1 π 2 being either the π + π − or π 0 π 0 pair, is completely described by the matrix element of the corresponding threebody decay amplitude. The η → π 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 η rest frame, and Q = T η + T π1 + T π2 = m η − 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 η as a consequence of the Bose-Einstein symmetry of the π 0 π 0 wave function, and for η → π + π − η as a consequence of charge-parity conservation. For the η → π 0 π 0 η Dalitz plot filled only with |X| 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 |1. − cX| and |1. + cX| 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 η 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 η → ππη 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 η 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 η mesons were produced in the charge-exchange reaction π − p → η n, with 5.4 · 10 3 η → π 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, η → π + π − η, was measured by the CLEO [30], 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 η → π + π − η decays, respectively. Compared to the other experiments, the measurement by CLEO [30] only determined 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 [31] for both the decay modes, which are based on the same data sample that was used to analyze the η → πππ decay [32], 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 η → ππη 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 η and η → π + π − η 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- Experimental results and theoretical calculations for the η → ππη matrix-element parameters, with statistical and systematic uncertainties for experimental data points plotted by the blue error bars and boxes, respectively. The η → π 0 π 0 η results from this work and by GAMS-4π [24] are shown with red boxes and the η → π + π − η 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 − NC and RChT for large-NC ChPT and Resonance Chiral Theory, respectively. The calculation from Ref. [29] is denoted as UGLσ.
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 η → ππη, 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 |a| ∼ |d| |b| ∼ |κ 21 | ∼ |κ 40 |, which was not confirmed by existing data, giving |a| ∼ |b| ∼ |d|. 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 η → ππη 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 [29]. The four nonets are two quark-antiquark and two four-quark 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 η and η → π + π − η decays, which is caused by strong final-state interactions, could especially be visible in the η → π 0 π 0 η spectrum below the π + π − threshold, where a pronounced cusp is expected. The contributions from the finalstate interactions for the η → ππη decays were calculated within the framework of non-relativistic effectivefield 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] (was 8% in the original work [1]). To check the reliability of such a prediction, a new measurement of η → π 0 π 0 η with good statistical accuracy and experimental resolution is needed.
The most recent dispersive analysis of the η → ππη 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 η amplitude based on an analysis of ∼ 1.2 · 10 5 η → π 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 Y X 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.

II. EXPERIMENTAL SETUP
The process γp → η p → π 0 π 0 ηp → 6γp was measured with the A2 experimental setup using the Crystal Ball (CB) [36] as a central calorimeter and TAPS [37,38] as a forward calorimeter. These detectors were installed in the energy-tagged bremsstrahlung photon beam of the Mainz Microtron (MAMI) [39,40]. The photon energies are determined by the Mainz end-point tagger (EPT) [41].
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 focal-plane detectors of the EPT magnetic spectrometer. It was especially built to conduct η 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 ±1.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 η threshold, E γ ≈ 1447 MeV.
Because the EPT experiments were mainly dedicated to studying η 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 η mass threshold.

III. EVENT SELECTION
The process γp → η 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 offline 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 pro-tons from the reaction γp → η 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 Timeof-Flight (TOF) method, based on the information from TAPS 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 least-squares 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 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 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 just by 3%.
To search for candidates to η → π 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, CL, (i.e., with a probability greater than 1%) were selected for further analysis. The pairing combination with the largest CL 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 CL or adding a constraint on the η mass.
The determination of the experimental acceptance for the process γp → η 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 (CL) 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 η → 6γ peak in the analysis of the experimental data and the corresponding MC simulation.
Because the experimental acceptance and the invariant-mass resolution depend on the η 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 → η 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 → η p differential cross sections obtained from the analysis of η → π 0 π 0 η → 6γ and η → γγ 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 η decay was generated close to the actual density distribution of the experimental plot, initially obtained by using the MC simulation with the η → π 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 finalstate 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 CL, 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 CL 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 ran-dom 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 CL > 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 CL was selected for further analysis. Similarly, for the γp → π 0 π 0 π 0 p → 6γp hypothesis, the combination with the best CL, from 15 possible pairings of the six photons to three π 0 mesons, was selected. The probability (or CL) distribution for experimental events selected as γp → π 0 π 0 ηp → 6γp with CL > 0.01 is shown in Fig. 3(a), after partially suppressing the γp → π 0 π 0 π 0 p → 6γp background by requiring CL(π 0 π 0 π 0 ) < 0.08 . This experimental distribution is described well by the sum of similar probability distributions for the MC simulations of γp → η 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 CL(π 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 η peak. Because fitting the γp → η p → π 0 π 0 ηp → 6γp hypothesis with an additional constraint on the η mass does not eliminate the background under the η peak, the experimental number of η mesons in every individual bin of various distributions was then determined by fitting the m(π 0 π 0 η) spectra corresponding to those bins.
This fitting procedure is illustrated in Fig. 4 for two different bins of the η → π 0 π 0 η Dalitz plot. To determine the number of measured η mesons, the experimental m(π 0 π 0 η) distribution was fitted with the sum of the signal line shape, fixed from a high-statistics γp → η p → π 0 π 0 ηp → 6γp MC sample, and a polynomial of order 4 for describing the background. The measured number of η mesons in a given bin was determined from the number of events in the corresponding MC distribution for η → π 0 π 0 η decays times the weight factor obtained from the fit. The uncertainty in the number of measured η 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 η 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 apllied 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 η decays was obtained by subtracting the background from the experimental spectra, based on the fit results (as shown in Fig. 4) for polynomial describing the background distribution.
The acceptance-corrected Dalitz plot can be obtained by correcting the measured number of η decays in each bin by the ratio of the reconstructed and generated events from the γp → η 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 η decays in each bin. Based on the analysis of the MC simulation for η → π 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 η 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 η mass. To determine the η → π 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 η 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 η Dalitz plot by using the η mass as a constraint in the kinematic-fit. The main advantage of using the η mass constraint is in the possibility of creating a sample of η → π 0 π 0 η decays, which allows fitting on an eventby-event basis, and in improving resolution in measured observables. Also, all events are reconstructed within the η → π 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 background could be suppressed by tightening selection criteria, this would also decrease the number of η → π 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 → η p differential cross sections by using η → π 0 π 0 η → 6γ and η → γγ decays with the method similar to Analysis I, i.e, by fitting the η signal above 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 η Dalitz plot, the γp → η 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 η mass) was also tested. The experimental sample of ∼ 1.23 · 10 5 η → π 0 π 0 η → 6γ decays, reconstructed by the fit with the four additional constraints, was then selected by requiring CL(γp → η p → π 0 π 0 ηp → 6γp) > 0.04 for the  3. (a) Probability (or CL) 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 → η 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 CL(γp → π 0 π 0 π 0 p → 6γp) < 0.08 . (b) Invariant-mass m(π 0 π 0 η) distributions for the same samples after passing an additional requirement CL(γp → π 0 π 0 ηp → 6γp) > 0.07, shown by the vertical line in (a). main hypothesis and CL(γ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 determines as 0.07 and 0.06, respectively. The background-estimation procedure of Analysis II is illustrated in Fig. 5. Invariant-mass distributions, m(6γ), for events reconstructed by testing γp → 6γp hypothesis, without any additional constraints on invariant masses, and selected by requiring the corresponding CL>1% are shown in Fig. 5(a) for the experimental data, the MC simulations of the signal channel γp → η 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 η peak is significantly larger than the η 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 CL>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 η → 6γ decays. The result of suppressing background contributions with CL(γp → π 0 π 0 ηp → 6γp) > 0.04 for the main hypothesis 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 → η 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 CL>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 CL>1%; (c) events reconstructed as in (b), but selected by requiring CL(γp → π 0 π 0 ηp → 6γp) > 4% and CL(γp → 3π 0 p → 6γp) < 0.75%; events for which CL(γp → η p → π 0 π 0 ηp → 6γp) < 4% from the fit with the fourth additional constraint on the η 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 CL(γp → η p → π 0 π 0 ηp → 6γp) > 4%, with the γp → η p → π 0 π 0 ηp → 6γp MC simulation shown by red points. and CL(γp → π 0 π 0 π 0 p) < 0.0075 for the background hypothesis is shown in Fig. 5(c). As seen, the background level under the η peak became sufficiently small not to cause much impact on the analysis of η → π 0 π 0 η decays. However, the background suppression resulted also in discarding 20% of good η → π 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 η mass and with the corresponding CL(γp → η 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 η mass. Then the events from Fig. 5(c) for which CL(γp → η 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 η → 6γ decays. The events shown in Fig. 5(d) also represent the final experimental and MC samples of η → π 0 π 0 η → 6γ decays, but with the kinematics reconstructed from the results of testing the γp → η p → π 0 π 0 ηp → 6γp hypothesis, involving the fourth additional constraint on the η 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 η 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 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 phasespace MC simulation of η → π 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 η 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 η 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 η 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 background-free 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 statisti- cal uncertainties. The magnitudes of these systematic uncertainties were determined from the change in the ratios depending on the kinematic-fit CL 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 η phase-space behavior and cannot mimic any narrow structures when divided by the corresponding phase-space spectra from the η 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 to this paper.

IV. RESULTS AND DISCUSSION
In Analysis I, to determine the η → π 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 η mesons, i is the corresponding detection efficiency, and f i (X, Y ) ∼ |M i (X, Y )| 2 is the theoretical function used in the fit and integrated over the given bin. The uncertainty σ i includes both the uncertainties 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 η 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 η signal in every bin and to the experimental resolution, which smeared η 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 CL selection criteria, the Dalitz-plot 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 into the first category. Test #1 accumulates the tests made to check the procedure that measures the η 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 η decays was obtained by subtracting the background from the experimental spectra, based on the fit results for polynomial of order 4, describing the background contribution. In test #2, the bin width in the m(π 0 π 0 η) spectra, used to measure the η signal, varied from 2 to 8 MeV/c 2 . In test #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 matrixelement 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 [47] between the parameter values obtained from fitting to the final, x f ± σ f , and test, x t ± σ 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 #4), and with the bins enlarged to 0.15 (test #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 #6 involved changes in the selection criteria based on the kinematic-fit CL for the γp → π 0 π 0 ηp → 6γp and γp → π 0 π 0 π 0 p → 6γp hypotheses. Variation of CL 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 #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 η matrixelement parameters: For convenience, these results are also plotted in Fig. 1, which compares the existing results and calculations for the standard matrix-element parametrization. For similarity with previous measurements, the correlation matrix is provided for fit #1, the results of which were taken as the main results in (8): 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 η 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 kinematicfit CL for both the γp → η p → π 0 π 0 ηp → 6γ and γ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 η decays with poorer resolution by tightening the CL(γp → η 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 η Dalitz plot are listed in Table I (#2), demonstrating good agreement with Analysis I (#1) within the uncertainties. Similarly 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 #2 and, in such 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 acceptancecorrected 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 #2, but still is in good agreement with it and fit #1 within the uncertainties.
The present results of Analysis I (shown also in Fig. 1) and II are consistent with the previous η → π 0 π 0 η measurement by GAMS-4π, a = −0.067 (16) pared 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,29]. 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 [29]. 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 [48]. 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 Y X 2 and κ 40 X 4 , the magnitudes of which were evaluated in Ref. [15]. The fit results after adding only the Y X 2 term are listed in Table I (#3, 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 η → ππη data with much higher statistical accuracy are needed for more reliable estimation of possible Y X 2 and X 4 contributions.
As discussed in Section I, a negative value of parameter b excludes a linear Y dependence of the η → ππη decay amplitude, with the matrix element parametrized according to Eq. (4). To allow a comparison to earlier measurements [24,26,30], this parametrization was also tested with the present data. The numerical results for Analyses I and II are listed in Table I (#5, 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 η 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 η → ππη decay amplitude parametrized within the NREFT framework [1]. In this model, the decay amplitude is decomposed up to two loops, A(η → ππη) = 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 twoloop 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 η → ππη Lagrangian by matching the standard Dalitz-plot parametrization, |M (X, Y )| 2 ∼ 1 + a Y + b Y 2 + d X 2 + ..., the tree amplitude can be parametrized as where a , b , and d are the tree-amplitude parameters that describe the η → ππη dynamics before the contributions from final-state interactions. Those parameters are also involved in the calculation of the one-and twoloop 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 , as πη final-state interactions turned out to be very small and were neglected in the NREFT approach. Because the tree amplitude is the same for both the η → π 0 π 0 η and η → π + π − η 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 ± 0.0051, where a 0 = 0.220 ± 0.005 and a 2 = −0.0444 ± 0.0010 [1,49]. 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 a 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 (#7), and the fit is shown in Fig. 7 by the black dash-dotted lines. The results obtained here for the tree-amplitude 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 (#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 and b , reflecting the decay dynamics before final-state interaction, look now closer to the linear Y dependence, which is typically expected for the tree amplitude [15,27]. As expected, the magnitude of parameter d is very close to d from the fit with the standard parametrization. The parameters obtained for the tree amplitude in the fits with the NREFT amplitude (#7-#11 in Table I) could be compared to the corresponding calculations for the leading-order terms of the η → π 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 (#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 tree-amplitude parameters according to the fit made with fixed a 0 − a 2 (#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 (#9). It demonstrates good consistency with the known value a 0 − a 2 = 0.2644 ± 0.0051 [49], though with much larger uncertainties, compared to it.
The results of fitting the η → π 0 π 0 η Dalitz plot from Analysis II with the NREFT amplitude, which are also listed in Table I (#10, 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 #7 are also compared to the dispersive analysis of the η → ππη 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 agreement within the uncertainties with this most recent calculation.

V. SUMMARY AND CONCLUSIONS
An experimental study of the η → π 0 π 0 η → 6γ decay has been conducted with the best up-to-date statistical accuracy, by measuring η mesons produced in the γp → η 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 η 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 η → ππη decays, but have smaller uncertainties. It was tested that including higher-order terms does not improve the description of the η → π 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 to the paper.  Fig. 7(b), including the NREFT fit #7, but compared to the dispersive analysis of the η → ππη 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.
No. LRU-FTI-123-2014) and the MSE Program "Nauka" (Project No. 3.825.2014/K). We thank the undergradu-ate students of Mount Allison University and The George Washington University for their assistance.