Observation of $\eta_c\to\omega\omega$ in $J/\psi\to\gamma\omega\omega$

Using a sample of $(1310.6\pm7.0)\times10^6$ $J/\psi$ events recorded with the BESIII detector at the symmetric electron positron collider BEPCII, we report the observation of the decay of the $(1^1 S_0)$ charmonium state $\eta_c$ into a pair of $\omega$ mesons in the process $J/\psi\to\gamma\omega\omega$. The branching fraction is measured for the first time to be $\mathcal{B}(\eta_c\to\omega\omega)= (2.88\pm0.10\pm0.46\pm0.68)\times10^{-3}$, where the first uncertainty is statistical, the second systematic and the third is from the uncertainty of $\mathcal{B}(J/\psi\to\gamma\eta_c)$. The mass and width of the $\eta_c$ are determined as $M=(2985.9\pm0.7\pm2.1)\,$MeV/$c^2$ and $\Gamma=(33.8\pm1.6\pm4.1)\,$MeV.

Using a sample of (1310.6 ± 7.0) × 10 6 J/ψ events recorded with the BESIII detector at the symmetric electron positron collider BEPCII, we report the observation of the decay of the (1 1 S0) charmonium state ηc into a pair of ω mesons in the process J/ψ → γωω. The branching fraction is measured for the first time to be B(ηc → ωω) = (2.88 ± 0.10 ± 0.46 ± 0.68) × 10 −3 , where the first uncertainty is statistical, the second systematic and the third is from the uncertainty of B(J/ψ → γηc). The mass and width of the ηc are determined as M = (2985.9 ± 0.7 ± 2.1) MeV/c 2 and Γ = (33.8 ± 1.6 ± 4.1) MeV.

I. INTRODUCTION
Although the η c was discovered already in 1980 [1], the properties of the lowest lying S -wave spin singlet charmonium state are still under investigation. Especially when considering the available data on the branching fractions (BFs) of different decay modes of the η c , it becomes obvious that this resonance is not fully understood yet. Several BFs are only measured very roughly or with large uncertainties, and the observed BFs sum up to only about 57%. Several peculiarities also arise, when the resonance parameters of this meson are studied in detail: The observed mass and decay width seem to vary by a large fraction from experiment to experiment, and also seem to be dependent on the production, and/or decay process in which they are observed. While the decay of the η c into a pair of φ mesons has been observed before (see e.g. Refs. [2], [3]), only an upper limit for the decay into two ω mesons has been set [4]. Apart from these measurements, the Belle experiment was able to determine the product BF B(γγ → η c ) × B(η c → ωω) [5]. The decay η c → 2(π + π − π 0 ), which should also contain a large fraction of the ωω channel, has been determined to be one of the strongest decay modes of the η c [6]. Predictions for the BFs of the η c into a pair of vector mesons have been recently published [7]. However, the predicted BFs for the decay modes η c → φφ and η c → ρρ are much smaller than those observed experimentally. The predictions are based on Next-to-Leading order (NLO) perturbative Quantum Chromodynamics (QCD) calculations and for the first time also include so-called highertwist contributions. It was found that these contributions do have a major impact on the BFs and lead to much larger values than expected from pure perturbative QCD. However, the effect is not strong enough to explain the experimentally determined BFs for the φφ and ρρ channels. The predictions for the BF of the η c → ωω process range from 9.1 × 10 −5 to 1.3 × 10 −4 , while the most sensitive experimental determination yielded an upper limit of < 3.1 × 10 −3 at the 90% confidence level [4].
In this paper we present the first measurement of the BF for the decay η c → ωω, where the η c is observed in the invariant mass of two ω mesons produced in the radiative decay J/ψ → γωω. The data set used for this analysis contains a total of (1310.6±7.0)×10 6 J/ψ events [8] produced in direct e + e − annihilations and recorded with the Beijing Spectrometer III (BESIII) detector. The mass, the width, and the yield of the η c signal are determined by means of a partial wave analysis (PWA) in the η c signal region to properly account for interference effects with other contributions to the ωω system.

II. DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector [9] is located at the electron positron collider BEPCII [10] at the Institute for High Energy Physics (IHEP), Beijing, China. The symmetric double-ring collider BEPCII provides a peak luminosity of 10 33 cm −2 s −1 at a center-of-mass energy of 3.77 GeV. The detector consists of four main components: A smallcell gas drift chamber with 43 layers directly surrounds the beam pipe. This main drift chamber (MDC) is filled with a 60% He, 40% C 3 H 8 gas mixture. It provides an average single-hit position resolution of 135 µm as well as a charged particle momentum resolution of 0.5% (0.6%) at 1 GeV/c in a 1 T (2009) or 0.9 T (2012) magnetic field, which is generated by a superconducting solenoid magnet. The dE/dx resolution of the MDC is 6% for electrons from Bhabha scattering. Surrounding the drift chamber, a plastic scintillator based time-of-flight system (TOF) for particle identification followed by a CsI(Tl)based electromagnetic calorimeter (EMC) is mounted. The EMC consists of 6240 crystals arranged in a cylindrical, barrel-shaped part and two end caps. The calorimeter provides an energy resolution of 2.5% (5%) for 1 GeV photons as well as a position resolution of 6 mm (9 mm) in the barrel (end caps). The time-of-flight system consists of 176 scintillator bars with a length of 2.4 m, arranged in a two-layer, barrel-shaped geometry and 96 fan-shaped scintillators in the end caps. All plastic scintillators of the time-of-flight system have a thickness of 5 cm. The system provides a K/π separation of 2σ for momenta up to ∼ 1 GeV/c with a time resolution of 80 ps (110 ps) in the barrel (end caps). The iron return yoke of the solenoid magnet is instrumented with 9 (8) layers of resistive plate chambers in the barrel (end cap) regions, yielding in total about 1272 m 2 of active area. The signals from these chambers can be used for muon identification with a position resolution of 2 cm.
Phase-space distributed Monte Carlo (MC) data sets of the signal channel are generated for optimizations of the event selection over the complete phase-space (26M events) as well as the minimization in the PWA containing only events in the η c mass range (2M events). The simulations are carried out using a Geant4-based simulation software, which includes a precise description of the BESIII geometry and material, the detector response and digitization models, as well as the detector running conditions and performance. The production of the J/ψ resonance is simulated by the MC generator KKMC [11]. The subsequent decay of the J/ψ into a radiative photon and a pair of ω mesons, as well as the three-body decays of the ω mesons into π + π − π 0 are generated using Be-sEvtGen [12], which is based on the EvtGen package [13].

III. EVENT SELECTION
We perform an exclusive reconstruction of the decay J/ψ → γωω, where both ω mesons are reconstructed in their decay into π + π − π 0 . Both π 0 mesons decay further into a pair of photons, thus yielding the final state π + π − π + π − 5γ. Candidate events are required to contain two pairs of oppositely charged tracks and at least five photon candidates.
Tracks of charged particles are reconstructed using the hit information from the MDC. A track is accepted as a charged particle candidate if the distance between the point of closest approach and the interaction point is smaller than 1 cm in the plane perpendicular to the beam and smaller than 10 cm in the beam direction. Furthermore, each track is required to be within the angular acceptance of the MDC, fulfilling the requirement on the polar angle | cos θ| < 0.93.
Pion candidates are selected from all good charged tracks, by exploiting the capabilities of particle identi-fication of the different subdetector systems. Using the information on the energy loss dE/dx measured with the MDC, as well as the information from the time-of-flight system, a likelihood is calculated under the hypotheses that the particle candidate under investigation is a pion (L(π)), kaon (L(K)) or proton (L(p)). Only candidates fulfilling the criteria L(π) > L(K) and L(π) > L(p) are accepted and retained for further analysis.
Photon candidates are showers detected with the EMC exceeding an energy of 25 MeV in the barrel (| cos θ| < 0.8) and 50 MeV in the end cap regions (0.86 < | cos θ| < 0.92), respectively. To reject photons originating from split-off effects, each photon candidate must lie outside a cone with an opening angle of 20 • around the impact point in the calorimeter of any charged track. Furthermore, photon candidates are only accepted if their hit time is within 700 ns of the event start time to suppress electronic noise and showers that are unrelated to the event.
To improve the momentum resolution of the ω candidates, suppress background and determine the correct combination of photons to form π 0 candidates, all events are kinematically fitted under the J/ψ → γπ + π − π 0 π + π − π 0 hypothesis for all possible combinations of photons. The fit is performed using six kinematic constraints, which are the energy and the three linear momentum components of the initial e + e − system, as well as the masses of the two π 0 candidates. The combination which yields the smallest χ 2 6C value for the kinematic fit is chosen and the event is kept for further analysis, if χ 2 6C < 25. This effectively reduces photon mis-combination. Finally the correct combination of two sets of three pions to form the two ω candidates must be found. The three pions are assigned to the ω candidate, for which they exhibit the closest Euclidean distance r from the nominal mass of the ω meson, given by Here, m(ω) indicates the nominal mass of the ω meson as listed in Ref. [14]. Figure 1 shows the 3π versus 3π invariant mass for all events retained after the selection procedure described above.
Two bands originating from the process J/ψ → γω3π, located at the nominal ω mass, are clearly visible in Fig. 1. Additionally, a flat, homogeneous background corresponding to J/ψ → γ6π events is visible. Events from both of these processes are also present under the clearly visible enhancement at the intersection of the two ω bands. To remove this type of background, an eventbased background subtraction method is used, which is described in the following section. After application of the background subtraction, a strict selection requirement around the intersection of the two bands is introduced. Distribution of the invariant masses of both three-pion systems appearing in the decay J/ψ → γ(π + π − π 0 )1(π + π − π 0 )2 for the chosen best combination of each event. The bands correspond to the mass of the ω meson; a clear enhancement at the intersection of the two bands is visible. The red circle indicates the signal region which is selected after application of the background subtraction method described in section IV.

