Observation of OZI suppressed decays $\chi_{cJ}\to\omega\phi$

Decays $\chi_{cJ}~(J=0,1,2)\to\omega\phi$ are studied using $(448.1\pm2.9)\times 10^{6} ~\psi(3686)$ events collected with the BESIII detector in 2009 and 2012. In addition to the previously established $\chi_{c0}\to\omega\phi$, first observation of $\chi_{c1} \to \omega \phi$ is reported in this paper. The measured product branching fractions are ${\cal{B}}(\psi(3686)\to\gamma\chi_{c0})\times{\cal{B}}(\chi_{c0}\to\omega\phi)=(13.83\pm 0.70\pm 1.01)\times10^{-6}$ and ${\cal{B}}(\psi(3686)\to\gamma\chi_{c1})\times{\cal{B}}(\chi_{c1}\to\omega\phi)=(2.67\pm 0.31\pm 0.27)\times10^{-6}$, and the absolute branching fractions are ${\cal{B}}(\chi_{c0}\to\omega\phi)=(13.84\pm 0.70\pm 1.08)\times10^{-5}$ and ${\cal{B}}(\chi_{c1}\to\omega\phi)=(2.80\pm 0.32\pm 0.30)\times10^{-5}$. We also find a strong evidence for $\chi_{c2}\to\omega\phi$ with a statistical significance of 4.8$\sigma$, and the corresponding product and absolute branching fractions are measured to be ${\cal{B}}(\psi(3686)\to\gamma\chi_{c2})\times{\cal{B}}(\chi_{c2}\to\omega\phi)=(0.91\pm0.23\pm0.12)\times10^{-6} $ and ${\cal{B}}(\chi_{c2}\to\omega\phi)=(1.00\pm0.25\pm0.14)\times10^{-5}$. Here, the first errors are statistical and the second ones systematic.


I. INTRODUCTION
The lowest triplet P -wave states of charmonium (the cc bound state), χ cJ (1P ), with quantum numbers I G J P C = 0 + J ++ and J = 0, 1, and 2, can be found abundantly in the electromagnetic decays ψ(3686) → γχ cJ with an approximate branching fraction of 30% [1]. The ψ(3686) meson can be directly produced at the e + e − colliders, such as the BEPCII [2], where the χ cJ meson are easily accessible by the electromagnetic decays ψ(3686) → γχ cJ .
The hadronic χ cJ decays are important probes of the strong force dynamics. Firstly, as is well known, the mass of the charm quark, m c ∼ 1.5 GeV/c 2 , is just between the perturbative and nonperturbative QCD domains in theoretical calculations. Due to the complexity and entanglement of the long-and short-distance contributions, there exist vary large theoretical uncertainties of branching ratios for the χ cJ → V V decays [3][4][5][6][7][8][9]. (In this paper, the symbol of V denotes the ω and φ mesons). The hadronic χ cJ decays provide a prospective laboratory to limit theoretical parameters and test various phenomenological models. Secondly, the χ cJ mesons have the same quantum numbers J P C as some glueballs and hybrids, although none of the glueball and hybrid states is seen experimentally until now. The hadronic χ cJ → V V decays are ideal objects to exploit the glueball-qq mixing and the quark-gluon coupling of the strong interactions at the relatively low energies. Thirdly, the χ cJ mesons are below the open-charm threshold. Most of the hadronic χ cJ decay modes are suppressed by the Okubo-Zweig-Iizuka (OZI) rule [10]. It is shown in the previous theoretical researches that the contributions from the intermediate glueballs or hadronic loops can scuttle the OZI rule in the χ cJ → V V decays [11][12][13][14], and avoid the so-called helicity selection (HS) rule (also called the "naturalness" which is defined as σ = (−1) S P by Ref. [15], where S and P are respectively the spin and parity of the particle.) in the χ c1 → V V decays [8,9].
The 'ideal mixing' between the φ and ω mesons is almost true in practice. So, the χ cJ , φ and ω mesons differ from each other in their components according to the quark model assignments. This fact causes the χ cJ → ωφ decay modes to be doubly OZI (DOZI) suppressed, and results in the branching fractions for the χ cJ → ωφ decays much less than those for the singly OZI-suppressed χ cJ → ωω, φφ decays [1,16]. The DOZI-suppressed χ cJ → ωφ decays have been observed based on the 106×10 6 ψ(3686) events accumulated with the BESIII detector in 2009, with significances of 10 σ, 4.1 σ and 1.5 σ for the χ c0 , χ c1 and χ c2 decays, respectively [16].
The precise measurements of χ cJ → ωφ decays will be helpful in clarifying the influence of the long-distance contributions in this energy region, understanding the vexatious theoretical dilemma surrounding the OZI and HS rules, and providing some useful information on the ω-φ mixing. In this paper, the χ cJ → ωφ decays will be reinvestigated via the radiative transitions ψ(3686) → γχ cJ with combined experimental data, i.e., (448.1±2.9)×10 6 ψ(3686) events collected with the BESIII detector during 2009 and 2012 [17].

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector operating at the BEPCII collider is described in detail in Ref. [2]. The detector is cylindrically symmetric and covers 93% of 4π solid angle. It consists of the following four sub-detectors: a 43-layer main drift chamber (MDC), which is used to determine momentum of the charged tracks with a resolution of 0.5% at 1 GeV/c in the axial magnetic field of 1 T; a plastic scintillator time-of-flight system (TOF), with a time resolution of 80 ps (110 ps) in the barrel (endcaps); an electromagnetic calorimeter (EMC) consisting of 6240 CsI(Tl) crystals, with photon energy resolution at 1 GeV of 2.5% (5%) in the barrel (endcaps); and a muon counter consisting of 9 (8) layers of resistive plate chambers in the barrel (endcaps), with position resolution of 2 cm.
The geant4-based [18,19] Monte Carlo (MC) simulation software boost [20] includes the geometry and material description of the BESIII detectors, the detector response and digitization models, as well as a database that keeps track of of the running conditions and the detector performance. MC samples are used to optimize the selection criteria, evaluate the signal efficiency, and estimate physics backgrounds. An inclusive MC sample of 5.1 × 10 8 ψ(3686) events is used for the background studies. The ψ(3686) resonance is produced by the event generator kkmc [21], where the initial state radiation is included, and the decays are simulated by evtgen [22] with known branching fractions taken from Ref. [1], while the unmeasured decays are generated according to lundcharm [23]. The signal is simulated with the decay ψ(3686) → γχ cJ generated assuming an electric-pole (E1) transition. The decay χ cJ → ωφ is generated using helamp [22], the helicity amplitude model where the angular correlation between ω decay and φ decay has been considered. Ref. [16] shows that the model describes the experimental angular distribution well. We assume χ cJ → ωφ and χ cJ → φφ have the same helicity amplitudes with the same helamp parameters. In addition, χ cJ states are simulated using a relativistic Breit-Wigner incorporated within the helicity amplitudes in the evtgen package [22]. The background decays χ cJ → ωK + K − , φπ + π − π 0 , and the nonresonant decay χ cJ → K + K − π + π − π 0 are generated using the phase space model.

