High-statistics measurement of the eta->3pi^0 decay at the Mainz Microtron

The largest, at the moment, statistics of 7x10^6 eta->3pi^0 decays, based on 6.2x10^7 eta mesons produced in the gamma p ->eta p reaction, has been accumulated by the A2 Collaboration at the Mainz Microtron, MAMI. It allowed a detailed study of the eta->3pi^0 dynamics beyond its conventional parametrization with just the quadratic slope parameter alpha and enabled, for the first time, a measurement of the second-order term and a better understanding of the cusp structure in the neutral decay. The present data are also compared to recent theoretical calculations that predict a nonlinear dependence along the quadratic distance from the Dalitz-plot center.


I. INTRODUCTION
For decades, the η → 3π decay has attracted much attention from theoretical and experimental studies as it gives access to fundamental physical constants. This decay, which is forbidden by isospin symmetry, mostly occurs due to the difference in the mass of the u and d quarks, with Γ(η → 3π) ∼ (m d − m u ) 2 [1]. Therefore, a precision measurement of this decay can be used as a sensitive test for the magnitude of isospin breaking in the * Electronic address: prakhov@ucla.edu Quantum Chromodynamics (QCD) part of the Standard Model (SM) Lagrangian. At the same time, the actual η → 3π dynamics involve a strong impact from ππ finalstate interactions, and the m d −m u magnitude cannot be approached without a precise experimental measurement of the η → 3π Dalitz plots, the density of which provides the information needed. Theoretical calculations of strong-interaction processes at low energy, which could typically be performed by using Chiral Perturbation Theory (χPTh) [1][2][3], were not very successful at describing the η → 3π density distributions observed experimentally. The main reason was in the final-state rescattering effects, the calculation of which turned out to be more reliable with dispersion relations [4,5], but still insuf-ficient to describe the experimental data. Meanwhile, the experimental progress in both the precise determination of the ππ phase shifts [6][7][8] and high-statistics data on the η → 3π 0 and η → π + π − π 0 decays [9][10][11][12][13][14] renewed the interest in theoretical studies of the η → 3π decay [15][16][17][18][19][20][21][22], which also included the extraction of the quark-mass ratio, Q 2 = (m 2 s −m 2 ud )/(m 2 d − m 2 u ) with m ud = (m u + m d )/2, from the data.
Due to the low energies of the decay pions, π 0 π 0 rescattering in η → 3π 0 is expected to be dominated by S waves. Such an assumption leads to the conventional leading-order parametrization |A(z)| 2 ∼ 1 + 2αz [23] of the η → 3π 0 amplitude, with only the quadratic slope parameter α, which was used in all previous measurements. Rather than fitting two-dimensional Dalitz plots, those measurements were based on the deviation of measured z distributions from the corresponding distributions obtained from the phase-space simulation of the η → 3π 0 decay, which is illustrated for both the Dalitz plot and z distribution in Fig. 1.
The result with the best accuracy, (α = −0.0322 ± 0.0012 stat ± 0.0022 syst , obtained by the A2 Collaboration at MAMI, was based on 3 × 10 6 observed η → 3π 0 decays [11]. Significant attention in that work was dedicated to a search for a possible cusp structure in the spectra below the π + π − threshold. Based on the ππ scattering length combination a 0 − a 2 , extracted from the analysis of K → 3π decays [32], and calculations within the framework of nonrelativistic effective field theory (NREFT) [15], the cusp effect was expected to be visible in the m(π 0 π 0 ) spectrum, reaching ∼ 1% at the 2π 0 threshold with respect to the spectrum in the case of the cusp absence. This calculation used the η → π + π − π 0 results from KLOE [12] to describe the charge-decay amplitude, assuming the isospin limit to connect it to the neutral decay. In principle, the predicted cusp magnitude should not change much even in the case of isospin breaking. However, the expected cusp structure was not confirmed experimentally in Ref. [11]. At the same time, the statistical accuracy of data points in the measured z distribution made it possible to indicate that the conventional leading-order parametrization |A(z)| 2 ∼ 1 + 2αz was not sufficient for the proper description of the η → 3π 0 decay amplitude. This indicates that the contributions from the higher-order terms in Eq. (3) need to be checked as well. The cusp structure cannot be described by polynomial expansion but, similar to the NREFT, the cusp range can be parametrized in the density function as ρ(s) = Re (1 − s/4m 2 π ± ), which results in ρ(s) = 0 for s ≥ 4m 2 π ± [33]. Then the density function is given by where the factor 2 in front of the cusp term is added for the consistency with the other terms. A better determination of the η → 3π decay parameters, needed for a precise determination of light-quark mass ratios, was recently the focus of many theoretical works. In Ref. [16], a detailed study of the η → 3π decays within the framework of the modified NREFT, in which final-state interactions were analyzed beyond one loop including isospin-breaking corrections, resulted in the extraction of the Dalitz-plot parameters for both the charged and neutral decays. The values obtained for the parametrization of the neutral decay with Eq. (3), α = −0.0246(49), β = −0.0042 (7), and γ = 0.0013(4), indicated nonzero contributions for the higher-order terms. Other η → 3π 0 calculations, involving parameter β, used a unitary dispersive model [18,19], in which substraction constants were fixed by fitting recent high-statistics η → π + π − π 0 data from WASA-at-COSY (1.74 × 10 5 decays) [13] and KLOE (4.7 × 10 6 decays) [14]. In contrast to Ref. [16], the latter calculations predicted a value of β consistent with zero. Another recent dispersive analysis [21] of the η → 3π decay amplitudes, in which the latest η → π + π − π 0 data from KLOE [14] were also fitted to determine subtraction constants, predicted a nonlinear z dependence for η → 3π 0 , which turned out to be in good agreement within the uncertainties with the measured z dependence from Ref. [11]. However, no numerical predictions were provided for the higher-order terms of Eq. (3). The most recent η → 3π calculation, which used the extended chiral Khuri-Treiman dispersive formalism [22], showed that the effect from the two light resonances f 0 (980) and a 0 (980) in the low energy region of the η → 3π decay is not negligible, especially for the neutral mode, and improves the description of the density variation over the Dalitz plot. The η → 3π 0 parameters obtained in Ref. [22] from their fitted amplitude, α = −0.0337 (12) and β = −0.0054(1), also predict nonzero contributions for the 2βz 3/2 sin(3φ) term.
Obviously, a better comparison of the experimental data with the recent η → 3π 0 calculations, going beyond the leading-order parametrization, should now be based on describing the two-dimensional density distribution of measured Dalitz plots, rather than on one-dimensional z distributions. To obtain reliable experimental results for the parametrization with Eq. (4), a new measurement of the η → 3π 0 Dalitz plot, with even higher statistical accuracy, is very important.
In this paper, we report on a new high-statistics measurement of the η → 3π 0 Dalitz plot, which is based on 7 × 10 6 detected decays. The A2 data used in the present analysis were taken in 2007 (Run I) and 2009 (Run II). Compared to the previous analysis of Run I reported in Ref. [11], the present analysis was made with an improved cluster algorithm, which increased the number of η → 3π 0 decays reconstructed in Run I from 3 × 10 6 to 3.5 × 10 6 . The γp → ηp → 3π 0 p → 6γp data from Run I and Run II used in this work were previously used to measure the γp → ηp differential cross sections, the analysis of which was recently reported in Ref. [35]. The new η → 3π 0 results were obtained with the parametrization involving the higher-order terms of the Dalitz-plot density function and the cusp term. The NREFT framework from Ref. [34] was also used to check whether the present η → 3π 0 data can be described together with the KLOE η → π + π − π 0 data [14], assuming the isospin limit. The experimental spectra are also compared to recent theo-retical calculations that predict a nonlinear dependence along the quadratic distance from the Dalitz-plot center.

II. EXPERIMENTAL SETUP
An experimental study of the η → 3π 0 decay was conducted via measuring the process γp → ηp → 3π 0 p → 6γp with 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 were determined by the Glasgow tagging spectrometer [41][42][43].
The CB detector is a sphere consisting of 672 optically isolated NaI(Tl) crystals, shaped as truncated triangular pyramids, which point toward the center of the sphere. The crystals are arranged in two hemispheres that cover 93% of 4π, sitting outside a central spherical cavity with a radius of 25 cm, which holds the target and inner detectors. In this experiment, TAPS was initially arranged in a plane consisting of 384 BaF 2 counters of hexagonal cross section. It was installed 1.5 m downstream of the CB center and covered 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, allowing for a higher count rate in the crystals near the photon-beam line. More details on the energy and angular resolution of the CB and TAPS are given in Refs. [11,44].
The present measurement used electron beams with energies of 1508 and 1557 MeV from the Mainz Microtron, MAMI-C [40]. The data with the 1508-MeV beam were taken in 2007 (Run I) and those with the 1557-MeV beam in 2009 (Run II). Bremsstrahlung photons, produced by the beam electrons in a 10-µm Cu radiator and collimated by a 4-mm-diameter Pb collimator, were incident on a liquid hydrogen (LH 2 ) target located in the center of the CB. The LH 2 target was 5 cm and 10 cm long in Run I and Run II, respectively. The total amount of material around the LH 2 target, including the Kapton cell and the 1-mm-thick carbon-fiber beamline, was equivalent to 0.8% of a radiation length X 0 , which was essential to keep the material budget as low as possible to minimize the conversion of final-state photons.
The target was surrounded by a Particle IDentification (PID) detector [45] used to distinguish between charged and neutral particles. The PID consists of 24 scintillator bars (50 cm long, 4 mm thick) arranged as a cylinder with the middle radius of 12 cm.
In Run I, the energies of the incident photons were analyzed up to 1402 MeV by detecting the postbremsstrahlung electrons in the Glasgow tagged-photon spectrometer (Glasgow tagger) [41][42][43], and up to 1448 MeV in Run II. The uncertainty in the energy of the tagged photons is mainly determined by the segmentation of the tagger focal-plane detector in combination with the energy of the MAMI electron beam used in the experiments. Increasing the MAMI energy increases the energy range covered by the spectrometer and also has the corresponding effect on the uncertainty in E γ . For both the MAMI energy settings of 1508 and 1557 MeV, this uncertainty was about ±2 MeV. More details on the tagger energy calibration and uncertainties in the energies can be found in Ref. [43].
The experimental trigger in Run I required the total energy deposited in the CB to exceed ∼320 MeV and the number of so-called hardware clusters in the CB (multiplicity trigger) to be two or more. In the trigger, a hardware cluster in the CB was a block of 16 adjacent crystals in which at least one crystal had an energy deposit larger than 30 MeV. Depending on the data-taking period, events with a cluster multiplicity of two were prescaled with different rates. TAPS was not included in the multiplicity trigger for these experiments. In Run II, the trigger threshold on the total energy in the CB was increased to ∼340 MeV, and the multiplicity trigger required three or more hardware clusters in the CB.

III. DATA ANALYSIS
The η → 3π 0 decays were measured via the process γp → ηp → 3π 0 p → 6γp from events having six or seven clusters reconstructed by a software analysis in the CB and TAPS together. Seven-cluster events were analyzed by assuming that all final-state particles were detected, and six-cluster events by assuming that only the six photons were detected, with the recoil proton going undetected. The offline cluster algorithm [46] was optimized for finding a group of adjacent crystals in which the energy was deposited by a single-photon electromagnetic (e/m) shower. This algorithm also works well for recoil protons. The software threshold for the cluster energy was chosen to be 12 MeV. Compared to the previous η → 3π 0 analysis of Run I [11], the cluster algorithm was improved for a better separation of e/m showers partially overlapping in the calorimeters, which is especially important for processes with large photon multiplicity in the final state and for conditions of the forward energy boost of the outgoing photons in the laboratory system. At the same time, the cluster algorithm has also to be efficient for reconstructing one photon splitting into two nearby e/m showers. The new optimization of the cluster algorithm was needed to improve its efficiency for higher energies of MAMI-C. Particularly for the process γp → ηp → 3π 0 p → 6γp, its reconstruction efficiency was improved by ∼ 17%, compared to the previous analysis [11].
The event identification was based on a kinematic fit, the details of which, including the parametrization of the detector information and resolutions were given in Ref. [11]. Many other details of the event selection in the present work are also very similar to the previous analysis. To test the γp → ηp → 3π 0 p → 6γp hypothesis, 15 combinations are possible to pair six photons into three neutral pions. To reduce the number of combinations tested with the kinematic fit, invariant masses of cluster pairs for each combination were tested prior to fitting. For seven-cluster events, where seven combinations are possible to select the proton cluster, this number was reduced by a cut on the cluster polar angle, the value of which is limited by the recoil-proton kinematics in the laboratory system. 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 into three pions was found to be quite small and could be further reduced by tightening a selection criterion on the kinematic-fit CL. Misidentification of the proton cluster with the photons was found to be negligibly small for seven-cluster events. The sixcluster sample, which includes ∼ 20% from all detected η → 3π 0 decays, had a small contamination from events in which one of the photons, instead of the proton, was undetected. Because such misidentification mostly occurred for clusters in TAPS, those events were successfully removed, based on the cluster's time-of-flight information, which provides good separation of the γp → ηp recoil protons from photons in the present energy range.
To minimize systematic uncertainties in the determination of experimental acceptance, Monte Carlo (MC) simulations of the production reaction γp → ηp were based on the actual spectra measured with the same data sets [35]. The η → 3π 0 decay was generated according to phase space (i.e., with the slope parameter α = 0). The simulated events were propagated through a GEANT (version 3.21) simulation of the experimental setup. To reproduce the resolutions observed in the experimental data, the GEANT output (energy and timing) was subject to additional smearing, thus allowing both the simulated and experimental data to be analyzed in the same way. Matching the energy resolution between the experimental and MC events was achieved by adjusting the invariant-mass resolutions, the kinematic-fit stretch functions (or pulls), and probability distributions. Such an adjustment was based on the analysis of the same data sets for reactions that could be selected with the kinematic fit practically without background from other reactions (namely, γp → π 0 p, γp → ηp → γγp, and γp → ηp → 3π 0 p were used). The simulated events were also tested to check whether they passed the trigger requirements.
For η → 3π 0 decays, physical background can only come from the γp → 3π 0 p events that are not produced from η decays. As shown in Ref. [47], those 3π 0 events are mostly produced via baryon decay chains, with a smaller fraction from γp → K 0 S Σ + → 3π 0 p. For selected γp → ηp → 3π 0 p events, this background is negligibly small near the η production threshold, and reaches ∼ 4% near beam energy E γ = 1.4 GeV. Because of the complicated dynamics of these background processes, they cannot be reproduced precisely with the MC simulation in order to be used for the background subtraction, and additional selection criteria have to be applied instead to reduce the remaining background to a level ≤ 1%. The initial level of the direct 3π 0 background under the η → 3π 0 peak can be seen in the m(3π 0 ) invariant-mass distributions for events selected at CL> 1% by testing the γp → 3π 0 p → 6γp hypothesis, which has no constraint on the η mass. These distributions are shown in Fig. 3. It was checked that the level of the direct 3π 0 background ≤ 1% in the η → 3π 0 data sample could be reached by requiring CL> 1.5% for the γp → ηp → 3π 0 p → 6γp hypothesis along with rejecting events having E γ > 1.3 GeV.
There are two more sources of background remaining in the selected γp → ηp → 3π 0 p → 6γp events and which could directly be subtracted from the experimental spectra. The first background is due to interactions of the bremsstrahlung photons in the windows of the target cell. The evaluation of this background is based on the analysis of data samples that were taken with the target cell emptied of liquid hydrogen. The weight for the subtraction of empty-target spectra is usually taken as a ratio of the photon-beam fluxes for the data samples with the full and the empty target. Because, in the present experiments, the amount of empty-target data were much smaller than with the full target, the subtraction of this background would cause larger statistical uncertainties. It was checked that, for the selection criteria used, the fraction of the empty-target background is ≤ 1%, and this background mostly contains actual η → 3π 0 decays that were just produced in interactions with the targetcell material. Thus, the subtraction of the empty-target background was neglected in the present analysis.
The second background was caused by random coincidences of the tagger counts with the experimental trigger. It mostly includes γp → ηp → 3π 0 p → 6γp events reconstructed with random E γ , resulting in poorer χ 2 and resolution after kinematic fitting. The subtraction of this background was carried out by using event samples for which all coincidences were random (see Ref. [11] for more details). The fraction of random background was 6.7% for Run I, and 6.9% for Run II. The actual background samples included much more events to diminish the impact from statistical fluctuations in the distributions used for the subtraction.  3); (c) the η → π + π − π 0 plot (without boundary bins) from the KLOE analysis of ∼ 4.7 × 10 6 decays [14].

