Observation of η c → ωω in J = ψ → γωω

M. Ablikim, M. N. Achasov, S. Ahmed, M. Albrecht , M. Alekseev, A. Amoroso, F. F. An, Q. An, Y. Bai, O. Bakina, R. Baldini Ferroli, Y. Ban, K. Begzsuren, D.W. Bennett, J. V. Bennett, N. Berger, M. Bertani, D. Bettoni, F. Bianchi, E. Boger, I. Boyko, R. A. Briere, H. Cai, X. Cai, A. Calcaterra, G. F. Cao, S. A. Cetin, J. Chai, J. F. Chang, W. L. Chang, G. Chelkov, G. Chen, H. S. Chen, J. C. Chen, M. L. Chen, S. J. Chen, X. R. Chen, Y. B. Chen, W. Cheng, X. K. Chu, G. Cibinetto, F. Cossio, H. L. Dai, J. P. Dai, A. Dbeyssi, D. Dedovich, Z. Y. Deng, A. Denig, I. Denysenko, M. Destefanis, F. De Mori, Y. Ding, C. Dong, J. Dong, L. Y. Dong, M. Y. Dong, Z. L. Dou, S. X. Du, J. Z. Fan, J. Fang, S. S. Fang, Y. Fang, R. Farinelli, L. Fava, F. Feldbauer, G. Felici, C. Q. Feng, M. Fritsch, C. D. Fu, Q. Gao, X. L. Gao, Y. Gao, Y. G. Gao, Z. Gao, B. Garillon, I. Garzia, A. Gilman, K. Goetzen, L. Gong, W. X. Gong, W. Gradl, M. Greco, L. M. Gu, M. H. Gu, Y. T. Gu, A. Q. Guo, L. B. Guo, R. P. Guo, Y. P. Guo, A. Guskov, S. Han, X. Q. Hao, F. A. Harris, K. L. He, F. H. Heinsius, T. Held, Y. K. Heng, Z. L. Hou, H. M. Hu, J. F. Hu, T. Hu, Y. Hu, G. S. Huang, J. S. Huang, X. T. Huang, X. Z. Huang, Z. L. Huang, T. Hussain, N. Hüsken, W. Ikegami Andersson, W. Imoehl, M. Irshad, Q. Ji, Q. P. Ji, X. B. Ji, X. L. Ji, H. L. Jiang, X. S. Jiang, X. Y. Jiang, J. B. Jiao, Z. Jiao, D. P. Jin, S. Jin, Y. Jin, T. Johansson, N. Kalantar-Nayestanaki, X. S. Kang, M. Kavatsyuk, B. C. Ke, I. K. Keshk, T. Khan, A. Khoukaz, P. Kiese, R. Kiuchi, R. Kliemt, L. Koch, O. B. Kolcu, B. Kopf, M. Kuemmel, M. Kuessner, A. Kupsc, M. Kurth, W. Kühn, J. S. Lange, P. Larin, L. Lavezzi, S. Leiber, H. Leithoff, C. Li, Cheng Li, D. M. Li, F. Li, F. Y. Li, G. Li, H. B. Li, H. J. Li, J. C. Li, J. W. Li, K. J. Li, Kang Li, Ke Li, L. K. Li, Lei Li, P. L. Li, P. R. Li, Q. Y. Li, W. D. Li, W. G. Li, X. L. Li, X. N. Li, X. Q. Li, X. L. Li, Z. B. Li, H. Liang, Y. F. Liang, Y. T. Liang, G. R. Liao, L. Z. Liao, J. Libby, C. X. Lin, D. X. Lin, B. Liu, B. J. Liu, C. X. Liu, D. Liu, D. Y. Liu, F. H. Liu, Fang Liu, Feng Liu, H. B. Liu, H. L. Liu, H. M. Liu, Huanhuan Liu, Huihui Liu, J. B. Liu, J. Y. Liu, K. Y. Liu, Ke Liu, L. D. Liu, Q. Liu, S. B. Liu, X. Liu, Y. B. Liu, Z. A. Liu, Zhiqing Liu, Y. F. Long, X. C. Lou, H. J. Lu, J. G. Lu, Y. Lu, Y. P. Lu, C. L. Luo, M. X. Luo, P. W. Luo, T. Luo, X. L. Luo, S. Lusso, X. R. Lyu, F. C. Ma, H. L. Ma, L. L. Ma, M.M. Ma, Q. M. Ma, X. N. Ma, X. Y. Ma, Y. M. Ma, F. E. Maas, M. Maggiora, S. Maldaner, Q. A. Malik, A. Mangoni, Y. J. Mao, Z. P. Mao, S. Marcello, Z. X. Meng, J. G. Messchendorp, G. Mezzadri, J. Min, T. J. Min, R. E. Mitchell, X. H. Mo, Y. J. Mo, C. Morales Morales, N. Yu. Muchnoi, H. Muramatsu, A. Mustafa, S. Nakhoul, Y. Nefedov, F. Nerling, I. B. Nikolaev, Z. Ning, S. Nisar, S. L. Niu, X. Y. Niu, S. L. Olsen, Q. Ouyang, S. Pacetti, Y. Pan, M. Papenbrock, P. Patteri, M. Pelizaeus, J. Pellegrino, H. P. Peng, Z. Y. Peng, K. Peters, J. Pettersson, J. L. Ping, R. G. Ping, A. Pitka, R. Poling, V. Prasad, M. Qi, T. Y. Qi, S. Qian, C. F. Qiao, N. Qin, X. S. Qin, Z. H. Qin, J. F. Qiu, S. Q. Qu, K. H. Rashid, C. F. Redmer, M. Richter, M. Ripka, A. Rivetti, M. Rolo, G. Rong, Ch. Rosner, M. Rump, A. Sarantsev, M. Savrié, K. Schoenning, W. Shan, X. Y. Shan, M. Shao, C. P. Shen, P. X. Shen, X. Y. Shen, H. Y. Sheng, X. Shi, X. D. Shi, J. J. Song, Q. Q. Song, X. Y. Song, S. Sosio, C. Sowa, S. Spataro, F. F. Sui, G. X. Sun, J. F. Sun, L. Sun, S. S. Sun, X. H. Sun, Y. J. Sun, Y. K. Sun, Y. Z. Sun, Z. J. Sun, Z. T. Sun, Y. T. Tan, C. J. Tang, G. Y. Tang, X. Tang, B. Tsednee, I. Uman, B. Wang, B. L. Wang, C.W. Wang, D. Wang, D. Y. Wang, H. H. Wang, K. Wang, L. L. Wang, L. S. Wang, M. Wang, Meng Wang, P. Wang, P. L. Wang, W. P. Wang, X. F. Wang, Y. Wang, Y. F. Wang, Z. Wang, Z. G. Wang, Z. Y. Wang, Zongyuan Wang, T. Weber, D. H. Wei, P. Weidenkaff, S. P. Wen, U. Wiedner, M. Wolke, L. H. Wu, L. J. Wu, Z. Wu, L. Xia, Y. Xia, Y. J. Xiao, Z. J. Xiao, Y. G. Xie, Y. H. Xie, X. A. Xiong, Q. L. Xiu, G. F. Xu, J. J. Xu, L. Xu, Q. J. Xu, X. P. Xu, F. Yan, L. Yan, W. B. Yan, W. C. Yan, Y. H. Yan, H. J. Yang, H. X. Yang, L. Yang, R. X. Yang, S. L. Yang, Y. H. Yang, Y. X. Yang, Yifan Yang, Z. Q. Yang, M. Ye, M. H. Ye, J. H. Yin, Z. Y. You, B. X. Yu, C. X. Yu, J. S. Yu, C. Z. Yuan, Y. Yuan, A. Yuncu, A. A. Zafar, Y. Zeng, B. X. Zhang, B. Y. Zhang, C. C. Zhang, D. H. Zhang, H. H. Zhang, H. Y. Zhang, J. Zhang, J. L. Zhang, J. Q. Zhang, J. W. Zhang, J. Y. Zhang, J. Z. Zhang, K. Zhang, L. Zhang, S. F. Zhang, T. J. Zhang, X. Y. Zhang, Y. Zhang, Y. H. Zhang, Y. T. Zhang, Yang Zhang, Yao Zhang, Yu Zhang, Z. H. Zhang, Z. P. Zhang, Z. Y. Zhang, G. Zhao, J. W. Zhao, J. Y. Zhao, J. Z. Zhao, Lei Zhao, Ling Zhao, M. G. Zhao, Q. Zhao, S. J. Zhao, T. C. Zhao, Y. B. Zhao, Z. G. Zhao, A. Zhemchugov, B. Zheng, J. P. Zheng, Y. H. Zheng, B. Zhong, L. Zhou, Q. Zhou, X. Zhou, X. K. Zhou, X. R. Zhou, Xiaoyu Zhou, PHYSICAL REVIEW D 100, 052012 (2019)


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 roughly or with large uncertainties, and the observed BFs sum up to only about 57%. Also the observed mass and decay width show a large variation from experiment to experiment, and may depend 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] and [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]. Recently published predictions of BFs for the decay modes η c → ϕϕ and η c → ρρ are much smaller than those observed experimentally [7]. These predictions are based on next-to-leading order (NLO) perturbative QCD calculations and for the first time also include the so-called higher-twist contributions. The latter were found to 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 in Ref. [ Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP 3 . 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 dataset used for this analysis contains a total of ð1310.6 AE 7.0Þ × 10 6 J=ψ events [8] produced in direct e þ e − annihilations and recorded with the Beijing Spectrometer III (BESIII) detector. The mass, width, and 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 Beijing Electron Positron Collider II (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 small-cell gas drift chamber (MDC), a time-of-flight system (TOF), an electromagnetic calorimeter (EMC), as well as a system for muon identification.
The 43-layer MDC directly surrounds the beam pipe and 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, the plastic scintillator based TOF system for particle identification followed by the CsI(Tl)-based EMC is mounted. The time-of-flight system consists of 176 scintillator bars with a length of 2.4 m, arranged in a two-layer, barrelshaped geometry, and 96 fan-shaped scintillators in the end caps. All plastic scintillators of the TOF 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 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 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) datasets 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 [11], 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 [12,13]. 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 BESEVTGEN [14], which is based on the EVTGEN package [15].

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 j cos θj < 0.93.
Pion candidates are selected from all good charged tracks, by exploiting the capabilities of particle identification 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 (j cos θj < 0.8) and 50 MeV in the end cap regions (0.86 < j cos θj < 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 that 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 miscombination to a level less than 1%. 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. [16]. 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 event-based 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.

IV. BACKGROUND SUBTRACTION
A sophisticated event-based method for background subtraction proposed in Ref. [17] is applied to events for which both three-pion invariant masses are located within a range of AE80 MeV around the nominal ω mass. Simpler methods, such as a two-dimensional side band subtraction, mostly require the analysis of a binned dataset, 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-to-background ratio Q in a very small cell of the available phasespace around each event. Therefore, a distinct kinematic variable is needed, for which parametrizations 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 corresponding ω candidates' helicity frame is given as where T π denotes the kinetic energy of the corresponding pion [18] and c is the speed of light. The parameter λ 0 takes its maximum value λ 0 max for totally symmetric decays with an angle of 120°between any pion pair (see Ref. [18]). 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 twodimensional 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 FIG. 1. 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 Sec. IV. 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 subfigures (c) and (d), where a good agreement can be seen. Figure 3 shows the invariant 3π mass and the normalized λ distribution for all preselected 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 Q-factor 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 cross-check 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 [19] using the helicity formalism [20] 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 spacelike 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 spacelike 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. [16], and ϑ denotes the polar angle in the respective helicity frame. The d-matrices do not depend on the azimuthal angle φ in contrast to the usual Wigner D-matrices, so that only the dependence on the polar angle ϑ remains. 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. [21][22][23], where h…i denotes the Clebsch-Gordan coefficients and B L ðqÞ are the Blatt-Weisskopf barrier factors as defined in Ref. [24]. 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 PC , L, S values is performed using the transformation where S is the total spin of the ωω system [20]. 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 spacelike 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 M1-transition matrix element [25], 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 toward 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 [26].
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 λ ω μ ðφ n ; ϑ n ; 0Þ ·λ μ ; ð8Þ where only the case μ ¼ 0 is allowed for this decay [27]. 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.
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 fit contribution to obtain a set of independent parameters.

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 [27] 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 whichn appears, so that the mean weight of an MC event is approximately 1 after the likelihood has been maximized. The integrals appearing in then 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 Q-factor 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 [28], can be used.
The negative log-likelihood function, which is actually minimized, now reads as

C. Fit strategy
Since the composition of the nonresonant contribution is not known a priori, different hypotheses are fitted to the selected dataset. These contain the η c component and one up to a maximum of four different nonresonant components. The nonresonant components are assumed to have the J P quantum numbers 0 − , 0 þ , 1 þ , or 2 þ , so that the most simple hypothesis is given as fη c ; 0 − g, and the most complex one by fη c ; 0 − ; 0 þ ; 1 þ ; 2 þ g. 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 PC ¼ 1 −þ ) as well as pseudo tensor contributions (J PC ¼ 2 −þ ) are tested and found to be insignificant. In total, about 45 hypotheses with different combinations of contributing waves were tested.
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 [29]. 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 overfitting. Theoretical considerations show [29] 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. 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 AE 58 η c events, which is the number    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 Sec. 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 AE 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. [16] 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 spacelike 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 parametrization of this factor is used. For this test, the CLEO parametrization is exchanged by the function E 2 γ;0 = ðE γ E γ;0 þ ðE γ − E γ;0 Þ 2 Þ, 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 [30]. 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 AE50 MeV=c 2 to estimate the uncertainty connected to the choice of the mass requirement. The partial wave fit is reperformed 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 reevaluated and the differences to the nominal result are taken as systematic uncertainties. (i) η c resonance parameters: We also reperformed 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. [16] 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 samples. We use this deviation as an estimation for the systematic uncertainty due to the detector resolution.

VII. BRANCHING FRACTION
Using the obtained results of the best fit to the data and the systematic uncertainties discussed above, the product BF of the decay chain J=ψ → γη c → γωω is determined as BðJ=ψ → γη c Þ · Bðη c → ωωÞ ¼ N η c N J=ψ B 2 ðω → π þ π − π 0 ÞB 2 ðπ 0 → γγÞϵ where the BFs Bðω → π þ π − π 0 Þ and Bðπ 0 → γγÞ are taken from Ref. [16], N η c is the η c signal yield determined from the best PWA fit, ϵ ¼ 3.42% is the detection and reconstruction efficiency, and N J=ψ ¼ ð1310.6 AE 7.0Þ × 10 6 [8] is the number of J=ψ events. Taking into account the measured BF for the J=ψ → γη c decay, which has large uncertainties, the BF of the η c decay is given by Bðη c → ωωÞ ¼ ð2.88 AE 0.10 stat: AE 0.46 syst: AE 0.68 ext: Þ × 10 −3 : The last quoted uncertainty corresponds to the error of the J=ψ → γη c BF and is the dominant uncertainty of this measurement.

VIII. MASS AND WIDTH OF THE η c
The mass and width of the η c are left as free parameters in the PWA fits. The systematic uncertainty of the extracted values is estimated from alternative fits with different fit ranges, different fit hypothesis, and the usage of the alternative damping factor. All sources of systematic uncertainties are assumed to be independent, and thus their deviations from the nominal result are added in quadrature. The values are found to be where the first uncertainties are statistical and the second systematic. The mass and width are consistent with the world average values.

IX. SUMMARY AND DISCUSSION
Using a sample of ð1310.6 AE 7.0Þ × 10 6 J=ψ events accumulated with the BESIII detector, we report the first observation of the decay η c → ωω in the process J=ψ → γωω. By means of a PWA, the branching fraction of η c → ωω is measured to be Bðη c → ωωÞ ¼ ð2.88 AE 0.10 stat: AE 0.46 syst: AE 0.68 ext: Þ × 10 −3 , where the external uncertainty refers to that arising from the branching fraction of the decay J=ψ → γη c . The obtained value is about 1 order of magnitude larger than what is expected from NLO perturbative QCD calculations including higher twist contributions. The mass and width of the η c are determined to be M ¼ ð2985.9 AE 0.7 stat: AE 2.1 syst: Þ MeV=c 2 and Γ ¼ ð33.8 AE 1.6 stat: AE 4.1 syst: Þ MeV. The extracted values for the mass and width of the η c are in good agreement with the world average values. This measurement provides new insights into the decay characteristics of charmonium resonances.