III. EVENT SELECTION
In this analysis, the φ mesons are reconstructed by K + K − , while ω by π + π − π 0 . Event candidates are required to have four well-reconstructed tracks from charged particles with zero net charge, and at least three good photon candidates.
A charged track reconstructed from MDC hits should have the polar angle, θ, |cosθ| < 0.93 and pass within ±10 cm of the interaction point along the beam direction and within 1 cm in the plane perpendicular to the beam. To separate K ± from π ± , we require that at least one track is identified as a kaon using dE/dx and TOF information. If the identified kaon has a positive (negative) charge, the second kaon is found by searching for a combination that minimizes |M K + K − −M φ |, among all identified kaons and the negative (positive) charged tracks, where M K + K − is the invariant mass of the identified kaon and an unidentified track with kaon mass hypothesis, and M φ is the nominal φ mass [1]. The remaining two charged tracks are assumed to be pions.
The photon energy deposit is required to be at least 25 MeV in the barrel region of the EMC (|cosθ| < 0.80) or 50 MeV in the EMC endcaps (0.86 < |cosθ| < 0.92). To suppress electronic noise and energy deposits unrelated to the event, the EMC time t of the photon candidates must be in coincidence with collision events within the range 0 ≤ t ≤ 700 ns. At least three photons are required in an event.
In order to improve the mass resolution, a fourconstraint (4C) kinematic fit is performed by assuming energy-momentum conservation for the ψ(3686) → 3γK + K − π + π − process. If the number of photons is larger than three, then looping all 3γK + K − π + π − combinations and the one with the smallest χ 2 4C is chosen.
The event is kept for further anal- In addition, is applied to suppressed the background with an extra photon in the final state.
The π 0 candidates are selected from the three γγ combinations as the pair with the minimum |M γγ − M π 0 |, where M π 0 is the nominal π 0 mass [1]. Figure 1(a) shows the plot of the K + K − versus π + π − π 0 invariant mass for the selected events in the χ cJ signal region ([3.3, 3.6] GeV/c 2 ), and a clear accumulation at the ω and φ masses is observed. The bottom-central square GeV/c 2 is taken as the ωφ signal region (labeled as D), and the five squares around the signal region are taken as the sideband regions (labeled as A, B and C), where M π + π − π 0 (M K + K − ) is the invariant mass of π + π − π 0 (K + K − ). The M K + K − distribution with M π + π − π 0 in the ω signal region in Fig. 1(b) shows a clear φ peak. Correspondingly, the M π + π − π 0 distribution with M K + K − within the φ signal region, as shown in Fig. 1(c), indicates clear ω and φ peaks. The latter is from the decay χ cJ → φφ → K + K − π + π − π 0 . Figure 2 shows the invariant mass spectrum M K + K − π + π − π 0 for events in the ωφ sideband regions (sub-figures labeled A, B, and C) and the signal region (sub-figure labeled D) with clear χ cJ peaks in all plots.
Analysis of the ψ(3686) inclusive MC sample indicates that the peaking background in the χ cJ signal region can be described by the sideband events. The data collected at √ s = 3.65 GeV with an integrated luminosity of approximately 1/15 of the ψ(3686) data are used to investigate non-resonant continuum background. After the same event selection criteria are applied, only a few events survive, and they do not have any obvious enhancements in the χ cJ mass region.