IV. RESULTS AND DISCUSSION
The full Dalitz plot obtained from ∼ 7 × 10 6 η → 3π 0 decays of Run I and Run II is shown in Fig. 4(a). Because there are three identical particles in the final state, variables X and Y can be determined in six different ways, with the same value for variable z and different angle φ from Eq. (3). Each of these six combinations in X and Y goes into six different sextants, repeating the density structure every 60 degrees. The difference between those sextants is only in their different orientation with respect to each other and to the plot binning. Also, this Dalitz plot is symmetric with respect to the Y axis. In principle, one sextant is sufficient to analyze the Dalitz-plot shape and to obtain the corresponding results with proper statistical uncertainties. Such a sextant plot, obtained for the angle range 30 • < φ < 90 • , is shown in Fig. 4(b). As seen, this sextant plot has bins with limited physical coverage not only along the external edge but also along angle φ = 30 • . To avoid any dependence of the results on such an effect and on the sextant orientation with respect to the plot binning, one half of the Dalitz plot (X < 0 or X > 0) can be used to analyze its shape. Because half of the plot has three entries per event, the parameter errors from fitting to such a plot must be multiplied by the factor of √ 3 to reflect the actual experimental statistics.
To obtain the η → 3π 0 plots shown in Figs. 4(a) and 4(b), the plots with the measured decays from Runs I and II were divided by the corresponding plots obtained from the analysis of the γp → ηp → 3π 0 p MC simulations for those data sets. Because η → 3π 0 decays were generated as phase space, the ratio of the experimental and the MC plots provides both the acceptance correction for the full area and the cancellation of the phase-space factor coming from the limited physical coverage, which is typical for boundary bins. Then those boundary bins can be treated in the same way as the inner bins while fitting the acceptance-corrected Dalitz plots with density functions. The only difference from the inner bins is in using X and Y coordinates averaged inside the boundary bins over the available phase space, instead of taking the bin centers. To combine the acceptance-corrected plots from different data sets (namely from Runs I and II), their normalization should be done in the same way. In the present analysis, an identical normalization was made by taking the weight of the MC Dalitz plot as the ratio of the event numbers in the experimental and the MC plots.
As shown in Fig. 4(a), the largest density of events is accumulated in the center of the η → 3π 0 Dalitz plot, with a smooth decrease of a few percent toward the plot edge. To compare such a structure with the charged decay, the acceptance-corrected η → π + π − π 0 Dalitz plot from KLOE [14] (with excluded boundary bins) is illustrated in Fig. 4(c), showing a sharp decrease in its density from the smallest Y to the largest. In the present work, this η → π + π − π 0 plot was used to check whether it could be described together with the η → 3π 0 data within the NREFT framework [34], assuming the isospin limit.
The advantage of analyzing the η → π + π − π 0 decay is the fact that the X and Y variables can be defined uniquely. Then the experimental raw (i.e., uncorrected for the acceptance) Dalitz plot can be fitted with the corresponding plots of the phase-space MC events weighted with the density-function terms. Because the weights are calculated from the generated variables, but filling the MC plots is done according to the reconstructed variables, such a fit takes into account both the experimental acceptance and resolution. For the η → 3π 0 decay, the X and Y generated in one sextant could be reconstructed in another sextant, which allows proper fitting a sextant of the raw Dalitz plot with the density function dependent only on z (which is the same for all pairs of X and Y ) but not on φ. Therefore, all fits with the higher-order terms were made only for the acceptance-corrected Dalitz plots. The sensitivity of the results to the experimental resolution, which could be determined by comparing to the fits to the raw Dalitz plots, was only checked for the leading-order parametrization.
The traditional z distributions, which were used in all The NREFT calculation by the Bonn group [16] is shown in (a) by the black long-dash-dotted line. The prediction from the dispersive analysis by the Bern group [21] is shown in (a) by the magenta long-dashed line with an error band. The prediction based on the extended chiral Khuri-Treiman formalism [22] is shown in (a) by the black dotted line. The fit of the combined z distribution with the leading-order term (fit no. 2 in Table I) Table I, are shown in (b) by the yellow solid and the green dashed lines, respectively. The isospin-limit results from fitting both the present η → 3π 0 and KLOE's η → π + π − π 0 [14] data within the NREFT framework from Ref. [34] are shown by the blue dash-dotted line. The isospin-breaking results from fitting solely the η → 3π 0 data within the same NREFT framework are shown by the red dotted line.
previous measurements of the slope parameters α, were obtained individually for Run I and Run II. Similar to the individual Dalitz plots, their normalization was based on the ratio of the total number of events in the experimental and the MC distributions, which allows the proper combination of the two independent measurements. The individual z distributions from Run I and Run II are compared in Fig. 5(a) with each other and with the earlier A2 data from Ref. [11], demonstrating good agreement within their statistical uncertainties. The combined z distribution, shown in Fig. 5(b), has a statistical accuracy in its 30 data points that appears to be sufficient to reveal the deviation from a linear dependence.
The ratios of the experimental m(π 0 π 0 ) invariant-mass distributions to phase space, in which a cusp structure is expected to be seen, were obtained in the same way as the z distributions. The agreement of the individual m(π 0 π 0 ) distributions from Run I, Run II, and the earlier A2 data from Ref. [11], can be seen in Fig. 6(a). The combined m(π 0 π 0 ) distribution is shown in Fig. 6(b), significantly improving the statistical accuracy in the cusp region, compared to the previous measurement [11].
In addition to fitting the present η → 3π 0 data with the density function from Eq. (4), the NREFT framework from Ref. [34] was used to check whether the neutraldecay data can be fitted well together with the KLOE η → π + π − π 0 data [14] by assuming the isospin limit. Next, the solely η → 3π 0 data were fitted in the same framework by assuming isospin breaking. In Ref. [34], the decay amplitude is decomposed into up to two loops, A(η → 3π) = A tree +A 1−loop +A 2−loop , with the tree amplitude complemented by final-state interactions of one and two loops. The tree amplitudes are parametrized as is the kinetic energy of pion i in the η rest frame. For the conventional Dalitz plot variables, the tree amplitudes can be rewritten as A tree (η → 3π 0 ) = u 0 +u 1 z and A tree (η → π + π − π 0 ) = v 0 + v 1 Y + v 2 Y 2 + v 3 X 2 , where, at the tree level, the quadratic slope parameter is α = u 1 /u 0 , and the coefficients u i and v i are strictly connected to K i and L i , respectively. Note that the shape of the actual η → 3π 0 Dalitz plot is determined by the total amplitude; therefore, a measured α could be different from the ratio u 1 /u 0 of the tree-amplitude coefficients. The coefficients K i and L i (or u i and v i ) are also involved in the calculation of A 1−loop and A 2−loop for both the neutral and charged decays. The cusp structure below 2m π ± appears in A(η → 3π 0 ) 1−loop , and the cusp sign and magnitude is mostly determined by the scattering length combination a 2 − a 0 [32] and the η → π + π − π 0 tree-amplitude coefficients L i . In the isospin limit, the coefficients of the tree amplitude for the neutral decay can be rewritten via the coefficients of the charged decay: , with Q η = m η − 3m π 0 . The isospin-limit fit to both the η → 3π 0 and η → π + π − π 0 Dalitz plots has only five free parameters (L i=1,2,3 and two normalization parameters), with fixed L 0 = 1. The η → 3π 0 data can also be fitted independently of the η → π + π − π 0 decay by assuming isospin breaking, which requires the addition of K 0 and K 1 as free parameters, but leaves just one normalization parameter.
Consistency of the present results for z and m(π 0 π 0 )  6: Ratios of the experimental m(π 0 π 0 ) invariant-mass distributions to phase space obtained (a) individually from Run I (blue circles) and Run II (red triangles), and (b) the combined results (black triangles). The earlier A2 data from Ref. [11] are depicted in (a) by green open squares. The NREFT calculation by the Bonn group [16] is shown in (a) by the black long-dash-dotted line. The prediction from the dispersive analysis by the Bern group [21,33] is shown in (a) by the magenta long-dashed line. The prediction based on the extended chiral Khuri-Treiman formalism [22] is shown in (a) by the black dotted line. The combined m(π 0 π 0 ) distribution is compared in (b) to the results of fitting a sextant (30 • < φ < 90 • ) of the acceptance-corrected η → 3π 0 Dalitz plot with the density function of Eq. (4): fits no. 1 (cyan long-dashed line), no. 4 (yellow solid line), and no. 6 (green dashed line) in Table I. The isospin-limit results from fitting both the present η → 3π 0 and KLOE's η → π + π − π 0 [14] data within the NREFT framework from Ref. [34] are shown by the blue dash-dotted line. The isospin-breaking results from fitting solely the η → 3π 0 data within the same NREFT framework are shown by the red dotted line.
with theoretical calculations that predict a nonlinear z dependence [16,21,22] is illustrated in Figs. 5(a) and 6(a). The results of fits to the present data with various density functions, including the NREFT fits, are depicted in Figs. 5(b) and 6(b). The fit results with the density function from Eq. (4) are also listed in Table I for different combinations of the density-function terms involved in a particular fit.
Fit no. 1 in Table I was made to a sextant (30 • < φ < 90 • ) of the acceptance-corrected Dalitz plot with the density function including only the leading-order term. Fit no. 2 was similar, but to the acceptance-corrected z distribution as in all previous measurements. As shown, the values obtained there for α are practically the same and are in agreement within the fit errors with the RPP value α = −0.0318 ± 0.0015 [23]. The magnitudes of the fit χ 2 /ndf values indicate that the use of the leading-order term only may be insufficient for a good description of the η → 3π 0 decay. Fit no. 2 is shown in Fig. 5(b) and fit no. 1 in Fig. 6(b) by the cyan long-dashed lines, confirming that it is not sufficient to use only the leading-order term. Fit no. 3 in Table I was made to the same sextant of the raw Dalitz plot with the technique taking both the acceptance and the experimental resolution into account (see the text above). This fit results in a slightly better χ 2 /ndf value and a slightly larger quadratic slope, which was expected because of some smearing of the acceptance-corrected distributions by the experimental resolution. In the end, the difference between the α results for the acceptance-corrected and the raw distributions can be considered as the magnitude of its systematic uncertainty due to the limited experimental resolution. Table I, which also involves the next density-function term 2βz 3/2 sin(3φ), does improve the χ 2 /ndf value, whereas including the 2γz 2 term in fit no. 5 practically does not. In addition, the parameters α and γ in fit no. 5 become strongly correlated, which results in large fit errors for them. Fit no. 4, shown in Figs. 5(b) and 6(b) by the yellow solid line, demonstrates a quite decent description of the z and m(π 0 π 0 ) distributions, except in the region where the cusp is expected. As shown in the m(π 0 π 0 ) distribution, the 2βz 3/2 sin(3φ) term curves the spectrum up at the lowest masses, which is opposite to the effect expected from the cusp. In the z distribution, the same term causes a kink up at z ≈ 0.75, which again is opposite to the effect expected from the cusp [11,15]. As shown in Figs. 5(a) and 6(a), the calculation within the framework of the modified NREFT [16] predicts a behavior that is very similar to fit no. 4, but with a smaller general slope. This can be explained by a smaller quadratic slope, α = −0.0246(49), and positive γ = 0.0013(4) from Ref. [16]. However, because of the large uncertainty in the calculated α, it is still in agreement with the corresponding value from fit no. 4. In contrast to the calculation from Ref. [16], the prediction based on the extended chiral Khuri-Treiman formalism [22] lies below the experimental data points, which is mostly determined by the larger quadratic slope, α = −0.0337 (12). At the same time, the predictions for the 2βz 3/2 sin(3φ) term, β = −0.0042(7) [16] and β = −0.0054(1) [22], are both in decent agreement with the corresponding value from fit no. 4. The experimental value for γ cannot be determined reliably in order to be compared with the prediction from Ref. [16]. As seen from fit no. 6 in Table I, further improvement in the description of the η → 3π 0 data was reached by adding the 2δ 3 i=1 ρ(s i ) term, which allows a cusp parametrization to be included in the density function. Such a fit results in a slightly smaller quadratic slope, compared to fit no. 4, but also in a stronger 2βz 3/2 sin(3φ) term. In Figs. 5(b) and 6(b), fit no. 6, which is shown by the green dashed line, demonstrates good agreement with both z and m(π 0 π 0 ) distributions. Based on the results of fit no. 6, the contributions from the 2βz 3/2 sin(3φ) and the cusp terms partially cancel each other in the z and especially in the m(π 0 π 0 ) distribution. Though, according to the result of fit no. 6 for the cusp term, the magnitude of the cusp effect at m(π 0 π 0 ) = 2m π 0 is almost 1%, its visibility here is strongly diminished by the 2βz 3/2 sin(3φ) term. The understanding of such a feature became possible due to fitting the η → 3π 0 Dalitz plot based on high experimental statistics.

