Light meson spectroscopy from Dalitz plot analyses of $\eta_c$ decays to $\eta' K^+ K^-$, $\eta' \pi^+ \pi^-$, and $\eta \pi^+ \pi^-$ produced in two-photon interactions

We study the processes $\gamma \gamma \to \eta_c \to \eta' K^+ K^-$, $\eta' \pi^+ \pi^-$, and $\eta \pi^+ \pi^-$ using a data sample of 519 $fb^{-1}$ recorded with the BaBar detector operating at the SLAC PEP-II asymmetric-energy $e^+e^-$ collider at center-of-mass energies at and near the $\Upsilon(nS)$ ($n = 2,3,4$) resonances. This is the first observation of the decay $\eta_c \to \eta' K^+ K^-$ and we measure the branching fraction $\Gamma(\eta_c \to \eta' K^+ K^-)/(\Gamma(\eta_c \to \eta' \pi^+ \pi^-)=0.644\pm 0.039_{\rm stat}\pm 0.032_{\rm sys}$. Significant interference is observed between $\gamma \gamma \to \eta_c\to \eta \pi^+ \pi^-$ and the non-resonant two-photon process $\gamma \gamma \to \eta \pi^+ \pi^-$. A Dalitz plot analysis is performed of $\eta_c$ decays to $\eta' K^+ K^-$, $\eta' \pi^+ \pi^-$, and $\eta \pi^+ \pi^-$. Combined with our previous analysis of $\eta_c \to K \bar K \pi$, we measure the $K^*_0(1430)$ parameters and the ratio between its $\eta' K$ and $\pi K$ couplings. The decay $\eta_c \to \eta' \pi^+ \pi^-$ is dominated by the $f_0(2100)$ resonance, also observed in $J/\psi$ radiative decays. A new $a_0(1700) \to \eta \pi$ resonance is observed in the $\eta_c \to \eta \pi^+ \pi^-$ channel. We also compare $\eta_c$ decays to $\eta$ and $\eta'$ final states in association with scalar mesons as they relate to the identification of the scalar glueball.

the identification of the scalar glueball.

I. INTRODUCTION
Scalar mesons remain a puzzle in light meson spectroscopy: they have complex structure, and there are too many states to be accommodated within the quark model without difficulty [1]. In particular, the structure of the isospin I= 1 2 Kπ S-wave is still poorly understood, which limits the precision of measurements involving a Kπ system in the final state, including recent searches for CP violation in B meson decay [2], and studies of new exotic resonances [3] and charmed mesons [4].
Decays of the η c , the lightest pseudoscalar cc state, provide a window on light meson states. The BABAR experiment first performed a Dalitz plot analysis of η c →K + K − π 0 and η c →K + K − η using an isobar model [5]. The analysis reported the first observation of K * 0 (1430)→Kη, and observed that η c decays into three pseudoscalars are dominated by intermediate scalar mesons. This newly observed K * 0 (1430) decay mode was expected to be small and in fact was not observed in the study of K − p→K − ηp interactions [6]. More recently, the BABAR experiment performed a measurement of the I= 1 2 Kπ S-wave amplitude from a Dalitz plot analyses of η c →KKπ [7]. Further information on the properties of the K * 0 (1430) resonance has been obtained by the CLEO experiment in an analysis of the D + →K − π + π + decay [8], and by the BESIII experiment, which observed its decay to Kη using χ c1 decays to η K + K − [9].
The existence of gluonium states is still an open issue for Quantum Chromodynamics (QCD). Lattice QCD calculations predict the lightest gluonium states to have quantum numbers J P C = 0 ++ and 2 ++ and to be in the mass region below 2.5 GeV/c 2 [10]. In particular, the J P C = 0 ++ glueball is predicted to have a mass around 1.7 GeV/c 2 . Searches for these states have been performed using many supposed "gluon rich" reactions such as radiative decays of the heavy quarkonium states J/ψ [11,12] and Υ (1S) [13]. However, despite intense experimental searches, there has been no conclusive experimental observation [14,15]. The identification of the scalar glueball is further complicated by possible mixing with standard qq states. The broad f 0 (500), * [16], f 0 (1500) [17,18], f 0 (1710) [19,21] and possibly f 0 (2100) [20] have been suggested as scalar glueball candidates. In the BESIII partial wave analysis of the radiative J/ψ decay to ηη [20], the authors conclude that the production rates of f 0 (1710) and f 0 (2100) are both about one order of magnitude larger than that of the f 0 (1500) and no clear evidence is found for f 0 (1370). A feature of the scalar glueball is that its ss decay mode should be favored with respect to uū or dd [22,23].
In the present analysis, we consider the three-body η c decays to η K + K − , η π + π − , and ηπ + π − , using twophoton interactions, e + e − →e + e − γ * γ * →e + e − η c . If both of the virtual photons are quasi-real, the allowed J P C values of any produced resonances are 0 ±+ , 2 ±+ , 4 ±+ ... [24]. Angular momentum conservation, parity conservation, and charge conjugation invariance imply that these quantum numbers also apply to these final states. The possible presence of a gluonic component of the η meson, due to the so-called gluon anomaly, has been discussed in recent years [25,26]. A comparison of the η and η content of η c decays might yield information on the possible gluonic content of resonances decaying to π + π − or K + K − . The γγ→η π + π − process has been recently studied by the Belle experiment [27] but no Dalitz plot analysis was performed.
This article is organized as follows. In Sec. II, a brief description of the BABAR detector is given. Section III is devoted to the event reconstruction and data selection. In Sec. IV, we describe the efficiency and resolution studies, while in Sec. V we report the measurement of the η c branching fraction. In Sec. VI we describe the Dalitz plot analysis methodology, and in Secs. VII, VIII, and IX we analyze η c decays to η K + K − , η π + π − , and ηπ + π − , respectively. The results are summarized in Sec. X.

II. THE BABAR DETECTOR AND DATASET
The results presented here are based on the full data set collected with the BABAR detector at the PEP-II asymmetric-energy e + e − collider located at SLAC, and correspond to an integrated luminosity of 519 fb −1 [28] recorded at center-of-mass energies at and near the Υ (nS) (n = 2, 3, 4) resonances. The BABAR detector is described in detail in ref. [29]. Charged particles are detected, and their momenta are measured, by means of a five-layer, double-sided microstrip detector and a 40layer drift chamber, both operating in the 1.5 T magnetic field of a superconducting solenoid. Photons are measured and electrons are identified in a CsI(Tl) crystal electromagnetic calorimeter. Charged-particle identification is provided by the measurement of specific energy loss in the tracking devices, and by an internally reflecting, ring-imaging Cherenkov detector. The pions tracking efficiency increases from 98% to 100% in the momentum range 0.5-3 GeV/c while the average kaon identification efficiency is 84%. Muons and K 0 L mesons are detected in the instrumented flux return of the magnet. Monte Carlo (MC) simulated events [30], with reconstructed sample sizes of the order 10 3 times larger than the corresponding data samples, are used to evaluate the signal efficiency and to determine background features. Two-photon events are simulated using the GamGam MC generator [31]. In this article, the inclusion of chargeconjugate processes is implied, unless stated otherwise.

III. EVENT RECONSTRUCTION AND SELECTION
A. Reconstruction of the η h + h − final state We first study the reactions where h + h − indicates a π + π − or K + K − system. The selection criteria are optimized for the η c signal, as described below. The η is reconstructed in the two decay modes η →ρ 0 γ, ρ 0 →π + π − , and η →ηπ + π − , η→γγ. To reconstruct these final states we select events in which the e + and e − beam particles are scattered at small angles, and hence are undetected, ensuring that both virtual photons are quasi-real. We consider photon candidates with reconstructed energy in the electromagnetic calorimeter greater than 100 MeV. All pairs of photon candidates are combined, assuming they originate from the e + e − interaction region, and pairs with invariant-mass within ±20 MeV/c 2 (±150 MeV/c 2 ) of the neutral pion (η meson) mass are considered π 0 (η) candidates. We consider events with exactly 4 wellmeasured charged-particle tracks with transverse momentum greater than 0.1 GeV/c, and fit them to a common vertex, which must be within the e + e − interaction region and have a χ 2 fit probability greater than 0.1%. Tracks are identified as either charged kaons or pions using a high-efficiency algorithm that rejects more than half the background with negligible signal loss. A track can be identified as both kaon or pion (or neither) at this point. For the η →ρ 0 γ selection, we allow the presence of only two γ candidates, where π 0 candidates are excluded. For the η →ηπ + π − we require exactly one η candidate, no more than three additional background photon candidates, and no π 0 candidate in the event. These selections are optimized on the data using as reference the η c signal.
To reconstruct η →ρ 0 γ decays, we consider π + π − pairs in the mass region 0.620 < m(π + π − ) < 0.875 GeV/c 2 . Each of these ρ 0 candidates is combined with all γ candidates, and any combination with invariant-mass in the range 0.935 < m(ρ 0 γ) < 0.975 GeV/c 2 is considered an η candidate. We compute the angle θ γ , defined as the angle between the π + and the γ in the π + π − rest frame. The distribution of θ γ is expected to be proportional to sin 2 θ γ [32]. We thus scan the ρ 0 γ mass spectrum with varying selection on | cos θ γ | and obtain a small reduction of the combinatorial background by requiring | cos θ γ | < 0.85. The above selection reduces the η signal and background yields by 3% and 17%, respectively.
To improve the mass experimental resolution, the η four-momentum is constructed by adding the momenta of the π + , π − , and γ, and computing the η energy by assigning the Particle Data Group (PDG) [33] nominal mass. This method, tested on MC simulations, improves the resolution by ≈ 20%.
To reconstruct η →ηπ + π − decays, we perform a kinematic fit to the η candidate, and require the ηπ + π − mass to be within ±2σ of the fitted η mass (956.8 ± 0.5) MeV/c 2 , where σ = 2.9 MeV/c 2 is the width of the resolution function describing the η signal. Similarly, to improve the experimental resolution, the η four momentum is constructed by adding the momenta of the π + , π − , and η, and computing the η energy by assigning the PDG mass.
Background arises mainly from random combinations of particles from e + e − annihilation, from other twophoton processes, and from events with initial-state photon radiation (ISR). The ISR background is dominated by events with a single high-energy photon recoiling against the reconstructed hadronic system, which in the mass region of interest is typically a J P C = 1 −− resonance [34]. We discriminate against ISR events by requiring the recoil mass M 2 rec ≡ (p e + e − − p rec ) 2 > 2 GeV 2 /c 4 , where p e + e − is the four-momentum of the initial state e + e − and p rec is the reconstructed four-momentum of the candidate η (η)h + h − system.
We define p T as the magnitude of the transverse momentum of the η h + h − system, in the e + e − rest frame, with respect to the beam axis. Well reconstructed twophoton events with quasi-real photons are expected to have low values of p T . Substantial background arises from γγ→2h + 2h − events, combined with a background photon candidate. These are removed by requiring We retain events with p T below a maximum value that is optimized with respect to the η c signal for each decay mode. We produce η h + h − invariant-mass spectra with different maximum p T values, and fit them to extract the number of η c signal events (N s ) (defined as the 2.93-3.03 GeV/c 2 interval) and the number of background events underneath the η c signal (N b ). We then compute the purity, defined as P = N s /(N s + N b ), the figure of merit S = N s / √ N s + N b , and their product, P S.
1. Reconstruction of the η π + π − final state For the final selection of the η π + π − final state, we require all four charged tracks to be positively identified as pions, using an algorithm based on multivariate anal- ysis [35] that is more than 98% efficient for the tracks in the sample, while suppressing kaons by a factor of at least seven.
Figures 1(a) and 1(b) show the p T distributions for selected events in the charmonium region. This region is defined as reconstructed invariant-mass m(η (η)h + h − ) > 2.7 GeV/c 2 . In the case of η →ρ 0 γ an upper mass requirement m(η π + π − ) < 3.5 GeV/c 2 is applied because of the large number of combinations produced by the presence of the γ. The data are compared with expectations from η c signal MC simulations; a signal from two-photon production is observed in the data in both cases, and is particularly clean for η →ηπ + π − . In a scan of the S, P , and P S variables as functions of the maximum p T value, we observe a broad maximum of S starting at 0.05 GeV/c for the η →ρ 0 γ decay candidates, and a maximum of P S at 0.15 GeV/c for the η →ηπ + π − candidates. We require p T < 0.05 GeV/c and p T < 0.15 GeV/c, respectively, as indicated by the dashed lines in the figures.
Figures 2(a) and 2(b) show the ρ 0 γ and ηπ + π − invariant-mass distributions, respectively, for events sat-isfying all selection criteria except that on these masses. Clear η signals are visible, and the shaded regions indicate the selection windows, (0.935-0.975) GeV/c 2 for η →ρ 0 γ and (0.948-0.966) GeV/c 2 for η →ηπ + π − . Figures 3(a,b) show the η π + π − invariant-mass spectra for the selected events in the data. Prominent η c signals are observed, and there is some activity in the η c (2S) mass region.
If there are multiple candidates in the same event, we retain them all. The fraction of events having two combinations in the η c mass region is 3% (and 3.4% in η c signal MC simulations) for η c →η π + π − with η →ρ 0 γ. No multiple candidates are found for η →π + π − π 0 or any of the other final states discussed below. For the η K + K − final state, we require the two charged tracks assigned to the η decay to be positively identified as pions and the other two to be positively identified as kaons. The algorithm is more than 92% efficient We study the reaction where η→γγ and η→π + π − π 0 .

η→γγ
For reaction (2), where η→γγ, we again consider wellmeasured charged-particle tracks with transverse momenta greater than 0.1 GeV/c and photons with energy greater than 0.1 GeV, and each pair of γ's is kinematically fitted to the π 0 and η hypotheses. We require exactly two selected tracks, fit them to a common vertex, and require the fitted vertex to be within the interaction region and the χ 2 probability of the fit to be greater than 0.1%. We retain events having exactly one η candidate, no π 0 candidates, and no more than three background γ's.
The p T distribution for such events in the charmonium mass region is compared with η c signal MC simulation in Fig. 4(a), where a clear signal of the two-photon reaction is observed. Optimizing the η c figure of merit (S) and purity (P ), we require p T < 0.1 GeV/c. The resulting ηπ + π − invariant-mass spectrum is shown in Fig. 6(a), where the η c signal can be observed together with some weak activity in the η c (2S) mass region.
For reaction (2), where η→π + π − π 0 , we require exactly four well-measured charged-particle tracks with the vertex χ 2 fit probability greater than 0.1%. In order to have sensitivity to low momentum π 0 mesons, we consider photons with energy greater than 30 MeV/c 2 .
We allow no more than two kinematically fitted π 0 candidates and no more than five background γ's. Candidate γγ→2π + 2π − events are removed by requiring p T (2π + 2π − ) > 0.05 GeV/c. Background ISR events are removed by requiring All four charged tracks are required to be loosely identified as pions.
The η candidates are reconstructed by combining every pair of oppositely charged tracks with each of the π 0 candidates in the event. The resulting π + π − π 0 invariantmass spectrum is shown in Fig. 5. A clean η signal can be seen; we select candidates in the mass region 538 < m(π + π − π 0 ) < 557 MeV/c 2 . The η is then reconstructed by adding the momentum three-vectors of the three pions and computing the η energy using its nominal PDG mass.
The p T distribution for such events in the charmonium mass region is compared with η c signal MC simulation in Fig. 4(b), where a clear signal of the two-photon reaction is observed. In this case, a maximum of the P S figure of merit leads to the requirement p T < 0.1 GeV/c. The resulting ηπ + π − invariant-mass spectrum is shown in Fig. 6(b), where the η c signal can be observed together with some weak activity in the η c (2S) mass region.

IV. EFFICIENCY AND ηc INVARIANT-MASS RESOLUTION
To compute the reconstruction and selection efficiency, MC signal events are generated using a detailed detector simulation [30,31] in which the η c mesons decay uniformly in phase space. These simulated events are reconstructed and analyzed in the same manner as data. We define the helicity angle θ H as the angle formed by the h + (where h = π, K), in the h + h − rest frame, and the η (η) direction in the h + h − η (h + h − η) rest frame. For each final state, we compute the raw efficiency in 50×50 intervals of the invariant-mass, m(h + h − ), and cos θ H , as the ratio of reconstructed to generated events in that interval.
To smoothen statistical fluctuations, the efficiency maps are parameterized as follows. We first fit the efficiency as a function of cos θ H in each of the 100 MeV/c 2 wide intervals of m(h + h − ), using Legendre polynomials up to L = 12: where m denotes the h + h − invariant-mass. For a given value of m(h + h − ), the efficiency is interpolated linearly between adjacent mass intervals. Figure 7 shows the resulting efficiency maps (m, cos θ H ) for the four η h + h − final states, and Fig. 8 shows the maps for the two ηπ + π − final states. The small regions of very low efficiency near | cos θ H | ∼ 1 are the result of the difficulty of reconstructing K ± mesons with laboratory momentum less than ≈ 200 MeV/c, and π ± mesons with laboratory momentum less than ≈100 MeV/c, due to energy loss in the beam pipe and inner-detector material.

V. YIELDS AND BRANCHING FRACTIONS
In this section, we fit the invariant-mass distributions to obtain the numbers of selected η c events, N η K + K − , N η π + π − , and N ηπ + π − , for each η or η decay mode. We then use the η K + K − and η π + π − yields to compute the ratio of branching fractions for η c to the η K + K − and η π + π − final states. This ratio is computed as for each η decay mode, where η K + K − and η π + π − are the corresponding weighted efficiencies described in the following Sec. V B.

A. Fits to the invariant-mass spectra
We determine N K + K − η and N π + π − η from η c decays by performing binned χ 2 fits to the η K + K − and η π + π − invariant-mass spectra, in the 2.7-3.3 GeV/c 2 mass region, separately for the two η decay modes. In these fits, the η c signal contribution is described by a simple Breit-Wigner (BW) function convolved with a fixed resolution function described above, with η c parameters fixed to PDG values [33]. An additional BW function is used to describe the residual background from ISR J/ψ events, and the remaining background is parameterized by a 2 nd order polynomial. The fitted η h + h − invariant-mass spectra are shown in Fig. 9. The fits generally describe the data well, although the fit to the η K + K − invariantmass spectrum for η →ηπ + π − (Fig. 9(d)), which has low statistics, appears to the eye to have a somewhat distorted lineshape. For this fit, we add two additional parameters by leaving free the parameters of the Gaussian component of the resolution function. To minimize the dependence of the N 's on the fit quality, the η c signal yields are obtained by integrating the data over the η c signal region after subtracting the fitted backgrounds.
Statistical errors on the η c yields are evaluated by generating, from each invariant-mass spectrum, 500 new spectra by random Poisson fluctuations of the content of each bin. The generated mass spectra are fitted using the same model as for the original one and the resulting distributions of the η c subtracted yields are fitted using a Gaussian function, whose σ is taken as the statistical uncertainty. The resulting yields and χ 2 per degree of freedom for the fits, χ 2 /ndf are reported in Table I.
We test the fitting procedure by leaving free the η c parameters and find agreement, within the errors, with world averages. For the decay η c →ηπ + π − , however, the fits without interference do not describe the data well for either η decay mode. Leaving free the η c parameters, the fits return masses shifted down by ≈ 10 MeV/c 2 with respect to PDG averages. We test the possibility Parametrized detection efficiencies in the cos θH vs. m(π + π − ) plane for simulated ηc→ηπ + π − events with (a) η→γγ and (b) η→π + π − π 0 . The average value of the efficiency is shown in each interval.
of interference effects of the η c with each non-resonant two-photon process [37], modifying the fitting function by defining where A nres is the non-resonant amplitude with |A nres | 2 described by a 2 nd order polynomial; the coherence factor c is the fraction of the non-resonant events that are true two-photon production of the same final state; the resonant contribution is A ηc = α · BW (m) · exp(iφ), where BW (m) is a simple Breit-Wigner with parameters fixed to PDG values; and α, φ, and c are free parameters. The sum of f (m) and the J/ψ contribution is convolved with the experimental resolution. Fits with interference and fixed PDG parameters give values of χ 2 /ndf = 77/54 (p-value=2.2%) and The lines are the results from the fits described in the text.
46/54 (p-value=77%) for η→γγ and η→π + π − π 0 decay modes, respectively. The fitted relative phases are φ = 1.41±0.02 stat ±0.02 sys rad and φ = 1.26±0.03 stat ±0.02 sys rad. Systematic uncertainties are related to the use of η c fixed parameters and on errors in the background shape. The fits, on the other hand, show little sensitivity to the c parameter. The fitted invariant-mass spectra are shown in Fig. 10, where reasonable descriptions of the data are evident. As a comparison we also fit the two mass spectra with no interference and fixed η c parameters and obtain the dotted lines distributions shown in Fig. 10 with corresponding χ 2 /ndf = 160/55 and χ 2 /ndf = 139/55, respectively.
We find that the interference model does not produce significant improvements in the description of the data for final states that include an η . As a cross check, we reanalyze the data reported in Ref. [5], and find no evidence for such interference effects also for the η c →ηK + K − decay mode.
Systematic uncertainties on the yields due to the fitting procedure are estimated by varying the η c parameters according to the PDG uncertainties. An additional uncertainty of 4% is assigned to the yield for η c →η K + K − with η →ηπ + π − due to the variation of the resolution function. We also take the integral of each full function used to describe the η c as an estimate of the yield, and take the difference as the systematic uncertainty. The quadratic sums of these uncertainties are given in Table I.

B. Branching fractions
We estimate η K + K − and η π + π − for the η c signals using the 2-D raw efficiency functions described in Sec. IV. Each event is first weighted by 1/ (m, cos θ H ). Since the backgrounds below the η c signals have different distributions in the Dalitz plot, we perform a sideband subtraction by assigning an additional weight of +1 to events in the η c signal region, defined as the (2.93-3.03) GeV/c 2 mass region, and a weight −1 to events in the sideband regions, (2.77-2.87) GeV/c 2 and (3.09-3.19) GeV/c 2 . The weights in the sideband regions are scaled by a small amount to match the fitted η c signal/background ratio, and added to those in the signal region, to produce the weighted yields shown in Table I.
Systematic uncertainties on the efficiencies have been evaluated as follows. The uncertainty due to the limited MC statistics is computed by generating 500 new efficiency tables, obtained from the original tables by random variation, according to a Poisson distribution, of the generated and reconstructed MC yields in each cell. The distributions of the resulting weights are fitted using a Gaussian function whose σ values are taken as systematic uncertainties and are listed in Table I. To estimate an uncertainty on the method of sideband subtraction, we use the average weights in the signal region, and take the difference as an uncertainty. The quadratic sums of these uncertainties are given in Table I.
We label with R 1 (ρ 0 γ) and R 2 (ηπ + π − ) the measurements of the branching fraction for the two η decay modes. In each case, the numerator and denominator involve the same number of charged tracks and γ's, so the systematic uncertainties on their reconstruction efficiencies cancel in the ratio. The only difference is the presence of two kaons in the numerator and two pions in the denominator. The uncertainties in the particle identification efficiencies are correlated; we assign a systematic uncertainty of 1% to the identification of each kaon and 0.5% to each pion. Table II summarizes the largest systematic uncertainties on the branching fraction, which arise from MC statistics, the use of the full fitting function in extracting the yield (labelled full-BW), the sideband subtraction in the efficiencies (labelled no-sideband), and the kaon/pion identification (labelled PID).
Adding the systematic uncertainties in quadrature, we obtain the following values of the branching ratios: and an average value of

VI. DALITZ PLOT ANALYSES
We perform Dalitz plot analyses of the η π + π − , η K + K − , and ηπ + π − systems in the η c mass region using unbinned maximum likelihood fits. The likelihood function is written as where: • N is the number of events in the signal region; • f sig is the fraction of those events attributed to η c decays; • for the n-th event, x n = m 2 (η/η h + ), y n = m 2 (η/η h − ), and • (x n , y n ) is the efficiency, parameterized as a function of x n = m(h + h − ) and y n = cos θ H (see Sec. IV); • c i is the complex amplitude of the i−th signal component; the c i are free parameters of the fit; • for the n-th event, A i (x n , y n ) describe the i − th complex signal-amplitude contribution; • k i is the magnitude of the i−th background component; the k i parameters are obtained by fitting the sideband regions; • for the n-th event, B i (x n , y n ) is the probabilitydensity function of the i-th background contribution; we assume that interference between signal and background amplitudes can be ignored; I: Information for the evaluation of the branching fractions. The reported yields are obtained from the integration of the ηc signal after background subtraction in the ηc signal region. The first error is statistical, the second systematic.
Amplitudes are parameterized as described in Refs. [38] and [39]. They include a relativistic Breit-Wigner function having a variable width modulated by the Blatt-Weisskopf [42] spin form factors and the relevant spinangular information. Note that these factors are both one for scalar resonances.
The efficiency-corrected fractional contribution f i due to resonant or non-resonant contribution (NR) is defined as follows: The f i do not necessarily sum to 100% because of interference effects. The uncertainty for each f i is evaluated by propagating the full covariance matrix obtained from the fit. The search for the amplitudes contributing to the signal or background is performed by starting with the largest resonance observed in the mass projections, which is taken as the reference amplitude with c 1 = 1 and phase zero. We then add, one by one, possible processes that could contribute to the decay, testing for an increase in the likelihood value. Amplitudes are discarded if no significant improvement in the likelihood (∆(−2 log L) > 2) is obtained. Each excluded resonance is reiterated many times in combination with other possible resonant contributions. Where possible, resonance parameters are left free, for comparison with existing values; otherwise, they are fixed to PDG values. Table III summarizes the information on the structure of the samples used in the Dalitz analyses. Yields and purities are computed in the η c signal region, defined as the mass ranges (2.93-3.03) GeV/c 2 for η h + h − and (2.92-3.02) GeV/c 2 for ηπ + π − .
Each Dalitz plot analysis deals with two sets of data contributing to the given η c final state, with different efficiencies and purities: η →ρ 0 γ and η →ηπ + π − for η c →η h + h − , η→γγ and η→π + π − π 0 for η c →ηπ + π − . Therefore we use the sum of two different likelihood functions, which share the free parameters and fitting model. Due to the lack of statistics we do not separate the contributing backgrounds for the two sets of data.
VII. DALITZ PLOT ANALYSIS OF ηc→η K + K − Figure 11 shows the Dalitz plot for the selected η c →η K + K − candidates in the data, for the two η decay modes combined. Figure 12(a)-(b) shows the two squared mass projections.
We observe that this η c decay mode is dominated by a diagonal band on the low mass side of the Dalitz plot. The m 2 (K + K − ) spectrum shows a large structure in the region of the f 0 (1710) resonance. The combined m 2 (η K ± ) invariant-mass spectrum shows a structure at threshold due to the K * 0 (1430) accompanied by weaker resonant structures.
We first fit the two η c sidebands separately, using an incoherent sum of amplitudes, which includes contribu- The K * 0 (1430) is a relatively broad resonance decaying to Kπ, Kη, and Kη . The measured Kη relative branching fraction is B(K * 0 (1430)→Kη) B(K * 0 (1430)→Kπ) = 0.092 ± 0.025 +0.010 −0.025 [5], while the Kη has only been observed in Ref. [9]. To describe the K * 0 (1430) lineshape in the Kη projection, we model it using a simplified coupled-channel Breit-Wigner function, which ignores the small Kη contribution. We parameterize the K * 0 (1430) signal as where m 0 is the resonance mass, g Kπ and g Kη are the couplings to the Kπ and Kη final states, and ρ j (m) = 2P/m are the respective Lorentz-invariant phase-space factors, with P the decay particle momentum in the K * 0 (1430) rest frame. The ρ 2 (m) function becomes imaginary below the Kη threshold. The values of m 0 and the g Kj couplings cannot be derived from the Kη system only, and therefore we make use of the Kπ S-wave measurement from BABAR [7]. We average the reported quasi model-independent (QMI) measurements of the Kπ S-wave from η c →K 0 S Kπ and η c →K + K − π 0 decays, and obtain the modulus squared of the amplitude and the phase shown in Fig. 13.
We perform a simultaneous binned χ 2 fit to the Kπ S-wave amplitude and phase from threshold up to 1.72 GeV/c 2 . Above this mass, other resonant contributions are present, which make the amplitude and phase more complicated. We model the Kπ S-wave in this region as: where BW Kπ (m) is given by Eq. ( 11), B(m) is an empirical background term, parameterized as and c, φ, and α are free parameters. The results of the fit are shown in Fig. 13 as the solid (red) lines. We obtain a χ 2 /ndf = 55/31 (χ 2 /ndf = 25/31 with included systematic uncertainties) and the K * 0 (1430) parameters listed in Table IV. We note a large statistical error on g 2 Kη that is expected because of the weak sensitivity of the Kπ S-wave to the opening of the Kη threshold. We also note the presence of a very small background term. We attempt to replace the background term with a BW function with parameters fixed to the PDG averages for the κ/K * 0 (700) resonance, but obtain a poor description of the data. For comparison, the K * 0 (1430) parameters used by BESIII in the Dalitz plot analysis of χ c1 →η K + K − [9] are those measured by the CLEO D + →K − π + π + Dalitz plot analysis [8], m = 1471.2 MeV/c 2 , g 2 Kπ = 0.299 GeV 2 /c 4 , and g 2 Kη = 0.0529 GeV 2 /c 4 . We perform a Dalitz plot analysis of the η c →η K + K − decay channel by using the η f 0 (1710) intermediate state as the reference amplitude. If there are regions of the phase space not well described by the fit, we add postulated K + K * − 0 , η f 0,2 , or η a 0 intermediate states, and accept them if ∆(−2 log L) > 2. At each stage, we test for the presence of a non-resonant contribution.
We describe the K * 0 (1430) according to Eq. (11) first with m 0 and g 2 Kπ parameters fixed to the values from the fit to the Kπ S-wave and g 2 Kη free. We observe little sensitivity to the g 2 Kη parameter, expressed by the large error, and therefore we also fix the value of this parameter to that from the fit to the Kπ S-wave.
The projections of the fit result are shown in Fig. 12, along with the largest signal components. To test the fit quality, we generate a large number of phase-space MC-simulated events, which are weighted by the likelihood function obtained by the fit. These MC-simulated events are then normalized to the observed yield and are superimposed to the data. To test the fit quality we also project the fit on the (m(K + K − ), cos θ H ) plane and compare data and simulation in each cell of the plane. Labelling with ndf = N cells − N par , where N cells is the number of cells having at least two expected events and N par the number of free parameters in the Dalitz analysis, we obtain χ 2 /ndf = 285/264 = 1.1 corresponding to a p-value of 18%.  Table V, together with their fitted fractions and relative phases. We label this fit as solution (A). The non-resonant contribution is consistent with zero.
We measure the f 0 (1710) parameters, listed in Table IV. In addition to the strong f 0 (1710)η and K * 0 (1430) + K − contributions there is evidence for a signal of the K * 0 (1950) + K − decay mode. We measure the parameters of the K * 0 (1950) (see Table IV) for which there is only one previous measurement from the LASS collaboration [40]. There are smaller contributions from f 0 (980)η , f 2 (1270)η , and f 0 (1510)η . The latter is indistinguishable from an f 2 (1525)η contribution, but for simplicity, we report only the f 0 (1510)η , which gives a slighly larger likelihood improvement.
Statistical significances of resonances contributing to the decay are evaluated using the Wilks' theorem [41] from the difference in log likelihood between fits with and without the specific signal component, taking into account the difference of two free parameters. For f 0 (1710)η and K * 0 (1950) + K − we obtain ∆(−2 log L) = 135.9 and ∆(−2 log L) = 15.3, respectively. The corresponding significances are listed in Table IV. We evaluate systematic uncertainties on the fitted fractions, phases, and resonance parameters. For resonances having parameters fixed to PDG values, we vary these parameters according to their PDG uncertainties. We modify the purity of the η c signal according to its statistical uncertainty. We replace the fitted efficiency with the raw efficiency, defined in Sec. IV. The Blatt-Weisskopf [42] form factor present in the relativistic BW functions, nominally fixed at 1.5 GeV −1 , is varied between 0 and 3.0 GeV −1 . The background description is modified by varying each resonant fraction by its statistical uncertainties in the fits to the sidebands. All the contributions are added in quadrature.
(16) This ratio can be written as where I Kη and I Kπ are the integrals over the η c phase space of the coupled-channel Breit-Wigner function describing the K * 0 (1430) in the η c →η K + K − and η c →π 0 K + K − decay modes (Eq. (11)). Using Eq. (17), we obtain the ratio of the couplings This ratio depends weakly on the resonance parameters, varying between 1.40 to 1.67. Therefore, we fix where we have included the change from solution (A) in the systematic uncertainty, as an estimate of the model uncertainty. Similarly, we use the estimates of the K * 0 (1430) mass and g 2 Kπ from solution (B), along with the differences from solution (A) (see Table IV), to obtain The inconsistency between the VIII. DALITZ PLOT ANALYSIS OF ηc→η π + π − Figure 15 shows the Dalitz plot for the selected η c →η π + π − candidates in the data, in the η c signal region, for the two η decay modes combined, and Figs. 16(a)-(b) show two squared-mass projections. We observe several diagonal bands in the Dalitz plot, in particular at the lower-left edge. There are corresponding structures in the m 2 (π + π − ) spectrum, including peaks attributable to the f 0 (980) and f 2 (1270) resonances, and a large structure at high π + π − mass. In the m 2 (η π ± ) spectrum, a large structure is present; there is no known resonance decaying to η π in this mass region, but this could be a reflection of the structure in the high m 2 (π + π − ) region.
We fit the two η c sidebands using an incoherent sum of amplitudes, which includes contributions from the ρ 0 (770), f 2 (1270), f 0 (1370), and f 0 (2100) resonances. To model the background in the η c signal region we take a weighted average of the fitted fractional contributions, and normalize using the results from the fit to the η π + π − invariant-mass spectrum.   A candidate for the large structure in the high π + π − mass region is the f 0 (2100) resonance, observed in radiative J/ψ decay to γηη [20]. We take f 0 (2100)η as the reference contribution, and perform a Dalitz plot analysis as described in Sec. VI. Again, no non-resonant contribution is needed, and the list of the resonances contributing to this η c decay mode is given in Table VI, together with their fitted fractions and relative phases.
The f 0 (2100) parameters are first left free in the fit, and we obtain the values listed in Table IV, which are in agreement with BESIII measurement (m = 2081 ± 13 +24

−36
MeV/c 2 , Γ = 273 +27+70 −24−23 ) MeV [20]. We then fix them to the values listed in the PDG. We also leave free the f 0 (500) parameters and obtain the values listed in Table IV which give a good description of the data. Given the low statistics, we do not assign systematic uncertainties to the fitted f 0 (500) resonance parameters, which are within the range of other measurements [33]. The f 0 (980) is parameterized by a coupled-channel Breit-Wigner function with parameters fixed to the measurement from Ref. [44]. To describe the small enhancement around 1.43 GeV/c 2 , we test both spin-2 and spin-0 hypotheses with free resonance parameters; we obtain ∆(−2 log L) = 2.4 in favor of the spin-2 hypothesis, so we attribute this signal to the f 2 (1430) resonance, and report the fitted parameter values in Table IV. We test the significance of this signal by removing it from the list of the resonances, obtaining ∆(−2 log L) = 23.8 and a significance of 4.4σ. Replacing the f 2 (1430) resonance with f 0 (1500) or f 0 (1370), we obtain poor fits with fractions from these possible contributions consistent with zero. The f 0 (2100) statistical significance is 10σ.
The projections of the fit result are compared with the data in Fig. 16. To test the fit quality, we generate a large number of phase-space MC-simulated events, which are weighted by the likelihood function obtained from the fit. These MC-simulated events are then normalized to the observed yield and superimposed to the data. We also project the fit on the (m(π + π − ), cos θ H ) plane and compare data and simulation in each cell, obtaining χ 2 /ndf = 409/386 = 1.1. The systematic uncertainties on the fitted fractions, phases and resonance parameters are evaluated as in the previous section.
IX. DALITZ PLOT ANALYSIS OF ηc→ηπ + π − Figure 17 shows the Dalitz plot for the selected η c →ηπ + π − candidates in the data, in the η c signal region, for the two η decay modes combined, and Figs. 18(a)-(b) show two squared-mass projections. We observe that the Dalitz plot is dominated by horizontal and vertical bands due to the a 0 (980) and diagonal bands due to resonances in the π + π − final state. The squaredmass projections show signals of f 0 (500), f 0 (980), and f 2 (1270).
The η c sidebands are also rich in resonant structure, and are fitted using an incoherent sum of amplitudes, including contributions from the a 0 (980), f 2 (1270), a 2 (1310), and f 2 (1950) resonances. We take a weighted average of the fitted fractions in the two sidebands, normalized using the results from the fit to the ηπ + π − invariant-mass spectrum, to estimate the background in the signal region, shown as the shaded regions in Figs. 18(a)-(b).
We take a 0 (980) + π − as the reference contribution, and perform a Dalitz plot analysis as described above. The resulting list of contributions to this η c decay mode is given in Table VII, together with fitted fractions and relative phases.
We find little sensitivity to the parameters of the f 0 (500) resonance, and therefore we use the parameters from the η c →η π + π − Dalitz plot analysis, listed in Table IV. A new a 0 (1700) resonance is observed in the ηπ ± invariant-mass spectrum, with fitted parameters listed in  onance is excluded from the fit is ∆(−2 log L) = 72.3, corresponding to a significance greater than 8σ. Possible contributions from the a 2 (1710) and f 0 (2100) resonances have been tested, but both are found to be consistent with zero.
We note the presence of a very large non-resonant scalar contribution, and in Table VII, we list both the sum of resonant contributions and the sum including the non-resonant contribution. A similar effect has been observed in charmless B decays [45]. This effect could be correlated with the interference of the η c with the twophoton continuum described in Sec. V.
We test the fit quality as described above, with the comparison in the (m(π + π − ), cos θ H ) plane giving χ 2 /ndf = 419/382 = 1.1. We evaluate systematic uncertainties as described above but adding an additional uncertainty due to the possible interference between intermediate resonances from the η c decay and those present in the background. To obtain the order of magnitude of the effect we compare the fits to the ηπ + π − mass spectra described in Sec. VI.A with and without the interference and obtain an average difference in the η c yield of the order of 26%. Multiplying this factor by the sum of all the resonant fractions given in Table VII, we obtain an estimate of the uncertainty of the order of 15% which is added in quadrature to the other sources of systematic uncertanties. We also vary the η c signal region width from 100 MeV/c 2 to 60 MeV/c 2 and add in quadrature the resulting differences in amplitudes fractions and phases as an additional source of systematic uncertainties.  We study the processes γγ→η K + K − , γγ→η π + π − , and γγ→ηπ + π − using a data sample of 519 fb −1 recorded with the BABAR detector operating at the SLAC PEP-II asymmetric-energy e + e − collider at center-of-mass energies at and near the Υ (nS) (n = 2, 3, 4) resonances. We observe η c decays to all the above final states and perform Dalitz plot analyses to measure intermediate resonant fractions and relative phases. Significant interference effects of the η c with the two-photon background are observed only for the decay η c →ηπ + π − .
We observe an enhanced contribution of f 0 (1710) in η c decays to η and an enhanced contribution of f 0 (1500) in η c decays to η. This effect may point to an enhanced gluonium content in the f 0 (1710) meson. A similar conclu-sion is drawn in the study of J/ψ radiative decays [20]. In particular, Ref. [21] finds that the production rate of the pure gauge scalar glueball in J/ψ radiative decays predicted by lattice QCD is compatible with the production rate of J/ψ radiative decays to f 0 (1710) and this suggests that f 0 (1710) has a larger overlap with the glueball compared to other glueball candidates (e.g., f 0 (1500)). The observation of f 0 (2100) in both J/ψ radiative decays and in η c →η π + π − allows to add this state in the list of the candidates for the scalar glueball.

XI. ACKNOWLEDGMENTS
We are grateful for the extraordinary contributions of our PEP-II colleagues in achieving the excellent luminosity and machine conditions that have made this work possible. The success of this project also re-