Measurement of the matrix elements for the decays $\eta^{\prime}\rightarrow\eta\pi^+\pi^-$ and $\eta^{\prime}\rightarrow\eta\pi^0\pi^0$

Based on a sample of $1.31\times10^9$ $J/\psi$ events collected with the BESIII detector, the matrix elements for the decays $\eta^{\prime}\to\eta\pi^{+}\pi^{-}$ and $\eta^{\prime}\to\eta \pi^{0}\pi^{0}$ are determined using 351,016 $\eta^{\prime}\to(\eta\to\gamma\gamma)\pi^+\pi^-$ and 56,249 $\eta^{\prime}\to(\eta\to\gamma\gamma)\pi^0\pi^0$ events with background levels less than 1\%. Two commonly used representations are used to describe the Dalitz plot density. We find that an assumption of a linear amplitude does not describe the data well. A small deviation of the obtained matrix elements between $\eta^{\prime}\to\eta\pi^{+}\pi^{-}$ and $\eta^{\prime}\to\eta \pi^{0}\pi^{0}$ is probably caused by the mass difference between charged and neutral pions or radiative corrections. No cusp structure in $\eta^{\prime}\to\eta \pi^{0}\pi^{0}$ is observed.


I. INTRODUCTION
The η 0 meson is well established, and its main decay modes are fairly well known [1]. However, η 0 decay dynamics remains a subject of extensive theoretical studies aiming at extensions of the chiral perturbation theory (ChPT). The two dominant hadronic decays, η 0 → ηπ þ π − and η 0 → ηπ 0 π 0 (called charged decay mode and neutral decay mode throughout the text, respectively), are believed to be an ideal place to study ππ and ηπ scattering [2,3], which may lead to a variation in the density of the Dalitz plot. Several extensions of the ChPT framework [4][5][6][7] and dispersive analysis based on the fundamental principles of analyticity and unitarity [8] have been applied to investigate the matrix element of η 0 → ηππ.
In experimental analyses, the Dalitz plot for the charged decay mode is usually described by the following two variables: For the neutral decay mode, the Dalitz plot has a twofold symmetry due to the two π 0 s in the final state. Hence, the variable X is replaced by Here, T π and T η denote the kinetic energies of a pion and η in the η 0 rest frame, Q ¼ m η 0 − m η − 2m π , and m π , m η , and m 0 η are the masses of the pion, η, and η 0 , respectively. Generally, the decay amplitude squared is parametrized as which is the so-called general representation. Here a, b, c and d are free parameters, and N is a normalization factor. The terms with odd powers in X are forbidden due to the charge conjugation symmetry in η 0 → ηπ þ π − and the wave function symmetry in η 0 → ηπ 0 π 0 . By considering the isospin symmetry, the Dalitz plot parameters for the *  charged and neutral decay modes should be the same. However, a small discrepancy, observed in previous measurements [9][10][11][12][13], is expected due to the mass difference between the charged and neutral pion or due to radiative corrections for the η 0 → ηπ þ π − mode [6].
A second parametrization for the decay amplitude squared used by previous experiments assumes a linear amplitude in Y and keeps the polynomial expansion in X, the so-called linear representation, where α is a complex number. The real part of α gives the linear term in Y for the Dalitz plot density, a ¼ 2ℜðαÞ, and the quadratic term is b ¼ ℜðαÞ 2 þ ℑðαÞ 2 , where ℜðαÞ and ℑðαÞ are the real and imaginary parts of α, respectively. The two representations are equivalent if b > a 2 =4, i.e., b should be at least larger than zero. Therefore, a negative value for b demonstrates that the ansatz of Eq. (4) does not describe the data.
In addition, the Dalitz plot for η 0 → ηπ 0 π 0 is expected to be affected by a cusp due to the π þ π − mass threshold. The size of this effect is predicted to be about 6% [14] (8% in original work [6]) within the framework of nonrelativistic effective field theory (NREFT), which is confirmed in a dispersive analysis [8]. An analogous cusp has been observed in the K þ → π þ π 0 π 0 [15] decay and allows us to determine the ππ S-wave scattering lengths [16][17][18][19].
The dynamics of the decays η 0 → ηπ þ π − and η 0 → ηπ 0 π 0 are studied in this work using η 0 mesons produced in the J=ψ → γη 0 decay. The present data sample of 1.31 × 10 9 J=ψ events accumulated with the BESIII detector is about 5 times of that used in the previous BESIII analysis [11].

II. DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector is a general-purpose magnetic spectrometer with a geometrical acceptance of 93% of 4π and is described in detail in Ref. [20]. The detector is composed of a helium-based drift chamber (MDC), a plastic-scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field (0.9 T in 2012). The solenoid is supported by an octagonal fluxreturn yoke with resistive-plate counters interleaved with steel for muon identification (MUC).
The GEANT4 [21] based Monte Carlo (MC) simulation software package BOOST [22] describes the geometry and material of the BESIII detector, as well as the detector response. It is used to optimize the event selection criteria, estimate backgrounds and determine the detection efficiencies. The production of the J=ψ resonance is simulated with KKMC [23,24], while the decays are generated with EVTGEN [25,26] for established modes using world-average branching fractions [1], and by LUNDCHARM [27] for the remaining decays. An inclusive MC sample of 1.2 × 10 9 J=ψ events is used to study the potential background contributions. The analysis is performed in the framework of the BESIII off-line software system (BOSS) [28].
For the reconstruction of J=ψ → γη 0 with η 0 → ηπ þ π − and η → γγ, candidate events must contain two tracks with an opposite charge and at least three photons. Each charged track reconstructed from the MDC hits is required to have a polar angle in the range j cos θj < 0.93 and to pass the interaction point within AE10 cm along the beam direction and within AE1 cm in the plane perpendicular to the beam. Photon candidates are reconstructed using isolated clusters of energy deposited in the EMC and required to have a deposited energy larger than 25 MeV in the barrel region (j cos θj < 0.80) or 50 MeV in the end cap region (0.86 < j cos θj < 0.92). The energy deposited in nearby TOF counters is included to improve the reconstruction efficiency and energy resolution. To eliminate clusters associated with charged tracks, the angle between the photon candidate and the extrapolation of any charged track to the EMC must be larger than 10°. A requirement on the EMC cluster timing with respect to the event start time (0 ≤ T ≤ 700 ns) is used to suppress electronic noise and energy deposits unrelated to the event.
Since the radiative photon from the J=ψ decay is always more energetic than those from the η decay, the photon candidate with the maximum energy in an event is assumed to be the radiative one. For each π þ π − γγγ combination, a six-constraint (6C) kinematic fit is applied, and the χ 2 6C is required to be less than 100. The fit enforces energymomentum conservation and constrains the invariant masses of γγ and ηπ þ π − to the nominal η and η 0 masses, respectively. If there are more than three selected photons in an event, the combination with the smallest χ 2 6C is retained. To estimate the background contribution, an alternative data sample is selected without applying the η and η 0 mass constraints in the kinematic fit. The π þ π − γγ invariant mass spectrum is shown in Fig. 1 after requiring the γγ invariant mass within the η signal region, ð0.518; 0.578Þ GeV=c 2 . A clear η 0 signal is observed with a low background level. In addition, a sample of 1.2 × 10 9 inclusive MC J=ψ decays is used to investigate potential backgrounds. Using the same selection criteria for the MC sample, no peaking background remains around the η 0 signal region. From this MC sample, the background contamination is estimated to be about 0.3%. This is consistent with an estimation obtained from an unbinned maximum likelihood fit to the Mðηπ þ π − Þ distribution, where the signal is described by the MC simulated shape convoluted with a Gaussian function representing the resolution difference between the data and MC simulation, and the background contribution is described by a third-order polynomial function. We therefore neglect the background contribution in the determination of the Dalitz plot parameters.
After the above requirements, 351,016 η 0 → ηπ þ π − candidate events are selected, with an averaged efficiency of 31.2% and a background contribution of less than 0.3%. Figure 2 shows the Dalitz plot in the variables X and Y for the selected events. The corresponding projections on X and Y are shown as the dots with error bars in Figs. 3(a) and 3(b), respectively. The resolution on the variables X and Y over the entire kinematic region, determined from the MC simulation, are 0.03 and 0.02, respectively.
Unbinned maximum likelihood fits to the data are performed to determine the free parameters in the decay amplitude squared [Eqs. (3) and (4)]. To account for the resolution and detection efficiency, the amplitude squared is convoluted with a function σðX; YÞ parametrizing the resolution and multiplied by a function εðX; YÞ parametrizing the detection efficiency. Both functions are derived from MC simulations. Two double Gaussian functions are used for σðX; YÞ, while εðX; YÞ is estimated as the average efficiencies of local bins. With the normalization, one derives the probability density function PðX; YÞ, which is applied in the fit, The integral over the full Dalitz plot range (DP) gives the normalization factor in the denominator. The fit is done by minimizing the negative log-likelihood value where PðX i ; Y i Þ is evaluated for an event i, and the sum runs over all accepted events. Imposing charge conjugation invariance by setting the coefficient of odd powers in X (c) to zero in the general representation, the fit yields following parameters: Here, the uncertainties are statistical only. The corresponding correlation matrix of the fit parameters is Projections of the fit result on X and Y are illustrated as the solid histograms in Fig. 3.
To check for the existence of a charge conjugation violating term, an alternative fit with a free parameter c is performed. The resultant value, c ¼ ð2.7 AE 2.4Þ × 10 −3 , is consistent with zero. Compared with the nominal fit results, the parameters a, b and d are almost unchanged, and the statistical significance for a nonzero value of the parameter c is determined to be 0.7σ only.
Alternative fits including the extra terms fY 3 þ gX 2 Y or eXY þ hXY 2 þ lX 3 in the general representation are also performed, resulting in f ¼ −0.004 AE 0.012, g ¼ 0.008 AE 0.010 or e ¼ 0.005 AE 0.007, h ¼ 0.004 AE 0.006, and l ¼ 0.007 AE 0.013, respectively, while the other parameters are unchanged.
A fit based on the linear representation is also performed and yields the following values: ℜðαÞ ¼ −0.034 AE 0.002; The imaginary part of α is consistent with zero. This can be understood by the observation that the coefficient b in the general representation is negative. Subsequently we will consider the fit result with ℑðαÞ fixed at zero. The parameters ℜðαÞ and d and their uncertainties remain the same as in Eq. (9), and the correlation coefficient between ℜðαÞ and d is −0.137. The log-likelihood value is lower by 33.9 compared with the fit using the general representation, which indicates that the linear representation is less compatible with the data. Projections on X and Y based on this result are illustrated as the dashed histograms in Fig. 3. The presented residuals show that the fit is slightly worse to describe the data in Y projection comparing to the general one.
The potential charge conjugation violating is also checked in the linear representation by performing an alternative fit with a free parameter c. The resultant value, c ¼ ð2.7 AE 2.4Þ × 10 −3 , is also consistent with zero, while the parameters ℜðαÞ and d are almost unchanged compared with the nominal fit results.