Fit no. 4 in
The isospin-limit NREFT fit to the present η → 3π 0 data together with KLOE's η → π + π − π 0 Dalitz plot [14] is shown in Figs. 5(b) and 6(b) by the blue dash-dotted line. As shown in the m(π 0 π 0 ) distribution, the major deviation of this fit from the data is in the cusp region, which is much more prominent in the fit curve. The description of the z distribution deviates from the data as well. The cusp magnitude obtained at m(π 0 π 0 ) = 2m π 0 is close to 1%, which is similar to the corresponding result of fit no. 6 in Table I. The discrepancy seems to come from inability of the isospin-limit fit to describe properly the 2βz 3/2 sin(3φ) term. Though the isospin-limit NREFT fit results in a good description of the charged decay, with χ 2 /ndf=1.072, it gives χ 2 /ndf=1.290 for the neutral decay. The numerical results for L i were obtained as L 0 = 1(0), L 1 = −4.004(31), L 2 = −41.55 (31), and L 3 = 5.28 (14), with K i recalculated from L i as K 0 = −2.322(7) and K 1 = 25.71(73).
A comparison of the results from the two NREFT fits indicates a strong isospin breaking between the charged and the neutral η → 3π decays, unless the NREFT framework in Ref. [34] could be improved for a better simultaneous description of both decay modes. As illustrated in Figs. 5(a) and 6(a), a recent dispersive analysis by the Bern group [21,33], in which the η → π + π − π 0 data [14] were used to determine subtraction constants, did provide predictions that described the η → 3π 0 data well.
The results of this work provide a strong indication that the parametrization of the η → 3π 0 decay with only the leading-order term is insufficient, and the RPP value α = −0.0318 ± 0.0015 [23] reflects a combined effect from higher-order terms and the cusp structure. As the results listed in Table I show, the values obtained for the quadratic slope parameter become smaller when the higher-order terms and the cusp are added, and those values for α are also closer to recent calculations reported in Refs. [16,18,19] (see also Fig. 2).
The exact systematic uncertainties in the results for α and for the other parameters are difficult to estimate reliably because the results themselves depend on the number of density-function terms included in the fit. The systematic effect due to the limited experimental resolution was discussed above for a fit with the leading-order term only (no. 3 in Table I). The sensitivity of the results to the sextant orientation with respect to the plot binning and to additional boundary bins was checked with fits to other sextants and to half of the Dalitz plot. All those tests demonstrated practically identical results, after multiplying the half-plot errors by the factor of √ 3 to correct for three entries per event (fits nos. 8-10 in Table I). The magnitudes of systematic effects for all parameters could also be understood by comparing fits to the independent data of Run I and Run II, which were taken with different MAMI beam energy and current, target length (resulting in different angular resolution), DAQ trigger, energy resolution of the calorimeters, etc. Those fits are listed as nos. 11-16 in Table I. As shown, the largest differences between the results from Run I and Run II were observed for parameters γ and δ; however, all results obtained from the different data sets are in agreement within the fit errors. The magnitude for parameter γ cannot be determined reliably from the experimental data because of the large correlation with parameter α. Therefore, the value obtained for α with the 2γz 2 term omitted actually reflects the combined effect from those two terms.
According to the present analysis, the density function of Eq. (4) with only three parameters is sufficient for a good description of the experimental η → 3π 0 Dalitz plot. The values obtained for these three parameters are α = −0.0265(10 stat )(9 syst ), β = −0.0074(10 stat )(9 syst ), and δ = −0.018(7 stat )(7 syst ), where the main numbers come from fit no. 6 in Table I, and the systematic uncertainties are taken as half of the differences between the results of fits nos. 13 and 14. The new result for the quadratic slope parameter α strongly indicates that its absolute value is smaller by ≈ 20%, compared to the previous measurements using the leading-order term only. The magnitude of the 2βz 3/2 sin(3φ) term is found to be different from zero by ∼ 5.5 standard deviations. The cusp magnitude obtained at m(π 0 π 0 ) = 2m π 0 from the 2δ 3 i=1 ρ(s i ) term is close to 1%, but with an uncer-tainty greater than 50%. This result is consistent with the prediction for the η → 3π 0 cusp magnitude made within the NREFT model [15]. The data presented in this work are expected to serve as a valuable input for new refined analyses by theoretical groups, which are interested in a better understanding of η → 3π decays and extracting the quark-mass ratios from such data.