IV. BACKGROUND SUBTRACTION
A sophisticated event-based method for background subtraction proposed in Ref. [15] is applied to events for which both three-pion invariant masses are located within a range of ±80 MeV around the nominal ω mass. Simpler methods, such as a two-dimensional side band subtraction, mostly require the analysis of a binned data set, while the goal here is to perform a PWA and thus an event-based method is preferred.
The method is based on analyzing the signal-tobackground ratio Q in a very small cell of the available phase-space around each event. Therefore a distinct kinematic variable is needed, for which parameterizations of both the signal and background shape are known for the events in these small cells. The first step is to assign a number of N nearest neighbors for each event, denoted as seed event. In order to measure distances between events, a metric has to be defined using the kinematic observables that span the phase space for the reaction. For this analysis, in total nine coordinates are used for the metric: the polar angle of the radiative photon in the J/ψ rest frame, where the z-axis is defined by the direction of the incoming positron beam, the angle between the two ω candidates' decay planes in the J/ψ rest frame, the invariant mass of the 2(π + π − π 0 ) system, the azimuthal and polar decay angles of the two ω candidates in the helicity frame of the corresponding ω candidate, as well as the two normalized slope parametersλ of the ω candidates' decays. The parameterλ characterized by the cross product of the two pion momenta in the corre- sponding ω candidates' helicity frame is given as and , where T π denotes the kinetic energy of the corresponding pion [16] and c is the speed of light. The parameter λ takes its maximum value λ max for totally symmetric decays with an angle of 120 • between any pion pair (see Ref. [16]). The distance between two events is given by the Euclidean distance considering all coordinates listed above.
For this analysis the two-dimensional m(3π) 1 versus m(3π) 2 distribution was chosen as the distinct kinematic variable. The signal is described with a two-dimensional Voigtian function, which is defined as the convolution of a Gaussian with a Breit-Wigner function, while the background consists of two different contributions: A two-dimensional linear function with individual slope parameters for the two 3-pion invariant masses is used to describe the homogeneous background. Additionally, the ω bands are described with a Voigtian function for the one, and a linear function for the corresponding other 3π invariant mass. These functional dependencies are determined using signal MC samples. Figure 2 (a) shows the 3π versus 3π distribution for the N = 200 nearest neighboring events of a seed event, while Fig. 2 (b) shows the function fitted to this data. The value of N should be as small as possible to ensure that the phase space cell of all selected neighbors is small and the assumption that the background behaves smoothly within the cell is satisfied, yet it has to be large enough to ensure stable and reliable single-event fits. The value is determined based on dedicated MC studies for this analysis by increasing N until stable fits are achieved. The MC samples are generated using an amplitude model obtained from a PWA fit so that all angular and invariant mass distributions of the recorded data are reproduced. The signal-to-background ratio at the location of the seed event is extracted from each single-event fit and represents the Q-factor for this event.
To illustrate the quality of these fits, the projections of fit function and data from Fig. 2 (a) to each of the 3π axes is shown in the sub-figures (c) and (d), where a good agreement can be seen. Figure 3 shows the invariant 3π mass and the normalizedλ distribution for all pre-selected events, as well as the distributions weighted by Q and (1 − Q) (both diagrams contain two entries per event, one for each ω candidate). The Q-weighted diagrams show a background-free ω signal and a linearly increasingλ distribution, starting at the origin, as it is expected for a pure ω signal.
The (1−Q)-weighted distributions contain background due to events without any intermediate ω resonances (linear shape in 3π invariant mass, flat distribution of λ), as well as events that only contain one instead of two ω mesons. The latter create a peaking structure in the invariant 3π mass as well as a slight increase of the (1 − Q)-weightedλ distribution. After all single-event fits are performed, the initially very large mass window for the ω candidates, which is needed to be able to fit the background component underneath the ωω signal, is replaced with a tighter requirement of 26 MeV around the two nominal ω masses, as indicated by the red circle in Fig. 1. Figure 4 shows the invariant ωω mass for the finally selected events within this narrow signal region without any weight, Q-weighted and (1 − Q)-weighted. In total 5128 events are selected in the signal region defined as m(ωω) ≥ 2.65 GeV/c 2 and with all other selection criteria applied as discussed above. The sum of the obtained Q-factors for these events yields 4489.31, so that about 12.5% of the initially selected events originate from background sources and are weighted out by the Qfactor method. All further analysis steps are performed using this weighted data sample. A strong signal of the η c is observed in this mass distribution.
The performance of the background suppression method is checked by selecting events from side-band regions in the 3π versus 3π mass distribution. A very good agreement between expectations from the side bands and the (1 − Q)-weighted data is found. This underlines the applicability of the method. Additionally, as a crosscheck and for tuning parameters like the number of neighbors, input-output checks are performed using different MC samples generated with amplitude models obtained from rough fits to the signal and sideband regions. Using the Q-factor method, the generated signal and background samples can be identified clearly and the remaining deviation from the generated sample is taken as a systematic uncertainty of the method.