IV. MEASUREMENT OF THE MATRIX ELEMENT
FOR THE DECAY η 0 → ηπ 0 π 0 In the reconstruction of J=ψ → γη 0 with η 0 → ηπ 0 π 0 and η=π 0 → γγ, candidate events must have at least seven photons and no charged track. The selection criteria for photon candidates are the same as those for η 0 → ηπ þ π − , except that the requirement on the angle between photon candidates and any charged track is not used. A requirement of an EMC cluster timing with respect to the most energetic photon (−500 ≤ T ≤ 500 ns) is also used. The photon with the largest energy in the event is assumed to be the radiative photon originating from the J=ψ decay. For the remaining clusters, pairs of photons are combined into π 0 =η → γγ candidates, which are subjected to a oneconstraint (1C) kinematic fit by constraining the invariant mass of the photon pair to be the nominal π 0 or η mass. The χ 2 for this 1C kinematic fit is required to be less than 25. To suppress π 0 miscombinations, the π 0 decay angle θ decay , defined as the polar angle of one of the decay photons in the γγ rest frame with respect to the π 0 flight direction, is required to satisfy j cos θ decay j < 0.95. Then an eight-constraint (8C) kinematic fit is performed for the γηπ 0 π 0 combination enforcing energy-momentum conservation and constraining the invariant masses of the three photon pairs and the ηπ 0 π 0 combination to the nominal π 0 =η and η 0 masses. If more than one combination is found in an event, only the one with the smallest χ 2 8C is retained. Events with χ 2 8C < 100 are accepted for further analysis.
To estimate the backgrounds, an alternative selection is performed where η and η 0 mass constraints in the kinematic fit are removed. The resulting π 0 π 0 η invariant mass spectrum is shown in Fig. 4, after requiring the γγ invariant mass within the η signal region, ð0.518; 0.578Þ GeV=c 2 . The inclusive MC study shows that the surviving backgrounds mainly consist of the peaking background η 0 → π 0 π 0 π 0 and a flat contribution from J=ψ → ωη with ω → γπ 0 and η → π 0 π 0 π 0 . From this MC sample, the background contamination is estimated to be about 0.9%, which is consistent with the estimation obtained from a fit to Mðηπ 0 π 0 Þ and therefore neglected in the determination of the Dalitz plot parameters. In the fit, the signal is described by the MC simulated shape convoluted with a Gaussian function representing the difference of the mass resolution between the data and MC simulation. The shape and the yield of the peaking background η 0 → π 0 π 0 π 0 are fixed according to the dedicated MC simulation [29]. A third-order polynomial function is used to represent the smooth background contribution.
After the above requirements, 56,249 η 0 → ηπ 0 π 0 candidate events are selected, with an averaged efficiency of 9.6% and a 0.9% background level. The Dalitz plot of selected events is displayed in Fig. 5. The corresponding projections on X and Y are shown as the dots with error bars in Fig. 6. The resolution on X and Y over the entire kinematic region, determined from the MC simulation, are 0.05 and 0.04, respectively.
As in the analysis of the η 0 → ηπ þ π − , an unbinned maximum likelihood fit method is used to determine the Dalitz plot parameters. The resolution is described with two double Gaussian functions, and the detection efficiencies in different X and Y bins are obtained from the MC simulation. From a dedicated study with the control sample of J=ψ → π þ π − π 0 , we find that the reconstruction efficiency for the π 0 candidate differs significantly between data and the MC simulation at low π 0 momenta. Thus, to describe the detection efficiency more accurately, an efficiency correction depending on the π 0 momentum is carried out, and the error of this correction will be considered in systematic uncertainty.
Considering the strict constraint from the symmetry of the wave function, only the fits without odd powers of X are performed. The fit based on the general representation yields the coefficient (with statistical uncertainties only) and the corresponding correlation matrix, Similarly to the case η 0 → ηπ þ π − , the fit gives a negative value of the coefficient b. The projections of the fit results on X and Y are shown as the solid histograms in Fig. 6. Extra terms fY 3 and gX 2 Y in the general representation are also added in an alternative fit, resulting in f ¼ −0.023 AE 0.028 and g ¼ 0.024 AE 0.025. The significance for nonzero values of f and g is 0.6σ.
In the fit based on the linear representation, the imaginary part of α (fitted to be 0.000 AE 0.038) does not contribute to the fit quality, as in the η 0 → ηπ þ π − case. Thus, the nominal fit omitting ℑðαÞ gives the results (statistical uncertainties only), The correlation coefficient between the two parameters is −0.170. Compared to the fit based on the general representation, the log-likelihood value is reduced by 13.7.
Projections on X and Y are illustrated as the dashed histograms in Fig. 6. Again, the fits based on the two different representations give similar results for the X projections, but slightly worse for Y in the linear case.