V. SUMMARY AND CONCLUSIONS
The largest, at the moment, statistics of 7 × 10 6 η → 3π 0 decays, based on 6.2×10 7 η mesons produced in the γp → ηp reaction, has been accumulated by the A2 Collaboration at the Mainz Microtron, MAMI. The results of this work provide a strong indication that the parametrization of the η → 3π 0 decay with only the leading-order term is insufficient, and the RPP value for α reflects the combined effect from higher-order terms and the cusp structure, whereas the actual quadratic slope is smaller by ≈ 20%. According to the analysis of the η → 3π 0 Dalitz plot, the cusp magnitude at m(π 0 π 0 ) = 2m π 0 is about 1%, but its visibility is strongly diminished by the second-order term of the density function, the magnitude of which is found to be different from zero by ∼ 5.5 standard deviations. The fits to the present η → 3π 0 and KLOE's η → π + π − π 0 data within the NREFT framework indicate a strong isospin breaking between the charged and the neutral decay modes. At the same time, the predictions based on the most recent dispersive analysis by the Bern group, in which the η → π + π − π 0 data were used to determine subtraction constants, were found to be in good agreement with the present η → 3π 0 data. The data points from the experimental Dalitz plot and the ratios of the z and m(π 0 π 0 ) distributions to phase space are provided as supplemental material to the paperl [48].