V. DATA ANALYSIS
We use a PWA to determine the number of η c candidates and the selection efficiency respecting all dimensions of the phase space simultaneously for the reaction under investigation. The amplitudes are constructed in our software [17] using the helicity formalism by describing the complete decay chain from the initial J/ψ state to the final state pions and photons. We assume that there are no other resonances nearby and thus the selected γωω events are described either as originating from the decay of the η c , or as phase space-like contributions with different J P quantum numbers of the ωω system, to consider tails of resonances that are located far away from the region of interest. For the amplitudes that describe the radiative decay of the J/ψ, an expansion into the electromagnetic multipoles of the radiative photon is applied. The decay of the η c as well as the phase space-like contributions are described using an expansion of the corresponding helicity amplitudes into the LS-scheme, where L denotes the orbital angular momentum between the two decay products and S their total spin.

A. Amplitude Model
The differential cross section of the reaction under study is expressed in terms of the transition amplitudes for the production and decay of all intermediate states and is given as Here, dΩ denotes an infinitesimally small element of the phase-space, and the function w is the transition probability from the initial to the final state. The outer (incoherent) sum runs over the helicity of the radiative photon, λ γ , as well as the z-component of the spin of the J/ψ, denoted with M . Furthermore, for all intermediate states X, a coherent summation over the helicity of the state (λ X ) as well as its daughter particles (λ ω1 , λ ω2 ) is performed. In this expression, X denotes the phase spacelike contributions with spin-parity J P , as well as the resonant η c component. The amplitudes for the J/ψ → γX process are given by where d denotes the Wigner d-matrices as defined in Ref. [14]. The d-matrices do not depend on the azimuthal angle ϕ in contrast to the usual Wigner D-matrices. The ϕ dependence vanishes for the J/ψ decay amplitudes, since both the electron and the positron beams are unpolarized. F represents the complex helicity amplitude, which is then expanded into radiative multipoles related to the corresponding final state photon using the transformation as given in Refs. [18] [19] [20], where ... denotes the Clebsch-Gordan coefficients and B L (q) are the Blatt-Weisskopf barrier factors as defined in Ref. [21]. Here, q is the linear momentum of one of the decay products in the J/ψ rest frame. q 0 is chosen as the breakup momentum for the X system and to coincide with the ωω mass threshold. Since the orbital angular momentum L between the decay products is not defined in the multipole basis, we use the minimal value L min depending on the spin-parity of X, which is expected to represent the dominant contribution. Due to this transformation, the helicities are replaced by a description based on the angular momentum J γ carried by the radiative photon. This way, the single terms of the expansion can be identified with electric or magnetic dipole, quadrupole and octupole transitions. The decay amplitudesÃ are given bỹ For these amplitudes an expansion into states with defined sets of J P C , L, S values is performed using the transformation where S is the total spin of the ωω system [22]. Also here, the normalized Blatt-Weisskopf factors are included as defined above. For the η c component, the break-up momentum q 0 is chosen to coincide with the nominal mass of the η c , while for all other contributions the ωω mass threshold is used. Since we assume that no resonances apart from the η c are nearby, the description of the dynamical part of the amplitudes for the phase space-like components (e.g. Breit-Wigner function) is omitted. For the line shape of the η c a modified relativistic Breit-Wigner function is used that takes the distortion due to the pure magnetic dipole transition J/ψ → γη c into account. The amplitude is modified by a factor E 3/2 γ , which originates from the M 1-transition matrix element [23] and corresponds to the expected E 3 γ dependency of the observed line shape. Since this factor leads to a good description around the pole mass but also introduces a diverging tail towards larger energies of the radiative photon (smaller invariant ωω masses), the amplitude for the η c is further modified using an empirical damping factor exp − E 2 γ 16β 2 with β = 0.065 GeV, in accordance with the factor used by the CLEO collaboration [24].
The decay amplitudes A of the ω resonances are directly proportional to the parameterλ introduced in Eq. (2). The normal vector n to the ω decay plane spanned by the three daughter particles in its helicity frame is described in terms of the Euler angles ϑ n , ϕ n and γ n = 0. With µ = J ω · n being the projection of the ω mesons spin to the direction of n, the amplitude reads as A Jω λω (ω → π + π − π 0 ) = 3 4π · D 1 * λωµ (ϕ n , ϑ n , 0) ·λ µ , (8) where only the case µ = 0 is allowed for this decay [25]. The free parameters varied in the minimization are the complex values a Jγ and α J X LS , as well as the mass and width of the η c . Symmetries arising from parity conservation and the appearance of two identical particles (ωω) are respected and lead to a reduction of free parameters in the fit.