V. SYSTEMATIC UNCERTAINTIES
Various sources of systematic uncertainties on the measured Dalitz plot parameters have been investigated, including tracking efficiency, kinematic fit, efficiency correction, and resolution. For the decay η 0 → ηπ 0 π 0 , additional uncertainties associated with photon miscombination, π 0 and η reconstruction are also considered.
Differences between the data and MC simulation for the tracking efficiency of charged pions are investigated using the control sample J=ψ → ppπ þ π − . A momentum dependent correction on the detection efficiency is obtained by comparing the efficiency between the data and MC simulation. Similarly, a momentum dependent correction for the η reconstruction efficiency is obtained with the control sample of J=ψ → γηπ þ π − . Then alternative fits are performed by incorporating the efficiency corrections for charged pions or η. Changes of the Dalitz plot parameters with respect to the nominal results are assigned as the systematic uncertainties. A momentum-dependent π 0 reconstruction efficiency correction has been applied in the nominal fit; the associated systematic uncertainties are estimated by changing the correction factor by one of its standard deviation and repeating the fit. In comparison with Mðπ 0 π 0 Þ without the π 0 reconstruction efficiency correction, it is found that this correction has little impact on the cusp region.
The possible miscombination of photons in signal MC samples has been studied by matching the generated photon pairs to the selected π 0 or η candidates. The fraction of events with wrong combinations is determined to be 2.7% for η 0 → ηπ 0 π 0 . Alternative fits are performed to the MC simulated sample with only truth-tagged events and the ones including miscombinations, individually. The difference between those two results are taken as the systematic uncertainties.
To estimate the uncertainties associated with the kinematic fitting procedure, the fit results are compared using a 4C (6C) instead of a 6C (8C) kinematic fit for η 0 → ηπ þ π − (η 0 → ηπ 0 π 0 ), and the corresponding changes in the fit parameters are taken as systematic uncertainties.
To estimate the uncertainties associated with the efficiency correction in Eq. (5), we change the Dalitz plot variables X and Y to the so-called square Dalitz plot variables MðηπÞ 2 and cos θ, where θ is the angle between the two pions in the rest frame of ηπ. Alternative fits are performed with the efficiency correction based on the newly defined Dalitz plot variable, and the resultant changes of the Dalitz plot parameters with respect to the nominal results are assigned as systematic uncertainties.
To estimate the impact from the nonflat resolution in the X-Y plane, the biases from input/output checks are taken as the systematic uncertainties. The impact from different resolutions of the Dalitz plot variables between data and the MC simulation is estimated by alternative fits varying the resolutions by AE10%. It is found that the change of the results is negligible. The effect of neglecting the residual background is checked by alternative fits including MC simulated backgrounds and found to be insignificant.  All of the above uncertainties are summarized in Table I. Assuming all the sources of systematic uncertainty are independent, the total systematic uncertainties for the Dalitz plot parameters are obtained by adding the individual values in quadrature, shown in the last row of Table I.
After the event selection criteria presented in Secs. III and IV, clean η 0 → ηπ þ π − and η 0 → ηπ 0 π 0 samples are selected. A comparison between the charged and neutral decay modes could be performed by dividing the acceptance corrected experimental distributions with the corresponding phase space distributions on variables X (absolute value for η 0 → ηπ þ π − ), Y, MðππÞ, and MðηπÞ, which are shown in Fig. 7, together with the Dalitz plot fit results based on the general representation. Although the statistical errors are large, the trends of the experimental distributions on Y and MðππÞ between the charged and neutral mode are obviously different, as the high statistical simulation based on the fit results on the general representation shows. At the same time, the difference on X and MðηπÞ are smaller. The observed differences are likely to be related to the ππ and ηπ final interaction.
The ratio between experimental and phase space distributions on Y and MðππÞ, Figs. 7(b) and 7(c), also provide the possibility to check the cusp effect. Overlaid on Fig. 7(b) is the prediction for η 0 → ηπ 0 π 0 in Ref. [8] based on the previous BESIII fit result for η 0 → ηπ þ π − [11], which are consistent with the experimental distribution within statistical errors. However, with current statistics, it is difficult to establish the structure (cusp effect) near the π þ π − mass threshold.