IV. SIGNAL EXTRACTION
The number of the χ cJ → ωφ events is determined by fitting the M K + K − π + π − π 0 distributions within the ωφ signal region (labeled as D in Fig. 1(a)). The signal is described by the MC simulated shape convolved with a Gaussian function, which is used to account for the difference in the χ cJ mass and resolution between data and MC simulation. The parameters of the Gaussian function are obtained using the sample ψ(3686) → γχ cJ → γφφ → γπ + π − π 0 K + K − .
The peaking backgrounds from the χ cJ → ωK + K − , φπ + π − π 0 , and the non-resonant K + K − π + π − π 0 background are estimated using the sideband regions labeled A, B, and C in Fig. 1(a). The total peaking background contribution, N bkg , is the sum calculated as: where N ωK + K − bkg , N φπ + π − π 0 bkg , and N n−r bkg are numbers of the aforementioned peaking background contributions. The contributions are determined using the following equations: where N A , N B , and N C are the numbers of the fitted χ cJ events in the A, B, and C regions, respectively; f C→A , f C→B , f C→D , f A→D , and f B→D are the relative scaling factors for the different regions. The factors are estimated using the corresponding MC simulation of χ cJ → K + K − π + π − π 0 , ωK + K − , and φπ + π − π 0 . For example, f A→D is the ratio of the χ cJ → ωK + K − yields between the D and A regions.
We perform a simultaneous unbinned maximum likelihood fit to the M K + K − π + π − π 0 distributions in the signal and sideband regions. The result of the fit is shown in Fig. 2. The parameters of the Gaussian functions accounting for the difference between data and MC simulation are assumed to be the same for the signal and sidebands. The shape of the distributions outside the χ cJ peaks is described by a polynomial function. The statistical significance of the χ c1 (χ c2 ) signal is determined by comparing the −2lnL value with the one from the fit without the χ c1 (χ c2 ) signal component, and considering the change in the number of degrees of freedom. The results are 12.3σ and 4.8σ for χ c1 and χ c2 , respectively. The extracted numbers of the χ cJ → ωφ events are given in Table I  The product branching fractions, B(ψ(3686) → γχ cJ ) ×B(χ cJ → ωφ) = B 1 × B 2 , are calculated as (5) where N ψ(3686) is the number of ψ(3686) events, B(ω), B(φ), and B(π 0 ) are the branching fractions of ω → π + π − π 0 , φ → K + K − , and π 0 → γγ, respectively [1]. The corresponding detection efficiencies, ǫ, are obtained from the MC simulations. The results for the product branching fractions are listed in Table I.
By using the world average values of B(ψ(3686) → γχ cJ ), the absolute branching fractions of χ cJ → ωφ are determined and also listed in Table I.