B. Fit procedure
Unbinned maximum likelihood fits are performed for all hypotheses, in which the probability function w is fitted to the selected data by varying the free parameters given by the complex amplitudes as well as the masses and widths, if applicable. Each amplitude can be expressed by a real magnitude and a phase, yielding two distinct fit parameters per amplitude. The likelihood function is given by [25] where N denotes the number of data events, n is defined as Ω is a vector of the phase-space coordinates and α of the complex fit parameters. The function w( Ω, α) is the transition probability function given in Eq. (3) and ( Ω) is the acceptance and reconstruction efficiency at the position Ω.
The function w is interpreted as a probability density function and the corresponding probabilities for all events are multiplied to obtain the total probability. A normalization of the extended likelihood function is achieved due to the exponential term in which n appears, so that the mean weight of an MC event is approximately 1 after the likelihood has been maximized. The integrals appearing in the n term as well as the denominator in the product in Eq. (9) are approximated using reconstructed, phase space distributed MC events. The events of the MC sample are propagated through the BESIII detector, reconstructed and selected with the same cuts as the data sample to account for the geometrical acceptance and selection efficiency in all dimensions of the phase-space.
The best description of the data sample is reached upon maximization of the likelihood L. Equation (9) is logarithmized so that the product is transformed into a sum. Finally, the event weights Q i obtained from the Qfactor method are also included in the likelihood function and a negative sign is added to the logarithmized function, so that commonly used minimizers and algorithms, in this case Minuit2 [26], can be used.
The negative log-likelihood function, which is actually minimized, now reads as (11)