VII. SUMMARY
With a sample of 1.31 × 10 9 J=ψ events collected with the BESIII detector, clean samples of 351,016 η 0 → ηπ þ π − events and 56,249 η 0 → ηπ 0 π 0 events are selected from J=ψ radiative decays. Then the most precise measurements of the matrix element for the η 0 → ηπ þ π − and η 0 → ηπ 0 π 0 decays as well as a search for the cusp effect in η 0 → ηπ 0 π 0 are performed. central value in Ref [8] fit uncertainty in Ref [8] phase uncertainty in Ref [8] (b) Both the general and the linear representations are used to determine the Dalitz plot parameters, and the corresponding results are summarized in Table II including the systematic uncertainties. The Dalitz plot parameters for both decays are in reasonable agreement and more precise than the previous measurements [9][10][11][12]. The results for η 0 → ηπ þ π − supersede the previous BESIII measurement [11], which used a subsample of the present data. As reported in Ref. [11], the discrepancy of the parameter a for η 0 → ηπ þ π − with respect to the VES value [10] is evident, which, at present, stands at about 3.8 standard deviations. The values of the parameter c in η 0 → ηπ þ π − are all consistent with zero within one standard deviation in both representations, in agreement with the charge conjugation conservation in the strong interaction. In addition, a discrepancy of 2.6 standard deviations for the parameter a is observed between η 0 → ηπ þ π − and η 0 → ηπ 0 π 0 processes, indicating an isospin violation. However, the result is not statistically significant enough to firmly establish such a violation, and additional effects, e.g., radiative corrections [6], should be considered in future experimental and theoretical studies.
A comparison between the results obtained from the general representation and the theoretical predictions within the framework of U(3) chiral effective field theory (EFT) incorporating with a relativistic coupled-channels approach [5] is given in Table II. In general, our results are compatible with the theoretical expectations. However, the theoretical prediction for the parameter a from η 0 → ηπ þ π − is about 2 times larger than our result, and the discrepancies on the parameter d for both η 0 → ηπ þ π − and η 0 → ηπ 0 π 0 are about four standard deviations. Table II also provides the predictions obtained in the frameworks of large-N C ChPT and resonance chiral theory (RChT) with the parameters a fixed according to the boundaries measured in Refs. [10,12]. The expected values are consistent with our results within two standard deviations in both decay modes, except that the parameter d in η 0 → ηπ þ π − is 3.1 standard deviations from the large-N C ChPT, and the parameter b in η 0 → ηπ 0 π 0 is 2.7 standard deviations from the RChT.
As previously mentioned, the linear and general representations are equivalent for the case of b > a 2 =4. However, the coefficients b for the Y 2 term are negative with 5.8 and 4.9 standard deviations to zero for η 0 → ηπ þ π − and η 0 → ηπ 0 π 0 , respectively, which implies that these two representations can not provide an identical description of data. In case of the linear representation, the results are in agreement with previous measurements and also provide a reasonable description on the X projection for both decay modes. However, the goodness of fit on the Y projections are worse than the general one. This is consistent with the conclusion reported by the VES Collaboration [10] that the linear representation can not describe the data well.
We also attempt to search for the cusp effect in the decay η 0 → ηπ 0 π 0 . Inspection of the π 0 π 0 mass spectrum around the π þ π − mass threshold does not show evidence of a cusp with current statistics.