V. SYSTEMATIC UNCERTAINTIES
The contribution of systematic effects on the product branching fractions from various sources is described in the following: 1. The tracking efficiency for π and K is investigated using control samples of J/ψ → ρπ [16] and ψ(3686) → π + π − K + K − , respectively. The difference in the efficiency for the track reconstruction between data and MC simulation is 1.0% per pion and per kaon. Assuming they are all correlated, the uncertainty due to tracking efficiency is 4%.
2. The particle identification (PID) efficiency for kaons is investigated with control samples of J/ψ → K * (892) 0 K 0 S + c.c., and the systematic uncertainty is determined to be 1% per kaon track [16]. In this analysis, only one of the two charged tracks is required to be identified as kaon. The bias on one of the two tracks being a kaon track will be much smaller than 1%. Therefore, the uncertainty due to PID efficiency is negligible.
3. The uncertainty of the photon reconstruction efficiency is studied using J/ψ → ρπ [24]. The difference between data and MC simulation is found to be 1.0% per photon, and the value 3% is taken as the systematic uncertainty.
The uncertainty of assuming ψ(3686) → γχ c1,2 as pure E1 transition is studied by taking the higherorder multipile amplitudes contribution [25] into account in the MC simulation. The resulting efficiency difference of 0.9% for χ c1 , and 0.5% for χ c2 , are taken as this systematic uncertainty. The uncertainty of modeling χ cJ → ωφ is studied by changing the model from helamp to a pure phase , detection efficiency (ǫ), the product branching fraction B1 × B2 = B(ψ(3686) → γχcJ ) × B(χcJ → ωφ), and the absolute branching fraction B(χcJ → ωφ). Here, the first uncertainty is statistical and the second systematic.
The total systematic uncertainties of MC generator are obtained to be 4.1% for χ c0 , 5.7% for χ c1 , and 1.4% for χ c2 by summing all individual contributions in quadrature, assuming two sources to be independent.
5. The uncertainty related to the π 0 mass window is studied by fitting the π 0 mass distribution of data and signal MC for the control sample ψ(3686) → π + π − π 0 . We obtain the π 0 detection efficiency, which is the ratio of the number of π 0 events selected with and without the π 0 mass window requirement, determined by integrating the fitted signal shape. The difference in the efficiency between data and MC simulation is 0.8%.
6. The uncertainty related to the M recoil π + π − mass window requirement is studied with the control sample ψ(3686) → π + π − J/ψ, J/ψ → µ + µ − . We obtain the M recoil π + π − detection efficiency, which is the ratio of the number of events with and without the M recoil π + π − mass window. The difference in the efficiency between data and MC simulation is 0.4%. 7. The systematic uncertainties associated with the 4C kinematic fit are studied with the track helix parameter correction method, as described in Ref. [26]. In the standard analysis, these corrections are applied. The difference of the MC signal efficiencies with the uncorrected track parameters are 1.5%, 1.9%, and 2.3% for χ c0 , χ c1 , and χ c2 decays, respectively. These values are taken as the uncertainties associated with the 4C kinematic fit.
8. The uncertainty related to the fitting comes from the fit range, resolution, ω and φ mass windows, sideband regions, and remaining backgrounds.
(a) The uncertainty due to the fit range is estimated by changing the range by ±5 MeV/c 2 in the mass spectrum. The largest differences for the branching ratios are 1.0% for χ c0 , 3.0% for χ c1 , and 10.0% for χ c2 . These numbers are assigned as the corresponding systematic uncertainties.
(b) The systematic effects from the detector resolution difference between data and MC simulation are studied with the control sample ψ(3686) → γχ cJ → φφ → K + K − π + π − π 0 . We change the difference by one standard deviation. No changes are found for the χ c0,1,2 signal yields and these systematic uncertainties are neglected.
(d) The uncertainty due to background estimates using the sidebands can be divided in two groups. One is due to the sideband ranges, the other is due to contributions of various intermediate states in χ cJ → K + K − π + π − π 0 in the MC simulation used to extract the scaling factors. The former can be estimated by changing the sideband range. By changing the mass region of M π + π − π 0 from [0. .072] GeV/c 2 , the differences of χ c0,1,2 signal yields are 0.6%, 4.6%, and 5.1%, respectively. For the non-resonant χ cJ → K + K − π + π − π 0 , a phase space process was used. The experimental distributions indicate the contribution of the intermediate states involving K * (892): χ cJ → K * 0K * 0 π 0 and χ cJ → K * + K − π + π − + c.c.. The corresponding MC distributions are mixed with the phase space model according to the ratios estimated from the fits to data to recalculate the scale factors related to the region C. The differences of the χ c0,1,2 signal yields are 0.2%, 2.0%, and 4.2%, respectively. The resulting differences due to the two preceding effects are found to be 0.6%, 5.0%, and 6.6% for χ c0,1,2 , respectively.
(e) The uncertainty from the non-χ cJ background is estimated by changing the polynomial from first to second order in fitting M K + K − π + π − π 0 mass spectrum. The differences in the final results are 0.4%, 0.5%, and 0.2%, respectively.
10. The number of ψ(3686) events is estimated from the number of inclusive hadronic events, as described in Ref. [17]. The uncertainty of the total number of ψ(3686) events is 0.6%. Table II summarizes the systematic uncertainties and their sources for the product branching fractions. The total systematic uncertainties are obtained by summing all individual contributions in quadrature, assuming all sources to be independent. For the uncertainties of absolute branching fractions χ cJ → ωφ, the uncertainty arising from ψ(3686) → γχ cJ transition rate is added.