C. Fit strategy
Since the composition of the non-resonant contribution is not known a priori, different hypotheses are fitted to the selected data set, which contain the η c component and one up to a maximum of four different non-resonant components. These non-resonant components are assumed to have the J P quantum numbers 0 − , 0 + , 1 + or 2 + , so that the most simple hypothesis is given as {η c , 0 − }, and the most complex one by {η c , 0 − , 0 + , 1 + , 2 + }. We also perform fits including higher spin contributions (J P = 4 + ) and the contribution of a spin-4 component is found to be not significant. Similarly, fits with contributions carrying exotic quantum numbers (e.g. J P C = 1 −+ ) as well as pseudo tensor contributions (J P C = 2 −+ ) are tested and found to be insignificant.
In order to be able to compare the quality of fits with different, generally not nested, hypotheses with different numbers of free parameters, two information criteria from model selection theory are utilized. The Bayesian Information Criterion (BIC) depends on the maximized value of the likelihood L, the number of free parameters k as well as the number of data points n, which is given by the sum of the Q-factors. It is defined as The BIC is based on the assumption that the number of data points n is much larger than the number of free parameters k [27]. This assumption is fulfilled for all fits performed here. The second criterion is the Akaike Information Criterion (AIC), which provides a different penalty factor compared to the BIC. It is defined as thus it is independent from the sample size n. In comparison to the BIC, the penalty term is much weaker, which increases the probability of over-fitting. Theoretical considerations show [27] that in general AIC should be preferred over BIC due to reasons of accurateness as well as practical performance.
As for the likelihood, also for BIC and AIC a more negative value indicates a better fit. The results for the five best hypotheses are listed in Table I. The overall best hypothesis is determined to be for which 21 parameters are free in the fit. The fit parameters are composed of the complex decay amplitudes a Jγ for the process J/ψ → γX in the radiative multipole schema, as well as the X → ωω decay amplitudes α J X LS after the transformation to the LS-scheme and the mass and width of the η c . Each complex decay amplitude yields two independent fit parameters (magnitude and phase), whereas the phase parameter for the J/ψ → γη c amplitude is fixed to zero as a global reference. Additionally, one magnitude and one phase parameter are fixed for the X → ωω decay amplitudes for each of the four fit contributions to obtain a set of independent parameters. A projection of this fit to the ωω invariant mass and other kinematically relevant variables is shown in Figs. 5 and 6. These figures also show efficiency-corrected versions of all mass spectra and angular distributions. The correction is performed using the PWA software and is therefore done in all dimensions of the phase-space simultaneously. The fit yields a total of 1705 ± 58 η c events, which is the number used for the calculation of the BF. The yields of all components are listed in Table II.
To estimate the overall goodness-of-fit, a global χ 2 value is calculated by comparing the histograms for data and fit projections in all relevant kinematic variables as defined for the metric used for the Q-factor background subtraction method (see Section IV). The global reduced χ 2 is calculated as where N data ij and N fit ij are the contents of the jth bin in the ith kinematic variable for data and fit histograms, respectively. The bin contents themselves are given by the sum of weights of the events for data (Q-weights) as well as fit (weights from the PWA fit) histograms. Accordingly, σ data ij and σ fit ij represent the corresponding sum of squared weights to account for the bin error in the weighted histograms. N bins is the sum of all bins considered and N params is the number of free parameters in the PWA fit. Bins with less than 10 effective events are merged with neighboring bins. For the best fit hypothesis H 0 , a value of χ 2 /ndf = 640/(609 − 21) = 1.09 is obtained, which indicates a good quality of the fit.

VI. SYSTEMATIC UNCERTAINTIES
Various sources of systematic uncertainties for the determination of the BF, the mass and the width of the η c are considered. The uncertainties arise from the reconstruction and fit procedure, background subtraction method, external BFs, kinematic fit, parameterization of the η c line shape and the number of J/ψ events in our data sample.
a. Number of J/ψ events Inclusive decays of the J/ψ are used to calculate the number of J/ψ events in the data sample used for this analysis. The sample contains (1310.6 ± 7.0) × 10 6 J/ψ decays, where the uncertainty is systematic only and the statistical uncertainty is negligible [8]. The uncertainty propagates to a systematic uncertainty on the η c → ωω BF of 0.5%.
b. Photon detection The detection efficiency for photons is studied using the well understood process J/ψ → π + π − π 0 . A systematic uncertainty introduced by the photon reconstruction efficiency of < 1% per photon is found. The systematic uncertainty for the reconstruction of the five signal photons in this analysis thus is conservatively taken to be 5%. c. Track reconstruction For the estimation of the systematic uncertainty arising from the reconstruction of charged tracks and the identification of pions, a detailed study of the process J/ψ → ppπ + π − is performed. It is found that a systematic uncertainty of 1% per pion is a reasonable estimation and thus the corresponding systematic uncertainty for the four charged pions in this analysis is set to 4%.
d. External branching fractions The uncertainties of the BFs entering this analysis, namely those of the decays J/ψ → γη c and ω → π + π − π 0 are taken from the world average values published in Ref. [14] and treated as systematic uncertainties. The uncertainty of B(π 0 → γγ) is negligible and is therefore excluded from Table III. It should be noted here that the uncertainty on the BF J/ψ → γη c is the dominant uncertainty in this analysis. e. Kinematic fit To estimate the systematic uncertainty of the kinematic fit, the charged track helix parameters in simulated data are smeared with a Gaussian function so that their distributions in MC and data match. The difference in efficiency between applying and not applying this correction for the given requirement on the χ 2 6C value of the kinematic fit is found to be 1.2% and is taken as the systematic uncertainty.
f. Q-factor method To estimate the systematic uncertainty introduced by the Q-factor method, tests with different dedicated MC samples are performed. Background and signal MC samples of different compositions are generated and subjected to the Q-factor method. The largest deviation between the number of generated signal events and the sum of the obtained Q-factors is obtained using a background sample that contains a peaking background contribution at the mass of the η c among other phase space-like contributions. The deviation is determined to be 0.9%, which is taken as the systematic uncertainty of the method.
g. η c damping factor To estimate the uncertainty due to the η c damping factor, an alternative parameterization of this factor is used. For this test the CLEO parameterization is exchanged by the function , where E γ denotes the energy of the radiative photon and E γ,0 is the most probable photon energy, corresponding to the mass of the η c [28]. The number of η c events and the efficiency are extracted from this fit and the difference between the resulting BF and the nominal result is measured to be 14.2%, which is assigned as a systematic uncertainty. The mass and width of the η c are left floating in this fit and their differences to the nominal result are considered as systematic uncertainties for the measurement of the resonance parameters.
h. Fit range While for the nominal result only events in the region m(ωω) > 2.65 GeV/c 2 are used, this lower mass limit is varied by ±50 MeV/c 2 to estimate the uncertainty connected to the choice of the mass requirement. The partial wave fit is re-performed for both scenarios and the largest deviation in the yield of the η c candidates is found to be 1.4 %. This value is taken as the systematic uncertainty due to the choice of the fitting mass range. Similarly, also the mass and width of the η c are re-evaluated and the differences to the nominal result are taken as systematic uncertainties.
i. η c resonance parameters We also re-performed the fit using fixed values for the resonance parameters of the η c . For this study, mass and width are set to their world average values published in Ref. [14] and a deviation of 1.0 % for the obtained yield of the η c signal is found, which is taken as a systematic uncertainty for the BF discussed in this paper.
j. Selection of fit hypothesis The results for the yield, mass and width of the η c are additionally evaluated for the second best hypothesis to estimate the uncertainty due to the choice of the hypothesis. The difference in the obtained number of observed η c events has a negligible effect on the extracted BF. The deviation of the mass is determined to be 0.6 MeV/c 2 while the width differs by 0.3 MeV, which are taken as systematic uncertainties.
k. Detector resolution To estimate the effect of the detector resolution we perform a dedicated MC study. Using all parameters obtained from the best PWA fit to data, we generate an MC sample and propagate the events through the BESIII detector simulation and reconstruction using the same criteria as for beam data. After performing a PWA fit to the reconstructed and selected MC sample we obtain a difference of 2.0 MeV/c 2 for the mass and 3.6 MeV for the width of the η c between the generated and reconstructed data sample. We use this deviation as an estimation for the systematic uncertainty due to the detector resolution.
The last quoted uncertainty corresponds to the error of the J/ψ → γη c BF and is the dominant uncertainty of this measurement.